0

就像在大多数物理问题中一样，我的情况也有边界，因此我想根据截断的高斯分布生成（使用 R）随机数。

这个想法是这些数字的平均值不应该取决于边界。我已经找到了 truncnorm 包，但它没有完成这项工作：

例如，这里是平均值为 0.1 和宽度为 0.1 的高斯的情况，但限制在 0 和 1 之间：

install.packages("truncnorm")
library(truncnorm)
vec=rtruncnorm(n=100000,a=0,b=1,mean=0.1,sd=0.1)
hist(vec,breaks=100)
mean(vec)
[1] 0.1289061

如您所见，最终平均值不是作为输入给出的平均值，通过使用标准 rnorm 函数并对结果进行子集化，我可以获得相同的结果。

我不想重新发明轮子，所以欢迎任何关于进一步包装的想法或建议！谢谢！

4

2 回答 2

5

当您截断分布时，这不是您所期望的吗？

#Example
x <- rnorm( 1e7 , mean = 0.1 , sd = 0.1 )
mean( x[ ! ( x < 0 | x > 1 ) ] )
#[1] 0.128814

#Visualising
hist( x , breaks = 100 , xlim = c(-1,1) )
#limits (red)
abline( v = 0 , col = "red" , lwd = 1 , lty = 2 )
abline( v = 1 , col = "red" , lwd = 1 , lty = 2 )
#truncated mean (green)
abline( v = mean( x[ !(x<0|x>1)] ) , col = "green" , lty = 2 , lwd = 1 )
#true mean (blue)
abline( v = 0.1 , col = "blue" , lty = 1 , lwd = 1 )

在此处输入图像描述

于 2014-04-24T09:28:53.777 回答
1

因此，我们可能必须区分截断之前和之后的平均值，并且您显然打算控制截断样本可能会收敛到的可观察平均值，尽管rnorm()（并且可能rtruncnorm()，我不知道）期望“之前”-means；虽然 stats.stackexchange.com 上的一些统计学家可能会为您提供更无懈可击的分析解决方案，但也许一些有趣的优化也可以帮助您找到合适的“之前”参数（您可能需要根据“之前”是否调整此代码 - sd 参数也应该修改）：

myrtruncnorm <- function(n,a,b,mean=0,sd=1) 
    qnorm(runif(n,pnorm(a,mean=mean,sd=sd),pnorm(b,mean=mean,sd=sd)),mean=mean,sd=sd)
set.seed(1)
optim(list(mean=.1,sd=.1), function(x)
    abs(mean(myrtruncnorm(n=100000,a=0,b=1,mean=x[[1]],sd=x[[2]]))-.1))
# returns mean=0.07785390 and sd=0.07777597, let's test that: 
x1 <- myrtruncnorm(100000,0,1,0.07785390,0.07777597)
hist(x1); mean(x1) # Is "mean=0.1003832" sufficiently close?
于 2014-04-24T10:18:23.807 回答