In: Statistics and Probability
I want to know why I can not get N(should contain 1000 number)
N<-vector()
for (k in 0:999){
y=vector()
for(i in 0:999){
r1=runif(1)
r2=runif(1)
x1=-log(r1)
x2=-log(r2)
k=(x1-1)^2/2
if (x2 >= k){
r<-runif(1)
if (r>0.5){
y[i]= x1
}else{
y[i]= -x1
}
}
y.1<- y[-which(is.na(y))]
number<-length(y.1)
}
N[k]<-number
}
N
There is certain problem in your code.
1. You have defined k in the for loop than again you defined the k=(x1-1)^2/2 in the program line. so this chage the value of k. First k should take the value 0 to 999 but here after running the k=(x1-1)^2/2 value changes. So you should assign different vactor. I have choosen K instead of k.
2. You have defined the array from 0. this cause problem since y[0] not defined. So either you start your loop by 1 to 1000 or you just change y[i] to y[ i+1] and N[k] to N[k+1].
3. You should close the second loop before line y.1<- y[-which(is.na(y))] . because till here we have the vector y and our work has done .
The correct code should be
N<-vector()
for (k in 0:999){
y=vector()
for(i in 0:999){
r1=runif(1)
r2=runif(1)
x1=-log(r1)
x2=-log(r2)
K=(x1-1)^2/2 ### Just change the k to K
if (x2 >= K){
r<-runif(1)
if (r>0.5){
y[i+1]= x1
}else{
y[i+1]= -x1
}
}
}
y.1<- y[-which(is.na(y))]
number<-length(y.1)
number
N[k+1]<-number
}
N