关于利用均匀分布随机变量产生任意分布变量的实现

　　首先交代下背景，课题需要产生服从高斯分布的随机变量，这个要求对于python，Matlab而言，也就是一个函数调用的事（其实C++的库里面也有，无奈之前不知道，(⊙o⊙)…），假如不调用，我们自己应当如何实现呢？或者再延伸下，如果我们需要产生任意分布，这下没函数调用了吧，那么我们应该怎么办呢？这就是一个比较有意思的问题了......

　　首先来看看我们有啥可以用的，一般而言我们是可以获取随机数的（其实这个随机数也是伪随机数，可以通过线性同余法来产生，这是后话了），我们可以把这个数看做服从均匀分布的随机变量，那么如何将一个均匀分布的变量映射到另一个分布的变量呢？

　　答案如下：

　　x为均匀分布的变量，y为另一分布的变量，利用其累计函数作为中间关系，来完成映射。那么问题又来了，为什么可以这样呢？如何来理解这一问题呢？

来两张图感性认识下：

  

前者是0到100的均匀分布累计函数，后者是N(50,25)的累计函数，可以发现：

1.二者均为递增函数(可以看做严格单调，尽管不明显)

2.函数值范围为[0,1]

映射的原理就是在由均匀分布变量的值找到其累计函数值，对应到图二中相同函数值所对应变量值，由上面两幅图可以发现，图一中几乎所有值都对应到了图二中的40~60中，这在直观上和正态分布的效果是相符的。

接下来就看看如何来程序实现，

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
//---------------------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;
}

 

posted @   Rainlin  阅读(6158)  评论(0)  编辑  收藏  举报
编辑推荐：
· 对象命名为何需要避免'-er'和'-or'后缀
· SQL Server如何跟踪自动统计信息更新?
· AI与.NET技术实操系列：使用Catalyst进行自然语言处理
· 分享一个我遇到过的“量子力学”级别的BUG。
· Linux系列：如何调试 malloc 的底层源码
阅读排行：
· C# 中比较实用的关键字，基础高频面试题！
· .NET 10 Preview 2 增强了 Blazor 和.NET MAUI
· 为什么AI教师难以实现
· 如何让低于1B参数的小型语言模型实现 100% 的准确率
· AI Agent爆火后，MCP协议为什么如此重要！
点击右上角即可分享
微信分享提示