# Draw independent samples from conditional distribution at a fast rate

I would like to draw samples from the following conditional distribution of the following (don't ask why):

fRandomFunction[x_, y_] :=
PDF[MultinormalDistribution[{3, 3}, {{0.5, 0}, {0, 0.5}}], {x, y}] +
PDF[MultinormalDistribution[{6, 6}, {{0.6, -0.5}, {-0.5, 0.6}}], {x,
y}] + PDF[MultinormalDistribution[{6, 6}, {{1, 0.5}, {0.5, 1}}], {x, y}] +
PDF[MultinormalDistribution[{6, 4}, {{0.6, 0.5}, {0.5, 0.6}}], {x,
y}] + PDF[MultinormalDistribution[{2, 8}, {{0.2, 0}, {0, 0.2}}], {x, y}]


At the moment I am using Mathematica's inbuilt ProbabilityDistribution function to draw individual independent samples:

fGenerateSample[aX_] :=
RandomVariate[
ProbabilityDistribution[fRandomFunction[aX, y], {y, 0, Infinity}], {1}
][[1]]


If I map this onto a random list of X inputs, this takes quite a while to output an independent sample:

lInputs = RandomVariate[NormalDistribution[5, 1], {1}];
lSamples = fGenerateSample /@ lInputs // Timing;


This takes about 1/20 second to execute on my computer. I would like a way of speeding this process up. Does anyone know of a way I could generate a single sample at a faster rate from the conditional density here?

Note this problem is quite different from my actual problem, where I need to essentially do a NestList with my fGenerateSample function, up to 100,000 times. That is why speed is of the essence.

• You said not to ask "why" but why do you call fRandomFunction a conditional distribution when it is the sum of 5 multivariate normal densities. Also, fRandomFunction is not a proper density until you divide by 5. And finally why does y go from 0 to positive infinity when fRandomFunction deals with values of y between minus and plus infinity? Or am I totally misunderstanding what you've written. (My running of your code takes about a minute rather than 1/20 of a second so I also need a new machine.)
– JimB
Apr 4, 2016 at 20:56
• I must agree with @JimBaldwin: "...don't ask why.", "... this problem is quite different from my actual problem..." - then why aren't you posting the actual problem, and for that matter, what is so Double-Naught Spy secret about the actual goal? There are many here that are not only Mathematica experts, but mathematicians & probabilists as well. Knowing what the end goal is for those readers may well result in a far superior way to meet the needs.
– ciao
Apr 5, 2016 at 6:03
• @ciao Well, at least the "actual" problem here is acknowledged. Apr 5, 2016 at 13:50

I agree with the comments. I am uncertain what the aim is, so I apologize if this is unhelpful. I only post to motivate clarification.

If the aim is a marginal of mixture of binormals:

m1 = MultinormalDistribution[{3, 3}, {{0.5, 0}, {0, 0.5}}];
m2 = MultinormalDistribution[{6, 6}, {{0.6, -0.5}, {-0.5, 0.6}}];
m3 = MultinormalDistribution[{6, 6}, {{1, 0.5}, {0.5, 1}}];
m4 = MultinormalDistribution[{6, 4}, {{0.6, 0.5}, {0.5, 0.6}}];
m5 = MultinormalDistribution[{2, 8}, {{0.2, 0}, {0, 0.2}}];
mix = MixtureDistribution[{1, 1, 1, 1, 1}, {m1, m2, m3, m4, m5}];
marg = MarginalDistribution[mix, 2];
Histogram[RandomVariate[marg, 100000]]
Histogram3D[RandomVariate[mix, 10000]]
Histogram3D[RandomVariate[#, 1000] & /@ {m1, m2, m3, m4, m5}]


If the aim is sum of binormals:

td = TransformedDistribution[{x1, y1} + {x2, y2} + {x3, y3} + {x4,
y4} + {x5, y5}, {{x1, y1} \[Distributed]
m1, {x2, y2} \[Distributed] m2, {x3, y3} \[Distributed]
m3, {x4, y4} \[Distributed] m4, {x5, y5} \[Distributed] m5}];
ListPlot[RandomVariate[td, 100000]]


Again apologies if this is unhelpful.

• I'd take a wild guess and imagine, that the 5 distributions can be arranged into one for a vector of length 10, possibly leading to a performance boost. Unless of course MMA does that under the hood. Apr 5, 2016 at 13:54