如何生成截尾正态分布随机数（思路）

如题，现在需要生成介于一定区间的正态随机数，但是不知道如何下手。

比如，生成介于10到20之间的正态随机数，要求均值为15，标准差为5
之所以有这样的需求，是因为如果单单生成均值为15，标准差为5的正态随机数，随机值当中会出现负值的情况，但是这跟现实不符，比如年龄是不可能出现负值的，尽管从统计上来讲是可能的
所以就需要截尾处理，使正态分布介于一定的区间以内

在VBA里可以用NORMINV函数生成正态随机数，但是不知道如何限定区间

曾经试过一些方法，比如先通过循环逐一生成均值为15，标准差为5的正态随机数，函数用NORMINV(RAND(),15,5)
然后循环内判断是否介于10-20之间，如果不在这个区间内，重新生成正态随机数
以此类推，但是结果是可以保证正态和区间，不能保证均值为10，标准差为5

也试过线性变化的方法来将正态随机数转换为一定区间内的分布，但也不能保证均值和标准差

目前找到一点R语言写的算法（可能是红色部分，也可能是全部），从来没用过R，高手能否解释一下这个代码的大致思路，我对VBA比较熟悉，感谢
#R program for dichotomous probit model
#read the data
x=as.matrix(read.table("c:\\bookdead.dat")[,3:9])
y=as.matrix(read.table("c:\\bookdead.dat")[,10])
#create variables, set values, and write out starting values
b=matrix(0,7); vb=solve(t(x)%*%x); ch<-chol(vb)
write(c(i,t(b)), file="c:\\dprob_gibbs.out", append=T, ncol=8)
#begin MCMC simulation
for(i in 2:5000){
#simulate latent data from truncated normal distributions
u=as.matrix(runif(length(y),min=0,max=1))
xb=as.matrix(x%*%b)
ystar=qnorm(y*u + u*(-1)^y*pnorm(0,mean=xb,sd=x[,1]) +y*pnorm(0,mean=xb,sd=1), mean=xb, sd=x[,1])
#simulate beta vector from appropriate mvn distribution
b<-vb%*%(t(x)%*%ystar) + t((rnorm(7,mean=0,sd=1))%*%ch)
write(c(i,t(b))), file="c:\\dprob_gibbs.out", append=T, ncolumns=8)
if(i%%10==0){print(c(i,t(b))),digits=2)}

?rnorm

这个方差肯定要比原来的小
可以有公式计算出方差，但是一般都比比原来的小
想要得到和原来一样方差那是不可能的， 除非区间相对于原来的方差大得多

不求原委的话，可以用bayesm包中的rtrun函数。
比较基本的想法是利用accept-reject algorithm, candidate function 可以取为一般的正态分布N(mu, sigma^2).

hugebear 发表于 2013-5-11 11:57
不求原委的话，可以用bayesm包中的rtrun函数。
比较基本的想法是利用accept-reject algorithm, candidate  ...

刚才在一本英文教材中找到了R语言的截尾正态分布随机数的生成代码，不过我不是很懂啥意思，请问您能否解释一下这个算法的大致原理，详见帖子内容，我已经重新编辑，并贴上R代码了

ntsean 发表于 2013-5-11 09:32
这个方差肯定要比原来的小
可以有公式计算出方差，但是一般都比比原来的小
想要得到和原来一样方差那是不 ...

谢谢，我贴了一个R代码，不过不是很了解，如果您熟悉的话，能否解释一下，感谢

你贴的程序是用Gibbs sampler做的，更简洁的代码如下(参见Robert, Casella: Introducing Monte Carlo Methods with R 7.4节）
rtnorm=function(n=1,mu=0,lo=-Inf,up=Inf，sigma = 1){
qnorm(runif(n,min=pnorm(lo,mean=mu,sd=sigma),
max=pnorm(up,mean=mu,sd=sigma)),
mean=mu,sd=sigma)}

hugebear 发表于 2013-5-11 13:35
你贴的程序是用Gibbs sampler做的，更简洁的代码如下(参见Robert, Casella: Introducing Monte Carlo Metho ...

您所说的我看不大懂，从来没有接触过R，我自己也在VBA里试了以下代码，主要思路是通过限制概率分布来达到限制边界的目的，比如根据公式(x-mean)/SD将概率分布限定在0.025-0.975之间，从而达到限定正态随机分布边界的目的。但是这个方法的缺陷是，当设定标准差较大的时候，生成的正态随机数的标准差会小于设定标准差。不知道您的R代码是不是这个意思，还是有更优的算法，请大侠指教。
Sub Test()
Randomize
Dim i&
Dim m, n As Single
    Max=15
    Min=5
    Mean=10
    SD=5
    m = 2 * Application.NormSDist((Max - Mean) / SD) - 1
    n = 1 - Application.NormSDist((Mean - Min) / SD)
For i = 1 To 1000
    Cells(i + 1, 1) = Application.NormInv(Rnd() * m + n, 10, 5)
Next
EndSub
当Mean=10，SD=1，边界为[5,15]的时候，随机数均值为10.01，标准差为1.00
当Mean=10，SD=2，边界为[5,15]的时候，随机数均值为9.99，标准差为1.90
当Mean=10，SD=3，边界为[5,15]的时候，随机数均值为9.98，标准差为2.38
当Mean=10，SD=4，边界为[5,15]的时候，随机数均值为9.98，标准差为2.59
当Mean=10，SD=5，边界为[5,15]的时候，随机数均值为9.99，标准差为2.71

思路是：
1. 在(a,b)上截取某个正态N(u, s^2)的去掉两个尾巴的分布（截尾正态），有四个参数a, b, u, s. 如果这个分布的MEAN=m(a,b,u,s), VAR=v(a,b,u,s),那么你的问题就变成已知MEAN（你题中是15）和VAR（你题中是5^2=25）求u,s(因为a,b已知，即你题中的10,20）,其实这就是个两个未知数的两个方程组成的方程组。即：求解15=m(10,20,u,s), 25=v(10,20,u,s)
2.当然你会问，m(),v()的函数形式不知，如何解方程组？很简单，从WIKIPEDIA，http://en.wikipedia.org/wiki/Truncated_normal_distribution 可以找到这两个函数的形式，太长，我不COPY了，那文中"MOMENTS"的地方写出前两个公式就是要找的函数。
3.不幸，你举的例子，方程是无解的，因为你的VAR（25）太大了，在10-20区间上MEAN=15的截尾正态最大VAR是8.3333333（那是10-20之间均匀分布的VAR，在此上的截尾正态之VAR不会超过它），所以我把你的例子改成MEAN=15，VAR=5（即标准差是SQRT（5）），这样方程组solution: u=15, s=2.5825
4. 后面就很容易了，算法如下：要得到(10,20)上截取某个正态N(15, 2.5825^2)的去掉两个尾巴的分布（截尾正态），用EXCEL写就是在：
A1格输入=NORMINV(RAND(),15,2.5825), B1输入=if(and(a1<20,a1>10),a1,-1)，拉A1，B1的公式充满A列和B列，B列中不是-1的数，就是你要求的截尾正态随机数（当然VAR=5，不是你题中要求的25，因为VAR=25是做不到的）。

TaskShare 发表于 2013-7-6 15:04
思路是：
1. 在(a,b)上截取某个正态N(u, s^2)的去掉两个尾巴的分布（截尾正态），有四个参数a, b, u, s. 如 ...

谢谢你的回复，你说的思路我之前也有尝试，之前感觉可能还存在更优的算法。后来我是直接限制[0,1]均匀随机数的概率大小，进而限制其取值范围