incomplete gamma function in R (conditional mean of Weibull to the power of N)

I am trying to calculate:

$$E(w^n | \underline{w} < w < \bar{w})$$

where $$w$$ follows a 2 parameter Weibull distribution $$w \sim W(\lambda,k)$$

From a previous question, I know the following formula for the expected value:

$$E(w^n | \underline{w} < w < \bar{w}) = \lambda^n \frac{\gamma(n/k+1,\bar{w}^{k/n})- \gamma(n/k+1,\underline{w}^{k/n})}{exp(-\underline{w}^{k/n}\lambda^{-k}) -exp(-\bar{w}^{k/n}\lambda^{-k}) }$$

But I'm trying to simulate this in R using the gammainc function and I'm getting a strange result. Consider the following variables:

w_u<-15
w_l<-4
shape<-3
scale<-8


We know $$E(w^n)=\lambda^n\Gamma(1+\frac{n}{k})$$

For $$n=3$$

scale^n*gamma(1+n/shape) = 512

Now, to calculate the conditional expectation I do:

n<-3
upterm1<- gammainc((w_u^(shape/n))*(scale^(-shape)),(n/shape)+1)[1]
upterm2<- gammainc((w_l^(shape/n))*(scale^(-shape)),(n/shape)+1)[1]
lowterm<- exp((-w_l^{shape/n})*scale^{-shape})-exp((-w_u^{shape/n})*scale^{-shape})

expected_wind_speed_n_cond<-(scale^n)*(upterm1-upterm2)/lowterm



Which returns a value close to 8.5, which does not make much sense, given the boundaries for $w$. Note that in R, the gammainc function seem to input parameters the other way around it is typically noted.

I get similar results when using:

upterm1<-pgamma(wr^(shape/n)*(scale^(-shape)),n/shape+1)*gamma(n/shape+1)
upterm2<-pgamma(win^(shape/n)*(scale^(-shape)),n/shape+1)*gamma(n/shape+1)


Any help would be appreciated.

• Abramowitz & Stegun didn't get the order of the arguments wrong, they got the upper vs lower tail wrong. That is, the first argument is the exponent on $t$ in the integral and the second is the limit of integration, both in the R package following Abramowitz & Stegun and in Wikipedia, but the sources disagree about whether you get the integral from 0 to $x$ or $x$ to $\infty$. Your code seems to have swapped the argument order but not swapped which tail you're using. Oct 5 '20 at 5:55
• @ThomasLumley: that sounds like the beginning of a good answer. Do you have the time and inclination to fill in a few details and post it? Oct 5 '20 at 8:38
• Your snippet "scale^3*shape(1+3/k) = 457.2" does not compute. Do you mean nn <- 3; scale^nn*gamma(1+nn/scale) (which in my R console yields a slightly different result, 455.12)? Oct 5 '20 at 8:39
• What is gammainc supposed to do? This function is not part of base R. Base R does supply a (standardized) incomplete gamma function: it is called pgamma.
– whuber
Oct 5 '20 at 14:50
• @StephanKolassa, my bad. fixed a typo. it's scale^n*gamma(1+n/shape) which gives 512. Oct 5 '20 at 15:18

I'll use a slightly different form for the conditional expectation of a Weibull random variable: $$E[W \ | \ a where the upper incomplete gamma function $$P \left( \alpha,x \right)$$ is the cumulative distribution function of a normalized gamma (scale=1) random variable defined as $$P \left( \alpha,x \right)=\frac{1}{\Gamma \left( \alpha \right) } \int_0^x t^{\alpha-1} e^{-t}dt$$

Based on whuber's demonstration in the linked question from the poster, we know that $$W^3$$ is also Weibull, with a shape parameter of $$k/3$$ and a scale parameter of $$\lambda^3$$. Using the lower and upper limiting points given above, we have

$$E[W^3 \ | \ 4^3 < W^3 < 15^3]=\frac{\lambda^3 \Gamma \left( \frac{3}{k}+1 \right) \left[ P \left( \frac{3}{k}+1,\left(\frac{15^3}{\lambda^3} \right)^{k/3} \right) - P \left( \frac{3}{k}+1,\left(\frac{4^3}{\lambda^3} \right)^{k/3} \right) \right] }{e^{-\left( 4^3 / \lambda^3 \right)^{k/3}}- \ e^{- \left( 15^3 / \lambda^3 \right)^{k/3}}} \$$

In R, using the pgamma function, we get

k <- 3
lambda <- 8
n <- 3

w_l <- 4
w_u <- 15

p_u <- pgamma((w_u^n/lambda^n)^(k/n),shape=n/k+1,scale=1)
p_l <- pgamma((w_l^n/lambda^n)^(k/n),shape=n/k+1,scale=1)

exp_l <- exp(-(w_l^n/lambda^n)^(k/n))
exp_u <- exp(-(w_u^n/lambda^n)^(k/n))


unifs <- runif(30000000)