constrained random numbers sampling using python (Monte-Carlo, Markov chains, pymc)

Question

I am trying to sample random numbers with constraints using Python and pymc library.

Here mins and maxes are arrays of minimums and maximums for each of 22 variables. It works fine in this case.

from pymc import *
X = Uniform('X', mins, maxes)
@potential
def s(X=X): 
 cons1 = X[0] < 225
 cons2 = X[0] > 405
 if cons1 or cons2:
     return -inf
 else:
     return 0.0
mc = MCMC([X, s])
mc.sample(10000)

However when I try to add one more constraint: that sum of variables must be equal to 1000, everything breaks:

from pymc import *
X = Uniform('X', mins, maxes)
@potential
def s(X=X): 
 cons1 = X[0] < 225
 cons2 = X[0] > 405
 cons3 = X.sum() >= 1000
 if cons1 or cons2 or cons3:
     return -inf
 else:
     return 0.0
mc = MCMC([X, s])
mc.sample(10000)

Error is: 'ZeroProbability: Potential s forbids its parents' current values'

Could you please recommend some solutions to this problem?

It looks like your constraints might be incompatible, that is no variable values exist to satisfy them. What are these mins and maxes? — Andrei, Dec 05 '17 at 11:35
Andrei, mins and maxes are intervals for each variable. Constraints are compatible: I have at least one solution calculated by hand. (but there are endless variety of them) However I need to obtain many more uniformly distributed samples. — dark980, Dec 05 '17 at 11:41
Piinthesky, I need to generate random numbers withing grouped constraints i.e. x1+x2+x3 <= 600; x3+x4+x5 >= 300; x1+x2+x3 +x4+...+xn == 1000; — dark980, Dec 05 '17 at 11:42
Generate random `x3` and with this constraint randomly the other `xi`? Or generate random n-tuple and discard, if they violate any constraint? — Mr. T, Dec 05 '17 at 11:51
Piinthesky, it was my first though, too. However with endless constraints this way of solving the problem is not working. (as well as many others that I tried) In the end, I come to Markov chains. Now I want to understand why this error happens and how to solve it. — dark980, Dec 05 '17 at 11:56
I think this answer explains what is going on here: https://stats.stackexchange.com/a/225958 — Andrei, Dec 05 '17 at 12:14
Ah, I didn't realize that your constraints are just examples of a more elaborate set of rules. Did you try, if `cons3` alone reproduces the error? — Mr. T, Dec 05 '17 at 12:15
@dark980 `x1+x2+x3 <= 600; x3+x4+x5 >= 300; x1+x2+x3 +x4+...+xn == 1000;` those are all conditions or you have more? How much is `n`? — Severin Pappadeux, Dec 06 '17 at 01:25

score 1 · Answer 1 · answered Dec 05 '17 at 17:33

Consider sampling from just two constrained uniform random variables on the standard interval, as pictured here.

The brown area represents the constrained sampling area based on some inequalities. The diagonal line represent the sum.

Notice: In 2-space the brown area has an area (otherwise known as a probability or measure) but the diagonal line has zero area, which is to say zero probability or measure.

The same holds true for higher dimensional spaces.

constrained random numbers sampling using python (Monte-Carlo, Markov chains, pymc)

1 Answers1