2
$\begingroup$

I would like to draw random numbers in R from a conditional Weibull distribution which is given by P(x<=t+T| x>T) = 1-(1-F(t+T))/(1-F(T)), where F is the Weibull-CDF and T is the time that has elapsed without event. I know the function rweibull for drawing random numbers, but is there a similar function for the conditional distribution?

$\endgroup$
  • $\begingroup$ There is the truncdist package which implements the inversion method mentioned by @jcken $\endgroup$ – Jarle Tufto Jun 3 at 8:42
  • $\begingroup$ @Jarle Tufto, thanks. Will look into this. $\endgroup$ – otwtm Jun 3 at 12:30
3
$\begingroup$

Assume $ T, \lambda, \kappa $ are fixed. We can use the Probability Integral Transform (PIT). In essense, the PIT says that the CDF of a random variable follows a $U(0,1)$ distribution, so drawing $U(0,1)$ samples and running them through $F^{-1}(x)$ gives draws of $X$. If we can draw samples of $x$, we can then estimate the $P(x < t + T | X > T)$ via monte carlo. Nb I use the form of the Weibull that is used on wikipedia [ https://en.wikipedia.org/wiki/Weibull_distribution ].

To apply the inversion method, we want to sample values of $x = t + T$.

We have \begin{align} p = P(x < t+T | X>T) &= 1 - \frac{1 - F(t+T)}{1 - F(T)} \\ & = 1 - \frac{1 - \left( 1 - \exp\left\{ -\lambda^{-\kappa}(t+T)^{\kappa} \right\} \right)}{1 - \left( 1 - \exp\left\{ -\lambda^{-\kappa}(T)^{\kappa} \right\} \right)}\\ & = 1 - \exp \left\{ \lambda^{-\kappa}T^\kappa - \lambda^{-\kappa}(t + T)^\kappa \right\} \end{align}

Some algebra leads to $t + T = \left\{ T^\kappa - \lambda^\kappa \log(1-p) \right\}^{1/\kappa}$

The following R code can then be used to draw samples of $X$.

draw.x <- function(n, lambda, kappa, T0){
  p <- runif(n)

  ( T0^kappa - lambda^kappa *log(1-p))^(1/kappa) 

}

## estimate P(X < 5+1 = 6 | X > 5)
n.draws <- 10^4

samples <- draw.x(n.draws, 1.5, 1, 5)

prob.estimate <- sum(samples < 6)/n.draws
prob.estimate
## if kappa = 1=> exponential distribution rate = 1/lambda
## compare to exponential 1/lambda
pexp(6-5, 1/1.5); prob.estimate


hist(samples); abline(v = 6,col = 2, lwd = 2)

```
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Amazing, thank you! I was aware of the Inverse Method (what you call PIT), but didn't realize that we can manually derive the inverse in this case. Helped me a lot. $\endgroup$ – otwtm Jun 3 at 12:32
0
$\begingroup$

In Mathematica you can do it easily,

 f[a_, b_, T_] = 
  ProbabilityDistribution[
   1 - (1 - CDF[WeibullDistribution[a, b], t + T])/
    CDF[WeibullDistribution[a, b], T], {t, T, \[Infinity]}];
data = RandomVariate[f[3.5, 1, 2], 25]

{2.32932, 2.63388, 2.29273, 3.60003, 3.64601, 3.27736, 2.89225, \ 2.9748, 3.18756, 2.49238, 2.24663, 3.58997, 3.282, 3.28332, 2.54555, \ 3.33817, 2.6722, 2.94851, 3.07248, 2.23742, 2.49302, 3.13033, \ 3.23821, 3.47035, 3.50443}

I am not sure about t and T, you can adjust them yourself.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Good to know, thanks for the answer! I was just looking for an R solution in this case. $\endgroup$ – otwtm Jun 3 at 12:33

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.