# Compute expectation numerically

I want to compute this expectation numerically with Mathematica. I could not figure out how to solve it numerically. Could you please help me?

$$E(x|p)=\sum _{i=0}^m \sum _{j=0}^m \binom{m}{i}\binom{m}{j}\frac{(j+1)}{(m+2)}p^i(1-p)^{m-i}\frac{Beta(i+j+1,2m-i-j+1)}{Beta(i+1,m-i+1)}$$

m is integer and p is probability. I want to compute this expectation numerically. Analytically, I obtained the solution.

The Mathematica format of the expression is below

Sum[
(Binomial[m, i])*(Binomial[m, j])*(Beta[i + j + 1, m + m - i - j + 1])*
((Beta[i + 1, m - i + 1])^-1)*(p^i)*((1 - p)^(m - i))*(j + 1)/(m + 2),
{i, 0, m}, {j, 0, m}]

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Apr 4 '15 at 23:25
• Please include this expression in Mathematica format, so that readers can work with it more easily. – bbgodfrey Apr 4 '15 at 23:27

I believe the symbolic evaluation is not correct but correct simplified expression can be found. I await insight from others wrt reasons. Here, as integer arguments of Beta have just changed to factorial.

func[a_, b_] := Factorial[a - 1] Factorial[b - 1]/Factorial[a + b - 1];
f[m_] := FullSimplify@
Sum[Binomial[m, i] Binomial[m, j] (j + 1) p^
i (1 - p)^(m - i) func[i + j + 1,
2 m \[Minus] i \[Minus] j +
1]/((m + 2) func[i + 1, m \[Minus] i + 1]), {i, 0, m}, {j, 0,
m}]
FindSequenceFunction[(f[#] /. p -> u) & /@ Range[10], x]


Using the found function:

sf[x_, u_] := (2 + 2 x + u x^2)/(2 + x)^2


Note: sf[14,6/10] yields 369/640 (0.576563).

Testing (not proof):

Row[Grid[#, Frame -> All] & /@
Partition[Table[{j, f[j], Simplify[sf[j, p]]}, {j, 1, 30}], 10]]


and for fun:

  Manipulate[
Plot[{Evaluate[sf[n, p]], p}, {p, 0, 1}, Frame -> True,
FrameLabel -> {"p", "E[x|p]"}, PlotLegends -> "Expressions"], {n, 2,
100, 1}]


• Thanks a lot for your help,Its really great work! – Jimmy Dur Apr 5 '15 at 7:09
• Nice work as always, plus to you. – ciao Apr 5 '15 at 8:38

This?

FullSimplify[Sum[
Binomial[m, i] Binomial[m, j] (j + 1)/(m + 2) p^i (1 - p)^(m - i)
Beta[i+j+1, 2m-i-j+1]/Beta[i+1,m-i+1], {i, 0, m}, {j, 0, m}]]

(* gives (2 (1+m)(1-p)^m)/(2+m)^2 *)


BUT your latex doesn't match your Mathematica and I don't know which to trust.

• Thanks for correction.I made a typo.and what you found is analytical solution I need to solve it numerically. – Jimmy Dur Apr 5 '15 at 0:16
• Your analytical solution is correct I made a mistake .I have already did some computation as you said for fixing m and p.but I am not sure its the numerical solution of this expectation.I was thinking if I can use some numerical methods and get a result after some iterations. – Jimmy Dur Apr 5 '15 at 0:40
• another thing is that when we do numerical solution as you wrote down, for p=0.6; m=14; Sum[Binomial[m, i] Binomial[m, j] (j+1)/(m+2) p^i (1-p)^(m-i) Beta[i+j+1, 2m-i-j+1]/Beta[i+1, m-i+1], {i, 0, m}, {j, 0, m}]] which gives us 0.576563 .But if I plug m=14 and p=0.6 in the analytical solution (2 (1+m)(1-p)^m)/(2+m)^2 , it gives 3.14573*10^-7. I dont know why but these two solutions do not agree.Thats kind of weird. – Jimmy Dur Apr 5 '15 at 0:51
• That is odd. My first guess is round-off error. So change 0.6 to 6/10, restart and try again. Still not the same. You might dig into this and see if you can figure this out. You might learn something from doing that. You might have found an "unexpected behavior" in MMA, but see if you can find an explanation before jumping to any conclusion. – Bill Apr 5 '15 at 1:53