關於利用均勻分布隨機變量產生任意分布變量的實現

　　首先交代下背景，課題需要產生服從高斯分布的隨機變量，這個要求對於python，Matlab而言，也就是一個函數調用的事（其實C++的庫里面也有，無奈之前不知道，(⊙o⊙)…），假如不調用，我們自己應當如何實現呢？或者再延伸下，如果我們需要產生任意分布，這下沒函數調用了吧，那么我們應該怎么辦呢？這就是一個比較有意思的問題了......

　　首先來看看我們有啥可以用的，一般而言我們是可以獲取隨機數的（其實這個隨機數也是偽隨機數，可以通過線性同余法來產生，這是后話了），我們可以把這個數看做服從均勻分布的隨機變量，那么如何將一個均勻分布的變量映射到另一個分布的變量呢？

　　答案如下：

　　x為均勻分布的變量，y為另一分布的變量，利用其累計函數作為中間關系，來完成映射。那么問題又來了，為什么可以這樣呢？如何來理解這一問題呢？

來兩張圖感性認識下：

  

前者是0到100的均勻分布累計函數，后者是N(50,25)的累計函數，可以發現：

1.二者均為遞增函數(可以看做嚴格單調，盡管不明顯)

2.函數值范圍為[0,1]

映射的原理就是在由均勻分布變量的值找到其累計函數值，對應到圖二中相同函數值所對應變量值，由上面兩幅圖可以發現，圖一中幾乎所有值都對應到了圖二中的40~60中，這在直觀上和正態分布的效果是相符的。

接下來就看看如何來程序實現，

//---------------------GaussianDistribution.h------------------------
#pragma once
//#include "Interface.h"
#define SP 1000//設置精確度，取整個正態分布的區間分為SP等份
namespace RAN
{
	class GaussianDistribution 	{
	public:
		GaussianDistribution(void);
		~GaussianDistribution(void);
		double u,d,step;
		double mp[SP],mf[SP];//mp為正態分布區間每等份所對應的概率密度，mf為累計累計概率密度
    

		void SetPara(double u,double d);
		double GetOne();
};
}


//---------------------GaussianDistribution.cpp------------------------
#include "GaussianDistribution.h"
#include <cmath>
#define   PI 3.1415926
#include <time.h>
GaussianDistribution::GaussianDistribution(void)
{
}


GaussianDistribution::~GaussianDistribution(void)
{
}

template<typename T>
int BinarySearch(T *array,T key)//二分查找，返回值所在數組中的位置，若無則返回其左邊值
{  
	int aSize=SP;
	if ( array == NULL || aSize == 0 )  
		return -1;  
	int low = 0;  
	int high = aSize - 1;  
	int mid = 0;  

	while ( low <= high )  
	{  
		mid = (low + high )/2;  

		if ( array[mid] < key)  
			low = mid + 1;  
		else if ( array[mid] > key )     
			high = mid - 1;  
		else  
			return mid;  
	}  
	return low;
} 

void GaussianDistribution::SetPara(double u,double d)
{
	this->u=u;
	this->d=d;
	step=10*d/SP;//取正態分布整個區間中以u為中心，5d為半徑的區間
	mf[0]=0;
	mp[0]=0;
	for (int i=1;i<SP;i++)
	{
		double x=u+(step*(i-SP/2));//計算pdf中的x值
		mp[i]=1/(sqrt(2*PI)*d)*exp(-(x-u)*(x-u)/(2*d*d));//x對應概率密度
		mf[i]=mf[i-1]+mp[i]*step;//累計概率密度
	}
	srand(time(0));
}
double GaussianDistribution::GetOne()
{
	//srand(time(0)); //use current time as seed for random generator
	int random= rand();
	double p= (random+0.0)/RAND_MAX;//均勻分布累計概率密度
	int k=BinarySearch(mf,p);//找到正態分布的對應
	double gg=u+(k-SP/2)*step;
	return gg;
}