Nothing
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
## ----events-------------------------------------------------------------------
str(eventDates <- c(as.Date(
c('1962-10-16', '1983-09-26')), Sys.Date()))
(daysBetween <- difftime(tail(eventDates, -1),
head(eventDates, -1), units='days'))
(yearsBetween <- as.numeric(daysBetween)/365.24)
names(yearsBetween) <- c('observed', 'censored')
str(yearsBetween)
## ----lik----------------------------------------------------------------------
Lik <- function(lambda, Times=yearsBetween){
Lik <- (exp(-sum(Times)/lambda) /
(lambda^(length(Times)-1)))
Lik
}
## ----logLik-------------------------------------------------------------------
logLk <- function(lambda, Times=yearsBetween){
logL <- (-sum(Times)/lambda -
(length(Times)-1)*log(lambda))
logL
}
## ----mle----------------------------------------------------------------------
(lambdaHat <- (sum(yearsBetween) /
(length(yearsBetween)-1)))
## ----chisq--------------------------------------------------------------------
#(chisq2 <- qchisq(
# c(.8, .95, .99, .999, 1-1e-6), 1))
(chisq2 <- qchisq(c(.2, .05, .01, .001, 1e-6),
1, lower.tail=FALSE))
## ----CI-----------------------------------------------------------------------
lambda <- lambdaHat+seq(-50, 1000, 50)
(logLR2 <- 2*(logLk(lambdaHat) - logLk(lambda)))
## ----CIlog--------------------------------------------------------------------
l_lam <- log(lambdaHat)+seq(-2, 6, 1)
(logLR2_u <- 2*(logLk(lambdaHat) - logLk(exp(l_lam))))
## ----CItheta------------------------------------------------------------------
theta <- (1/lambdaHat + seq(-.016, .06, .004))
(logLR2_th <- 2*(logLk(lambdaHat) - logLk(1/theta)))
## ----plot_u-------------------------------------------------------------------
makePlots <- FALSE
library(grDevices)
outType <- ''
#outType = 'png'
switch(outType,
svg=svg('yrs2Armageddon.svg'),
# need png(..., 960, 960), because the default 480
# is not sufficiently clear to easily read the labels
png=png('yrs2Armageddon.png', 960, 960)
)
op <- par(mar=c(6, 4, 4, 2)+.1)
# Experiment with the range of "seq" here until
# head and tail of logLR2_u. are just over 6.63:
u. <- log(lambdaHat)+seq(-1.86, 4.36, .02)
lam. <- exp(u.)
logLR2_u. <- 2*(logLk(lambdaHat) - logLk(lam.))
head(logLR2_u., 1)
tail(logLR2_u., 1)
if(makePlots){
plot(lam., logLR2_u., type='l', bty='n', log='x',
xlab='', ylab='', las=1, axes=FALSE, lwd=2)
axis(1, padj=-1)
axis(2, las=1)
# xlab = \lambda:
# Greek letters did not render in GIMP 2.10.8 on 2018-12-30,
# so don't use svg until this is fixed.
switch(outType,
# cex doesn't work properly with svg > GIMP
# Therefore, I can NOT use svg
svg={cex2 <- 2; mtext('lambda', 1, 1.6, cex=cex2)},
png={cex2 <- 2; mtext(expression(lambda), 1,
1.6, cex=cex2)},
{cex2 <- 1.3; mtext(expression(lambda), 1,
1.6, cex=cex2)}
)
lamTicks <- axTicks(1)
thTicks <- 1/lamTicks
axis(1, lamTicks, thTicks, line=3, padj=-1)
switch(outType,
svg=mtext('theta == 1/lambda', 1, 4.9, cex=cex2),
mtext(expression(theta == 1/lambda), 1, 4.9, cex=cex2)
)
abline(h=chisq2, col='red', lty=c('dotted', 'dashed'),
lwd=2)
(CI.8 <- range(lam.[logLR2_u. <= chisq2[1]]))
text(lambdaHat, chisq2[1],
paste0('80% CI =\n(',
paste(round(CI.8), collapse=', '), ')'),
cex=cex2)
(CI.95 <- range(lam.[logLR2_u. <= chisq2[2]]))
text(lambdaHat, chisq2[2],
paste0('95% CI =\n(',
paste(round(CI.95), collapse=', '), ')'),
cex=cex2)
abline(v=CI.8, col='red', lty='dotted', lwd=2)
abline(v=CI.95, col='red', lty='dashed', lwd=2)
(CI.99 <- range(lam.[logLR2_u. <= chisq2[3]]))
text(lambdaHat, chisq2[3],
paste0('99% CI =\n(',
paste(round(CI.99), collapse=', '), ')'),
cex=cex2)
abline(v=CI.8, col='red', lty='dotted', lwd=2)
abline(v=CI.95, col='red', lty='dashed', lwd=2)
abline(v=CI.99, col='red', lty='dashed', lwd=2)
if(outType != '')dev.off()
par(op)
}
## ----plot_lin-----------------------------------------------------------------
# copy the code from the last snippet
# and delete "log='x'", then adjust the placement
# of CI text
switch(outType,
svg=svg('yrs2Armageddon_lin.svg'),
# need png(..., 960, 960), because the default 480
# is not sufficiently clear to easily read the labels
png=png('yrs2Armageddon_lin.png', 960, 960)
)
op <- par(mar=c(6, 4, 4, 2)+.1)
u. <- log(lambdaHat)+seq(-1.86, 4.36, .02)
lam. <- exp(u.)
logLR2_u. <- 2*(logLk(lambdaHat) - logLk(lam.))
head(logLR2_u., 1)
tail(logLR2_u., 1)
if(makePlots){
plot(lam., logLR2_u., type='l', bty='n',
xlab='', ylab='', las=1, axes=FALSE, lwd=2)
axis(1, padj=-1)
axis(2, las=1)
# xlab = \lambda:
# Greek letters did not render in GIMP 2.10.8 on 2018-12-30,
# so don't use svg until this is fixed.
switch(outType,
# cex doesn't work properly with svg > GIMP
# Therefore, I can NOT use svg
svg={cex2 <- 2; mtext('lambda', 1, 1.6, cex=cex2)},
png={cex2 <- 2; mtext(expression(lambda), 1,
1.6, cex=cex2)},
{cex2 <- 1.3; mtext(expression(lambda), 1,
1.6, cex=cex2)}
)
lamTicks <- axTicks(1)
thTicks <- 1/lamTicks
axis(1, lamTicks, thTicks, line=3, padj=-1)
switch(outType,
svg=mtext('theta == 1/lambda', 1, 4.9, cex=cex2),
mtext(expression(theta == 1/lambda), 1, 4.9, cex=cex2)
)
abline(h=chisq2, col='red', lty=c('dotted', 'dashed'),
lwd=2)
(CI.8 <- range(lam.[logLR2_u. <= chisq2[1]]))
#text(lambdaHat, chisq2[1],
text(400, chisq2[1],
paste0('80% CI =\n(',
paste(round(CI.8), collapse=', '), ')'),
cex=cex2)
(CI.95 <- range(lam.[logLR2_u. <= chisq2[2]]))
#text(lambdaHat, chisq2[2],
text(800, chisq2[2],
paste0('95% CI =\n(',
paste(round(CI.95), collapse=', '), ')'),
cex=cex2)
abline(v=CI.8, col='red', lty='dotted', lwd=2)
abline(v=CI.95, col='red', lty='dashed', lwd=2)
(CI.99 <- range(lam.[logLR2_u. <= chisq2[3]]))
#text(lambdaHat, chisq2[3],
text(3000, chisq2[3],
paste0('99% CI =\n(',
paste(round(CI.99), collapse=', '), ')'),
cex=cex2)
abline(v=CI.8, col='red', lty='dotted', lwd=2)
abline(v=CI.95, col='red', lty='dashed', lwd=2)
abline(v=CI.99, col='red', lty='dashed', lwd=2)
if(outType != '')dev.off()
}
par(op)
## ----plot_inverse-------------------------------------------------------------
switch(outType,
svg=svg('yrs2Armageddon_inverse.svg'),
png=png('yrs2Armageddon_inverse.png', 960, 960)
)
op <- par(mar=c(6, 4, 4, 2)+.1)
# This will require more changes than just deleting log='x':
u. <- log(lambdaHat)+seq(-1.86, 4.36, .02)
lam. <- exp(u.)
logLR2_u. <- 2*(logLk(lambdaHat) - logLk(lam.))
head(logLR2_u., 1)
tail(logLR2_u., 1)
if(makePlots){
plot(-1/lam., logLR2_u., type='l', bty='n',
xlab='', ylab='', las=1, axes=FALSE, lwd=2)
thTicks <- (-axTicks(1))
axis(1, -thTicks, abs(1/thTicks), padj=-1)
axis(2, las=1)
# xlab = \lambda:
# Greek letters did not render in GIMP 2.10.8 on 2018-12-30,
# so don't use svg until this is fixed.
switch(outType,
# cex doesn't work properly with svg > GIMP
# Therefore, I can NOT use svg
svg={cex2 <- 2; mtext('lambda', 1, 1.6, cex=cex2)},
png={cex2 <- 2; mtext(expression(lambda), 1,
1.6, cex=cex2)},
{cex2 <- 1.3; mtext(expression(lambda), 1,
1.6, cex=cex2)}
)
axis(1, -thTicks, thTicks, line=3, padj=-1)
switch(outType,
svg=mtext('theta == 1/lambda', 1, 4.9, cex=cex2),
mtext(expression(theta == 1/lambda), 1, 4.9, cex=cex2)
)
abline(h=chisq2, col='red', lty=c('dotted', 'dashed'),
lwd=2)
(CI.8 <- range(lam.[logLR2_u. <= chisq2[1]]))
#text(lambdaHat, chisq2[1],
text(-.02, chisq2[1],
paste0('80% CI =\n(',
paste(round(CI.8), collapse=', '), ')'),
cex=cex2)
(CI.95 <- range(lam.[logLR2_u. <= chisq2[2]]))
#text(lambdaHat, chisq2[2],
text(-.04, chisq2[2],
paste0('95% CI =\n(',
paste(round(CI.95), collapse=', '), ')'),
cex=cex2)
abline(v=CI.8, col='red', lty='dotted', lwd=2)
abline(v=CI.95, col='red', lty='dashed', lwd=2)
(CI.99 <- range(lam.[logLR2_u. <= chisq2[3]]))
#text(lambdaHat, chisq2[3],
text(-.06, chisq2[3],
paste0('99% CI =\n(',
paste(round(CI.99), collapse=', '), ')'),
cex=cex2)
abline(v=-1/CI.8, col='red', lty='dotted', lwd=2)
abline(v=-1/CI.95, col='red', lty='dashed', lwd=2)
abline(v=-1/CI.99, col='red', lty='dashed', lwd=2)
if(outType != '')dev.off()
}
par(op)
## ----simExp-------------------------------------------------------------------
set.seed(1)
simExp <- rexp(1000)
if(makePlots){
qqnorm(simExp, datax=TRUE)
qqnorm(simExp, datax=TRUE, log='x')
qqnorm(1/simExp, datax=TRUE)
}
## ----rlam---------------------------------------------------------------------
library(invgamma)
set.seed(123)
rlambda2 <- function(n, sumTimes=lambdaHat){
# -sumTimes/log(runif(n))
# k <- 2; rinvgamma(n k-1, scale=sumTimes)
rinvgamma(n, 1, rate=sumTimes)
}
simLam <- rlambda2(1e4)
quantile(simLam)
mean(simLam)
## ----betaSSQ------------------------------------------------------------------
Dev2 <- function(shapes, p=c(.3, .6), q=c(.1, .9)){
devs <- (q - pbeta(p, shapes[1], shapes[2]))
sum(devs^2)
}
# test
Dev2(c(1, 1))
## ----betaSolve----------------------------------------------------------------
(betaSolve <-optim(c(1,1), Dev2,
method="L-BFGS-B", lower=c(0,0)))
## ----meanBeta-----------------------------------------------------------------
a.b <- sum(betaSolve$par)
(meanBeta <- betaSolve$par[1]/a.b)
(varBeta <- with(betaSolve, par[1]*par[2] /
(a.b^2 * (a.b+1))))
## ----mcArmageddon-------------------------------------------------------------
mcArmageddon <- function(N,
betapars=betaSolve$par,
sumTimes=lambdaHat){
# 1. Start time
start <- proc.time()
# 2. Q ~ B(\alpha, \beta)
Q <- rbeta(N, betapars[1], betapars[2])
# 3. K ~ NB(Q, 1)
K <- (1+rnbinom(N, 1, Q))
# 4. Time[i] <- sum(rlambda2(K[i]))
Time <- numeric(N)
for(i in 1:N){
Time[i] <- sum(rlambda2(K[i],
sumTimes=sumTimes))
}
attr(Time, 'Qbar') <- mean(Q)
attr(Time, 'quantileQ') <- quantile(Q)
attr(Time, 'Kbar') <- mean(K)
attr(Time, 'quantileK') <- quantile(K)
# 5. quantiles
cat('meanTime = ', mean(Time), '\n')
print(quantile(Time))
# 6. elapsed.time
et <- (proc.time()-start)
# 7. Return et as an attribute
attr(Time, 'elapsed.time') <- et
cat('et = ', et, '\n')
Time
}
set.seed(1)
(mcArm1 <- mcArmageddon(10))
## ----mcArm2-------------------------------------------------------------------
set.seed(2)
mcArm2 <- mcArmageddon(100)
attributes(mcArm2)
## ----mcArm3-------------------------------------------------------------------
set.seed(3)
mcArm3 <- mcArmageddon(1000)
attributes(mcArm3)
## ----mcArm4-------------------------------------------------------------------
set.seed(4)
mcArm4 <- mcArmageddon(1e4)
attributes(mcArm4)
## ----mcArm5-------------------------------------------------------------------
set.seed(5)
mcArm5 <- mcArmageddon(1e5)
attributes(mcArm5)
## ----million------------------------------------------------------------------
set.seed(6)
mcArm6 <- mcArmageddon(1e6)
attributes(mcArm6)
## ----MCsumStats---------------------------------------------------------------
mean(mcArm6)
mean(mcArm6<40)
mean(mcArm6<60)
quantile(mcArm6, c(1e-6, 1e-3))
## ----1e7----------------------------------------------------------------------
if(!fda::CRAN()){
# Don't run this with CRAN tests,
# because it takes too long
set.seed(7)
mcArm7 <- mcArmageddon(1e7)
print(attributes(mcArm7))
print(mean(mcArm7))
print(quantile(mcArm7, c(1e-6, 1e-3)))
}
## ----sort7--------------------------------------------------------------------
if(fda::CRAN()){
mcArm. <- mcArm6
} else mcArm. <- mcArm7
mcArm.s <- sort(mcArm.)
quantile(mcArm.s)
## ----qqnorm7------------------------------------------------------------------
str(qq7 <-as.data.frame(qqnorm(mcArm.s,
plot.it=FALSE)))
## ----indx---------------------------------------------------------------------
N. <- length(mcArm.s)
#index.5 <- c(1:10, seq(20, 100, 10),
# seq(200, 1000, 100),
# seq(3000, (N./2)-2000, 1000))
index.5 <- c(1:1000,
seq(2000, (N./2)-2000, 1000))
index <- c(index.5, N.+1-rev(index.5))
tail(index, 30)
length(index)
# yes: I think I did this right.
## ----plot7--------------------------------------------------------------------
switch(outType,
svg=svg('yrs2ArmageddonQQ.svg'),
# need png(..., 960, 960),
# because the default 480 is not sufficiently
# clear to easily read the labels
png=png('yrs2ArmageddonQQ.png', 960, 960)
)
op <- par(mar=c(5, 5, 4, 5)+.1)
if(makePlots){
with(qq7, plot(y[index], x[index], type='l',
log='x', las=1, bty='n', lwd=2,
xlab='', ylab='',
cex.lab=2, axes=FALSE) )
# xlab='years to Armageddon',
# ylab='standard normal scores',
axis(1, cex.axis=2)
axis(2, cex.axis=2, las=1)
probs <- c(.001, .01, .1, .25, .5, .75,
.9, .99, .999)
z <- qnorm(probs)
if(outType==''){
cex.txt <- 1.5
cex.ax4 <- 1.3
} else {
cex.txt <- 3
cex.ax4 <- 2
}
axis(4, z, probs, cex.axis=cex.ax4,
las=1, line=-.5)
p40 <- mean(mcArm.s<40)
p60 <- mean(mcArm.s<60)
z40.60 <- qnorm(c(p40, p60))
max7 <- tail(mcArm.s, 1)
lines(c(rep(40, 2), max7),
c(-5, rep(z40.60[1], 2)),
lty='dotted', lwd=2, col='red')
lines(c(rep(60, 2), max7),
c(-5, rep(z40.60[2], 2)),
lty='dashed', lwd=2, col='purple')
text(15, -5, '40', col='red', cex=cex.txt)
text(200, -2.5, '60', col='purple', cex=cex.txt)
text(.2*max7, z40.60[1]-.6,
paste0(round(100*p40), "%"), cex=cex.txt,
col='red')
text(.2*max7, z40.60[2]+.6,
paste0(round(100*p60), "%"), cex=cex.txt,
col='purple')
}
par(op)
if(outType != '')dev.off()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.