蒙特卡洛正态分布java,R：蒙特卡洛模拟和正态分布问题

I am trying to solve the following exercise:

Let Z_n be maximum of n standard normal observations. Estimate what n should be so that P(Z_n>4)=0.25

I have tried following code and I know the answer is about n=9000 because it returns aproximately 0.25.

I should change my code so that n is the output and not the input.

n=9000

x1

length(x1[x1>4])/length(x1)

How can I do that?

Thanks for helping!

解决方案

Well, you could select appropriate range and then just do binary search. Just remember, result will depend on number of samples and RNG seed.

Zn

max(rnorm(n))

}

Sample

set.seed(312345) # sample same sequence of numbers

x

sum( x > 4.0 )/N

}

P

BinarySearch

lo

hi

s_lo

s_hi

if (s_lo > P)

return(list(-1, 0.0, 0.0)) # wrong low end of interval

if (s_hi < P)

return(list(-2, 0.0, 0.0)) # wrong high end of interval

while (hi-lo > 1) {

me

s_me

if (s_me >= P)

hi

else

lo

cat("hi = ", hi, "lo = ", lo, "S = ", s_me, "\n")

}

list(hi, Sample(N, hi-1), Sample(N, hi))

}

q

print(q[1]) # n at which we have P(Zn(n)>4)>=0.25

print(q[2]) # P(Zn(n-1)>4)

print(q[3]) # P(Zn(n)>4)

As a result, I've got

9089

0.24984

0.25015

which looks reasonable. It is quite slow though...