Stima della densità di probabilità della somma di variabili casuali uniformi in python

Question

Ho due variabili aleatorie X e Y, che sono uniformemente distribuiti sul simplex:

voglio valutare la densità della loro somma:

Dopo aver valutato quanto sopra integrale, il mio obiettivo finale è quello di calcolare il seguente integrale:

per calcolare ° e primo integrale, sto generando punti uniformemente distribuiti in simplex e poi controllando se appartengono alla regione desiderata nell'integrale di cui sopra e prendendo la frazione di punti per valutare la densità di cui sopra.

Una volta calcolata la densità di cui sopra, sto seguendo una procedura simile per calcolare l'integrale logaritmo di cui sopra per calcolare il suo valore. Tuttavia, questo è stato estremamente inefficiente e richiedendo molto tempo per 3-4 ore. Qualcuno mi può suggerire un modo efficace per risolvere questo in Python? Sto usando il pacchetto Numpy.

Ecco il codice

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

Non sono sicuro che sia una risposta alla tua domanda, ma cominciamo

In primo luogo, ecco alcune esempi di codici e la discussione come campionare correttamente da Dirichlet (n) (aka simplex), via gammavariate() o via -log(U) come avete fatto, ma con una corretta impugnatura per il potenziale caso d'angolo, link

un problema con il codice che posso vedere è che, per esempio, per le dimensioni di campionamento sion = 2 simplex stai ottenendo tre (!) numeri uniformi, ma saltandone uno quando si esegue la comprensione degli elenchi per x. Questo è sbagliato. Per campionare Dirichlet n-dimension devi ottenere esattamente n U (0,1) e trasformarli (o campioni n da gammavariate).

Ma, la soluzione migliore potrebbe essere solo utilizzare numpy.random.dirichlet(), è scritto in C e potrebbe essere il più veloce di tutti, vedere link.

Ultimo, a mio modesto avviso, non stai valutando correttamente log(PDF(X+Z)). Ok, ci sono alcuni, ma cosa è PDF(X+Z) a questo punto?

Questo

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

assomiglia PDF? Come sei riuscito ad averlo?

Test semplice: integrazione di PDF(X+Z) per l'intera area X+Z. Ha prodotto 1?

UPDATE

Sembra che potremmo avere idee diverse quello che noi chiamiamo simplex, Dirichlet vertici ecc Sono praticamente insieme a this definition, dove nel d spazio dim abbiamo d punti e d-1 simplex è convesso di collegamento . La dimensione Simplex è sempre uno spazio inferiore a causa della relazione tra le coordinate.Nel caso più semplice, d = 2, 1-simplex è un segmento che collega i punti (1,0) e (0,1), e dalla distribuzione di Dirichlet Ho l'immagine

Nel caso di d = 3 e 2-simplex abbiamo punti triangolo collegamento (1,0,0), (0,1,0) e (0,0,1)

codice, pitone

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