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

 本文的目的是利用（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。