Python/numpy - conditional sampling of variables, distribution of subsequent value is based on result of previous value

Question

I am trying to generate a random sample of multiple variables which are loosely related to each other. Meaning that "allowed" values of some variables depend on the value which is set for another variable.

For simplicity let's imagine that I have just two variables - A and B and let's say that both of them have uniform or gaussian distribution (we don't really care which exact distribution they follow and can accept both). For discussion let's assume both have uniform distribution.

Let's say that variable A can take any value between 0 and 100. We can easily sample from this distribution, say, 1000 data points.

Now, we also want to generate values for variable B, which can take any value between, say, 50 and 150. The catch here is that there is a constraint in resulting sample - sum of values A and B must be between 60 and 160. Constraints generally can be similar - with limiting boundaries of acceptable values in combination of variables - but there is number of them. Sometimes also variables depend on each other - for instance some variable C again can be random but needs to be, for instance such that A<=C<=300. Just to give an idea.

Final catch is that each time we run the sampling process precise boundaries of sampling are changing (for example in one case A can be between 0 and 100 as above, next day it needs to be between -10 and 75 etc). Basically from day to day precise boundaries of sampling are evolving.

Right now we do it in a very inefficient way - generate completely random grid of A and B values independently, than eliminate all of the A and B combinations which don't satisfy constraints which we specify and than use them in subsequent steps. For example such grid could look like:

However, as you guess it is super-inefficient. In reality we have a lot of variables (30+) and large set of constraints we apply. Completely random generation of grid leads to instances where after applying all constraints we end up with no points satisfying all constraints if we don't use large enough sample size - and to ensure we always have at least some points we need to generate grid with millions points. Beyond that each time we re-run the sampling procedure we get different resulting dataset - sometimes all points are getting eliminated, sometimes we get 10 points as result and sometimes - 1000.

So my question is - is there a way to do it more efficiently in a "statistically correct way", ideally in a way which will allow us to specify how many sample points satisfying all constraints we want to get in the end of the day. Any guidance or pointers to some code examples will be much appreciated.

Answer 1

From what I understand, what you want to do is to sample two numbers $A$ and $B$, that need to pass two constraints:

• they have fixed ranges $ A \in [A_L, A_U]$ and $B \in [B_L, B_U]$,
• their sum falls into a given range $ A+B = C \in [C_L, C_U]$.

If those are the only constraints, it's fairly simple: you can sample conditionally by sampling $A$ from the appropriate uniform distribution, then subtract it from the range for the sum $[C_L-A, C_U-B]$, and then sample $B$ from the intersection of this range and the range for $B$. You can use a similar procedure but starting with $B$, the only thing you need is the interval arithmetics. You can also start with sampling $C$, then $B$ from the interval that makes $C$ possible, then $A=C-B$ if it fits the range of acceptable values.

Answer 2

I'll describe what could be done with the simple example you gave, but it isn't a sure thing that they generalize.

What is known (uniform case): $$\begin{aligned} a | b=60 &\sim \mathrm{Uniform}(0, 100) \\ b | a=10 &\sim \mathrm{Uniform}(50, 150) \\ c | a=0 &\sim \mathrm{Uniform}(0, 300) \\ a + b &\in [60, 160] \\ a &\in [0, c] \end{aligned}$$

This seems like it could be sampled from by calculating the marginal cdfs and then using inverse transform sampling. Sample in the order such that you get the values you need to know exact cdfs of the others. Complement with rejection sampling as you've done previously when necessary.

If your expressions are to complicated, then Metropolis-Hastings or some other Markov Chain Monte Carlo technique would probably be a better bet. The python packages emcee and PyMC3 are okay in my experience.