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:


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:

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})


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:


Any help would be appreciated.

  • 2
    $\begingroup$ 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. $\endgroup$ Oct 5 '20 at 5:55
  • $\begingroup$ @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? $\endgroup$ Oct 5 '20 at 8:38
  • $\begingroup$ 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)? $\endgroup$ Oct 5 '20 at 8:39
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Oct 5 '20 at 14:50
  • $\begingroup$ @StephanKolassa, my bad. fixed a typo. it's scale^n*gamma(1+n/shape) which gives 512. $\endgroup$
    – Raimundo
    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<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 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))

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

[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))

[1] 570.8276

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.