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

发布时间:2023年12月22日

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



前言

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


一、beta分布理论知识

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

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

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

Beta分布的概率密度函数:
在这里插入图片描述
其中, Γ ( x ) \Gamma(x) Γ(x)gamma函数 B ( α , β ) B(\alpha,\beta) B(α,β)beta函数
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

  • Beta分布的均值
    α α + β \frac{\alpha}{\alpha+\beta} α+βα?

  • Beta分布的方差
    α β ( α + β ) 2 ( α + β + 1 ) \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)} (α+β)2(α+β+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;
}


文章来源:https://blog.csdn.net/weixin_62355917/article/details/135152752
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。