4
$\begingroup$

What's the best way to model situations with dependencies? What I mean is well-defined arrangements like:

If I choose a 6-sided die with probability $p$ and a 20-sided one otherwise, and then roll it twice, what's the probability distribution on the sum of the two rolls?

Here's my current attempt:

p = 0.3
sides = TransformedDistribution[6+(1-d)*14, d \[Distributed]BernoulliDistribution[p]]    
rolls[dist_] := TransformedDistribution[x+y, {x,y} \[Distributed] ProductDistribution[{dist, 2}]]
sumDist = ParameterMixtureDistribution[rolls[DiscreteUniformDistribution[{1,n}]], n \[Distributed] sides]
RandomVariate[sumDist, 10]

I get this formula for the random variate--it looks like $n$ isn't properly being passed through:

RandomVariate[
ParameterMixtureDistribution[
  TransformedDistribution[\[FormalX]1 + \[FormalX]2, {\[FormalX]1, \
\[FormalX]2} \[Distributed] 
    ProductDistribution[{DiscreteUniformDistribution[{1, n}], 2}]], 
  n \[Distributed] 
   TransformedDistribution[
    6 + 14 (1 - \[FormalX]), \[FormalX] \[Distributed] 
     BernoulliDistribution[0.3]]], 10]

My follow-on question will be to compute the probability of which die was chosen based on seeing a sum of e.g 11, so I'm trying to keep that initial choice in the model explicitly.

$\endgroup$
2
  • 2
    $\begingroup$ Hint: you're over-complicating this... step back and think about it. $\endgroup$
    – ciao
    Commented Aug 21, 2015 at 10:13
  • $\begingroup$ Conditioned[] will be useful here. $\endgroup$ Commented Aug 21, 2015 at 11:44

2 Answers 2

6
$\begingroup$

Here's a start for you. p is probability of choosing die 1, f1/f2 are number of faces (starting at 1) for die 1/2:

p1 = 3/10
f1 = 6
f2 = 20
d = MixtureDistribution[{p1, 1 - p1}, {DiscreteUniformDistribution[{1, f1}], 
                                       DiscreteUniformDistribution[{1, f2}]}];

d2 = TransformedDistribution[a + b, {a, b} \[Distributed] ProductDistribution[{d, 2}]];

(* PMF of sums *)
Probability[x == #, x \[Distributed] d2] & /@ Range[2, 40] // Short

(* {289/40000,289/20000,867/40000,289/10000,<<32>>,147/40000,49/20000,49/40000} *)

You can also just use PDF[d2,z] to get PMF, will be slower...

Random variates from the distribution is just:

RandomVariate[d2,10]
(* {8, 25, 9, 22, 7, 20, 21, 15, 23, 20} *)
$\endgroup$
4
$\begingroup$

Just building on @ciao (@rasher) to deal with second part:

c66 = TransformedDistribution[
   a + b, {a, b} \[Distributed] 
    DiscreteUniformDistribution[{{1, 6}, {1, 6}}]];
p66 = Probability[x == 11, x \[Distributed] c66];
c620 = TransformedDistribution[
   a + b, {a, b} \[Distributed] 
    DiscreteUniformDistribution[{{1, 6}, {1, 20}}]];
p620 = Probability[x == 11, x \[Distributed] c620];
c2020 = TransformedDistribution[
   a + b, {a, b} \[Distributed] 
    DiscreteUniformDistribution[{{1, 20}, {1, 20}}]];
p2020 = Probability[x == 11, x \[Distributed] c2020];
pdice[p_] := 
 With[{r = {p^2, 2 p (1 - p), (1 - p)^2} {p66, p620, p2020}}, {Total@
    r, r/Total@r}]

pdice will produce probability sum of pair of dice is 11 and the probability it came from {6,6}, {6,20}||{20,6}, {20,20} given total is 11 using $P(D={6,6}|Z=11)P(Z=11)=P(Z=11|D={6,6})P(D={6,6})$ etc.

So

Manipulate[
 PieChart[pdice[probability][[2]], 
  ChartLegends -> {"6-6", "6-20", "20-20"}, 
  LabelingFunction -> (Placed[
      Panel[NumberForm[#1, 2], FrameMargins -> 0], 
      "RadialCallout"] &)], {probability, 0, 1, 
  Appearance -> "Labeled"}]

enter image description here

and "reality checking" cf @ciao:

pdice[3/10]
Probability[x == 11, x \[Distributed] d2]

yield:{153/4000, {20/153, 28/51, 49/153}} and 153/4000 respectively.

$\endgroup$
2
  • $\begingroup$ Oh, +1 on the cool-as-*$*% graphics manipulate... I always enjoy your examples! $\endgroup$
    – ciao
    Commented Aug 22, 2015 at 6:38
  • $\begingroup$ @ciao thanks...could have done more tersely and general...alas lots to do :) $\endgroup$
    – ubpdqn
    Commented Aug 22, 2015 at 6:42

Your Answer

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

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