Practical — Distribution Theory (STS-107)

Topics Covered

Objectives

1. Knowing the manual procedures and also their implementation using R.
2. Generation of random samples from any distribution.
3. Identifying an appropriate probability distribution for the given data.
4. Fitting and testing the probability distribution.
5. Drawing the probability distribution curves and stating the nature of the distributional curve and its properties for the given data sets.

The Ten Practicals

Practical 1 — Uniform Random Numbers

1. Problem

(a) Generate five uniform random numbers by hand with the linear congruential generator \(a = 5\), \(c = 3\), \(m = 16\) and seed \(x_0 = 7\). (b) Generate 1000 uniform numbers in R and test whether they are a sample from \(U(0,1)\).

2. Aim

To generate pseudo-random uniform numbers by a linear congruential generator, and to check a large sample against U(0, 1).

3. Formula

Starting from a seed \(x_0\), the \(u_i\) lie in \([0,1)\) and, for well chosen \(a, c, m\), pass the usual tests of uniformity and independence. \(D\) is the Kolmogorov–Smirnov distance between the sample's distribution function and that of \(U(0,1)\).

Applying it:

1. Multiply the last \(x\) by \(a\), add \(c\), and keep the remainder on division by \(m\).
2. Divide by \(m\) to get \(u\); repeat.
3. For the large sample, compare the mean and variance with \(\tfrac12\) and \(\tfrac1{12}\), and test with Kolmogorov–Smirnov.

4. Calculation

lcg <- function(n, a = 5, c = 3, m = 16, x0 = 7) {
  x <- numeric(n); x[1] <- (a * x0 + c) %% m
  for (i in 2:n) x[i] <- (a * x[i - 1] + c) %% m
  x / m
}
lcg(5)                         # the hand working above
set.seed(2025)                 # reproducibility: same seed, same sample
u <- runif(1000)
c(mean = mean(u), var = var(u))   # compare with 0.5 and 1/12 = 0.08333
ks.test(u, "punif")            # formal test of uniformity
hist(u, breaks = 20, freq = FALSE, main = "1000 uniform variates")
abline(h = 1, col = "red", lwd = 2)      # the U(0,1) density

[1] 0.3750 0.0625 0.5000 0.6875 0.6250
      mean        var 
0.50822464 0.08337176 

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  u
D = 0.022013, p-value = 0.7177
alternative hypothesis: two-sided

5. Result

(a) The five numbers are \(0.3750, 0.0625, 0.5000, 0.6875, 0.6250\). The modulus bounds the period at \(m = 16\), so this toy generator repeats quickly; real ones use \(m\) of the order of \(2^{31}\) or more. (b) The 1000 values have mean 0.508 and variance 0.0834, against 0.5 and 0.0833, and \(D = 0.022\) with \(p = 0.72\): uniformity is not rejected.

Interpretation: the numbers are not random at all — they are a completely determined sequence that behaves like one, which is why the term is pseudo-random and why the seed reproduces a run exactly.

Practical 2 — Discrete Variates by Inverse Transform

1. Problem

Using the five uniform numbers of Practical 1, generate a sample from the Poisson distribution with \(\lambda = 1.5\) by the inverse transform. Then generate samples from the binomial, Poisson, geometric and negative binomial distributions in R.

2. Aim

To generate discrete random variates by inverting the cumulative distribution function, and with R's built-in generators.

3. Formula

The value \(x\) is returned with probability \(F(x) - F(x-1) = P(X = x)\), which is exactly what is wanted. For the geometric (failures before the first success) the inverse is in closed form, \(x = \left\lceil \dfrac{\ln(1-u)}{\ln(1-p)} \right\rceil - 1\), so no table is needed. A binomial Bin\((n,p)\) variate can also be generated as the count of successes in \(n\) independent Bernoulli trials: draw \(n\) uniforms and count how many fall below \(p\). This is slower than the inverse transform for large \(n\) but needs no cumulative table.

Applying it:

1. Tabulate \(P(X = x)\) and \(F(x)\).
2. For each \(u\), find the smallest \(x\) with \(u \le F(x)\).

4. Calculation

Then read each \(u\) against the last column:

# Inverse transform written out, to show the method rather than hide it
inv_poisson <- function(u, lambda) {
  x <- 0
  cdf <- exp(-lambda)
  while (u > cdf) {
    x <- x + 1
    cdf <- cdf + exp(-lambda) * lambda^x / factorial(x)
  }
  x
}
sapply(c(0.3750, 0.0625, 0.5000, 0.6875, 0.6250), inv_poisson, lambda = 1.5)
# Built-in generators
set.seed(2025)
rbinom(10, size = 5, prob = 0.4)
rpois(10, lambda = 1.5)
rgeom(10, prob = 0.3)          # counts FAILURES before the first success
rnbinom(10, size = 4, prob = 0.6)

[1] 1 0 1 2 2
 [1] 3 2 2 2 3 2 3 1 2 2
 [1] 1 4 2 0 1 3 1 1 2 0
 [1] 4 8 5 0 0 1 4 1 6 2
 [1] 0 1 2 3 0 0 6 1 4 1

5. Result

The Poisson(1.5) sample is \(1, 0, 1, 2, 2\), and R's inverse transform returns the same values exactly. R's built-in generators give samples of ten from each of the four distributions.

Practical 3 — Continuous Variates

1. Problem

(a) From the uniforms \(0.1273, 0.4820, 0.7391, 0.9052, 0.0361\), generate a sample from the exponential distribution with \(\theta = 0.5\). (b) From the pairs \((u_1, u_2) = (0.2740, 0.6130)\) and \((0.8821, 0.1094)\), generate standard normal variates by the Box–Muller method. (c) Generate samples from the normal, exponential, gamma, beta and Cauchy distributions in R.

2. Aim

To generate continuous random variates by the inverse transform, by the Box–Muller transformation, and by composition.

3. Formula

Exponential. \(F(x) = 1 - e^{-\theta x}\). Setting \(u = F(x)\) and solving, \( u = 1 - e^{-\theta x} \Rightarrow e^{-\theta x} = 1 - u \Rightarrow x = -\ln(1-u)/\theta \).

Box–Muller. From two independent uniforms, \(z_1\) and \(z_2\) are two independent \(N(0,1)\) variates. It is the polar form of the two-dimensional normal: \(-2\ln u_1\) is an Exponential\((\tfrac12)\), which is \(\chi^{2}_{2}\), giving the squared radius, and \(2\pi u_2\) is the angle, uniform on the circle.

Others by composition. Gamma\((k, \theta)\) with integer \(k\) is the sum of \(k\) independent Exponential\((\theta)\) variates. Beta\((p,q)\) is \(G_1/(G_1 + G_2)\) with \(G_1 \sim\) Gamma\((p,1)\) and \(G_2 \sim\) Gamma\((q,1)\) independent. Cauchy is \(\tan\!\left[\pi(u - \tfrac12)\right]\), or equivalently the ratio of two independent standard normals.

Applying it:

1. Exponential: compute \(-\ln(1-u)\) for each \(u\) and divide by \(\theta\).
2. Normal: for each pair, compute the radius \(\sqrt{-2\ln u_1}\) and the angle \(2\pi u_2\); multiply the radius by the cosine and the sine.

4. Calculation

Sample mean \(= 1.812126\), against the theoretical \(1/\theta = 2\).

Box–Muller with \(u_1 = 0.2740\), \(u_2 = 0.6130\):

A second pair, \(u_1 = 0.8821\), \(u_2 = 0.1094\): \(\sqrt{-2\ln u_1} = 0.500899\), \(\cos = 0.772911\), \(\sin = 0.634515\), giving \(z_1 = 0.387150\) and \(z_2 = 0.317828\).

u <- c(0.1273, 0.4820, 0.7391, 0.9052, 0.0361)
x <- -log(1 - u) / 0.5         # exponential, theta = 0.5, by inverse transform
round(x, 6); mean(x)
# Box-Muller written out, to check against the hand working above
box_muller <- function(u1, u2) {
  r <- sqrt(-2 * log(u1))
  c(r * cos(2 * pi * u2), r * sin(2 * pi * u2))
}
round(box_muller(0.2740, 0.6130), 6)
round(box_muller(0.8821, 0.1094), 6)
# Built-in generators: the mean of 1000 from each
set.seed(2025)
round(c(exp = mean(rexp(1000, rate = 0.5)), norm = mean(rnorm(1000)),
        gamma = mean(rgamma(1000, shape = 3, rate = 2)),
        beta = mean(rbeta(1000, shape1 = 2, shape2 = 5))), 4)
rcauchy(5, location = 0, scale = 1)   # no mean to compare: the Cauchy has none

[1] 0.272327 1.315560 2.687236 4.711972 0.073535
[1] 1.812126
[1] -1.220292 -1.048876
[1] 0.387150 0.317828
   exp   norm  gamma   beta 
2.0296 0.0158 1.4833 0.2904 
[1] -3.4108950  0.7219658  2.1882630 -2.4570693 -0.0309465

5. Result

(a) The exponential sample is \(0.2723, 1.3156, 2.6872, 4.7120, 0.0735\), with mean 1.81 against \(1/\theta = 2\). Five observations is far too few to expect agreement; the point of the check is that the value is of the right order, not that it matches. (b) The normal variates are \(-1.2203, -1.0489\) and \(0.3872, 0.3178\). (c) With 1000 values, the sample means are close to the theoretical means: 2 for the exponential, 0 for the normal, \(3/2\) for the gamma and \(2/7 = 0.286\) for the beta.

Practical 4 — Fitting a Discrete Distribution

1. Problem

(a) Five coins were tossed 100 times and the number of heads recorded. Fit a binomial distribution and test the fit.

(b) Accidents per day at a junction over 200 days are given below. Fit a Poisson distribution and test the fit.

2. Aim

To fit binomial and Poisson distributions by the method of moments and to test each fit by the chi-square test of goodness of fit.

3. Formula

Degrees of freedom are the number of cells, less one for the total, less one for each parameter estimated. For the Poisson the method of moments and maximum likelihood agree, both giving the sample mean.

Applying it:

1. Estimate the parameter from the sample mean.
2. Compute the expected frequencies; the last cell takes the remainder so the totals agree.
3. Form \(\chi^{2}\), find its degrees of freedom, and compare with the 5% point (or find the p-value).

4. Calculation

(a) Binomial. \(\sum f x = 0(2) + 1(14) + 2(20) + 3(34) + 4(22) + 5(8) = 0 + 14 + 40 + 102 + 88 + 40 = 284\), so \(\bar x = 284/100 = 2.84\), \(\hat p = 2.84/5 = 0.568\), \(\hat q = 0.432\), and \(E_x = 100 \binom{5}{x} (0.568)^{x}(0.432)^{5-x}\):

The expected frequencies sum to \(100.000000\), which is the arithmetic check that \(\hat p\) and the binomial probabilities were computed correctly. Degrees of freedom are \(6 - 1 - 1 = 4\), so \(\chi^{2} = 4.007575\) on 4 d.f., \(P(\chi^{2}_{4} > 4.007575) = 0.404982\); the tabulated 5% point is 9.488.

(b) Poisson. \(\sum f x = 0 + 70 + 80 + 60 + 32 + 10 = 252\), so \(\hat\lambda = 252/200 = 1.26\), and \(E_x = 200\, e^{-1.26}(1.26)^{x}/x!\) for \(x = 0, \ldots, 4\), the last cell taking the remainder:

\(6 - 1 - 1 = 4\) degrees of freedom: \(\chi^{2} = 1.551074\), \(P(\chi^{2}_{4} > 1.551074) = 0.817558\).

A caution on the last cell. Its expected frequency is \(1.8838\), below the usual rule of thumb that no expected frequency should fall under 5. Strictly the \(4\) and \(5+\) cells should be pooled, giving \(E = 7.8416\) against \(O = 10\) and \(3\) degrees of freedom. That pooled test gives \(\chi^{2} = 1.438000\) on 3 d.f. with \(P(\chi^{2}_{3} > 1.438000) = 0.696652\).

f <- c(2, 14, 20, 34, 22, 8); x <- 0:5
N <- sum(f); phat <- sum(f * x) / (N * 5)
phat
E <- N * dbinom(x, size = 5, prob = phat)
chisq <- sum((f - E)^2 / E)
round(c(chisq = chisq, df = 4, p = pchisq(chisq, df = 4, lower.tail = FALSE)), 6)
# Poisson
g <- c(60, 70, 40, 20, 8, 2); y <- 0:5
lam <- sum(g * y) / sum(g); lam
Eg <- c(sum(g) * dpois(0:4, lam), sum(g) * ppois(4, lam, lower.tail = FALSE))
round(c(chisq = sum((g - Eg)^2 / Eg), p = pchisq(sum((g - Eg)^2 / Eg), 4, lower.tail = FALSE)), 6)
# the 4 and 5+ cells pooled: 3 d.f.
O4 <- c(g[1:4], sum(g[5:6])); E4 <- c(Eg[1:4], sum(Eg[5:6]))
round(c(chisq = sum((O4 - E4)^2 / E4), p = pchisq(sum((O4 - E4)^2 / E4), 3, lower.tail = FALSE)), 6)

[1] 0.568
   chisq       df        p 
4.007575 4.000000 0.404982 
[1] 1.26
   chisq        p 
1.551074 0.817558 
   chisq        p 
1.438000 0.696652 

5. Result

(a) \(\chi^{2} = 4.008 < 9.488\) (\(p = 0.405\)), so the binomial model is not rejected. Interpretation. \(\hat p = 0.568\), not \(0.5\), so the coins are not quite fair, or 100 trials is too few to tell. The goodness-of-fit test answers a different question from "is \(p = 0.5\)": it asks whether the shape is binomial, and it is.

(b) \(\chi^{2} = 1.551\) (\(p = 0.818\)): an excellent fit, with no evidence against the Poisson. Pooling the small last cell, as should be done and stated rather than passed over, leaves the conclusion unchanged (\(p = 0.697\)).

Choosing a Continuous Family

Practical 5 — Fitting a Continuous Distribution

1. Problem

R's built-in data set precip gives the average annual precipitation, in inches, of 70 cities of the United States (from 7 to 67 inches). Fit a normal distribution by maximum likelihood and test the fit by the chi-square test, with the classes under 20, 20–30, 30–40, 40–50 and 50 or more.

2. Aim

To fit a continuous distribution by maximum likelihood and to test the fit by chi-square on grouped frequencies.

3. Formula

The maximum likelihood estimates of the normal are the sample mean and the variance with divisor \(n\). Each class \((a_j, b_j]\) gets the probability the fitted normal gives it; the degrees of freedom are the number of classes, less one, less the two parameters estimated.

Applying it:

1. Compute \(\hat\mu\) and \(\hat\sigma\).
2. Standardise the class boundaries and read \(\Phi\) at each.
3. Expected frequency of each class = \(n\) × its probability; form \(\chi^{2}\) on \(5 - 1 - 2 = 2\) d.f.

4. Calculation

\(\sum x = 2442\), so \(\hat\mu = 2442/70 = 34.8857\); \(\sum(x - \bar x)^{2} = 12963.19\), so \(\hat\sigma^{2} = 185.1884\) and \(\hat\sigma = 13.6084\).

Boundaries 20, 30, 40, 50 standardise to \(z = -1.0939, -0.3590, 0.3758, 1.1107\), with \(\Phi(z) = 0.1370, 0.3598, 0.6465, 0.8666\).

\(\chi^{2} = 12.84\) on \(5 - 1 - 2 = 2\) d.f.; the 5% point is \(\chi^{2}_{0.05,\,2} = 5.991\), and \(P(\chi^{2}_{2} > 12.84) = 0.0016\).

x <- precip; n <- length(x); n
mu <- mean(x); s <- sqrt(mean((x - mu)^2))    # maximum likelihood: divisor n
round(c(mu = mu, sigma = s), 4)
br <- c(-Inf, 20, 30, 40, 50, Inf)
O <- table(cut(x, br))
P <- diff(pnorm(br, mu, s)); E <- n * P
cbind(O = as.vector(O), E = round(E, 4), contrib = round((as.vector(O) - E)^2 / E, 6))
chisq <- sum((O - E)^2 / E)
round(c(chisq = chisq, df = 5 - 1 - 2, p = pchisq(chisq, 2, lower.tail = FALSE)), 6)
hist(x, breaks = c(0, 20, 30, 40, 50, 70), main = "precip"); curve(dnorm(x, mu, s), add = TRUE)   # unequal classes: hist() draws densities

[1] 70
     mu   sigma 
34.8857 13.6084 
      O       E  contrib
[1,] 13  9.5905 1.212078
[2,]  5 15.5947 7.197824
[3,] 25 20.0679 1.212145
[4,] 21 15.4118 2.026251
[5,]  6  9.3350 1.191471
    chisq        df         p 
12.839770  2.000000  0.001629 

5. Result

\(\chi^{2} = 12.84 > 5.991\) (\(p = 0.0016\)): the normal distribution is rejected for these data. Most of the discrepancy is in one class: only 5 cities have 20–30 inches where the fitted normal expects 15.6, while 13 have under 20 where it expects 9.6. The very dry cities form a group of their own, which a single bell-shaped curve cannot describe; a fitted distribution is only as good as the test of its shape says it is.

Practical 6 — Goodness of Fit of a Cauchy Distribution

1. Problem

Twenty observations are given below. Fit a Cauchy distribution and test the goodness of fit at 5%.

−91.62, 5.09, −138.90, 5.19, 6.13, −17.55, 7.98, 0.17, 7.01, 5.38, 6.42, 5.59, 5.88, −18.36, 9.85, 3.35, 2.31, 2.87, 4.84, 8.66

2. Aim

To fit a Cauchy distribution, which has no mean and no variance, by its quantiles, and to test the fit by the Kolmogorov–Smirnov test.

3. Formula

The Cauchy has no mean and no variance, so no moment-based method can be used at all, and a chi-square test is awkward because the tails are so heavy that cells must be pooled aggressively. Fit by quantiles instead: the quartiles of the Cauchy are at \(\theta \pm \lambda\). Then test with Kolmogorov–Smirnov, which uses the whole distribution function and needs no cells — and which is legitimate here by the Glivenko–Cantelli lemma of Probability Theory, Unit 3.

Applying it:

1. Sort the data; find the median and the quartiles.
2. Estimate \(\theta\) and \(\lambda\).
3. At each ordered value compute \(F(x_{(i)})\) and the two gaps to the empirical distribution function; \(D\) is the largest.
4. Compare \(D\) with the 5% critical value for \(n = 20\), 0.294.

4. Calculation

Sorted, the quartiles are \(Q_1 = 1.775\), median \(= 5.14\), \(Q_3 = 6.2025\), so \(\hat\theta = 5.14\) and \(\hat\lambda = (6.2025 - 1.775)/2 = 2.2138\).

The largest gap is \(D = 0.1690\), at \(x_{(4)} = -17.55\).

set.seed(107)                  # how the data were drawn
x <- round(rcauchy(20, location = 5, scale = 2), 2)
q <- quantile(x, c(0.25, 0.5, 0.75)); q
theta <- unname(q[2]); lambda <- unname((q[3] - q[1]) / 2)
round(c(theta = theta, lambda = lambda), 4)
ks.test(x, "pcauchy", location = theta, scale = lambda)

   25%    50%    75% 
1.7750 5.1400 6.2025 
 theta lambda 
5.1400 2.2138 

	Exact one-sample Kolmogorov-Smirnov test

data:  x
D = 0.16904, p-value = 0.5604
alternative hypothesis: two-sided

5. Result

\(D = 0.169 < 0.294\) (\(p = 0.56\)): the Cauchy distribution with \(\hat\theta = 5.14\) and \(\hat\lambda = 2.21\) fits the data. The estimates are close to the values the data were drawn from, 5 and 2, although three of the twenty observations lie below −17: the median and quartiles are not moved by such values, where a mean would be ruined. Because \(\theta\) and \(\lambda\) were estimated from the same data, the stated \(p\) is approximate, and on the safe side.

Practical 7 — Fitting a Two-Parameter Gamma

1. Problem

Ten times to failure, in hours, are \(31,\ 58,\ 87,\ 119,\ 153,\ 192,\ 238,\ 296,\ 378,\ 520\). Fit a two-parameter gamma distribution: start from the method of moments and refine by the method of scoring.

2. Aim

To fit a two-parameter gamma distribution by maximum likelihood, starting from the moment estimates and iterating by the method of scoring.

3. Formula

The moments give starting values from \(\bar x = k/\theta\) and \(s^{2} = k/\theta^{2}\) (divisor \(n\)). The maximum likelihood equations for the gamma have no closed form: \(\hat\theta = k/\bar x\), and \(k\) solves \(S(k) = 0\), where \(\psi\) is the digamma function and \(\psi'\) the trigamma. Scoring divides the score by the information and repeats.

Applying it:

1. Compute \(\bar x\), \(s^{2}\), and the moment estimates.
2. Compute \(\ln \bar x\) and the mean of the \(\ln x\).
3. Apply the scoring step until \(k\) stops changing; then \(\hat\theta = \hat k/\bar x\).

4. Calculation

\(\bar x = 207.2\), \(s^{2} = 21455.36\), so \(\hat k_0 = 207.2^{2}/21455.36 = 2.000984\) and \(\hat\theta_0 = 207.2/21455.36 = 0.009657\).

\(\ln \bar x = 5.333685\) and \(\overline{\ln x} = 5.037875\). At \(k_0\): \(\ln k_0 = 0.693639\), \(\psi(k_0) = 0.423419\), \(\psi'(k_0) = 0.644537\), so

So \(\hat k = 1.839405\) and \(\hat\theta = 1.839405/207.2 = 0.008877\).

t <- c(31, 58, 87, 119, 153, 192, 238, 296, 378, 520); n <- length(t)
xbar <- mean(t); s2 <- mean((t - xbar)^2)
k0 <- xbar^2 / s2; th0 <- xbar / s2
round(c(xbar = xbar, s2 = s2, k0 = k0, theta0 = th0), 6)
k <- k0; lg <- mean(log(t))        # scoring on the shape k
for (it in 1:4) {
  score <- n * (log(k) - log(xbar) + lg - digamma(k))
  info  <- n * (trigamma(k) - 1 / k)
  k <- k + score / info
  cat(sprintf("iteration %d: k = %.6f\n", it, k))
}
round(c(k = k, theta = k / xbar, mean = k / (k / xbar)), 6)

        xbar           s2           k0       theta0 
  207.200000 21455.360000     2.000984     0.009657 
iteration 1: k = 1.824240
iteration 2: k = 1.839272
iteration 3: k = 1.839405
iteration 4: k = 1.839405
         k      theta       mean 
  1.839405   0.008877 207.200000 

5. Result

The fitted gamma has shape \(\hat k = 1.839\) and rate \(\hat\theta = 0.00888\) per hour; its mean, \(\hat k/\hat\theta = 207.2\) hours, is the sample mean, as maximum likelihood for the gamma always gives. Scoring settled in three steps from the moment estimate \(\hat k_0 = 2.001\). Since \(\hat k > 1\), the hazard rate rises with age, the same conclusion the Weibull of Practical 9 reaches for these data.

Practicals 8–10 — One Dataset, Three Families

Practical 8 — Fitting a Lognormal

1. Problem

Fit a two-parameter lognormal distribution to the ten failure times \(31,\ 58,\ 87,\ 119,\ 153,\ 192,\ 238,\ 296,\ 378,\ 520\) by maximum likelihood, and find the fitted median, the fitted mean and \(P(X > 200)\).

2. Aim

To fit a lognormal distribution by maximum likelihood, which is the normal fit on the log scale.

3. Formula

Since \(\ln X \sim N(\mu, \sigma^{2})\), maximum likelihood is the normal fit to the logs. (The divisor is \(n\), not \(n-1\): maximum likelihood, not the unbiased estimator.)

Applying it:

1. Take the natural log of each time.
2. Compute \(\hat\mu\) and \(\hat\sigma^{2}\).
3. Transform back for the median and mean; \(P(X > 200) = 1 - \Phi\{(\ln 200 - \hat\mu)/\hat\sigma\}\).

4. Calculation

The logs are 3.433987, 4.060443, 4.465908, 4.779123, 5.030438, 5.257495, 5.472271, 5.690359, 5.934894, 6.253829, whose sum is \(50.378748\), so

t <- c(31, 58, 87, 119, 153, 192, 238, 296, 378, 520)
mu <- mean(log(t)); s2 <- mean((log(t) - mu)^2)
round(c(mu = mu, sigma2 = s2, sigma = sqrt(s2)), 6)
round(c(median = exp(mu), mean = exp(mu + s2 / 2),
        P200 = 1 - pnorm((log(200) - mu) / sqrt(s2))), 6)

      mu   sigma2    sigma 
5.037875 0.686784 0.828724 
    median       mean       P200 
154.142088 217.297367   0.376658 

5. Result

The fitted lognormal has \(\hat\mu = 5.0379\) and \(\hat\sigma = 0.8287\): median life 154 hours, mean life 217 hours, and \(P(X > 200) = 0.377\).

Interpretation. The lognormal fitted mean, \(217.30\), is within \(1\%\) of the Weibull's \(215.40\) (Practical 9) — the two families agree closely in the middle. They part company in the tail: the lognormal gives \(P(X > 200) = 0.3767\) against the Weibull's \(0.4377\), a relative difference of \(14\%\). Choosing between them on fit statistics alone is unreliable with ten observations, so the choice should rest on what is known about the failure mechanism.

Practical 9 — Fitting a Weibull by the Weibull Plot

1. Problem

Fit a two-parameter Weibull distribution to the ten failure times \(31,\ 58,\ 87,\ 119,\ 153,\ 192,\ 238,\ 296,\ 378,\ 520\) by the Weibull plot, and find the fitted mean life and the reliability at 200 hours.

2. Aim

To fit a Weibull distribution by least squares on its linearised distribution function, and to judge the fit by the straightness of the plot.

3. Formula

Plotting \(y = \ln\{-\ln[1-F]\}\) against \(x = \ln t\) gives a straight line of slope \(k\) and intercept \(-k\ln\lambda\). \(F\) at each point is estimated by the median rank \((i - 0.3)/(n + 0.4)\), the standard small-sample plotting position.

Applying it:

1. Rank the times; compute the median ranks \(F_i\).
2. Compute \(x = \ln t\) and \(y = \ln\{-\ln(1-F)\}\).
3. Fit the least-squares line; its slope is \(\hat k\), and \(\hat\lambda\) comes from the intercept.
4. Judge the fit by the correlation of the plotted points; use the fitted model.

4. Calculation

Numerator: \(-175.943100 + 263.537931 = 87.594831\). Denominator: \(2606.696660 - 2538.018250 = 68.678410\). So \(\hat k = 87.594831/68.678410 = 1.275435\).

The correlation of the plotted points is \(r = 0.999215\). Then

t <- c(31, 58, 87, 119, 153, 192, 238, 296, 378, 520); n <- length(t)
i <- 1:n; Fi <- (i - 0.3) / (n + 0.4)
x <- log(t); y <- log(-log(1 - Fi))
fit <- lm(y ~ x)
k <- unname(coef(fit)[2]); lambda <- unname(exp(-coef(fit)[1] / k))
round(c(k = k, lambda = lambda, r = cor(x, y)), 6)
round(c(mean = lambda * gamma(1 + 1 / k), R200 = exp(-(200 / lambda)^k)), 6)
plot(x, y, pch = 19, xlab = "ln t", ylab = "ln(-ln(1-F))"); abline(fit, col = "red")

         k     lambda          r 
  1.275435 232.297728   0.999215 
      mean       R200 
215.397136   0.437716 

5. Result

\(\hat k = 1.275\) and \(\hat\lambda = 232.3\) hours; mean life 215.4 hours and \(R(200) = 0.438\). A Weibull plot as straight as \(r = 0.9992\) is as strong an endorsement of the family as this sample size can give.

Interpretation. \(\hat k = 1.275 > 1\), so the hazard rate rises with age — these components wear out, and an exponential model would be wrong. But \(1.275\) is only moderately above 1, so the wear-out is mild; with ten observations the standard error on \(\hat k\) is large enough that \(k = 1\) could not be firmly ruled out.

Practical 10 — Fitting a Pareto, and Reading the Warning

1. Problem

Fit a two-parameter Pareto distribution to the ten failure times \(31,\ 58,\ 87,\ 119,\ 153,\ 192,\ 238,\ 296,\ 378,\ 520\) by maximum likelihood, and say what the fitted model can and cannot report.

2. Aim

To fit a Pareto distribution by maximum likelihood, and to check that the quantities to be reported exist under the fitted parameters.

3. Formula

By Unit 1, section 4, the Pareto mean requires \(\alpha > 1\) and the variance \(\alpha > 2\).

Applying it:

1. Take \(\hat x_m\) as the smallest observation.
2. Compute \(\sum\ln(t_i/\hat x_m)\) and \(\hat\alpha\).
3. Check what exists before reporting anything; report only what does.

4. Calculation

\(\hat x_m = 31\); \(\sum_i \ln(t_i/31) = 16.038876\), so \(\hat\alpha = 10/16.038876 = 0.623485\).

Here \(\hat\alpha = 0.623 < 1\), so under the fitted model neither the mean nor the variance exists. Substituting into \(\alpha x_m/(\alpha - 1)\) returns \(-51.33\), a negative mean life, which is the arithmetic announcing that the formula has been used outside its range of validity. Probabilities remain well defined even when moments do not:

t <- c(31, 58, 87, 119, 153, 192, 238, 296, 378, 520); n <- length(t)
xm <- min(t); alpha <- n / sum(log(t / xm))
round(c(sum_log = sum(log(t / xm)), alpha = alpha, P200 = (xm / 200)^alpha), 6)
if (alpha <= 1) message("alpha <= 1: the fitted Pareto has no mean")

  sum_log     alpha      P200 
16.038876  0.623485  0.312740 
alpha <= 1: the fitted Pareto has no mean

5. Result

\(\hat\alpha = 0.623\), \(\hat x_m = 31\): the fitted Pareto has no mean and no variance, and gives \(P(X > 200) = 0.313\).

Interpretation. The Pareto is the wrong family for this dataset, and the fitted \(\hat\alpha\) says so plainly. The data have a definite scale — failures cluster between 30 and 520 hours — whereas the Pareto describes something with no characteristic scale at all. A negative fitted mean is not a numerical accident to be rounded away; it is the model reporting its own inapplicability, and the right response is to keep the Weibull. This is the single most useful habit in distribution fitting: check that the fitted parameters lie in the region where the quantities you intend to report actually exist.

What the Practical Record Should Contain