【问题标题】:Collect random dots under normal distribution using R使用 R 收集正态分布下的随机点
【发布时间】:2023-04-03 23:09:01
【问题描述】:

我正在尝试通过首先绘制标准正态分布曲线然后在图上均匀生成点来编写 R 脚本。我想收集标准普通 PDF 下的所有点数。但是当我计算金额时,它与图上显示的不一样。

这是我的代码,

N = 20
normpdf = function(x) {
    f = (1/sqrt(2*pi))*exp(x^2/(-2))
    return(f)
}        

x = seq(-5,5,length=N)
normprobdist = normpdf(x)    
plot(x,normprobdist,type="l")
par(new=TRUE) 
u = NULL
dots_undercurve = NULL
for (i in 1:length(x)) { 
     u[i] = runif(x,0,1)
     if (u[i]<=normprobdist[i]) {
                dots_undercurve[i] = TRUE
     } 
     else {
              dots_undercurve[i] = FALSE
     }

}

sum(dots_undercurve)
plot(x,u,xaxt='n', yaxt = 'n', ann=FALSE)

sum(dots_undercurve) 显示曲线下的点数与图上显示的不同。我的编程知识很差,所以如果有人能找到我的脚本有什么问题吗? 谢谢

【问题讨论】:

  • 为什么将x 传递给runif?它期望绘制的随机值的数量作为第一个参数。
  • 我的错,我当然应该通过1

标签: r


【解决方案1】:

您不需要 for 循环或自己定义 pdf:

N <- 1e6
set.seed(42) #for reproducibility of random numbers
x <- runif(N, -5, 5)
y <- runif(N, 0, 1)

yunder <- y < dnorm(x) #which dots are under curve
sum(yunder)/N * (10 * 1) #the expected value is about 1
#[1] 0.99914

curve(dnorm, -5, 5, ylim=c(0,1))
points(x, y, col=c("red", "green")[yunder+1], pch=".")

【讨论】:

  • 嘿,谢谢!但在我的任务中,我不允许使用内置分发功能 dnorm :)
  • 你也可以使用x &lt;- seq(-5, 5, length.out=N)。没关系。
  • 不允许你使用dnorm 是愚蠢的。接下来,他们禁止你使用exp
  • 嗯,这有点意思,他们希望我们学习一些编程。所以我想要做的是计算说 P(R
  • 好吧,您可以用您的normpdf 替换我的代码中的dnorm。他们应该尝试教你 R 编程,这意味着尽可能使用矢量化而不是 for 循环。
【解决方案2】:

这里有很多错误: 我认为你想要的是这样的:

rm(list=ls(all=TRUE))
options(warn=-1)
N = 20
normpdf = function(x) {
  f = (1/sqrt(2*pi))*exp(x^2/(-2))
  return(f)
}
x = seq(-5,5,length=N)
normprobdist = normpdf(x)    
par(new=TRUE) 
u = runif(N,0,1)
dots_undercurve = NULL
for (i in 1:length(x)) { 
  if (u[i]<=normprobdist[i]) {
  dots_undercurve[i] = TRUE
} 
else {
  dots_undercurve[i] = FALSE
} 
}
sum(dots_undercurve)
plot(x,u)
lines(x,normprobdist,type="l")

【讨论】:

  • 我想指出,这段代码效率很低,不能很好地适应大型N
  • 问题是 sum(dots_undercurve) 给出的点数与图中所见的点数不同
猜你喜欢
  • 1970-01-01
  • 2014-02-15
  • 2015-03-05
  • 2011-03-17
  • 2021-09-03
  • 1970-01-01
  • 2017-07-24
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多