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.ndarrayan 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_simulationsX = []for i in np.arange(self.a): #number 3 of plan of attackuniform_samples = uniform().rvs(size=n) #number 1 and 2 of plan of attackX.append((-1*np.log(1-uniform_samples)))Y = []for i in np.arange(self.b): #number 4 of plan of attackuniform_samples = uniform(0, 1).rvs(size=n) #number 1 and 2 of plan of attackY.append(-1*np.log(1-uniform_samples))#5 of plan of attack - create Gamma(a,1) from Expos by summing all      #the a uniforms samplesX = 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 samplesY = 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;
}