You can use uniroot(...)
for this.
[Note: If the point of this exercise is to implement your own version of a Newton Raphson technique, let me know and I'll delete the answer.]
If I'm understanding this correctly, you want to generate random samples from a distribution with probability density function f
and cumulative density F
where
f = x*exp(-x)
F = 1 - (1+x)*exp(-x)
As you imply, this can be done by generating a random sample from U[0,1]
and transforming that according to the inverse CDF of F
. The procedure is very similar to the ones posted here and here, except that you already have an expression for the CDF.
f <- function(x) x*exp(-x)
F <- function(x) 1-(1+x)*exp(-x)
F.inv <- function(y){uniroot(function(x){F(x)-y},interval=c(0,100))$root}
F.inv <- Vectorize(F.inv)
x <- seq(0,10,length.out=1000)
y <- seq(0,1,length.out=1000)
par(mfrow=c(1,3))
plot(x,f(x),type="l",main="f(x)")
plot(x,F(x),type="l",main="CDF of f(x)")
plot(y,F.inv(y),type="l",main="Inverse CDF of f(x)")
Then, generate X ~ U[0,1]
and Z = F.inv(X)
.
set.seed(1)
X <- runif(1000,0,1) # random sample from U[0,1]
Z <- F.inv(X)
par(mfrow=c(1,1))
hist(Z, freq=FALSE, breaks=c(seq(0,10,length=30),Inf), xlim=c(0,10))
lines(x,f(x),type="l",main="Density function", col="red",lty=2)