0
$\begingroup$

In a paper[1] on generating survival data from a Cox proportional hazard model, they specify following cumulative distribution function:

\begin{equation} F(t^0) = \frac{1}{5}\sum_{x=0}^{4} \exp(-\exp(\beta x) t^0). \end{equation}

with $\beta = 0.5\log(2)$

They generate random numbers from this distribution but don't specify how. I don't see a way to invert the distribution and inverse transform sampling.
How would you sample from this kind of distribution?

[1]: Mackenzie, Todd, and Michal Abrahamowicz. "Marginal and hazard ratio specific random data generation: applications to semi-parametric bootstrapping." Statistics and Computing 12.3 (2002): 245-252.

$\endgroup$
  • 2
    $\begingroup$ This is a mixture distribution. An expression like this could, for instance, arise as a kernel density estimate using an Exponential kernel (even though this is not a KDE: it's a scale mixture rather than a location mixture). In this fashion your question is completely answered at stats.stackexchange.com/questions/321542. $\endgroup$ – whuber Jan 16 '18 at 16:52
2
$\begingroup$

You can think of this as a mixture of exponential distributions with five components. Each mixture component has a probability of $1/5$, and the parameters of the individual components are $\exp(\beta x), x = 0, \dots, 4$.

Sampling from this distribution can be done in two steps, sample code in R. First, randomly select the mixture component to sample from:

x <- sample(5,1) - 1

and second, generate the r.v. from the appropriate component:

t0 <- rexp(1, exp(beta*x))

which can be made much more efficient:

beta <- 0.5*log(2)
N <- 1000
t0 <- rexp(N, exp(beta*(sample(5, N, replace=TRUE)-1)))

We can compare the results to calculated values of the density function as follows:

x <- seq(0.1,7,0.1)
pdf <- rep(0, length(x))
for (i in 1:length(x)) {
  pdf[i] <- mean(dexp(x[i], exp(beta*(0:4))))
}

hist(t0, freq=FALSE)
lines(pdf~x, lwd=2)

which generates the following plot:

enter image description here

which in turn indicates that we likely haven't messed up anywhere.

$\endgroup$

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.