R:蒙特卡洛模拟和正态分布问题 [英] R: Problem with MonteCarlo Simulation and Normal Distribution

问题描述

我正在尝试解决以下问题:

I am trying to solve the following exercise:

令Z_n为n个标准正常观测值的最大值.估计n应该是什么，以便P(Z_n> 4)= 0.25

我尝试了下面的代码，我知道答案大约是n = 9000，因为它的返回值约为0.25. 我应该更改代码，以便n是输出而不是输入.

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 <- sapply(1:n, function(i){max(rnorm(n=n,0,1))})
length(x1[x1>4])/length(x1)

我该怎么做?

感谢您的帮助！

推荐答案

好吧，您可以选择适当的范围，然后只进行二进制搜索.请记住，结果将取决于样本数量和RNG种子.

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 <- function(n) {
    max(rnorm(n))
}

Sample <- function(N, n) {
    set.seed(312345) # sample same sequence of numbers
    x <- replicate(N, Zn(n))
    sum( x > 4.0 )/N
}

P <- 0.25

BinarySearch <- function(n_start, n_end, N) {
    lo <- n_start
    hi <- n_end

    s_lo <- Sample(N, lo)
    s_hi <- Sample(N, 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 <- (hi+lo) %/% 2
        s_me <- Sample(N, me)
        if (s_me >= P)
            hi <- me
        else
            lo <- me

        cat("hi = ", hi, "lo = ", lo, "S = ", s_me, "\n")
    }
    list(hi, Sample(N, hi-1), Sample(N, hi)) 
}    

q <- BinarySearch(9000, 10000, 100000) # range [9000...10000] with 100K points sampled

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)

结果，我得到了

9089
0.24984
0.25015

看起来很合理.虽然很慢...

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

这篇关于R:蒙特卡洛模拟和正态分布问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！