MixtureDistribution in Mathematica 7

Question

Is it possible to use a mixture of several distributions in Mathematica 7.0? Consider for example the following code :

arms = Table[
         NestWhileList[
            {#[[1]] + RandomInteger[{-2, 1}], # + RandomReal[{-1, 1}] & /@ #[[2]]} &, 
            {RandomInteger[{25, 30}], {0, 0, 0}}, #[[1]] > 15 &],
            {10}
       ];

points = Table[# + RandomReal[
                    MixtureDistribution[
                       {1, 1}, 
                       {NormalDistribution[-0.1, 0.3], NormalDistribution[0.1,0.3]}
                    ]] & /@ #[[2]], {#[[1]]}] &/@Flatten[arms, 1];

ListPointPlot3D[points, BoxRatios -> {1, 1, 1}, ImageSize -> 800, 
                PlotStyle -> {Blue, PointSize[Small]}]

The command "MixtureDistribution" isn't recognized in Mathematica 7.0.

The question is : Is there a way to mix distributions in Mathematica 7 ? I never found any documentation on this, so I suspect that the answer is simply "No. Upgrade to v8".

Sometimes, in the software world, there are options and commands that aren't well documented. So this is why I'm asking the question anyway.

Answer 1

Although version 7 does not support creating arbitrary distributions, implementing a MixtureDistribution for your needs is not all that hard. The following shows you how, and implemented so that the syntax and usage is the same as in v8+.

I'll start with a CDF for MixtureDistribution so that you get the idea:

Clear@MixtureDistribution
MixtureDistribution /: CDF[MixtureDistribution[wts_List, dist_List]] :=
    With[{normWts = Normalize[wts, Total], cdfs = CDF[#, \[FormalX]] & /@ dist}, 
        Evaluate[Total[normWts cdfs] /. \[FormalX] -> #] &
    ]
MixtureDistribution /: CDF[m : MixtureDistribution[wts_List, dist_List], x_] := CDF[m]@x

Now let's try to plot the CDF for a custom mixture distribution:

ℳ = MixtureDistribution[{1, 1}, {NormalDistribution[], NormalDistribution[4, 1/2]}];
Plot[CDF[ℳ]@x, {x, -3, 6}, Filling -> Axis]

You can also add an up-value (look up UpValues) for RandomReal to generate random numbers from this distribution. (Note that this is a quick and dirty solution to help you progress with whatever v8 code you have).

MixtureDistribution /: RandomReal[MixtureDistribution[wts_List, dist_List], n_Integer] := 
    RandomReal @@@ Tally@RandomChoice[Normalize[wts, Total] -> dist, n] // Flatten // RandomSample

MixtureDistribution /: RandomReal[m : MixtureDistribution[wts_List, dist_List]] := 
    First@RandomReal[m, 1]

Now you can use this as you would any other distribution:

RandomReal[ℳ]
(* 1.0797 *)

RandomReal[ℳ, 10]
(* {-0.777609, 0.196815, 4.14167, 3.51495, -0.465554, 3.56857, 3.9301,  4.40071, 3.47119, -0.244056} *)

Histogram[RandomReal[ℳ, 5000], {-3, 6, 0.1}]

Adding additional definitions that you need (for PDF, multiple dimensions in RandomReal, etc.) is left as an exercise for you :)

Answer 2

For the limited purpose of generating random variates, you can implement MixtureDistribution yourself:

RandomReal[mixtureDistribution[p_List, f_List], opts___] ^:= 
  Block[{q = p / Plus @@ p, g, i, j, x},
   i = RandomChoice[q -> Range[Length[q]], opts];
   x = RandomReal[#, opts] & /@ f;
   j = RotateRight[Range[Depth[i]], 1];
   MapThread[Part, {Transpose[x, j], i}, Depth[i] - 1]
   ];

This generates a set of indexes i into the mixture according to the specified probability distribution p (a list of positive numbers which will be normalized to sum to unity), in parallel with a similarly-structured set of random values x. At the deepest level, elements of i are indexes and elements of x at that same level are lists of random values obtained from the list of distributions in f. The final line matches each index in i with its corresponding list of values in x, using that index to select the appropriate value.

For example,

f = mixtureDistribution[{2, 1}, {UniformDistribution[{0, 1}], UniformDistribution[{-1, 0}]}];
RandomReal[f]

$0.643672$

RandomReal[f, {2, 3, 1}]

$\{\{\{0.107656\},\{-0.498868\},\{0.577273\}\},\{\{0.162742\},\{0.0406792\},\{0.270292\}\}\}$

Timing indicates some 90% of the effort occurs in the final MapThread operation used to match the indexes in i with the values in x: perhaps someone knowledgeable about optimizing Mathematica structural operations could propose a way to speed this up.

Answer 3

MixtureDistribution was introduced in version 8. You can find information like that on the bottom of the functions' documentation page.