估计python中均匀随机变量之和的概率密度

Posted

技术标签:

【中文标题】估计python中均匀随机变量之和的概率密度【英文标题】:Estimating the probability density of sum of uniform random variables in python 【发布时间】:2016-05-25 11:32:55 【问题描述】:

我有两个随机变量 X 和 Y,它们均匀分布在单纯形上:

我想评估它们总和的密度:

在评估上述积分之后,我的最终目标是计算以下积分:

为了计算第一个积分,我在单纯形中生成均匀分布的点,然后检查它们是否属于上述积分中的所需区域,并取点的分数来评估上述密度。

一旦我计算出上述密度,我将按照类似的程序来计算上述对数积分来计算其值。然而,这非常低效并且需要花费很多时间,例如 3-4 个小时。谁能建议我用 Python 解决这个问题的有效方法?我正在使用 Numpy 包。

这里是代码

import numpy as np
import math
import random
import numpy.random as nprnd
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
#This function checks if the point x lies the simplex and the negative simplex shifted by z
def InreqSumSimplex(x,z):
    dim=len(x)
    testShiftSimpl= all(z[i]-1 <= x[i] <= z[i] for i in range(0,dim)) and (sum(x) >= sum(z)-1)
    return int(testShiftSimpl)

def InreqDiffSimplex(x,z):
    dim=len(x)
    testShiftSimpl= all(z[i] <= x[i] <= z[i]+1 for i in range(0,dim)) and (sum(x) <= sum(z)+1)
    return int(testShiftSimpl)
#This is for the density X+Y
def DensityEvalSum(z,UniformCube):
    dim=len(z)
    Sum=0
    for gen in UniformCube:
        Exponential=[-math.log(i) for i in gen] #This is exponentially distributed
        x=[i/sum(Exponential) for i in Exponential[0:dim]] #x is now uniformly distributed on simplex

        Sum+=InreqSumSimplex(x,z)

    Sum=Sum/numsample

    FunVal=(math.factorial(dim))*Sum;
    if FunVal<0.00001:
        return 0.0
    else:
        return -math.log(FunVal)
#This is for the density X-Y
def DensityEvalDiff(z,UniformCube):
    dim=len(z)
    Sum=0
    for gen in UniformCube:
        Exponential=[-math.log(i) for i in gen]
        x=[i/sum(Exponential) for i in Exponential[0:dim]]

    Sum+=InreqDiffSimplex(x,z)

    Sum=Sum/numsample

    FunVal=(math.factorial(dim))*Sum;
    if FunVal<0.00001:
        return 0.0
    else:
        return -math.log(FunVal)
def EntropyRatio(dim):    
    UniformCube1=np.random.random((numsample,dim+1)); 
    UniformCube2=np.random.random((numsample,dim+1))

    IntegralSum=0; IntegralDiff=0

    for gen1,gen2 in zip(UniformCube1,UniformCube2):

        Expo1=[-math.log(i) for i in gen1];        Expo2=[-math.log(i) for i in gen2]

        Sumz=[ (i/sum(Expo1)) + j/sum(Expo2) for i,j in zip(Expo1[0:dim],Expo2[0:dim])] #Sumz is now disbtributed as X+Y

        Diffz=[ (i/sum(Expo1)) - j/sum(Expo2) for i,j in zip(Expo1[0:dim],Expo2[0:dim])] #Diffz is now distributed as X-Y

    UniformCube=np.random.random((numsample,dim+1))

    IntegralSum+=DensityEvalSum(Sumz,UniformCube) ; IntegralDiff+=DensityEvalDiff(Diffz,UniformCube)

    IntegralSum= IntegralSum/numsample; IntegralDiff=IntegralDiff/numsample

    return ( (IntegralDiff +math.log(math.factorial(dim)))/ ((IntegralSum +math.log(math.factorial(dim)))) )

Maxdim=11
dimlist=range(2,Maxdim)
Ratio=len(dimlist)*[0]
numsample=10000

for i in range(len(dimlist)):
    Ratio[i]=EntropyRatio(dimlist[i])

【问题讨论】:

你能告诉你当前的代码吗? 您对n 的哪种值感兴趣? @MarkDickinson:我实际上对更高的 n 值感兴趣,比如高达 100,200 等。但我需要绘制从 n=2 到 200 的所有值。这就是我想让它高效的原因. @MaxNoe:大约 100 行 Python 代码。如何上传代码? 您是否分析了代码?实际上需要这么长时间?您可以为此使用 profilehooks 模块。 【参考方案1】:

不确定这是否是您问题的答案,但让我们开始

首先,这里是一些代码示例和讨论如何通过gammavariate() 或通过-log(U) 正确地从 Dirichlet(n)(又名单工)采样,但对潜在的极端情况有适当的处理,link

我所看到的代码的问题是,例如,对于采样维度 = 2 单纯形 你得到三个(!)统一数字,但在对x 进行列表理解时跳过了一个。这是错误的。要对 n 维 Dirichlet 进行采样,您应该得到准确的 n U(0,1) 并进行变换(或 n 来自 gammavariate 的样本)。

但是,最好的解决方案可能是使用 numpy.random.dirichlet(),它是用 C 编写的,可能是最快的,请参阅 link。

最后一个,以我的拙见,你没有正确估计log(PDF(X+Z))。好的,你发现有些是,但此时 PDF(X+Z) 是什么?

这样做

testShiftSimpl= all(z[i]-1 <= x[i] <= z[i] for i in range(0,dim)) and (sum(x) >= sum(z)-1)
return int(testShiftSimpl)

看起来像 PDF?你是怎么得到它的?

简单测试:在整个X+Z 区域上整合PDF(X+Z)。它产生了 1 吗?

更新

看起来我们可能有不同的想法,我们称之为单纯形,Dirichlet 等。我非常赞同 this definition,在 d 暗空间中,我们有 d 点和 d-1 单纯形是凸包连接顶点。单纯形维数总是 由于坐标之间的关系,比空间小一。在最简单的情况下,d=2,1-simplex 是连接点 (1,0) 和 (0,1) 的线段,从 Dirichlet 分布中我得到了图片

d=3 和 2-单纯形的情况下,我们有三角形连接点 (1,0,0)、(0,1,0) 和 (0,0,1)

代码,Python

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

import math
import random

def simplex_sampling(d):
    """
    Sample one d-dim point from Dirichet distribution
    """
    r = []
    sum = 0.0

    for k in range(0, d):
        x = random.random()
        if x == 0.0:
            return make_corner_sample(d, k)

        t = -math.log(x)
        r.append(t)
        sum += t

    norm = 1.0 / sum

    for k in range(0, d):
        r[k] *= norm

    return r

def make_corner_sample(d, k):
    """
    U(0,1) number k is zero, it is a corner point in simplex
    """
    r = []
    for i in range(0, d):
        if i == k:
            r.append(1.0)
        else:
            r.append(0.0)

    return r

N = 500 # numer of points to plot
d = 3   # dimension of the space, 2 or 3

x = []
y = []
z = []

for k in range(0, N):
    pt = simplex_sampling(d)

    x.append(pt[0])
    y.append(pt[1])
    if d > 2:
        z.append(pt[2])

if d == 2:
    plt.scatter(x, y, alpha=0.1)
else:
    fig = plt.figure()
    ax  = fig.add_subplot(111, projection='3d')
    ax.scatter(x, y, z, alpha=0.1)

    ax.set_xlabel('X Label')
    ax.set_ylabel('Y Label')
    ax.set_zlabel('Z Label')

plt.show()

【讨论】:

上述条件确保 z-x 位于单纯形区域,这是我们进行密度评估所需的区域。所以我正在计算满足上述条件的单纯形中点的分数,这是对 pdf 的估计。 对于单纯形内点的生成,我没有使用您指出的 Dirichlet 分布过程。但我的程序是,如果 U1,...,U_n+1 以速率 1 呈指数分布,则 (U1/U_1+..U_n+1,....., U_n/U_1+....+U_n+1 ) 在单纯形上是均匀的。这就是我在列表理解过程中跳过一个的原因。

以上是关于估计python中均匀随机变量之和的概率密度的主要内容,如果未能解决你的问题,请参考以下文章

人工智能数学基础--概率与统计14:连续随机变量的指数分布威布尔分布和均匀分布

人工智能数学基础--概率与统计14:连续随机变量的指数分布威布尔分布和均匀分布

已知分布函数如下,求概率密度,请写出具体步骤

概率论与数理统计

概率论与数理统计 Chapter2. 随机变量及概率分布

机器学习中的概率模型和概率密度估计方法及VAE生成式模型详解之四(第2章)