2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Is this supposed to be real? $\endgroup$ – acl Jun 27 '14 at 13:23
  • $\begingroup$ Looking at the table of values provided by the authors of the reference, yes. $\endgroup$ – dr.blochwave Jun 27 '14 at 13:23
  • 1
    $\begingroup$ But what happens for values $z<-(3/8+\sigma^2)$? $\endgroup$ – acl Jun 27 '14 at 13:28
  • 1
    $\begingroup$ 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. $\endgroup$ – acl Jun 27 '14 at 13:29
  • 1
    $\begingroup$ Then add this information to the question! OK, so this is answered then? $\endgroup$ – acl Jun 27 '14 at 13:31
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – dr.blochwave Jun 27 '14 at 13:36
2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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