Schätzung der Wahrscheinlichkeitsdichte der Summe einheitlicher Zufallsvariablen in Python

Question

Ich habe zwei Zufallsvariablen X und Y, die auf dem Simplex gleichmäßig verteilt sind:

Ich möchte die Dichte ihrer Summe bewerten:

nach der obigen Integral Auswertung mein Endziel ist es, das folgende Integral zu berechnen:

th zu berechnen Als erstes Integral erzeuge ich gleichmäßig verteilte Punkte in Simplex und überprüfe dann, ob sie zu der gewünschten Region im obigen Integral gehören und nehme den Bruchteil von Punkten, um die obige Dichte zu bewerten.

Sobald ich die obige Dichte berechnet habe, befolge ich ein ähnliches Verfahren, um das obige Logarithmusintegral zu berechnen, um seinen Wert zu berechnen. Dies war jedoch extrem ineffizient und dauerte so 3-4 Stunden. Kann mir jemand einen effizienten Weg vorschlagen, dies in Python zu lösen? Ich benutze Numpy-Paket. Hier

ist der Code

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

nicht sicher, es ist eine Antwort auf Ihre Frage, aber sie

Ersten Start, hier sind einige Code-Beispiele und Diskussion, wie man richtig von Dirichlet Probe (n) (simplex aka), über gammavariate() oder über -log(U) wie du aber mit der richtigen Griff für mögliche Ecke Fall link

Problem mit Ihrem Code, wie ich sehen kann, ist, dass, sagen wir, für die Probenahme dimen sion = 2 simplex Sie erhalten drei (!) einheitliche Zahlen, überspringen aber eins, wenn Sie Listenverstehen für x machen. Das ist falsch. Um n-dimensionales Dirichlet zu testen, sollten Sie genau n U (0,1) erhalten und dann (oder n Proben von gammavariate) transformieren.

Aber, beste Lösung könnte nur numpy.random.dirichlet() verwenden, ist es in C geschrieben und könnte von allen am schnellsten sein, siehe link.

Letzte, meiner bescheidenen Meinung nach, sind Sie nicht richtig schätzen log(PDF(X+Z)). Ok, du findest einige, aber was ist PDF(X+Z) an dieser Stelle?

Ist diese

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

wie PDF aussieht? Wie hast du es geschafft?

Einfacher Test: Integration von PDF(X+Z) über ganze X+Z Bereich. Hat es 1 produziert?

UPDATE

Sieht aus wie wir verschiedene Ideen haben könnten, was wir simplex nennen, Dirichlet usw. bin ich ziemlich viel zusammen mit this definition, wo in d dim Raum, den wir d Punkte und d-1 simplex haben konvex Rumpf verbindet Scheitel . Simplex-Dimension ist immer eins weniger als Platz aufgrund der Beziehung zwischen Koordinaten.Im einfachsten Fall d = 2 ist, ist 1-simplex-Segment ein Verbindungspunkte (1,0) und (0,1) und von Dirichlet Verteilung habe ich das Bild erhielt

Bei von d = 3 und 2-simplex wir haben Dreieck Verbindungspunkte (1,0,0), (0,1,0) und (0,0,1)

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