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.
 A: I'll use a slightly different form for the conditional expectation of a Weibull random variable: $$E[W \ | \ a<W<b]=\frac{\lambda \Gamma \left( \frac{1}{k}+1 \right) \left[ P \left( \frac{1}{k}+1,\left(\frac{b}{\lambda} \right)^k \right) - P \left( \frac{1}{k}+1,\left(\frac{a}{\lambda} \right)^k \right) \right] }{e^{-\left( a / \lambda \right)^k}-e^{- \left( b / \lambda \right)^k}} \ ,$$ where the lower 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))

answer <- (lambda^n*gamma(n/k+1)*(p_u-p_l))/(exp_l-exp_u)

answer 
[1] 570.846

Here is simulation code to confirm:
unifs <- runif(30000000)
w3 <- lambda^n*(-log(exp(-(w_l^n/lambda^n)^(n/k))*(1-unifs)+unifs*(exp(-(w_u^n/lambda^n)^(n/k))))^(n/k))

mean(w3)
[1] 570.8276

