Stack Exchange
Numerical integration in R. unirroot Error

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

  roots = uniroot.all(G,c(0,1))
  probExtinction = roots[2]
  integrand = probExtinction*exp((-9)*lambda)*lambda^(14)

  return(integrand)
}

ExpectedProb = integrate(fun,0,Inf)

The error I get when I run that code is:
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 
G = function(s){exp(lambda*(s-1))-s}
roots = uniroot.all(G,c(0,1))
probExtinction = roots[2])

Anyone knows how I can fix this?

Very annoying that this post was migrated here because now the LateX formats are all gone… How hard is it to include LateX support in this forum as well? Pretty important if you ask me, since math and programming are quite interwoven. Any suggestion on how I should make the math part more readable?

Yes it is. You can use this in the meantime to use latex https://codecogs.com/latex/eqneditor.php

How does that help me to display LaTex in this site?

An anchor tag is created at the bottom of the site you can use

What library is uniroot.all from? Is there supposed to be an extra -s in your Gs function? I'm not seeing it in the equations you gave.

@mickey the package is rootSolve I used. Yes there is, I understand why it confused you, I changed the name of the function Gs to just G. And G = exp(lambda*(s-1))-s. Apologies for the confusion but the code is now correct as it is.

I kept getting a similar issue as you. The problem is somewhere with the integrate function (which I have little experience using). Instead, I just sampled a bunch of gamma variables and evaluated that s(lambda) function at each one. Take the mean and presto:
library(rootSolve)

a = 15
b = 9 
n = 10000
lambda = rgamma(n, a, b)    # lambda ~ Gamma(15, 9)

# 's' is s(lambda)
s = sapply(lambda, function(l){
    G = function(s) exp(l*(s-1))-s  # solve G(s) = s
    min(uniroot.all(G, c(0, 1)))    # smallest root
    })  

mean(s)
# plot(cumsum(s) / 1:n)

The mean would be your estimate of the integral (not unlike doing mean(exp(rnorm(1000))) to estimate the mean of exponentiated normals). You can plot the cumulative mean to see the convergence.
