Lesson 5 - Monte Carlo First¶
Objectives¶
After this lesson, you should be able to do the following:
Explain the difference between a model, a simulation, and an observable.
Distinguish a deterministic view from a stochastic view of the same physical system.
Use a random number as a computational coin for a yes/no event.
Describe the basic Monte Carlo pattern: one history, one tally, many histories, one estimate.
Read and write the logic for a simple neutron source-detector Monte Carlo simulation.
Separate parameters that change the physical model from parameters that change the sampling noise.
Why Monte Carlo first?¶
When I took on this course, I knew we’d need to introduce probability and statistics to support later courses. Rather than treat these as a separate language (like many math courses, in my experience), we’ll consider an application for which probability and statistics offer an enabling toolbox.
But that starts with the history of one neutron.
This neutron is born at a source. It moves through a material. It might bounce off a nucleus. It might be absorbed. It might escape the volume of interest. It might reach a detector. One neutron has one life. A detector response is not the life of one neutron. It is a statement about many possible neutron lives.
That is the core reason Monte Carlo simulation is fundamental to nuclear engineering. Particle transport is naturally a story of individual histories. We can often write equations for the average behavior, and those equations are important. But the event-by-event view is also powerful because it mirrors the physical narrative:
source particle
-> random flight or step
-> possible interaction
-> possible termination
-> tally
-> repeat many times
-> estimate response
In later courses on shielding and nuclear reactor physics, those histories will become very, very complex. A neutron will have an energy, a direction, a position, material-dependent cross sections, scattering kernels, and detector response functions. For now, we will build the smallest version of the idea.
This lesson is called Monte Carlo first because we will use randomness before formalizing all of probability. We will use the language of events, tallies, histories, and estimates. In the next lessons, we will give more formal names to the probability objects that are already doing work here.
Models, simulations, and observables¶
A model is a representation of a physical system. It is the collection of assumptions, balance equations, probabilistic rules, parameters, data, and numerical approximations used to describe what matters and how those pieces relate.
A model is not the real system. It is a chosen description of the real system.
A simulation is the execution of a model. It is what happens when the model is run to produce predictions, histories, fields, tallies, or other outputs.
The model is the representation. The simulation is the run.
Before building either one, we should identify the observable. The observable is the quantity we want to predict, estimate, or understand. Examples include transmitted intensity, remaining nuclei, detector counts, neutron flux, absorbed dose, leakage probability, and probability of detection.
This distinction matters because the same physical system can be modeled in more than one way. A model is a choice about what is kept, what is averaged, what is ignored, and what is sampled.
For example, suppose a beam travels through an attenuating slab. If the observable is only the transmitted fraction, then one model may describe a smooth intensity field. Another model may follow individual particles. Both can aim at the same observable.
The observable keeps us honest. Without it, we can add details forever and still not know what question we are answering.
Deterministic and stochastic views¶
A deterministic model has no random branching. Once the inputs are fixed, the output is fixed. If you run the same deterministic model twice with the same parameters, initial conditions, boundary conditions, and numerical settings, you get the same result.
A deterministic model often describes an averaged quantity. Things like temperature, pressure, and concentration are commonly treated this way. Many microscopic events have already been averaged into the input parameters.
A stochastic model includes random events. The same inputs do not define one path through the model. They define a collection of possible paths and how likely each path is.
A Monte Carlo simulation is a stochastic simulation that estimates an observable by sampling many random histories. One history is one possible realization of the model. The answer comes from many histories summarized together.
Despite their differences, deterministic and Monte Carlo simulation are often two routes to the same physical quantity.
The deterministic route averages first, then solves for a smooth response. The Monte Carlo route samples individual events first, then averages afterward.
A first physical example: attenuation through a slab¶
The physical system¶
A narrow beam enters an attenuating material of thickness \(L\). The material might be tinted glass, fog, murky water, shielding, or any medium that removes some of the beam as it travels.
The observable is the transmitted fraction:
We will describe the same system in two ways.
First, we will use a deterministic continuum model. Then, we will use a stochastic particle model.
Deterministic continuum model¶
Let \(I(x)\) be the beam intensity after traveling distance \(x\) into the material. Suppose the material removes a fixed fraction of the remaining intensity per unit length. Then the intensity satisfies
where \(\mu\) is the attenuation coefficient; such coefficients are almost always evaluated based on a combination of measurements and modeling.
The minus sign says the intensity decreases as the beam travels. The factor \(I\) says the amount removed over a small distance is proportional to how much remains.
With \(I(0)=I_0\), the solution is
At the far side of the slab,
so the transmitted fraction is
This model does not follow individual particles. It treats intensity as a smooth, continuous quantity. It is a continuum model and gives the average response directly.
Stochastic particle model¶
Now imagine launching identical particles into the slab one at a time. Each particle has a history. It either reaches the far side or it is removed inside the material.
For one particle, define a yes/no outcome:
A single particle does not give the transmitted fraction. A single particle gives only one yes/no result.
To connect this event-by-event picture to the continuum model, divide the slab into many thin layers of width \(\Delta x\). Suppose a particle entering one layer has probability approximately \(\mu\Delta x\) of being removed in that layer and probability \(1-\mu\Delta x\) of surviving to the next layer.
If the slab contains
layers, then the probability of surviving all layers is approximately
As the layers become thinner,
Thus the stochastic particle model and the deterministic attenuation equation predict the same transmitted fraction in the thin-layer limit.
The deterministic model gives \(T=e^{-\mu L}\) directly. The stochastic model estimates that same value by simulating many particle histories.
For \(N\) simulated particles,
The symbol \(T\) denotes the model quantity. The symbol \(\hat{T}\) denotes the estimate from one finite simulation.
That distinction is central to Monte Carlo methods:
\(T\) is what the model predicts in the ideal infinite-sampling limit.
\(\hat{T}\) is what we observe with a finite number of histories.
Random numbers as computational coins¶
A computer can generate pseudorandom numbers between 0 and 1. In Python, we often get them from a NumPy random number generator:
[1]:
import numpy as np
rng = np.random.default_rng()
U = rng.random()
The value \(U\) behaves like a random draw from the interval \([0,1)\). We might think of \(U\) as a computational, \(n\)-sided die whose possible face values \(v = i/n, i = 0, 1, 2, \ldots, n-1\) uniformly fill that interval as \(n\) grows very, very large.
Suppose an event should occur with probability \(p\), where \(0\le p\le 1\). Draw \(U\) between 0 and 1 and use the rule
If \(U<p\), say the event happens. Otherwise, say the event does not happen.
The interval from 0 to \(p\) occupies a fraction \(p\) of the unit interval. Therefore a uniformly selected number lands in that interval with probability \(p\).
In pseudocode:
draw U
if U < p:
event happens
else:
event does not happen
This tiny pattern is enough to build many physical simulations.
The random number does not create physics by itself. The model supplies the probability rule. The random number executes that rule for one history.
Probability is not observed frequency¶
A probability belongs to the model. An observed frequency belongs to a finite simulation run.
If the model says an event has probability \(p=0.3\), that does not mean exactly 3 out of every 10 simulated histories will show the event. It means that, over many independent histories, the observed fraction tends to hover near 0.3, with less wiggle as the number of histories increases.
This is one of the first places where Monte Carlo forces careful language.
Bad language:
The probability was 0.3, so 30% happened.
Sharper language:
The model probability was p = 0.3. In this finite simulation, the observed fraction was 0.287.
The model probability is not noisy. The estimate is noisy because we sampled a finite number of histories.
This distinction will become the basis for uncertainty quantification. We will not develop that theory fully in this lesson, but we should begin using the language correctly now.
The Monte Carlo pattern¶
The simplest Monte Carlo simulation has four pieces.
One history¶
A history is one complete realization of the model.
A history begins when a particle is born and ends when that particle is detected, absorbed, leaves the volume of interest, or otherwise terminates.
One tally¶
A tally records the outcome of interest from each history.
For a yes/no observable, the tally can be written as
Many histories¶
One history is not the answer. Monte Carlo simulation repeats the random experiment many times. Each history produces one tally.
One estimate¶
After \(N\) histories, the Monte Carlo estimate is the tally average:
For a yes/no event, this average is simply the fraction of histories in which the event happened.
The entire pattern can be written as
for i in 1, 2, ..., N:
simulate one history
compute that history's outcome
update the tally with the score
compute averages
A brief historical note¶
The phrase Monte Carlo is a reference to games of chance. The modern computational method became famous in the mid-twentieth century, especially through work by Stanislaw Ulam, John von Neumann, Nicholas Metropolis, and others during and after the wartime development of large-scale scientific computing.
Monte Carlo is not merely a way to add randomness to a problem. It is a way to compute quantities that are difficult to obtain by direct analysis, especially when the system has many possible paths.
Nuclear engineering is full of such systems. Neutron histories branch by scattering, absorption, leakage, fission, energy change, and geometry. Monte Carlo lets us build those branches into the simulation one event at a time.
A toy source-detector problem¶
Consider a deliberately simplified one-dimensional geometry. A particle (let’s call it a neutron) begins at the source (\(x = 0\)). To the right is a material region and then a detector boundary (at \(x = L\)). To the left is a leakage boundary. The neutron moves one fixed step at a time.
At each step, the neutron chooses a direction. If it reaches the detector boundary, it is detected. If it crosses the leakage boundary, it leaks out. If it remains inside the material, it may be absorbed or may scatter and continue.
The observable is the probability that a source neutron reaches the detector.
We will estimate that probability by simulating many neutron histories and counting the fraction that end in detection.
For neutron history \(i\), define
After \(N\) histories, estimate the detector response by
Here \(\hat{p}_d\) is the simulated detector response. It is the fraction of source neutron histories that ended in detection.
This model is intentionally crude. It does not include continuous flight distances, angular scattering distributions, energy changes, cross-section libraries, fission, or realistic detector physics. Those omissions are not mistakes. They are how we isolate the algorithmic structure.
Rules for one neutron history¶
A single history is a loop with three possible terminal outcomes: detection, leakage, or absorption.
Step |
Rule |
Result |
|---|---|---|
Start |
Place one neutron at the source position. |
The history is active. |
Choose direction |
Draw \(U\). If \(U<0.5\), move right; otherwise, move left. |
A fair computational coin sets the direction. |
Move |
Update the neutron position by one fixed step. |
The neutron is at a new position. |
Check detector |
If the neutron reaches the detector boundary, stop. |
Return tally 1. |
Check leakage |
If the neutron crosses the leakage boundary, stop. |
Return tally 0. |
Sample interaction |
If the neutron is still inside the material, draw another random number. |
Absorb with probability \(p_{\text{absorb}}\); otherwise scatter and continue. |
The order matters: move first, check boundaries immediately, and sample a material interaction only if the neutron remains inside the material. In debugging code, a maximum number of steps is useful as a guard against accidental infinite loops; that limit is not part of the physical model.
Flowchart before code¶
Monte Carlo algorithms branch. A flowchart helps make the branches visible before they become code.
For one neutron history, the logical structure is
start at source
choose direction
move one step
if detector reached:
return 1
if leakage boundary crossed:
return 0
if absorbed:
return 0
otherwise:
scatter and continue
Notice the order. We move first, then check boundaries, then sample the material interaction if the neutron is still inside the material.
That ordering is a model choice. Other physical models may require a different order. For example, a continuous-flight transport model samples a distance to collision, moves to the next collision or boundary, and then samples the collision outcome. We will eventually move in that direction. For this first model, fixed steps are simpler.
A good flowchart should show three kinds of logic:
actions, such as starting at the source or moving one step;
decisions, such as checking a boundary or sampling absorption;
terminal outcomes, such as detected, leaked, or absorbed.
If a history can loop, the loop should be visible. If a history can terminate, every terminal path should return a score to be tallied.
Python implementation¶
Here is one direct implementation of the toy model.
[ ]:
import numpy as np
def run_one_history(n_steps, p_absorb, rng=None, max_steps=10_000):
"""Simulate one simplified neutron history.
Parameters
----------
n_steps : int
Position of the detector boundary in "steps."
p_absorb : float
Probability of absorption after a step.
rng : numpy.random.Generator
Random number generator.
max_steps : int
Upper limit to prevent accidental infinite loops.
Returns
-------
int
1 if the neutron is detected; 0 otherwise.
"""
if rng is None:
rng = np.random.default_rng()
x = 0
for _ in range(max_steps):
# Choose direction with a fair computational coin.
if rng.random() < 0.5:
direction = 1
else:
direction = -1
# Move one fixed step.
x += direction
# Boundary checks come immediately after motion.
if x >= n_steps:
return 1
if x < 0:
return 0
# If still inside, sample absorption.
if rng.random() < p_absorb:
return 0
# Otherwise the neutron scattered and continues.
# It is good practice to report unexpected issues!
print("Warning: history did not terminate before max_steps!")
The function returns a single score. It does not return the whole path. That is because the observable is detector response. If the observable were path length, number of collisions, or energy deposited, the tally would be different.
Now repeat many histories:
[ ]:
def run_simulation(n_histories, n_steps, p_absorb, seed=None):
"""Estimate the detector response from many neutron histories."""
rng = np.random.default_rng(seed)
detected = 0
for _ in range(n_histories):
detected += run_one_history(
n_steps=n_steps,
p_absorb=p_absorb,
rng=rng,
)
return detected / n_histories
An example run might look like this:
[4]:
response = run_simulation(
n_histories=10_000,
n_steps=10,
p_absorb=0.1,
seed=12345,
)
print(response)
0.0059
If you change the seed, the answer may change slightly. That is not a bug. It is sampling noise.
If you increase n_histories, the estimate will usually wiggle less from run to run. That does not change the physical model. It changes the quality of the estimate.
Reading the output¶
The output of the simulation is an estimate:
It is tempting to say, “the detector probability is \(\hat{p}_d\).” A sharper statement is: “For this model and this finite simulation run, the estimated detector response was \(\hat{p}_d\).”
The actual model probability \(p_d\) is the ideal quantity we are trying to estimate. We do not know it exactly just because we ran a finite simulation.
This is not a philosophical detail. It affects how we compare simulations.
Suppose two simulation runs give 0.041 and 0.044. Is the second model really more detectable? Or did two finite Monte Carlo estimates simply wiggle? We cannot answer that without thinking about sampling uncertainty. That topic belongs to the later statistics lesson, but the habit starts here.
What changes the model? What changes the run?¶
Monte Carlo results can change for different reasons. We should separate them.
Model parameters¶
Model parameters define the physical or mathematical system being represented.
In the attenuation example, \(\mu\) and \(L\) are model parameters.
In the neutron toy model, model choices include
geometry;
detector boundary location;
leakage boundary location;
step size;
probability of absorption;
probability of moving left or right;
boundary rules;
interaction rules.
Changing any of these changes the model. It asks a different physical or mathematical question.
Sampling parameters¶
Sampling parameters change how the model is executed. The number of histories \(N\) and the random seed are sampling parameters.
Increasing \(N\) usually reduces sampling noise. It does not change the underlying model.
Changing the seed changes the particular random realization. It should not systematically change the model prediction. If changing the seed creates a large persistent shift, either the simulation is too noisy or the code has a problem.
What we have already used¶
We have already used several probability and statistics ideas, even though we have not fully formalized them.
Events: detection, absorption, leakage, scattering, and transmission are possible outcomes of a random process.
Indicators: a yes/no tally equals 1 when the event of interest happens and 0 otherwise.
Probability rules: every random decision is controlled by a rule that assigns likelihoods to possible outcomes; later we will call these distributions.
Estimates: a Monte Carlo result approximates a model quantity, such as a probability, average response, detector count rate, flux, or other observable.
Sampling noise: finite numbers of random histories make estimates wiggle from run to run; later we will quantify that wiggle.
Deterministic and Monte Carlo views of the same problem¶
Return to attenuation through a slab.
The deterministic model solves
and obtains
The stochastic model launches many particles and estimates the same transmitted fraction from the number that reach the far side.
These are not unrelated tricks. They are two ways of arranging the averaging.
The deterministic view says: average the beam into a smooth intensity and solve for the average response.
The Monte Carlo view says: simulate individual histories, tally the outcome of each history, then average the tallies.
The neutron source-detector model follows the same pattern. A deterministic transport method would try to describe the average particle population or flux directly. A Monte Carlo method follows many particles and estimates the detector response from their histories.
This is why we start with Monte Carlo. It makes the physical story visible.
Looking ahead¶
The next lessons will formalize the objects we have already used.
In the next lesson, we will name random variables and distributions. We will describe what rng.random() is sampling from, and we will define the probability rules behind random decisions.
After that, we will learn how to transform uniform random numbers into samples from other distributions. That will let us replace fixed neutron steps with sampled flight distances.
Finally, we will quantify Monte Carlo uncertainty. We will ask how much a finite estimate should wiggle and how to report that wiggle responsibly.