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.

I am trying to evaluate the expression

Expectation[(B x (-6 b + 12 a n0 + B x + B^2 x^2))/(2 (6 a + B x + B^2 x^2)), 
  Distributed[x, BetaDistribution[α, β]]]

The evaluation runs for a few seconds and then Mathematica displays the unevaluated expression without any error messages or any other output. Is this because this random variable transformation is too complex for Mathematica to handle? If so, shouldn't some message be generated?

Edit

The same expression with an UniformDistribution

Expectation[(B x (-6 b + 12 a n0 + B x + B^2 x^2))/(2 (6 a + B x + B^2 x^2)), 
  Distributed[x, UniformDistribution[{α, β}]]]

completes and returns a formula.

share|improve this question
    
I think you will have to tell us how all the free variables in your Expectation expression are defined before we can help you. –  m_goldberg Jul 17 '13 at 7:31
    
All they are positive numbers. I have no more constrains, I am trying to find a general formula for the expectation with a,b,B and n0 as parameters. –  papirrin Jul 17 '13 at 7:42
    
@m_goldberg That would only be necessary if it concerned NExpectation. –  Sjoerd C. de Vries Jul 17 '13 at 9:37
1  
Returning unevaluated when no result could be obtained is a standard procedure in Mathematica . Try for instance Integrate[Sin[Sin[x]], x]. I believe it is in fact related to the core of Mathematica as a pattern matcher / rule user. It returns what it has got when it doesn't have any remaining rule to apply. –  Sjoerd C. de Vries Jul 17 '13 at 9:45
    
Note that Expectation effortlessly finds the expectation of a polynomial in x with non-numeric coefficients. I guess your expression is just a bit to complicated. –  Sjoerd C. de Vries Jul 17 '13 at 10:10
show 3 more comments

1 Answer

up vote 5 down vote accepted

Let's write out the expectation as an integral. We first write the rational function in a more canonical form: $$ R(x) = \frac{B}{2} x + \frac{3 B x}{B^2 x^2+B x+ 6 a}\left(2 a n_0 -a-b\right) $$ Clearly, then: $$ \mathbb{E}\left(R(X)\right) = \frac{B}{2} \frac{\alpha}{\alpha+\beta} + B \left( n_0 -\frac{a+b}{2 a}\right) \mathbb{E}\left( \frac{6 a X}{B^2 X^2+B X+ 6 a}\right) $$ The remaining rational function can be decomposed into a series using Chebyshev's polynomials: $$ \begin{align} \frac{6 a x}{B^2 x^2+B x+ 6 a} &= x \frac{1}{1 + 2 \sqrt{\frac{1}{24 a}} \left(\frac{B}{\sqrt{6a}} x\right) + \left( \frac{B}{\sqrt{6a}} x\right)^2 } \\ &= x \sum_{n=0}^\infty U_n\left( -\sqrt{\frac{1}{24 a}} \right) \left( \frac{B}{\sqrt{6a}} x\right)^n \end{align} $$ where $U_n(x)$ denotes Chebyshev polynomials of the second kind. Hence, in some open region of the parameter space, where the interchanging of the sum and the expectation is warranted, we have: $$ \mathbb{E}\left(R(X)\right) = \frac{B}{2} \frac{\alpha}{\alpha+\beta} + B \left( n_0 -\frac{a+b}{2 a}\right) \sum_{n=0}^\infty U_n\left(-\sqrt{\frac{1}{24 a}} \right) \left( \frac{B}{\sqrt{6a}}\right)^n \frac{ \left(\alpha\right)_{n+1}}{\left(\alpha+\beta\right)_{n+1}} $$ where $\left(a\right)_n$ denotes Pochhammer's symbol, as comes from $$ \mathbb{E}\left(X^{n+1}\right) = \frac{ \left(\alpha\right)_{n+1}}{\left(\alpha+\beta\right)_{n+1}} $$

Here is a numerical confirmation in Mathematica:

In[44]:= Block[{al = 2, be = 3, B = 1/3, a = 1, n0 = 1, b = 2/3}, 
 NExpectation[(B x (-6 b + 12 a n0 + B x + B^2 x^2))/(2 (6 a + B x + 
       B^2 x^2)), Distributed[x, BetaDistribution[al, be]]]]

Out[44]= 0.0881815

In[43]:= Block[{al = 2, be = 3, B = 1/3, a = 1`, n0 = 1, b = 2/3, 
   nmax = 500}, 
  B/2 al/(al + be) + 
   B (n0 - (a + b)/(2 a)) Sum[
     ChebyshevU[n, -Sqrt[(1/(24 a))]] (B/Sqrt[6 a])^
      n Pochhammer[al, n + 1]/Pochhammer[al + be, n + 1], {n, 0, 
      nmax}]] // N

Out[43]= 0.0881815

It is unlikely that this particular sum admits a closed form evaluation.

EDIT
It was pointed out to me that Mathematica can actually evaluate the sum:

In[9]:= Sum[
 ChebyshevU[n, -Sqrt[1/(24*a)]]*(B/Sqrt[6*a])^
   n*(Pochhammer[α, n + 1]/
        Pochhammer[α + β, n + 1]), {n, 0, Infinity},
   Assumptions -> α > 0 && β > 0 && a > 0 && B > 0]

Out[9]= (12*I*Sqrt[6]*a*Gamma[1 + α]*Gamma[α + β]*
    Hypergeometric2F1[1, 1 + α, 1 + α + β,
          -((2*B)/(Sqrt[a]*(Sqrt[1/a] + I*Sqrt[(-1 + 24*a)/a])))] -
      I*Sqrt[6]*Gamma[1 + α]*Gamma[α + β]*
    Hypergeometric2F1[1, 1 + α, 1 + α + β,
          -(((Sqrt[1/a] + I*Sqrt[(-1 + 24*a)/a])*B)/(12*Sqrt[a]))] +
      12*I*Sqrt[6]*a*Gamma[1 + α]*Gamma[α + β]*
    Hypergeometric2F1[1, 1 + α, 1 + α + β,
          -(((Sqrt[1/a] + I*Sqrt[(-1 + 24*a)/a])*B)/(12*Sqrt[a]))] +
      (Sqrt[6]*Sqrt[(-1 + 24*a)/a]*Gamma[1 + α]*
      Gamma[α + β]*
      Hypergeometric2F1[1, 1 + α,
       1 + α + β, -(((Sqrt[1/a] + I*Sqrt[(-1 + 24*a)/a])*
            B)/
                  (12*Sqrt[a]))])/Sqrt[1/a])/(Sqrt[
    12 - Sqrt[6]*Sqrt[1/a]]*
      Sqrt[12 + Sqrt[6]*Sqrt[1/a]]*
   a*(Sqrt[1/a] + I*Sqrt[(-1 + 24*a)/a])*Gamma[α]*
      Gamma[1 + α + β])
share|improve this answer
    
Lucid analysis indeed (+1). Thank you. –  Stefan Jul 18 '13 at 15:28
add comment

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.