Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
add comment

3 Answers 3

up vote 6 down vote accepted

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 :)

share|improve this answer
3  
P.S.: I felt compelled to answer with a workaround for version 7, since you've been a good sport (in the comments) at accepting the conventions of this community and show a willingness to learn and participate. Enjoy! –  rm -rf Mar 13 '13 at 18:18
    
Thanks a lot. This solution should be very usefull to me (and hopefully to others too). I was having difficulties with the lack of custom distributions in Mma 7, and my 3D distributions of points (for nebula models) were too simplistic. –  Cham Mar 13 '13 at 18:26
    
@rm-rf You ought to revisit your use of RandomSample, because the discretization it introduces may be too crude for some applications. At least it's worth noting that you are making an approximation there. Also, it appears there is an error in its implementation: sample output suggests you are not implementing a true mixture, in which the proportions of outcomes from each component are random: these proportions look fixed. Consequently, in simulations based on random variates, the variances will be too low and the empirical distributions too uniform. –  whuber Mar 13 '13 at 18:34
    
I just feel obliged to point out that this is not mere nit-picking: subtle errors in implementation of RNGs can create profound differences in applications. For instance, in any simulation based on a mixture in which the bottom 50% of the distributions have no overlap with the top 50%, the median might have no variation in your implementation, whereas the correct median could have arbitrarily high variation. This can lead to subtle, elusive bugs that might not even be detected (the worst kind: they just create wrong results). That's a lot to brush under the rug with a "Q&D" disclaimer! –  whuber Mar 13 '13 at 18:45
    
@whuber Please see my edit –  rm -rf Mar 13 '13 at 19:11
show 3 more comments

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.

share|improve this answer
    
I am putting a constructive interpretation on the question, taking it to ask "how do I implement MixtureDistribution so I can go on with my coding," rather than reading it literally in the sense of "what is wrong," which--especially because this problem pertains to an older version of MMA--seems unproductive. –  whuber Mar 13 '13 at 18:00
    
This is interesting. Thanks a lot for this solution. I'll try to work with it. –  Cham Mar 13 '13 at 18:19
add comment

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

share|improve this answer
3  
I should add that this seems to imply that the OP hasn't even bothered to look up the function he was asking questions about. The missing documentation in v7 should have provided a hint. –  Sjoerd C. de Vries Mar 13 '13 at 12:31
    
I did checked the help in Mma, of course. And it did responded that the command existed, but didn't gave any details. I was suspecting that the command wasn't implemented in that version. This is why I asked the question here, and also to know if there was an alternative to mix several distributions in Mma 7.0. So I don't understand the negative reactions here. There are other users who still use Mma 7.0, and it's usefull to know that we actually can't mix distributions in that version. I now feel discouraged to ask any question on this forum ! :-( –  Cham Mar 13 '13 at 15:31
    
You can answer questions as well, you know! :) –  cormullion Mar 13 '13 at 16:07
    
@cormullion : what do you mean ? –  Cham Mar 13 '13 at 16:11
    
@Cham I agree with Sjoerd. "I now feel discouraged to ask any question on this forum ! :-(" - you shouldn't feel discouraged. However, we must maintain the reasonably high standards for the questions asked here, if we want this place to be the world's primary place to ask professional questions about Mathematica, and not degrade into student forum or the like. With all due respect, the sort of information you asked for is trivial to locate using Help. –  Leonid Shifrin Mar 13 '13 at 16:23
show 4 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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