Note sure if PYMC is the best place for this question as I am not sure that it can be done by sampling, but here it is:

Suppose that I have the following:

- Three Dirichlet distributions of dimension k: distr1, distr2, distr3
- A vector A of dimension 3 and a vector B of dimension k.

Does anyone know about a way to find 3 probability vectors p1, p2 and p3 such that:

- p1, p2 and p3 are in an area of high probability for distr1, distr2, distr3 respectively
- They satisfy the linear constraint A[0]*p1 + A[1]*p2 + A[2]*p3 = B

(You can assume that the set of probability vectors p1, p2, p3 satisfying the linear constraint is not empty)

Intuitively the idea is to sample the probability vectors from each Dirichlet probability distributions, but shift the sampling in directions so that we get closer to the linear constraints being satisfied.

Any help or direction appreciated in case you have ever seen a similar case of constrained problem

thank you!