# How to use MMA to simulate the German tank problem

In the statistical theory of estimation, the German tank problem consists of estimating the maximum of a discrete uniform distribution from sampling without replacement. In simple terms, suppose there exists an unknown number of items which are sequentially numbered from 1 to N. A random sample of these items is taken and their sequence numbers observed; the problem is to estimate N from these observed numbers.

Assuming tanks are assigned sequential serial numbers starting with 1, suppose that four tanks are captured and that they have the serial numbers: 19, 40, 42 and 60. Let N equal the total number of tanks predicted to have been produced, 60 equal the highest serial number observed and 4 equal the number of tanks captured.

m = 60;(*m equal the highest serial number observed*)
k = 4;(*k equal the number of tanks captured*)
μ = (m - 1) (k - 1)/(k - 2)
σ = Sqrt[((k - 1) (m - 1) (m - k + 1))/((k - 3) (k - 2)^2)]
{μ - σ, μ + σ} // N


I want to use the Monte Carlo method to simulate the value of N with the help of MMA. If possible, find the 90% confidence interval for N.

n = 100000;
Table[Count[Table[RandomChoice[Range[m], 1], n] // Flatten,
60], {m, {60, 100, 140, 180, 220, 260}}]


Through the above method, it can be seen that when the number of tanks is 60, the probability of obtaining the maximum value of the observation is the largest, but this problem is still not well solved.

• Please show what you've tried so far. – JimB Jul 21 '20 at 0:33

## 2 Answers

You want to sample "without replacement" so it is RandomSample rather than RandomChoice.

Matching the usual notation is awkward with Mathematica because N is reserved so I'm taking the liberty to modify all of the notation.

Here is an example of how well the usual Frequentist estimator works for 60 tanks and a sample size of 4.

(* True number of tanks *)
nTanks = 60;
(* Sample size *)
sampleSize = 4;
(* Confidence level *)
confidenceLevel = 0.9;

(* Simulate n random samples of size sampleSize *)
n = 100000;
SeedRandom;
samples = RandomSample[Range[nTanks], sampleSize] & /@ Range[n];

(* Find point estimates and approximate confidence intervals *)
maxSample = Max[#] & /@ samples;
estimates = maxSample (1 + 1/sampleSize) - 1.;
lowerCLs = maxSample/((1 + confidenceLevel)/2)^(1/sampleSize);
upperCLs = maxSample/((1 - confidenceLevel)/2)^(1/sampleSize);

(* Estimated mean and standard deviation of estimator *)
Mean[estimates]
(* 59.9897 *)
StandardDeviation[estimates]
(* 11.9683 *)

(* Lower 5% and upper 95% quantiles:  central 90% of the observations *)
quantiles = Quantile[estimates, {0.05, 0.95}]
(* {36.5,74.} *)

(* What is the percent of confidence intervals that include the true number of tanks? *)
coverage =100*Length[Select[
Transpose[{lowerCLs, upperCLs}], #[] <= nTanks <= #[] &]]/n // N
(* 89.054 *)

• Thank you for your answer, but your result is obtained by applying the relevant calculation formula, I am more interested in the corresponding result obtained by using the violent simulation method. – A little mouse on the pampas Jul 21 '20 at 2:24
• Given the way you've phrased your request, I can't tell if there is a written language issue or that you don't understand sampling and what a confidence interval is. You can use the 100,000 values in estimates to look at the interval where the middle 90% of the simulation estimates are. But that interval is not a confidence interval. Is that the kind of interval you want? – JimB Jul 21 '20 at 3:29
• Yes, I need this code. I don’t require it to be consistent with the theoretical results. The interval of numerical simulation and the interval of theoretical calculation are similar to meet my needs. – A little mouse on the pampas Jul 21 '20 at 3:39
• Added in code for quantiles. – JimB Jul 21 '20 at 5:16
• You're complaining that the mean value Mean[estimates] with 100,000 simulations is 9.99 rather than 10 (the true number of tanks) ??? It would appear that the estimator is pretty much unbiased for the true number of tanks. – JimB Jul 21 '20 at 16:49

I've often thought it was more fun (and more pedagogical) to pose the Monte Carlo version of this problem in a pessimistic way: How often would this observation occur for some hypothetical "true" result? This is pessimistic, in the sense that it focuses on the probability of being wrong.

The first step is to simulate the possible observed maxima for some hypothetical number of total items:

simulatedObservations[hypotheticalNumber_, nObservations_, nTrials_] :=
Max /@ Table[
RandomSample[ Range[hypotheticalNumber], nObservations], {nTrials}]


Then use this function to compute the probability that the simulated hypothetical observations are consistent with your "real life" observation:

fractionOfMatchingObservations[observedMax_, nObservations_, nTrials_:10^4][hypotheticalNumber_] := With[
{observationsMatch =
Length@Select[LessEqualThan[observedMax]]@
simulatedObservations[hypotheticalNumber, nObservations, nTrials]},
observationsMatch/nTrials] // N


If the hypothetical true value is 60 tanks, then 100% of the simulated observations agree with your actual observation.

On the other hand, if the true value is 120 tanks, 5.5% of the simulated observations agree with your observed maximum of 60. (by fractionOfMatchingObservations[60, 4]). Do you still want to fight this battle?

One can of course generate a list of results by mapping over the hypothetical values to create a list of pairs:

results = {#, fractionOfMatchingObservations[60, 4][#]} & /@
Range[60, 1000];

ListLinePlot[results, PlotRange -> {{60, 200}, All}] This can be converted into an estimate of the expected value by first normalizing this into a CDF, and then computing the PDF by taking the derivative, and finally computing the expected value using the list of hypothetical true number & probability pairs

cdf[results_List] := With[
{hypotheticalNumber = results[[All, 1]],
normalization = Total@results[[All, 2]],
accumulated = Accumulate@results[[All, 2]]},
Transpose[{hypotheticalNumber, accumulated/normalization}]]

pdf[cdf_List] := With[
{hypotheticalNumber = Most@cdf[[All, 1]],
pdf = Differences[cdf[[All, 2]]]},
Transpose[{hypotheticalNumber, pdf}]]

expectedValue[pdf_List] := pdf[[All, 1]].pdf[[All, 2]]


Computing the result for our experiment:

expectedValue@pdf@cdf[results] (*84.3044*)


Note that this direct simulation converges very slowly with the maximum hypothetical number of tanks considered to the analytical expectation value result). For example, if one considers a hypothetical maximum of 200, the expectation value computed in this way is about 79. Showing this convergence is left as an exercise to the reader.

• But the theoretically calculated mean is 88.5, which cannot be reflected in your simulation. – A little mouse on the pampas Jul 21 '20 at 3:35
• While sampling with replacement (RandomChoice) is not "wrong", it is not how the problem is described in the Wiki page: "In the statistical theory of estimation, the German tank problem consists of estimating the maximum of a discrete uniform distribution from sampling without replacement." Also, generally the maximum is then modified (i.e., increased) to provide a "better" estimator. So the maximum isn't used as the estimate of the number of tanks. – JimB Jul 21 '20 at 3:37
• @Jim: Updated for the sampling without replacement – Joshua Schrier Jul 21 '20 at 13:51
• @Ordinaryusers68 Made expectation value explicit – Joshua Schrier Jul 21 '20 at 13:51