python/C 生成beta分布的随机数

python/C 生成beta分布的随机数

---

---

前言

想把一个算法用C语言实现，其中涉及到了beta分布取随机数，记录一下结果

---

一、beta分布理论知识

参考博客：
Beta分布及其应用

Beta分布是一个定义在[0,1]区间上的连续概率分布族，它有两个正值参数，称为形状参数，一般用α和β表示。在贝叶斯推断中，Beta分布是Bernoulli、二项分布、负二项分布和几何分布的共轭先验分布。Beta分布的概率密度函数形式如下：

• 定义域：[0,1]
• 参数：α , β 均为正值参数，又称为形状参数

Beta分布的概率密度函数:

其中，$\Gamma (x)$为gamma函数，$B(\alpha ,\beta )$为beta函数，

• Beta分布的均值
  $$\frac{\alpha }{\alpha +\beta }$$

• Beta分布的方差
  $$\frac{\alpha \beta }{(\alpha +\beta {)}^2(\alpha +\beta +1)}$$

• Beta分布的图形
  

从Beta分布的概率密度函数的图形我们可以看出，Beta分布有很多种形状，但都是在0-1区间内，因此Beta分布可以描述各种0-1区间内的形状（事件）。因此，它特别适合为某件事发生或者成功的概率建模。同时，当α=1，β=1的时候，它就是一个均匀分布。
beta分布主要有 α和 β两个参数，这两个参数决定了分布的形状，从上图及其均值和方差的公式可以看出：
1）α/(α+β)也就是均值，其越大，概率密度分布的中心位置越靠近1，依据此概率分布产生的随机数也多说都靠近1，反之则都靠近0。
2）α+β越大，则分布越窄，也就是集中度越高，这样产生的随机数更接近中心位置，从方差公式上也能看出来。

---

二、python 生成服从beta分布的随机数

参考文章

一个满足beta分布的随机变量可以被两个服从Gamma分布的独立随机变量X,Y推导得到：

P.S. python是有相应的库函数可以用来生成beta分布随机数，scipy.stats.beta，下面的代码是通过数值近似从均匀分布随机数给出beta分布随机数的代码

    def do_Beta_Universality_Uniform(self)->np.ndarray:
        """
        Creates an array with samples from the desired beta(a,b) distribution.

        Returns
        -------
        beta_rv : numpy.ndarray
            an array with samples from the desired beta(a,b) distribution.

        """
        # 1) Create "a" and "b" many Expos by drawing n samples
        # from the Uniform distr. and plugging them into F^(1) of the
        #Expo (-log(1-U))
        n = self.number_of_simulations
        X = []
        for i in np.arange(self.a): #number 3 of plan of attack
            uniform_samples = uniform().rvs(size=n) #number 1 and 2 of plan of attack
            X.append((-1*np.log(1-uniform_samples)))
        Y = []
        for i in np.arange(self.b): #number 4 of plan of attack
            uniform_samples = uniform(0, 1).rvs(size=n) #number 1 and 2 of plan of attack
            Y.append(-1*np.log(1-uniform_samples))

        #5 of plan of attack - create Gamma(a,1) from Expos by summing all      
        #the a uniforms samples
        X = np.array(X)
        X = X.sum(axis=0)
        
        #6 of plan of attack - create Gamma(b,1) from Expos by summing all  
        #the b uniform samples
        Y = np.array(Y)
        Y = Y.sum(axis=0)

        #7 of plan of attack -the final Beta r.v. is X/(X+Y)
        beta_rv = X/(X+Y)

        return beta_rv

---

三、C语言生成服从beta分布的随机数

算法是一样的，只不过C语言没有生成均匀分布的随机数的库函数，需要从rand()（生成正态分布随机数）出发：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

/* 
	函数功能：	产生（0,1）区间上均匀分布的随机数
	输入参数说明：
		0		给定区间的下限
		1		给定区间的上线
		seed	长整型指针变量， *seed 为伪随机数的种子
*/
double uniform_data(long int * seed)
{
	double t;
	*seed = 2045.0 * (*seed) + 1;
	*seed = *seed - (*seed / 1048576) * 1048576;
	t = (*seed) / 1048576.0;
	return t;
}

// 生成贝塔分布随机数
double rand_beta_dist(double alpha, double beta, long int* s) {

    // 生成服从 [0, 1] 均匀分布的随机数
    double x=0,y=0;
    double tmp1=0;
    for(int i=0; i<alpha; i++)
    {
        tmp1 = uniform_data(s);
        // printf("%f\n", tmp1);
        x = x + (-log(1-tmp1));
    }
    double tmp2=0;
    for(int i=0; i<beta; i++)
    {
        tmp2 = uniform_data(s);
        // printf("%f\n", tmp2);
        y = y + (-log(1-tmp1));
    }
    // 生成贝塔分布随机数
    return x / (x + y);
}




int main() {
    srand((unsigned int)time(NULL));

    double alpha = 2.0;  // 根据需要更改
    double beta = 200.0;   // 根据需要更改

	long int s;
    s = 1000;
    printf("%f\n", log(10));
    for (int i = 0; i < 10; ++i) {
        double sample = rand_beta_dist(alpha, beta, &s);
        printf("%f\n", sample);
    }

    return 0;
}