蓝色艺林

导航

以均匀随机分布生成任意斜率随机线性分布

 本文的目的是利用(0,1)上的均匀分布随机数生成器生成区间为(imin,imax),斜率为slope的任意归一化线性随机数生成器。

借助(0,1)上的均匀随机数生成器,可以通过反函数法生成任意分布的随机数生成器。

对于C++,生成(imin,imax)上随机数生成器的代码为

double uniform_dist(double imin,double imax)
{
    int temp;
    while((temp=rand())==RAND_MAX){
        ;
    }
    return temp/RAND_MAX*(imax-imin)+imin);
}

其中imin为随机变量最小值,imax为随机变量最大值。

在调用时,只需

x=uniform_dist(0,1);

就可以生成$(0,1)$上的任意随机变量。

 

对任累积意分布函数$F$,其概率密度为$f(x)$,则有

 

$F(x_1)=P(x\le x_1)=\int _{-\infty} ^{x_1} f(x') dx$

由累积分布函数的值范围$F \in (0,1),则若$F$为$(0,1)$上的随机变量,即可通过$F(x_1)$的反函数来求得与此对应的随机变量$x_1$。

 

因此要生成概率密度函数为$f(x)=slope*x+b$区间为$(imin,imax)$的随机变量,分为2步

1. 通过归一化求出b的值

$ \int _{-\infty} ^{\infty} f(x)dx=\int _{imin} ^{imax} f(x)dx=1$

可以求出

$b=\frac {1-\frac {1}{2} slope (imax^2-imin^2)}{imax-imin}$

2. 由累积分布函数

$F(x_1)= \int _{imin} ^{x_1} f(x)dx=F_1(x_1)-F_1(imin)$

其中$F_1$为概率密度函数$f(x)$的原函数。由$F(x_1)$的反函数$F^{-1} (x_1)$可以求得的随机变量$x_1 \in (imin,imax)$,并且服从概率密度函数$f(x)$。

带入方程得

$\frac {1}{2}a x_1^2+\frac {1-\frac {1}{2} slope (imax^2-imin^2)}{imax-imin}x_1 \\ -\frac {1}{2}\ slope\ imin^2-\frac {1-\frac {1}{2} slope (imax^2-imin^2)}{imax-imin} imin - F(x_1)=0$

$ a=\frac {1}{2}slope$

$ b=\frac {1-\frac {1}{2}\ slope \ (imax^2-imin^2)}{imax-imin}$

$ c=-\frac {1}{2}\ slope\ imin^2-\frac {1-\frac {1}{2} slope (imax^2-imin^2)}{imax-imin} imin - F(x_1)$

解得

$ x_1=\frac {-b + \sqrt {b^2-4ac}}{2a}$  

 

C++实现程序如下

#include<iostream>
#include<vector>
#include<time.h>
#include<math.h>
#include<stdlib.h>
using namespace std;

//*******************************************************************
//  uniform distribution
//*******************************************************************

double uniform_dist(double imin,double imax) { int temp; while((temp=rand())==RAND_MAX){ ; } return((double) temp/RAND_MAX*(imax-imin)+imin); }

//*******************************************************************
//  linear distribution for any slope
//*******************************************************************
double linear_dist(double a,double b,double imin,double imax,double slope) { double c; c=-0.5*slope*imin*imin-b*imin-uniform_dist(0,1); return ((-b+sqrt(b*b-4*a*c))/(2a)); }
int main() {
  //===========================================================
  //  sequence to be sampled
  //===========================================================
const N=1000000; vector< double> x(N); srand((unsigned)time(NULL));

  //===========================================================
  //  parameters of the linear distribution
  //===========================================================
double imin,imax,slope; imin=0.0; imax=3.0; slope=-0.2; double a,b,criterion; a=0.5*slope; b=(1-0.5*slope*(imax*imax-imin*imin))/(imax-imin);

  //*********************************************************
  //  judge the reasonability of the linear distribution
  //*********************************************************
if(slope>0){ criterion=slope*imin+b; if(criterion<0){ cerr<<"The probability of the variable must be larger than 0"<<endl; exit(1); } } else if(slope<0){ criterion=slope*imax+b; if(criterion<0){ cerr<<"The probability of the variable must be larger than 0"<<endl; exit(1); } }
  //**********************************************************
  //  sampling the random data
  //**********************************************************
for(int i=0;i!=x.size();++i){ x[i]=linear_dist(a,b,imin,imax,slope); }
  //**********************************************************
  //  output the data
  //**********************************************************
for(int i=0;i!=x.size();++i){ cout<<x[i]<<""; } return 0;

注意:

对斜率的选取和随机变量区间的选取需要考虑归一化的限制,也即是需要考虑每一个随机变量的概率要大于0。

 

posted on 2017-06-04 20:05  蓝色艺林  阅读(580)  评论(0编辑  收藏  举报