Application to Data

Each generation and yield pair in the first 24 readings is one point on each plot: the reading 4.87 at generation 51 is the highest point on the standard plot and the point $(51, \ln 4.87) \approx (51, 1.583)$ on the semi-log plot. R spells $\ln(\text{yield})$ as log(yield_KT), because log() in R is the natural logarithm.

The standard plot is not linear, because the readings sit near 0 through generation 45 and rise steeply after it. The semi-log plot is close to linear from generation 5 to generation 44, because the slopes between consecutive points there all lie between 0.814 and 0.832. Between generations 5 and 10, for example, the semi-log slope is

The readings through generation 44 therefore suit an exponential model:

The semi-log plot bends down over the last readings, where the slope between generations 50 and 51 falls to

Those last readings grow more slowly than the exponential model. Example 3 follows that slowdown into the later readings.

(b) The semi-log form of the exponential model is $\ln(\text{yield}) = m\,\text{generation} + b$, with $k = m$ and $C = e^{b}$, as derived on the Lesson 5 sheet, section 4. The two points are the first and last readings on the semi-log plot, $(1, \ln(6.66\times10^{-18}))$ and $(51, \ln 4.87)$. Substitute them into the slope formula of the two-point method (Lesson 6 sheet, section 2):

Here the model keeps no horizontal shift, so that its $C$ compares directly with the fitted $C$ of part (c). Solve for the intercept $b$, the height of the line at generation 0, from the first point:

R returns $k = 0.8226701$ and $C = 2.925453\times10^{-18}$, the hand values at full precision.

(c) Fit the semi-log line to the two vectors from part (a), then convert its slope and intercept to $k$ and $C$ by the method of the Lesson 6 sheet, sections 4 and 5:

R returns the slope $m = 0.824163$, the intercept $b = -39.880049$, and $C = 4.789773\times10^{-18}$. The output of coef() names each number by its letter in the formula, so cf[1] is the slope $m$ and cf[2] is the intercept $b$.

The fitted intercept is $-39.880049 - (-40.373) \approx 0.493$ above the hand intercept of part (b), because fitModel() uses all 24 readings and the hand line uses two. The fitted $C$ is therefore about $e^{0.493} \approx 1.637$ times the hand value. The best-fit model is

The purple curve is the model, drawn with col = "purple". It stays within 0.01 kilotons of every reading through generation 47 and then rises above the readings. The last line returns $8.604304 \approx 8.604$ kilotons at generation 51, against the reading of 4.87.

Rebuild the model of Example 1 and evaluate it at generation 61:

R returns a predicted yield of 32,659.54 kilotons at generation 61. The prediction is an extrapolation as defined on the Lesson 10 sheet, section 6, because generation 61 lies outside the inputs of the first 24 readings, generation 1 to 51.

Rebuild the model of Example 1 and evaluate it at the four generations:

R returns the four model values in the second column of the table, and the third column gives the readings from the Plutonium Record:

The exponential model does not fit the ten later readings. It is 2.5 times the reading at generation 52 and 1,608.8 times the reading at generation 61, because the later readings level off near 20 kilotons and the model grows without bound. The readings level off because the sphere expands until it is no longer critical and the chain reaction stops.

(b) Estimate the four parameters by the method of the Lesson 9 sheet, section 4. The figure marks each estimate on the readings from generation 40 on: the two levels $v$ and $v + A$, the halfway point at the mean $m$, and the spread $s$ as the input distance from $m$ to the five-sixths level.

• The vertical shift $v$ is the level on the left, near 0 where the readings start: $v = 0$.
• The amplitude $A$ is the whole rise from the level on the left to the level on the right. The last four readings, from 19.90 to 20.30, sit near the level 20: $A = 20$.

• The mean $m$ is the input at the halfway output $v + A/2 = 10$. The readings are 7.92 at generation 52 and 11.50 at generation 53. The output 10 lies $(10 - 7.92)/(11.50 - 7.92) \approx 0.58$ of the way between them:

  $$m \approx 52 + 0.58 \approx 52.6.$$

• The estimate of the spread $s$ is the input distance from $m$ to the five-sixths level $v + 5A/6 \approx 16.7$, which stands in for the pnorm(1) point (Lesson 9 sheet, section 2). The readings are 14.70 at generation 54 and 17.00 at generation 55. The output 16.7 lies $(16.7 - 14.70)/(17.00 - 14.70) \approx 0.87$ of the way between them, at about generation $54 + 0.87 \approx 54.9$:

  $$s \approx 54.9 - 52.6 = 2.3.$$

The last line prints $A = 20.15251639$, $m = 52.61711780$, $s = 2.37144147$, and $v = 0.02243469$. The best-fit model is

The red curve is the model, drawn with col = "red". It lies within 0.17 kilotons of every reading, at both levels and along the rise between them.

Rebuild the sigmoid model of Example 4 and evaluate it at generations 47 and 53:

R returns predictions of $0.2023269 \approx 0.202$ kilotons at generation 47 and $11.39113 \approx 11.391$ kilotons at generation 53. Both predictions are interpolations, because 47 and 53 lie inside generations 1 to 61, the inputs of all 34 readings. At generation 53 the reading is 11.50, and the model is within 1 percent of it. At generation 47 the reading is 0.323, and the model is 37 percent low.

Evaluate the model at generations 52 and 53, where the readings cross 10 kilotons:

R returns $8.029908 \approx 8.030$ kilotons at generation 52, below 10, and $11.39113 \approx 11.391$ kilotons at generation 53, above 10. The model is increasing, because pnorm increases with generation and $A \approx 20.153$ is positive. So the first generation with a yield of at least 10 kilotons is generation 53.