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.

Apologies for stating my problem poorly in the first instance, thank you for the help in narrowing down the issue.

Problem

For those interested, it's for the exact unbiased inverse of the generalized Anscombe Transform, as given in this paper: http://dx.doi.org/10.1109/TIP.2012.2202675.

I'm trying to integrate, for a given value of $y$ and $\sigma$, the function

$$ \int_{-\infty}^{+\infty} f_{\sigma}\left ( z \right ) \sum_{k=0}^{\infty}\left ( \frac{y^{k}\mathrm{e}^{-y}}{k!\sqrt{2\pi \sigma^{2}}} \mathrm{e}^{-\frac{\left ( z-k \right )^{2}}{2\sigma^{2}}} \right ) \mathrm{d}z $$

where

$$ f_{\sigma}\left ( z \right ) = \left\{\begin{matrix} 2\sqrt{z+\frac{3}{8}+\sigma^{2}}, & z > -\frac{3}{8}-\sigma^{2}\\ 0, & z \leq -\frac{3}{8}-\sigma^{2} \end{matrix}\right. $$

which gives the equation in the reference below:

$$ \int_{-\infty}^{+\infty} 2\sqrt{z+\frac{3}{8}+\sigma ^{2}}\sum_{k=0}^{\infty}\left ( \frac{y^{k}\mathrm{e}^{-y}}{k!\sqrt{2\pi \sigma^{2}}} \mathrm{e}^{-\frac{\left ( z-k \right )^{2}}{2\sigma^{2}}} \right ) \mathrm{d}z $$

In the reference, the values of the above function are tabulated, presumably in Matlab, for $\sigma \in \left \{ 0.01,...,50 \right \}$ and $y \in \left \{0,...,200 \right \}$. The authors have made this table available for download, but I want to do it myself.

However, I can't get it to work - Mathematica won't evaluate it:

func[z_, y_, sig_] := 
 Piecewise[{{2 Sqrt[z + 3/8 + sig^2],z > -3/8 - sig^2}},0.]*
  Sum[((y^k )* Exp[-y])/(k! * Sqrt[2 \[Pi] sig^2]) * 
    Exp[-((z - k)^2) / (2 sig^2)], {k, 0, Infinity}]
NIntegrate[func[z, 5, 2], {z, -Infinity, Infinity}]

Have I made a mistake somewhere? Or misunderstood the integral/problem as given?

Update

A suggestion is to replace the infinite sum with an upper limit of ~50 as it converges quickly. This does allow Mathematica to now evaluate the integral.

share|improve this question
    
Is this supposed to be real? –  acl Jun 27 at 13:23
    
Looking at the table of values provided by the authors of the reference, yes. –  blochwave Jun 27 at 13:23
1  
But what happens for values $z<-(3/8+\sigma^2)$? –  acl Jun 27 at 13:28
1  
Anyway try uplim = 50; arg = N[2 Sqrt[z + 3/8 + 2^2]* Total@Table[((5^k)*Exp[-5])/(k!*Sqrt[2 \[Pi] 2^2])* Exp[-((z - k)^2)/(2 2^2)], {k, 0, uplim}] ]; Plot[arg, {z, -10, 20}] NIntegrate[arg, {z, -\[Infinity], \[Infinity]}]. This constructs the series by hand. –  acl Jun 27 at 13:29
1  
Then add this information to the question! OK, so this is answered then? –  acl Jun 27 at 13:31

2 Answers 2

up vote 1 down vote accepted

It cannot be real, since zunder the radical goes to minus infinity. Anyway, if you only need a numerical table, why do not you do something like this:

    f[y_, sig_, n_] := 
 NIntegrate[
  2 Sqrt[z + 3/8 + sig^2]*
   Sum[((y^k)*Exp[-y])/(k!*Sqrt[2 \[Pi] sig^2])*
     Exp[-((z - k)^2)/(2 sig^2)], {k, 0, 
     n}], {z, -\[Infinity], \[Infinity]}]

Here nis the number of terms in the sum. It converges rather fast. To check it let us do the following:

 f[5, 2, #] & /@ {10, 100, 1000}

(*   {5.92884 + 0.000418281 I, 6.03812 + 0.000418544 I, 
 6.03812 + 0.000418544 I}   *)

The result is complex, which is to be expected, and not related to the number of terms in the sum.

share|improve this answer
    
Thanks for your answer - as @acl has pointed out I managed to miss of a somewhat crucial piece of information regarding $z$ (now added), for which I must apologise. –  blochwave Jun 27 at 13:36

Assuming that if the integrand is complex we should take it to be zero, you can do this:

uplim = 50;
arg = If[Im@# > 0, 0, #] &@
   N[2 Sqrt[z + 3/8 + 2^2]*
     Total@Table[((5^k)*Exp[-5])/(k!*Sqrt[2 \[Pi] 2^2])*
        Exp[-((z - k)^2)/(2 2^2)], {k, 0, uplim}]
    ];
Plot[arg, {z, -10, 20}]
NIntegrate[arg, {z, -\[Infinity], \[Infinity]}]

This constructs the series by hand, up to some upper value (I worked out that 50 is enough by looking at the plot

share|improve this answer

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.