Ước tính mật độ xác suất của tổng của các biến ngẫu nhiên thống nhất trong python

Question

Tôi có hai biến ngẫu nhiên X và Y được phân bố đều trên simplex:

tôi muốn đánh giá mật độ tổng của chúng:

sau khi đánh giá tích ở trên, mục tiêu cuối cùng của tôi là để tính tích phân sau:

để tính thứ e tích phân đầu tiên, tôi tạo ra các điểm phân bố đồng nhất trong simplex và sau đó kiểm tra nếu chúng thuộc về vùng mong muốn trong tích phân trên và lấy phần điểm để đánh giá mật độ trên.

Khi tôi tính toán mật độ trên, tôi đang làm theo một quy trình tương tự để tính tích phân lôgarit ở trên để tính toán giá trị của nó. Tuy nhiên, điều này đã vô cùng không hiệu quả và mất rất nhiều thời gian 3-4 giờ. Bất cứ ai có thể cho tôi một cách hiệu quả để giải quyết điều này trong Python? Tôi đang sử dụng gói Numpy.

Đây là mã

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]) 

Answer 1

Không chắc đó là một câu trả lời cho câu hỏi của bạn, nhưng chúng ta hãy bắt đầu

Thứ nhất, đây là một số mẫu mã và thảo luận làm thế nào để đúng cách lấy mẫu từ Dirichlet (n) (aka simplex), thông qua gammavariate() hoặc thông qua -log(U) như bạn đã làm nhưng với xử lý thích hợp cho trường hợp góc tiềm năng, link

Vấn đề với mã của bạn như tôi có thể thấy là, ví dụ để lấy mẫu làm mờ sion = 2 simplex bạn đang nhận được ba (!) số đồng phục, nhưng bỏ qua một khi thực hiện việc hiểu danh sách cho x. Cái này sai. Để lấy mẫu Dirichlet n-dimension bạn sẽ nhận được chính xác n U (0,1) và biến đổi sau đó (hoặc n mẫu từ gammavariate).

Nhưng, giải pháp tốt nhất có thể chỉ sử dụng numpy.random.dirichlet(), được viết bằng C và có thể là nhanh nhất, xem link.

Cuối cùng, theo ý kiến ​​khiêm tốn của tôi, bạn không ước tính đúng cách log(PDF(X+Z)). Ok, bạn tìm thấy một số, nhưng PDF(X+Z) vào thời điểm này là gì?

Liệu này

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

trông giống như PDF? Làm thế nào bạn quản lý để có được nó?

Thử nghiệm đơn giản: tích hợp PDF(X+Z) trên toàn bộ khu vực X+Z. Nó đã sản xuất 1?

CẬP NHẬT

Hình như chúng ta có thể có những ý tưởng khác nhau mà chúng ta gọi simplex, Dirichlet vv Tôi khá nhiều cùng với this definition, nơi ở d không gian mờ chúng tôi có d điểm và d-1 simplex được thân lồi kết nối đỉnh . Kích thước đơn giản luôn là một ít hơn không gian do mối quan hệ giữa các tọa độ.Trong trường hợp đơn giản nhất, d = 2, 1-simplex là một đoạn kết nối điểm (1,0) và (0,1), và từ phân phối Dirichlet Tôi đã có hình ảnh

Trong trường hợp của d = 3 và 2-simplex chúng tôi có hình tam giác kết nối điểm (1,0,0), (0,1,0) và (0,0,1)

Mã, 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()