Oszacowanie gęstości prawdopodobieństwa sumy jednolitych zmiennych losowych w Pythonie

Question

Mam dwie zmienne losowe X i Y, które są równomiernie rozmieszczone na simplex:

chcę oceniać gęstość ich sumy:

Po przeanalizowaniu powyższych integralną mój ostatecznym celem jest obliczenie następującej całki:

aby obliczyć th Pierwszą całką generuję równomiernie rozłożone punkty w simpleksie, a następnie sprawdzam, czy należą one do pożądanego regionu w powyższej całce i biorę ułamek punktów, aby ocenić powyższą gęstość.

Po obliczeniu powyższej gęstości postępuję w podobny sposób, aby obliczyć powyższy logarytm całkowy, aby obliczyć jego wartość. Jednak było to niezwykle nieefektywne i wymagało dużo czasu, na przykład 3-4 godzin. Czy ktoś może zaproponować mi skuteczny sposób rozwiązania tego w Pythonie? Używam pakietu Numpy.

Oto kod

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

Nie jestem pewien, że jest to odpowiedź na twoje pytanie, ale zacznijmy

Po pierwsze, oto niektóre przykłady kodu i dyskusja jak prawidłowo próbki z Dirichlet (n) (aka simplex), poprzez gammavariate() lub poprzez -log(U) jak ty, ale z odpowiednim uchwytem do potencjalnego przypadku narożnego, link

problem z kodem jak widzę to, że powiedzmy, dla DIMEN próbkowania sion = 2 simplex otrzymujesz trzy (!) jednolite numery, ale pomijasz jedną, gdy robisz spis treści ze zrozumieniem dla x. To jest źle. Aby wypróbować Dirichlet o wymiarze n, powinieneś uzyskać dokładnie n U (0,1) i następnie przekształcić (lub n próbek z gammavariate).

Ale najlepszym rozwiązaniem może być po prostu użycie numpy.random.dirichlet(), jest napisane w C i może być najszybsze ze wszystkich, zobacz link.

Ostatnia, moim skromnym zdaniem, nie oceniasz prawidłowo log(PDF(X+Z)). Ok, niektóre są, ale co to jest PDF(X+Z) w tym momencie?

Czy to

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

wygląda PDF? Jak udało ci się go zdobyć?

Prosty test: integracja PDF(X+Z) na całym obszarze X+Z. Czy to spowodowało 1?

UPDATE

Wygląda na to, możemy mieć różne pomysły, co nazywamy Simplex, Dirichleta wierzchołki itp mam dość dużo wraz z this definition, gdzie w d przyćmionym przestrzeni mamy d punktów i d-1 simplex jest wypukła łączące kadłub . Wymiar prostokąta jest zawsze o jeden mniejszy niż spacja z powodu relacji między współrzędnymi.W najprostszym przypadku, d = 2, 1-simplex jest segment łączący punkty (1,0) i (0,1), a od rozkładu Dirichleta Mam obraz

w przypadku od d = 3, 2-simplex mamy trójkąt łączący punkty (1,0,0), (0,1,0) i (0,0,1)

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