rm(list=ls())

n<-25000

u<-runif(n)
v<-rweibull(n,0.4,9)

y1<-ifelse(u<(1+exp(-0.981+0.1*pmin(v,24)))/
             (1+exp(-0.981+0.1*24)),v,NA)

y<-y1[!is.na(y1)][1:10000]

mean(y)

length(y)

summary(y1)
mean(!is.na(y1))

hist(y,nclass=20,freq=FALSE)
t<-c(1:1000)
ft<-12*(.4/9)*(t/9)^(.4-1)*exp(-(t/9)^.4)/(1+exp(-.981+.1*pmin(t,24))) 
lines(t,ft,col=2)
#####  The normalizing constant should have been 12*.4: I interchanged months and years in the write-up. 

plot(ecdf(y))
plot(ecdf(y),xlim=c(0,24))