0
$\begingroup$

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.

$\endgroup$
3
  • 4
    $\begingroup$ 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.) $\endgroup$
    – JimB
    Apr 4, 2016 at 20:56
  • 2
    $\begingroup$ 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. $\endgroup$
    – ciao
    Apr 5, 2016 at 6:03
  • $\begingroup$ @ciao Well, at least the "actual" problem here is acknowledged. $\endgroup$
    – LLlAMnYP
    Apr 5, 2016 at 13:50

1 Answer 1

3
$\begingroup$

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}]

enter image description here

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]]

enter image description here

Again apologies if this is unhelpful.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – LLlAMnYP
    Apr 5, 2016 at 13:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.