满足任意概率密度函数分布的随机变量生成算法

写在前面:
  又一次看到物理书上电子云示意图时,我就想可不可以由计算机模拟画一幅电子云图呢,它只给出了概率密度函数,我可以生成满足它的随机点吗?

概率密度的定义
  概率简单地说就是在空间某处附近一确定的范围内发生某事件的可能性,而概率密度则是概率与空间范围之比的极限。在讨论更复杂一些的空间之前,可以从一维的线性空间出发定义它们。
  在一维空间中,用数轴表示所有的点,定义Ф(x)为:随机变量i,当i<x时的概率,因此得Ф(-∞)=0,Ф(+∞)=1(因为全概率为1),这个条件是直接由定义得出,称为归一化条件。
  由此当x取数轴上的点时,随机变量i取到[x,x+△x](△x>0)中时的概率应该为Ф(x+△x)-Ф(x),则定义i取到这个区间的平均概率密度p^=(Ф(x+△x)-Ф(x))/ △x,此时如果极限:
  lim(Ф(x+△x)-Ф(x))/ △x=p(x) 存在
  △ x->0
  则称极限值为这一点的概率密度,由数学知识这正是导数的定义,于是有p(x)=dФ/dx,于是有△Ф=∫p(x)dx,积分区域为△,虽然这是从一维情况推出,这是概率的一般形式定义。于是归一化条件就相应为∫p(v)dv=1,积分区域为全空间。
  由此概率密度函数就唯一决定了随机变量的分布,反过来我们也说随机变量服从某种概率密度函数的分布。

如何由p函数模拟随机变量的生成
  首先我们要有一个这样的随机函数,C的原型如下:
double random();/*返回一个大于等于0,小于1之间的等概率随机双精度小数*/
  
-->规格化的随机变量<--
  这个算法的实现将在以后给出算法,事实上笔者最开始学习编程时就是学的这种随机函数,它是BASIC里的标准函数,这里只是为了讨论的需要,因为由它可以非常方便的得到任意范围内的随机数,包括浮点数或整数。
  然后我们再拿来这样一个任意的概率密度函数p(x),且满足归一化条件∫p(x)dx=1,积分区域为[0,1)。最后再给出这个函数在[0,1)上的上界(最大值)h。然后用下面的方法生成服从该分布的随机变量。
  首先得到一个等概率的随机的二维的点(x,y),它落在矩形区域[0,1;0,h]内,同时我们还要求这样一个条件y<p(x),如果得到的点不满足这个条件则重新取点,也就是说把那些不满足此条件的点‘去掉’。最后我说这里的x满足p(x)的分布。证明如下:

证明

  首先由算法过程将得到在y=0,x=0,x=1,y=p(x)四条曲线围成的二维区域内的等概率随机点。
  任取x属于[0,1),则在x的一个邻域内[x,x+△x](取△x,使这个邻域也在[0,1)内),上述随机变量i取到这个邻域的概率记为△Ф(x),由几何性质△Ф(x)=S'/S,其中S'为这一区域围成的曲顶梯形的面积,而S为整个y=p(x)曲线围成的面积,由归一化条件马上知S=1于是△Ф(x)就是那块曲顶梯形的面积,由中值定理,有△Ф(x)=p(x') △x,其中x'为上面那个邻域内的某一点,现令△x->0,则由夹逼定理x'->x,于是得p(x)=dФ(x)/dx,这正是p(x)的定义,于是产生的随机变量即满足p(x)的分布。证毕。

C++语言算法描述
double random_rickone(double (*p)(double),double h)
{
   double x=random(),y=random()*h;
   if(y<p(x))return x
   return random_rickone(p,h);
}

阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页