Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: computing expectations for probability

  • To: mathgroup at smc.vnet.net
  • Subject: [mg105300] Re: [mg105288] computing expectations for probability
  • From: danl at wolfram.com
  • Date: Fri, 27 Nov 2009 06:27:38 -0500 (EST)
  • References: <200911260404.XAA09561@smc.vnet.net>

> hi all,
>
> i am trying to see if there is a simplified expression for a
> particular random variable using mathematica. i have two independent
> binomial random variables, X and Y, and i'd like to compute the
> expectation of X/(X+Y). more formally, X is distributed Binomial(p1,
> n) and Y is distributed Binomial(p2, n).  note that p1 and p2 can be
> different, but the n's are the same. i would like to compute E(X/(X
> +Y)). since the ratio is undefined when X+Y = 0, i would like to
> assume that X+Y >= 1.
>
> is there a way to compute this expectation and see if it has a
> simplified analytic form subject to this constraint?
>
> i defined my function to compute the probability of this random
> variable as follows:
>
> myvar[x_, y_] :=
>  Binomial[n, x]*p1^x*(1 - p1)^(n - x)*Binomial[n, y]*
>   p2^y*(1 - p2)^(n - y)*(x/(x + y))
>
> now i just want to take the expectation of this. since the value of X
> ranges from 0 to N and the value of Y ranges from 0 to N, it should be
> possible to simply take a sum over X from 0 to N and sum of Y from 0
> to Y to see what this expression evaluates to and whether it has a
> simple analytic form.
>
> the trick is to take the sum with the constraint of X + Y >= 1 -- how
> can i express this in Mathematica?
>
> i am also open to computing the expectation of X/(X+Y) by assuming
> they are independent Poissons if that makes the math easier, i.e.
> using poisson approximation to binomial.
>
> thank you for your help.
>
>
>
>> > where X, Y are independent but p1 and p2 could be different. I am
>> > trying to compute the mean E(X/(X+Y)) but i am stuck. is there an easy
>> > form for this mean? i looked it up in several probability books but
>> > could not get it.
>>
>> > i am willing to make either a Poisson approximation or a Normal
>> > approximation to thebinomialfor my application if i could compute
>> > this mean. if i assume that X and Y are distributed Poisson, i know
>> > that E(X+Y) has a nice form as a function of the rates of the two
>> > poissons, but i still do not know how to get the distribution of X/(X
>> > +Y) so that i can compute the expectation E(X/(X+Y)).

The business about starting at 0 or 1 is probably not relevant. I'm not
able to get Sum to do the thing analytically.

Sum[j*PDF[BinomialDistribution[n, p1], j]*
  PDF[BinomialDistribution[n, p2], k]/(j + k), {j, n}, {k, n},
 Assumptions -> {0 < p1 < 1, 0 < p2 < 1}]

Sum[((1 - p1)^(-j + n) p1^j (1 - p2)^(-k + n) p2^
  k Binomial[n, j] Binomial[n, k])/(j + k), {j, n}, {k, n},
 Assumptions -> {0 < p1 < 1, 0 < p2 < 1}]

Nor could I get an analytic approximation using Integrate. But I think the
alteration below gives the Poisson approximation you have in mind.

In[111]:= approx =
 Sum[j*PDF[PoissonDistribution[mu1], j]*
   PDF[PoissonDistribution[mu2], k]/(j + k), {j, Infinity}, {k,
   Infinity}]
Out[111]= (E^(-mu1 -
  mu2) (-E^mu1 mu1 + E^(mu1 + mu2) mu1 + mu2 - E^mu1 mu2))/(mu1 + mu2)

In[119]:= approx2 = Simplify[approx /. {mu1 -> n*p1, mu2 -> n*p2}]
Out[119]= (E^(-n (p1 + p2)) (E^(n (p1 + p2)) p1 + p2 -
   E^(n p1) (p1 + p2)))/(p1 + p2)

A quick sanity check is to plug in explicit values for the parameters and
compare to explicit summing of the Binomial case.

In[127]:= eval1 = approx2 /. {n -> 14, p1 -> 1/2, p2 -> 1/3}
Out[127]= (6 (1/3 - (5 E^7)/6 + E^(35/3)/2))/(5 E^(35/3))

In[128]:= eval2 =
 With[{n = 14, p1 = 1/2, p2 = 1/3},
  Sum[j*PDF[BinomialDistribution[n, p1], j]*
    PDF[BinomialDistribution[n, p2], k]/(j + k), {j, n}, {k, n}]]
Out[128]= 1886578214434897919/3143703825373593600

In[129]:= N[{eval1, eval2}]
Out[129]= {0.5906, 0.600113}

With p1=1/10 and p2=4/5 I get for the numeric values
{0.1111, 0.105788} which again seem respectably close to one another.

Daniel Lichtblau
Wolfram Research




  • Prev by Date: Re: .NET/LINK USER GUIDE Comment: Calling
  • Next by Date: Re: ItemSize and ImageSize->Full , I do not understand what happens
  • Previous by thread: computing expectations for probability distributions in mathematica
  • Next by thread: Re: computing expectations for probability distributions in mathematica