Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Comparison of two vectors resulted after simulation

I would like to apply the Rejection sampling method to simulate a random vector Y=(Y_1, Y_2) of a uniform distribution from a unit disc D = { (X_1 , X_2) \in R^2: \sqrt{x^2_1 + x^2_2} ≤ 1} such that X = (X_1 , X_ 2) is random vector of a uniform distribution in the square S = [−1, 1]^2 and the joint density f(y_1,y_2) = \frac{1}{\pi} 1_{D(y_1,y_2)}.

In the rejection method, we accept a sample generally if f(x) \leq C * g(x). I am using the following code to :

x=runif(100,-1,1)
y=runif(100,-1,1)

d=data.frame(x=x,y=y)
disc_sample=d[(d$x^2+d$y^2)<1,]
plot(disc_sample)

I have two questions:

{Using the above code, logically, the size of d should be greater than the size of disc_sample but when I call both of them I see there are 100 elements in each one of them. How could this be possible. Why the sizes are the same.} THIS PART IS SOLVED, thanks to the comment below.

The question now

Also, how could I reformulate my code to give me the total number of samples needed to get 100 samples follow the condition. i.e to give me the number of samples rejected until I got the 100 needed sample?

Thanks to the answer of r2evans but I am looking to write something simpler, a while loop to store all possible samples inside a matrix or a data frame instead of a list then to call from that data frame just the samples follow the condition. I modified the code from the answer without the use of the lists and without sapply function but it is not giving the needed result, it yields only one row.

i=0
samps <- data.frame()
goods <- data.frame()
nr <- 0L
sampsize <- 100L
needs <- 100L
while (i < needs) {
  samps <- data.frame(x = runif(1, -1, 1), y = runif(1, -1, 1))
  goods <- samps[(samps$x^2+samps$y^2)<1, ]
i = i+1
}

and I also thought about this:

i=0
j=0
samps <- matrix()
goods <- matrix()
needs <- 100

while (j < needs) {
  samps[i,1] <- runif(1, -1, 1)
  samps[i,2] <- runif(1, -1, 1)
  if (( (samps[i,1])**2+(samps[i,2])**2)<1){
  goods[j,1] <- samps[i,1]
  goods[j,2] <- samps[i,2]
}
else{
i = i+1
}
}

but it is not working.

I would be very grateful for any help to modify the code.

like image 638
Sophie Allan Avatar asked Mar 03 '20 00:03

Sophie Allan


People also ask

How to compare two vectors in R?

You can use the following basic syntax to compare two vectors in R: #check if two vectors are identical identical (vector_1, vector_2) #display items that are in both vectors intersect (vector_1, vector_2) #display items that are only in first vector, but not in second vector setdiff (vector_1, vector_2)

What is the vector norms relationship for similarity?

If similarity score means any sensible measure of the length of | x - y |, the vector norms relationship is provide the results. Definition. Two norms ||⋅|| α and ||⋅|| β on a vector space X are called equivalent if there exist positive constants c 1 and c 2 such that c 1 || x || β ≤ | |x || α ≤ c 2 || x || β, ∀ x ∈X.

How to check if two STL vectors contain same elements?

Quickly check if two STL vectors contain same elements or not. Unlike normal C/C++ arrays, we don’t need to do element by element comparison to find if two given vectors contain same elements or not. In case of vectors, the operator “==” is overloaded to find the result quickly. Below is an example to demonstrate same.

Why is the value of false returned when intersecting two vectors?

The two vectors are not identical, so a value of FALSE is returned. The following code shows how to use the intersect () function to display the items that exist in both vectors:


1 Answers

As to your second question ... you cannot reformulate your code to know precisely how many it will take to get (at least) 100 resulting combinations. You can use a while loop and concatenate results until you have at least 100 such rows, and then truncate those over 100. Because using entropy piecewise (at scale) is "expensive", you might prefer to always over-estimate the rows you need and grab all at once.

(Edited to reduce "complexity" based on homework constraints.)

set.seed(42)
samps <- vector(mode = "list")
goods <- vector(mode = "list")
nr <- 0L
iter <- 0L
sampsize <- 100L
needs <- 100L
while (nr < needs && iter < 50) {
  iter <- iter + 1L
  samps[[iter]] <- data.frame(x = runif(sampsize, -1, 1), y = runif(sampsize, -1, 1))
  rows <- (samps[[iter]]$x^2 + samps[[iter]]$y^2) < 1
  goods[[iter]] <- samps[[iter]][rows, ]
  nr <- nr + sum(rows)
}
iter                # number of times we looped
# [1] 2
out <- head(do.call(rbind, goods), n = 100)
NROW(out)
# [1] 100
head(out) ; tail(out)
#            x          y
# 1  0.8296121  0.2524907
# 3 -0.4277209 -0.5668654
# 4  0.6608953 -0.2221099
# 5  0.2834910  0.8849114
# 6  0.0381919  0.9252160
# 7  0.4731766  0.4797106
#               x          y
# 221 -0.65673577 -0.2124462
# 231  0.08606199 -0.7161822
# 251 -0.37263236  0.1296444
# 271 -0.38589120 -0.2831997
# 28  -0.62909284  0.6840144
# 301 -0.50865171  0.5014720
like image 58
r2evans Avatar answered Oct 24 '22 09:10

r2evans