अजगर

Question

में एक समान यादृच्छिक परिवर्तनीय की राशि की प्रायिकता घनत्व आकलन मैं दो यादृच्छिक चर एक्स और वाई जो समान रूप से सिंप्लेक्स पर वितरित कर रहे हैं:

:

मैं उनका योग के घनत्व का मूल्यांकन करना चाहते

वीं की गणना करने के लिए:

ऊपर अभिन्न मूल्यांकन करने के बाद, मेरा अंतिम लक्ष्य अभिन्न निम्नलिखित गणना करने के लिए है ई पहला अभिन्न, मैं सरल में समान रूप से वितरित बिंदु उत्पन्न कर रहा हूं और फिर जांच कर रहा हूं कि वे उपर्युक्त अभिन्न अंग में वांछित क्षेत्र से संबंधित हैं और उपरोक्त घनत्व का मूल्यांकन करने के लिए अंक का अंश ले रहे हैं।

एक बार जब मैं उपरोक्त घनत्व की गणना करता हूं तो मैं इसके मूल्य की गणना करने के लिए उपरोक्त लॉगरिदम अभिन्न की गणना करने के लिए एक समान प्रक्रिया का पालन कर रहा हूं। हालांकि, यह बेहद अक्षम है और इस तरह के 3-4 घंटे बहुत समय ले रहा है। क्या कोई मुझे पाइथन में हल करने का एक प्रभावी तरीका सुझा सकता है? मैं 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]) 

Answer 1

सुनिश्चित नहीं हैं कि यह आपके प्रश्न का उत्तर है, लेकिन चलो शुरू

सबसे पहले जाने है, यहाँ कुछ कोड नमूने और चर्चा कैसे ठीक (एन) के Dirichlet से नमूने के लिए है अपने कोड के साथ (उर्फ सिंप्लेक्स), gammavariate() के माध्यम से या -log(U) तुमने किया था के रूप में के माध्यम से, लेकिन संभावित कोने मामले के लिए उचित संभाल के साथ, link

समस्या के रूप में मैं देख सकते हैं कि कहते हैं, है, नमूने Dimen के लिए sion = 2 simplex आपको तीन (!) वर्दी संख्याएं मिल रही हैं, लेकिन x के लिए सूची समझने पर एक को छोड़ना। ये गलत है। एन-आयाम नमूना करने के लिए Dirichlet आपको n यू (0,1) प्राप्त करना चाहिए और फिर (या n स्तनपायी से नमूने) को बदलना चाहिए।

लेकिन, सबसे अच्छा समाधान केवल numpy.random.dirichlet() का उपयोग किया जा सकता है, यह सी में लिखा गया है और सभी में से सबसे तेज़ हो सकता है, 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) 

पीडीएफ की तरह लग रहा है? आप इसे कैसे प्राप्त करने में कामयाब रहे?

सरल परीक्षण: पूरे X+Z क्षेत्र में PDF(X+Z) का एकीकरण। क्या यह 1 उत्पादन किया?

अद्यतन

ऐसा लगता है कि हम अलग अलग विचारों क्या हम सिंप्लेक्स फोन, Dirichlet आदि मैं काफी this definition, जहां d मंद अंतरिक्ष में हम d अंक है और d-1 सिंप्लेक्स के साथ उत्तल पतवार कनेक्ट कर रहा है कर रहा हूँ कोने हो सकता है । सिंपलक्स आयाम हमेशा समन्वय के बीच संबंध के कारण अंतरिक्ष से कम है।सामान्य स्थिति में, d = 2, 1-सिंप्लेक्स एक खंड जोड़ने अंक (1,0) और (0,1) है, और Dirichlet वितरण से मैं चित्र

मामले में मिल गया है d = 3 और 2-सिंप्लेक्स की हमारे पास त्रिकोण जोड़ने अंक (1,0,0), (0,1,0) और (0,0,1)

संहिता, अजगर

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