Htykuut

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 $lambdasimGamma(15,9)$. Let's for convention now use the function $s(Lambda)$, where $LambdasimGamma(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)cdotpi(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)dlambda=int_0^{infty}s(lambda)frac{9^{15}}{Gamma(15)}lambda^{15-1}e^{-9lambda}dlambda.$$

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

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

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?

The solution of
$$e^{lambda(s-1)}-s=0$$ is given in terms of Lambert function and it is
$$s=-frac{Wleft(-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$.