在Python中估算均匀随机变量之和的概率密度

Question

我有两个随机变量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])

- pikachuchameleon

---

@MarkDickinson：我实际上对n的更高值感兴趣，比如100、200等。但我需要绘制从n=2到200的所有值的图表。这就是为什么我想让它更有效率的原因。 - pikachuchameleon

你对代码进行了分析吗？实际上是什么导致了这么长的时间？你可以使用 profilehooks 模块来进行分析。 - MaxNoe

您可以在问题中发布代码，它不会变得太大。 - MaxNoe

@MaxNoe：我已经添加了代码。如果需要进一步的注释，请告诉我。 - pikachuchameleon

抱歉，我相信它。我的错。不错的方法。 - Mark Dickinson

Answer 1

“不确定这是否是你问题的答案，但让我们开始吧。首先，这里有一些代码示例和讨论如何正确地从Dirichlet（n）（又名simplex）中进行采样，通过gammavariate（）或通过像您所做的-log（U）以及对潜在角落情况的适当处理link。我发现你的代码存在问题，比如对于采样维度=2的simplex，你得到了三个！均匀数，但在列表推导式中跳过了一个。这是错误的。要从n维Dirichlet中采样，您应该恰好获得n个U（0,1），然后进行转换（或从gammavariate中取n个样本）。但是，最好的解决方案可能只是使用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 吗？”

“更新”

“看起来我们可能对称谓单纯形、狄利克雷等有不同的理解。我更倾向于 this definition 中的定义，在 d 维空间中，有 d 个点，d-1 单纯形是连接顶点的凸包。由于坐标之间的关系，单纯形维数总是比空间少一。在最简单的情况下，d=2，1 单纯形是连接点 (1,0) 和 (0,1) 的线段，从狄利克雷分布中可以得到这张图片。”

在 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()

- Severin Pappadeux

---

以上条件确保了z-x位于单纯形区域内，这是我们进行密度评估所需的。因此，我计算满足上述条件的单纯形中点的比例，这是概率密度函数的估计值。 - pikachuchameleon

此外，对于在单纯形内生成点的过程，我没有使用您指出的狄利克雷分布程序。但是我的程序是这样的：如果U1，...，U_n + 1的速率为1，则（U1 / U_1 +..U_n + 1，.....，U_n / U_1 +....+ U_n + 1）在单纯形上均匀分布。这就是为什么在列表理解期间我要跳过一个的原因。 - pikachuchameleon