Estimating the probability density of sum of uniform random variables in python

Not sure it is an answer to your question, but let's start

First, here is some code samples and discussion how to properly sample from Dirichlet(n) (a.k.a. simplex), via gammavariate() or via -log(U) as you did but with proper handle for potential corner case, link

Problem with your code as I can see is that, say, for sampling dimension=2 simplexyou're getting three (!) uniform numbers, but skipping one when doing list comprehension for x. This is wrong. To sample n-dimension Dirichlet you should get exactly n U(0,1) and transform then (or n samples from gammavariate).

But, best solution might be just use numpy.random.dirichlet(), it is written in C and might be fastest of all, see link.

Last one, in my humble opinion, you're not properly estimating log(PDF(X+Z)). Ok, you find some are, but what is PDF(X+Z) at this point?

Does this

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

looks like PDF? How did you managed to get it?

Simple test: integration of PDF(X+Z) over whole X+Z area. Did it produce 1?

UPDATE

Looks like we might have different ideas what we call simplex, Dirichlet etc. I'm pretty much along with this definition, where in d dim space we have d points and d-1 simplex is convex hull connecting vertices. Simplex dimension is alwaysone less than space due to relation between coordinates. In simplest case, d=2, 1-simplex is a segment connecting points (1,0) and (0,1), and from Dirichlet distribution I've got the picture

In the case of d=3 and 2-simplex we have triangle connecting points (1,0,0), (0,1,0) and (0,0,1)

Code, Python

from mpl_toolkits.mplot3d import Axes3Dimport matplotlib.pyplot as pltimport mathimport randomdef 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 rdef 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 rN = 500 # numer of points to plotd = 3   # dimension of the space, 2 or 3x = []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()