R:蒙特卡洛模拟和正态分布的问题
R: Problem with MonteCarlo Simulation and Normal Distribution
我正在尝试解决以下练习:
设 Z_n 为 n 个标准正态观测值的最大值。估计 n 应该是多少 P(Z_n>4)=0.25
我试过下面的代码,我知道答案大约是 n=9000,因为它 returns 大约是 0.25。
我应该更改我的代码,使 n 成为输出而不是输入。
n=9000
x1 <- sapply(1:n, function(i){max(rnorm(n=n,0,1))})
length(x1[x1>4])/length(x1)
我该怎么做?
感谢您的帮助!
好吧,您可以 select 适当的范围,然后只进行二进制搜索。请记住,结果将取决于样本数量和 RNG 种子。
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
这看起来很合理。虽然很慢...
我正在尝试解决以下练习:
设 Z_n 为 n 个标准正态观测值的最大值。估计 n 应该是多少 P(Z_n>4)=0.25
我试过下面的代码,我知道答案大约是 n=9000,因为它 returns 大约是 0.25。 我应该更改我的代码,使 n 成为输出而不是输入。
n=9000
x1 <- sapply(1:n, function(i){max(rnorm(n=n,0,1))})
length(x1[x1>4])/length(x1)
我该怎么做?
感谢您的帮助!
好吧,您可以 select 适当的范围,然后只进行二进制搜索。请记住,结果将取决于样本数量和 RNG 种子。
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
这看起来很合理。虽然很慢...