Generate Random Joint from Known Marginals

I hope you are doing well. I've got a question.

My ultimate goal is to model an arbitrary language based on a pre-defined set of marginals and joints.

Let us assume I want to model language for words $a...d$ and I have known marginals $p(a)$ to $p(d)$. We also know that $p(a)>p(b)>...$ and so on (Zipfs Law).

1) I want to generate a symmetric random joint probability matrix $p(a,b)$ with zero diagonal based on the marginals:

    N=4
p_of_w = [ 0.48,  0.24,  0.16,  0.12]

p_of_w_rev = list(reversed(p_of_w)) #Reverse ordered marginals

jpm = np.zeros([N,N]) # Zero joint prob mat

for i in range(N-1):
p_vector = np.random.dirichlet(np.ones(N-i-1)) * (p_of_w_rev[i] - np.sum(jpm[i]))
p_vector = np.insert(p_vector, 0,0)
jpm[i,i:N] = jpm[i:N,i] = p_vector


The np.sum(jpm[:,i]) and np.sum(jpm[i]) are all equal to the corresponding marginals, except for the last one. There's something wrong with my reasoning, however I am somewhat stuck.

2) Is there a more elegant solution to this problem?

3) Let us assume I want to introduce some sparsity in the joint which is inversely proportional to the corresponding marginal $p(a)$. I.e. higher: marginal, less probability of sparse entries in the corresponding row of the joint. Can you think of a solution?

Ok, so its technically possible under the relaxation that jmp[0,0] and jmp[-1,-1] are allowed to be nonzero.

    p_of_w_rev = list(reversed(p_of_w))
jpm = np.zeros([N,N])

for i in range(N):
if i > 0 and i != N-1:
p_vector = np.random.dirichlet(np.ones(N-i-1)) * (p_of_w_rev[i] - np.sum(jpm[i]))
p_vector = np.insert(p_vector, 0,0)
jpm[i,i:N] = jpm[i:N,i] = p_vector
elif i == 0:
p_vector = np.random.dirichlet(np.ones(N-i)) * (p_of_w_rev[i] - np.sum(jpm[i]))
jpm[i,i:N] = jpm[i:N,i] = p_vector
elif i == N-1:
jpm[i,i] = p_of_w_rev[i] - np.sum(jpm[i])