---
title: "Time to nuclear Armageddon"
author: "Spencer Graves"
date: "1/1/2019"
output: html_document
# On 2019-01-05 from a computer 
# running macOS 10.14.2, 
# I got printer errors 
# opening the html in a browser and 
# trying to print from there 
# or creating a PDF from there and 
# trying to print that.  
# I was successful after switching to 
# output: word_document

vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Time to nuclear Armageddon}
  %\VignetteKeyword{nuclear Armageddon}
  %\SweaveUTF8
  %\usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Introduction 

This is a companion to the [Wikiversity article on "Time to extinction of civilization"](https://en.wikiversity.org/wiki/Time_to_extinction_of_civilization).  That article assumes the [1962 Cuban Missile Crisis](https://en.wikipedia.org/wiki/Cuban_Missile_Crisis) and the [1983 Soviet nuclear false alarm incident](https://en.wikipedia.org/wiki/1983_Soviet_nuclear_false_alarm_incident) provide one observation on the time _between_ major nuclear crises, with a second time between such crises being _censored_ at the present.

What can we say about the distribution of the time between such major nuclear crises?  

With one observed time and a second censored, we can construct a [likelihood function](https://en.wikipedia.org/wiki/Likelihood_function), which we can then use to estimate the mean time between such crises and the uncertainty in that estimate.  With further estimates of the probability that such a crisis would lead to a nuclear war and nuclear winter, we can simulate such times and obtain plausible bounds on uncertainty in our estimates.  

This methodology could later be expanded to consider a larger list of nuclear crises with a broader range of probabilities for each crisis escalating to a nuclear war and winter.  The fact that no such nuclear war has occurred as of this writing puts an upper limit on such probabilities.  A rough lower limit can be estimated from comments from people like Robert McNamara and Daniel Ellsberg, both of whom have said that as long as there are large nuclear arsenals on earth, it is only a matter of time before a nuclear crises escalates to such a nuclear Armageddon.  McNamara was US Secretary of Defense during the 1962 Cuban Missile Crisis, and Ellsberg as a leading nuclear war planner advising McNamara and the rest of President Kennedy's team during that crisis.  For more on this, see the companion [Wikiversity article on "Time to extinction of civilization"](https://en.wikiversity.org/wiki/Time_to_extinction_of_civilization).  

We start by being explicit about the observed and censored times between major nuclear crises.  

## Times of major nuclear crises 

```{r 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)
```

## Likelihood of times between major nuclear crises 

[Appendix 1 of that Wikiversity  article](https://en.wikiversity.org/wiki/Time_to_extinction_of_civilization#Appendix_1._Estimating_the_mean_time_between_nuclear_crises) provides the following likelihood assuming we observe times between major nuclear crises, $T_1, ..., T_{k-1},$ plus one censoring time, $T_k$, and the times between such crises follow an [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution) that does not change over time:  

\[
L(\lambda | \mathbf{T}) = \exp[−S_k / \lambda ] / \lambda^{k-1}
\]

where $\mathbf{T}$ = the vector consisting of $T_1, ..., T_k$, and 

\[
S_k = \sum_{i=1}^k{T_i}.
\]

The exponential distribution is the simplest lifetime distribution.  It is widely used for applications like this and seems reasonable in this context.  

[For setting math in RMarkdown, we are following 
[Cosma Shalizi (2016) "Using R Markdown for Class Reports"](https://www.stat.cmu.edu/~cshalizi/rmarkdown/#math-in-r-markdown).]

We code this as follows:  

```{r lik}
Lik <- function(lambda, Times=yearsBetween){
  Lik <- (exp(-sum(Times)/lambda) / 
      (lambda^(length(Times)-1)))
  Lik
}
```

From this, we compute the log(likelihood) as follows:  

\[
l(\lambda | \mathbf{T}) = [(−S_k / \lambda) - (k-1)\log(\lambda)].
\]

We code this as follows:  

```{r logLik}
logLk <- function(lambda, Times=yearsBetween){
  logL <- (-sum(Times)/lambda - 
        (length(Times)-1)*log(lambda))
  logL
}
```

By differentiating $l$ with respect to $\lambda$ or $u = \log(\lambda)$ or $\theta = 1/\lambda$, we get a score function that is zero when $\lambda$ is $\sum T_i/(k-1)$, where the "-1" comes from assuming that only the last of Times is censored.  

The value of parameter(s) that maximize the likelihood is (are) called [maximum likelihood estimates (MLEs)](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation), and it is standard to distinguish an MLE with a circumflex (^).  We use this convention to write the following:  

\[
\hat\lambda = \sum T_i / (k-1).
\]

This is commonly read "lambda hat".  We code it as follows:  

```{r mle}
(lambdaHat <- (sum(yearsBetween) / 
                 (length(yearsBetween)-1)))
```

From [Wilks' theorem](https://en.wikipedia.org/wiki/Wilks%27_theorem), we know that 2*log(likelihood ratio) is approximately [Chi-squared](https://en.wikipedia.org/wiki/Chi-squared_distribution) with degrees of freedom equal to the number of parameters estimated, which is 1 in this case.  

```{r 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))
```

In reality, because of the questionable nature of our assumptions, we may wish to place less confidence in these numbers than what is implied by the stated confidence levels.  However, we will not change these numbers but leave it to the reader to downgrade them as seems appropriate.   

Also, in the following, we will mark the 80, 95, and 99 percent confidence intervals on the plots, leaving the more extreme tails for separate computations.  For now, we want to plot 2*log(likelihood) in a neighborhood of the MLE.  For this, we will focus on the region that is closer than `chisq2/2` of the maximum:  

```{r CI}
lambda <- lambdaHat+seq(-50, 1000, 50)
(logLR2 <- 2*(logLk(lambdaHat) - logLk(lambda)))
```

After several attempts at adjusting the parameters for the [`seq`](https://www.rdocumentation.org/packages/base/versions/3.5.2/topics/seq) function while including 0 in the sequence, I gave up trying to get a range with 0 inside and just over [`qchisq(0.99, 1)`](https://www.rdocumentation.org/packages/stats/versions/3.5.2/topics/Chisquare) = 6.63 on both ends: Obviously, I got it for $\lambda$ small.  However, it seemed infeasible to do this with only a relatively few points for $\lambda$ large: The MLE here is only the second of 22 evenly-spaced points on the $\lambda$ scale while the value for the `22nd` point is not close to the target 6.63.  Let's try $u = \log(\lambda)$:  

```{r CIlog}
l_lam <- log(lambdaHat)+seq(-2, 6, 1)
(logLR2_u <- 2*(logLk(lambdaHat) - logLk(exp(l_lam))))
```

This seems more sensible, though still somewhat skewed, with the MLE as the third of 9 evenly-spaced points on the $u$ scale.  What about $\theta = 1/\lambda$? 

```{r CItheta}
theta <- (1/lambdaHat + seq(-.016, .06, .004))
(logLR2_th <- 2*(logLk(lambdaHat) - logLk(1/theta)))
```

This looks worse:  With $\hat\theta = 1/\hat\lambda = 0.0178$, it seems infeasible to get a sequence of only a few equally-spaced positive numbers that include 0 and still produce numbers just over [`qchisq(0.99, 1)`](https://www.rdocumentation.org/packages/stats/versions/3.5.2/topics/Chisquare) = 6.63 on either ends, as we created on the $\lambda$ scale, let alone both as we have on the $u$ scale.  

This suggests we should parameterize our analysis in terms of $u = \log(\lambda)$.  To confirm this, let's create plots on all three scales, starting with $u$.  

However, to save space on CRAN, we will not plot them by default;  to see the plots, a user will need to manually set `makePlots <- TRUE`:  

```{r 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)
}
```

Let's produce this same plot without the log scale for $\lambda$:  

```{r 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)
```

The plot vs. $\log(\lambda)$ is obviously skewed, but this linear plot is vastly worse.  

What about linear in $\theta = 1/\lambda$?  

```{r 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)
```

Clearly, we don't want to mess with the $\theta = 1/\lambda$ scale, and $\u = \log(\lambda) seems the best for understanding what's happening here.  

## Monte Carlo the time between major nuclear crises

If we had a probability distribution for $\lambda$, $u = \log(\lambda)$, or $\theta = 1/\lambda$, we could simulate that.  To get such, we recall that one statement of [Bayes' theorem](https://en.wikipedia.org/wiki/Bayes'_theorem) is that the posterior is proportional to the likelihood times the prior.  

However, what should we use as a prior?  The Wikipedia article on the [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution) describes several more or less standard priors for that distribution. There's not just one, and they all seem more complicated than what we need here.  Instead, we will use the [improper prior](https://en.wikipedia.org/wiki/Prior_probability#Improper_priors) that is uniform in $u = \log(\lambda)$.  To support this, we note that the exponential distribution is closer to the lognormal than to a normal, and the distribution of the reciprocal of an exponential random variable is even farther from normal, as we see in the following simulation:  

```{r simExp}
set.seed(1)
simExp <- rexp(1000)
if(makePlots){
  qqnorm(simExp, datax=TRUE)
  qqnorm(simExp, datax=TRUE, log='x')
  qqnorm(1/simExp, datax=TRUE)
}
```

Let's rewrite the above likelihood in terms of $u$:

\[
L(u | \mathbf{T}) = \exp[−S_k e^{-u} - (k-1)u].  
\]

With an improper prior locally uniform in $u = \log(\lambda)$, we get the following:  

\[
P(a < \lambda \leq b | \mathbf{T}) \propto \int_{\log(a)}^{\log(b)}{\exp[-S_k e^{-u}] e^{-(k-1)u} du}
\]

Let's transform this back by replacing $u$ with $\lambda = e^u$:  

\[
P(a < \lambda \leq b | \mathbf{T}) \propto \int_a^b \exp[-S_k / \lambda] \lambda^{-k} d\lambda
\]

This says that the posterior for $\lambda$ follows an [inverse-gamma distribution](https://en.wikipedia.org/wiki/Inverse-gamma_distribution) with shape parameter $(k-1)$ and scale $S_k$.  The moments for this distribution are as follows:  

\[
\mathbb{E}(\lambda^r | \mathbf{T}) = S_k^r \Gamma(k-1-r) / \Gamma(k-1).
\]

If $(k-1-r)$ is an integer less than 1, this is infinite, which it is for $k$ = 2 and $r$ = 1, the case of most interest here.  This, in turn, means that the sample moments of real or Monte Carlo data will be highly erratic.  

This also elevates the priority for increasing $k$ by considering a larger list of nuclear crises, as previously mentioned.  

Functions to compute the density, cumulative probability distribution (CDF), quantiles, and random numbers for this distribution are available in the CRAN package [`invgamma`](https://www.rdocumentation.org/packages/invgamma/versions/1.1/topics/invgamma).  We will use those in the following.  

```{r 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)
```

This distribution is obviously highly skewed as expected, with a mode (MLE) at 56 years and a mean estimated here at 439  

## Probability that a major nuclear war might lead to the extinction of civilization

The analysis in the Wikiversity article on ["time to extinction of civilization"](https://en.wikiversity.org/wiki/Time_to_extinction_of_civilization#Probability_that_a_nuclear_crisis_would_generate_a_major_nuclear_war) includes estimated of the probability that a major nuclear crisis like the [1962 Cuban Missile Crisis](https://en.wikipedia.org/wiki/Cuban_Missile_Crisis) or the [1983 Soviet nuclear false alarm incident](https://en.wikipedia.org/wiki/1983_Soviet_nuclear_false_alarm_incident) would lead to a major nuclear war.  The numbers given there ranged from 0.3 to 0.6 with a typical number of 0.45.  For present purposes, we shall assume that (0.3, 0.6) represent an 80 percent, equal tail [confidence interval](https://en.wikipedia.org/wiki/Confidence_interval) of a [beta distribution](https://en.wikipedia.org/wiki/Beta_distribution) and estimate its two shape parameters, $\alpha$ and $\beta$.  

Doing this requires an iteration, because no simple formula exists for this. We want to shape1 = $\alpha$ and shape2 = $\beta$ to satisfy the following:  

> `0.1 = pbeta(0.3, shape1, shape2)`    
> `0.9 = pbeta(0.6, shape1, shape2)`

The most reliable way to solve equations like these is to convert this into a minimization problem and use something like [`optim`](https://www.rdocumentation.org/packages/stats/versions/3.5.2/topics/optim):  

```{r 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))
```

The beta distribution with parameters (1, 1) is just the uniform distribution.  Manual computation shows that this is the correct answer for this case.   

```{r betaSolve}
(betaSolve <-optim(c(1,1), Dev2, 
            method="L-BFGS-B", lower=c(0,0)))
```

What's the mean of this distribution?  

Recall that the mean of the [$B(\alpha, \beta)$ distribution](https://en.wikipedia.org/wiki/Beta_distribution) is $\alpha / (\alpha+\beta)$, and its variance is as follows:  

\[
\mathrm{var}(Q) = \alpha\beta / [(\alpha+\beta)^2 (\alpha+\beta+1)]
\]

```{r 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))))
```

That's quite close to the representative value of 0.45 discussed in the companion Wikiversity article on ["time to extinction of civilization"](https://en.wikiversity.org/wiki/Time_to_extinction_of_civilization).  

## Monte Carlo time to extinction of civilization

This section will start with a function to generate `N` random times to Armageddon as follows:   

1.  start timing
2.  Generate `N` random variates `Q` ~ [$B(\alpha, \beta)$](https://stat.ethz.ch/R-manual/R-patched/library/stats/html/Beta.html) indicating the probability that each simulated crisis in a sequence would produce a nuclear Armageddon.  
3.  From this, generate `N` random variables `K` ~ [$NB(Q, 1)$](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/NegBinomial.html) indicating the number of simulated crises in a series required to produce one nuclear Armageddon with `Q[i]` = probability of each being the last, for `i = 1, ..., N`. 
4.  For each `i`, compute `Time[i] <- sum(rlambda2(K[i]))`.  
5.  compute `elapsed.time` 
6.  Return (`Time`, `gammapars`, `gammaGOF`, 
`elapsed.time`)

First do this with `set.seed(1), N`=10 and time the result.  Then `set.seed(j), N`=10^j, j = 2, 3, ..., timing each one.  Save the results until the time gets too long to continue or we get to `N` = 1e7.  

```{r 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))
```
  
This all looks sensible.  Let's try larger sample sizes:   

```{r mcArm2}
set.seed(2)
mcArm2 <- mcArmageddon(100)
attributes(mcArm2)
```

```{r mcArm3}
set.seed(3)
mcArm3 <- mcArmageddon(1000)
attributes(mcArm3)
```

`N` = 1000 still takes only 0.009 seconds.  

```{r mcArm4}
set.seed(4)
mcArm4 <- mcArmageddon(1e4)
attributes(mcArm4)
```

The time was still only 0.074 seconds, so let's try `N=1e5`:  

```{r mcArm5}
set.seed(5)
mcArm5 <- mcArmageddon(1e5)
attributes(mcArm5)
```

The time was still only 0.535 -- well under 10 times `N= 1e4`.  

What about a million?  

```{r million}
set.seed(6)
mcArm6 <- mcArmageddon(1e6)
attributes(mcArm6)
```

This too just over 5 seconds.  

For the Wikiversity article on "Time to extinction of civilization", we'd like the percentages of these times that are less than 40 and 60 years, representing roughly the remaining lives of half the people currently alive today and the time remaining in the twenty-first century as of this writing, as well as the quantiles of one in a million and one in a thousand chances:  

```{r MCsumStats}
mean(mcArm6)
mean(mcArm6<40)
mean(mcArm6<60)
quantile(mcArm6, c(1e-6, 1e-3))
```

Let's see if we can generate 1e7 random times in, hopefully, just over 50 seconds:  

```{r 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)))
}
```

Let's make a normal probability plot of `log(mcArm7)`. With this many points, the standard [`qqnorm`](https://www.rdocumentation.org/packages/stats/versions/3.5.2/topics/qqnorm) function can take a long time creating the plot.  Let's start by sorting the points as a separate step:  

```{r sort7}
if(fda::CRAN()){
  mcArm. <- mcArm6
} else mcArm. <- mcArm7

mcArm.s <- sort(mcArm.)
quantile(mcArm.s)
```

Next, let's call `qqnorm` without plotting:  

```{r qqnorm7}
str(qq7 <-as.data.frame(qqnorm(mcArm.s,
                          plot.it=FALSE)))
```

Let's cut the data down to the first and last 10 plus 9 of the next 90 from each end plus 1 percent of the rest:  

```{r 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.  
```

```{r 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()
```