Compute integral for expected value in R.

Question

From a previous assignment, I found that the extinction probability of of a branching process is given by solving the equation

$$G(s)=s,$$

where $G(s)=e^{\lambda(s-1)}$ is the PGF and noting the smallest root in terms of $s$. The solution to the equation

$$e^{\lambda(s-1)}-s=0$$

will be a function of $\lambda$, namley $s(\lambda)$ which can be computed numerically for a fixed $\lambda.$

Now I'm asked to numerically find the expected value of this probability given that $\lambda\sim\Gamma(15,9)$. Let's for convention now use the function $s(\Lambda)$, where $\Lambda\sim\Gamma(15,9)$, I know that the expected value of a function of a random variable $X$ is given by

$$\mathbb{E}[g(X)]=\int_{\infty}^{\infty}g(x)\cdot\pi(x)dx,$$

where $\pi(x)$ is the pdf of $X$. So in my case I get

$$\mathbb{E}[s(\Lambda)]=\int_{0}^{\infty}s(\lambda)\pi(\lambda)d\lambda=\int_0^{\infty}s(\lambda)\frac{9^{15}}{\Gamma(15)}\lambda^{15-1}e^{-9\lambda}d\lambda.$$

Since $9^{15}/\Gamma(15)\approx 2361.7$ we have that the integral to be numerically computed in R is

$$2361.7\int_{0}^{\infty}s(\lambda)\lambda^{14}e^{-9\lambda}.$$

I figured, since we don't have a closed form expression for $s(\lambda)$, that I need nestled functions, first express the solution of $G(s)-s=0$ in one function depending on $s$, keeping $\lambda$ constant and then wrap this with another function only depending on $\lambda,$ like this:

fun = function(lambda){
    sLambda = function(s){exp(lambda*(s-1))-s}

    roots = uniroot.all(sLambda,c(0,1))
    probExtinction = roots[2]       # since this root is the smallest one.
    integrand = probExtinction*exp((-9)*lambda)*lambda^(14)

    return(integrand)
 }

    ExpectedProb = 2361.7*integrate(fun,0,Inf)$value

But here is where trouble starts. I get the error below. I understand what the error means, it's that if the function values at the endpoints are not oposite, then the function never crosses the x-axis, thus no roots will be found. However, plotting the function $f(x)=e^{a(x-1)}-x$ I see that the function dramatically changes behavior when $a$ passes $1$. I don't see how to fix my code in order to deal with this issue.

Error in uniroot(f, lower = xseq[i], upper = xseq[i + 1], ...) : 
  f() values at end points not of opposite sign
In addition: Warning messages:
1: In lambda * (s - 1) :
  longer object length is not a multiple of shorter object length
2: In if (is.na(f.lower)) stop("f.lower = f(lower) is NA") :
  the condition has length > 1 and only the first element will be used
3: In if (is.na(f.upper)) stop("f.upper = f(upper) is NA") :
  the condition has length > 1 and only the first element will be used

Whats weird is that when I only run this code I don't get the uniroot-error:

lambda = 2 
sLambda = function(s){exp(lambda*(s-1))-s}
roots = uniroot.all(sLambda,c(0,1))
probExtinction = roots[2])

Anyone knows how I can fix this?

Answer 1

The solution of $$e^{\lambda(s-1)}-s=0$$ is given in terms of Lambert function and it is $$s=-\frac{W\left(-\lambda\, e^{-\lambda } \right)}{\lambda }$$ This function is available in R (have a look here).

When plotted, the integrand looks very nice but you still need numerical integration. If I did it properly, the result is $\approx 0.4002$.