# Multivariate error modeling and uncertainty quantification using importance
(re-)weighting for Monte Carlo simulations in particle transport

Pia Stammer  Lucas Burigo  Oliver Jäkel  Martin Frank  Niklas Wahl  Karlsruhe
Institute of Technology, Steinbuch Centre for Computing, Karlsruhe, Germany
German Cancer Research Center - DKFZ, Department of Medical Physics in
Radiation Oncology, Heidelberg, Germany  HIDSS4Health - Helmholtz Information
and Data Science School for Health, Karlsruhe/Heidelberg, Germany  Heidelberg
Institute for Radiation Oncology (HIRO), Heidelberg, Germany  Heidelberg Ion
Beam Therapy Center — HIT, Department of Medical Physics in Radiation
Oncology, Heidelberg, Germany

###### Abstract

Fast and accurate predictions of uncertainties in the computed dose are
crucial for the determination of robust treatment plans in radiation therapy.
This requires the solution of particle transport problems with uncertain
parameters or initial conditions. Monte Carlo methods are often used to solve
transport problems especially for applications which require high accuracy. In
these cases, common non-intrusive solution strategies that involve repeated
simulations of the problem at different points in the parameter space quickly
become infeasible due to their long run-times. Intrusive methods however limit
the usability in combination with proprietary simulation engines. In [51], we
demonstrated the application of a new non-intrusive uncertainty quantification
approach for Monte Carlo simulations in proton dose calculations with normally
distributed errors on realistic patient data. In this paper, we introduce a
generalized formulation and focus on a more in-depth theoretical analysis of
this method concerning bias, error and convergence of the estimates. The
multivariate input model of the proposed approach further supports almost
arbitrary error correlation models. We demonstrate how this framework can be
used to model and efficiently quantify complex auto-correlated and time-
dependent errors.

##  1 Introduction

The linear Boltzmann equation has been used to model transport processes in a
large variety of applications [e.g. 44, 50, 10, 3]. Computational methods in
radiation therapy require a discretization of the seven-dimensional phase
space. These discretizations can be grouped into deterministic and stochastic
approaches. Deterministic approaches [3, 13], such as the spherical harmonics
($P_{N}$) [9] and the discrete ordinates ($S_{N}$) [56, 16] method for the
angular discretization, are computationally efficient. However, stochastic
Monte Carlo methods are increasingly preferred for radiation therapy due to
their accuracy and flexibility [58, 23]. They also have the advantage of a
known error estimate, an intuitive physical interpretation and being well-
suited for parallel implementation. Under the assumption of exact geometry,
this theoretically enables an arbitrarily exact approximation at the cost of
longer run-times. However, in radiation therapy, human error as well as lack
of information can cause deviations of the treatment reality from the model.
Therefore, in addition to the physical distribution of phase space parameters,
uncertainties also have a significant impact on the outcome of the simulations
[32, 27].

When using the Monte Carlo method to solve the transport problem, quantifying
the impact of these uncertainties can be challenging, as many UQ methods, such
as (quasi-) Monte Carlo [5] or stochastic collocation [61], require numerous
repeated runs of the solver at different points in the uncertain space.
Further, the incorporation and sampling of multivariate error correlation
models, which are essential for a realistic representation of certain error
types, can be extremely tedious in scenario-based methods. Most other
approaches either cannot provide the required accuracy [e.g. 14, 6], are only
aimed at deterministic solution methods [e.g. 1, 42] or are intrusive [e.g.
47], which is not realizable with the often proprietary softwares used for
Monte Carlo simulations.

We introduce a non-intrusive method, that generates estimates of the expected
value and variance of radiation therapy doses from a single MC simulation run.
The approach is based on the concept of importance sampling [24, 18] and
mimics the computation of scenarios in the parameter space by (re-)weighting
the results of a single realization. The computational overhead of a complex
simulation is thus reduced to that of the (re-)weighting step, effectively
making the uncertainty quantification problem a scoring problem independent
from the underlying physical simulation.

Similar strategies have been previously employed, e.g. to estimate the average
of samples from different distributions [2, 53], for more efficient updates in
optimization [48] and reinforcement learning [45, 17] or to artificially
increase the sample size in costly simulations [52]. However, to the best of
our knowledge there is no previous research on computing a conditional
variance by (re-)weighting a single original sample. Also, in the field of
medical physics and in particular dose calculation, we are not aware of any
previous use of importance (re-)weighting for the reuse of samples, neither
for dose estimates (which are averages) nor for the variance with respect to
input uncertainties. Note that a detailed application of the proposed approach
to intensity modulated proton therapy (IMPT) treatment plans and several
patient cases from the open CORT (common optimization for radiation therapy)
data set [8] is described in our previous publication [51].

In the following, we first introduce the problem of particle transport with
uncertain input in the context of radiation therapy (section 2), then show how
to use importance sampling to quantify uncertainties (section 4), discuss
mathematical properties of this approach (section 5) as well as different
models for the input uncertainties (section 6), and lastly show results for a
water phantom as well as a prostate patient (section 8).

##  2 Solution of the deterministic linear Boltzmann equation using Monte
Carlo

We are interested in particle transport for an application in radiation
therapy. The behavior of particles can be described by the Boltzmann transport
equation. Here we focus on protons, which are especially sensitive to
uncertainties due to the sharp peak in their depth-dose curve [32, 30].

Let $\mathcal{D}_{r}$ be a bounded region containing the voxelized patient
anatomy and $\Phi(\boldsymbol{r},\boldsymbol{\Omega},E)$ the proton fluence at
location $\boldsymbol{r}\in\mathcal{D}_{r}\subset\mathbb{R}^{3}$, with unit
velocity direction $\boldsymbol{\Omega}\in\mathbb{S}^{2}$ and energy
$E\in\mathbb{R}_{+}$. The linear Boltzmann equation reads

$\displaystyle\Omega\cdot\nabla\Phi(\boldsymbol{r},\boldsymbol{\Omega},E)+\sigma_{t}(\boldsymbol{r},E)\Phi(\boldsymbol{r},\boldsymbol{\Omega},E)$
$\displaystyle=Q^{Sca}(\boldsymbol{r},\boldsymbol{\Omega},E)+S(\boldsymbol{r},\boldsymbol{\Omega},E)\;,$
(1)

where $Q^{Sca}$ is the energy loss and scattering component and is given as

$Q^{Sca}(\boldsymbol{r},\boldsymbol{\Omega},E)=\int_{0}^{\infty}\int_{\boldsymbol{\Omega^{\prime}}\in
S^{2}}\sigma_{s}(\boldsymbol{r},\boldsymbol{\Omega}\cdot\boldsymbol{\Omega^{\prime}},E\rightarrow
E^{\prime})\Phi(\boldsymbol{r},\boldsymbol{\Omega^{\prime}},E^{\prime})d\boldsymbol{\Omega^{\prime}}dE^{\prime}\;.$
(2)

Note that the total and differential cross sections $\sigma_{t}$ and
$\sigma_{s}$ are material and particle dependent and
$S(\boldsymbol{r},\boldsymbol{\Omega},E)$ is the source function, which
characterizes the spatial, angular, and energy distribution of the radiation
source. The dose can then be determined from the fluence as

$D(\boldsymbol{r})=\frac{1}{\rho(\boldsymbol{r})}\int_{0}^{\infty}\int_{\boldsymbol{\Omega}\in
S^{2}}L(\boldsymbol{r},E)\Phi(\boldsymbol{r},\boldsymbol{\Omega},E)d\boldsymbol{\Omega}dE\;.$
(3)

Here, $L(\boldsymbol{r},E)$ is the stopping power, i.e., the rate at which
particles lose energy and $\rho(\boldsymbol{r})$ is the tissue density at
position $\boldsymbol{r}$.

###  2.1 The Monte Carlo method

The Monte Carlo method is a numerical integration method, which is based on
evaluating the integrand at random positions sampled from the probability
distribution of the integration variable. For realizations
$\boldsymbol{x_{i}}\,,i=1,...,N$ of the random variable
$\boldsymbol{X}\in\mathbb{R}^{d}$ with probability density function
$p:\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$, the estimator is defined as
follows:

$I(f)=E_{p}[f(\boldsymbol{X})]=\int_{\mathbb{R}^{d}}f(\boldsymbol{X})p(\boldsymbol{X})d\boldsymbol{X}\approx\frac{1}{N}\sum_{i=1}^{N}f(\boldsymbol{x_{i}}):=I^{N}(f)\;.$
(4)

Note, that with the Strong Law of Large Numbers, one can show

$P\left(\lim_{N\rightarrow\infty}I^{N}(f)\,=I(f)\right)=1\;,$ (5)

i.e., the Monte Carlo estimator converges to the true value in probability.
Further, the independence of the random variables $\boldsymbol{x_{i}}$ and the
Central Limit Theorem provide a formula for the variance of the estimator:

$Var(I^{N}(f))=\frac{v_{f}^{2}}{N}\;,$ (6)

which depends solely upon the variance $v_{f}^{2}$ of the integrand and the
number of samples $N$. The estimator thus converges with
$O(\frac{1}{\sqrt{N}})$ and its accuracy can be controlled through the number
of computed realizations.

###  2.2 Randomized quasi-Monte Carlo

Due to the slow convergence and therefore often large sample sizes required
for Monte Carlo estimates, numerous variations of the method exist which
reduce the variance of the estimator.
The quasi-Monte Carlo method replaces the randomly sampled points with a set
of deterministic low discrepancy points
$P_{N}:=\\{\boldsymbol{u_{1}},...,\boldsymbol{u_{N}}\\}$, which cover the
domain more uniformly [5]:

$I_{qMC}^{N}(f):=\frac{1}{N}\sum_{i=1}^{N}f(\boldsymbol{u_{i}})\;.$ (7)

This can for example be achieved using grid points or number theoretic
sequences such as the Sobol sequence [49], used later in this paper. These
typically result in a faster convergence, however it is less straightforward
to compute the error compared to the standard Monte Carlo method [29].
Randomized quasi-Monte Carlo therefore reintroduces randomness to the method,
while retaining the convergence properties of quasi-MC. The resulting point
set is a low discrepancy sequence, in which each individual point is uniformly
distributed [29]. In the case of a grid rule, this can be achieved e.g.
rotating the lattice or randomly drawing points within the grid cells from a
uniform distribution. For number theoretic sequences, different random
shuffling and scrambling111Shuffling refers to a randomization of the sample
index, while scrambling refers to the randomization of the sample value.
strategies have been developed [37, 35, 28, 26]. We use a Sobol sequence with
linear scrambling and a random digit shift according to [35]. For square-
integrable functions, the resulting randomized quasi-Monte Carlo approach
preserves the properties of Sobol sequences, ensuring a variance of
$O\left(\frac{\log(N)^{d}}{N}\right)$ for $N$ $d$-dimensional sample points,
while also generating an unbiased estimate and allowing for a statistical
estimation of the integration error [37, 38].

###  2.3 The Boltzmann equation in integral form

In order to solve equation 1 using a variant of the Monte Carlo method, it
needs to be given in its integral form. We first execute a change of
variables, where $\boldsymbol{r}$ is parametrized with respect to a reference
location $\boldsymbol{r_{0}}$

$\displaystyle\boldsymbol{r}=\boldsymbol{r_{0}}+s\boldsymbol{\Omega}\;.$ (8)

For all combinations of $\boldsymbol{r^{\prime}}$ and
$\boldsymbol{\Omega^{\prime}}$, such a reference location
$\boldsymbol{r_{0}^{\prime}}$ is determined, s.t.
$\boldsymbol{r_{0}^{\prime}}+s\boldsymbol{\Omega^{\prime}}$:

$\displaystyle
Q_{Sca}(\boldsymbol{r},\boldsymbol{\Omega},E)=\int_{0}^{\infty}\int_{\boldsymbol{\Omega^{\prime}}\in
S^{2}}\sigma_{s}(\boldsymbol{r_{0}^{\prime}}+s\boldsymbol{\Omega^{\prime}},\boldsymbol{\Omega}\cdot\boldsymbol{\Omega^{\prime}},E\rightarrow
E^{\prime})\Phi(\boldsymbol{r_{0}^{\prime}}+s\boldsymbol{\Omega^{\prime}},\boldsymbol{\Omega^{\prime}},E^{\prime})d\boldsymbol{\Omega^{\prime}}dE^{\prime}\;,$
(9)
$\displaystyle\frac{d}{ds}\Phi(\boldsymbol{r_{0}}+s\boldsymbol{\Omega},\boldsymbol{\Omega},E)+\sigma_{t}(\boldsymbol{r_{0}}+s\boldsymbol{\Omega},E)\Phi(\boldsymbol{r_{0}}+s\boldsymbol{\Omega},\boldsymbol{\Omega},E)=Q^{Sca}(\boldsymbol{r_{0}}+s\boldsymbol{\Omega},\boldsymbol{\Omega},E)+S(\boldsymbol{r_{0}}+s\boldsymbol{\Omega},\boldsymbol{\Omega},E)\;.$
(10)

For simplicity we now rename
$\boldsymbol{r_{0}}+s\boldsymbol{\Omega}\rightarrow s$,

$\frac{d}{ds}\Phi(s,\boldsymbol{\Omega},E)+\sigma_{t}(s,E)\Phi(s,\boldsymbol{\Omega},E)=Q^{Sca}(s,\boldsymbol{\Omega},E)+S(s,\boldsymbol{\Omega},E)\;.$
(11)

Lastly we multiply with the integrating factor
$\exp(-\int_{0}^{s}\sigma_{t}(s^{\prime},E)ds^{\prime})$ and integrate both
sides to obtain

$\displaystyle\Phi(s,\boldsymbol{\Omega},E)=$
$\displaystyle\;\Phi(\boldsymbol{r_{0}},\boldsymbol{\Omega},E)\exp\left(-\int_{0}^{s}\sigma_{t}(s^{\prime},E)ds^{\prime}\right)$
$\displaystyle+\int_{0}^{s}\left(Q^{Sca}(s^{\prime\prime},\boldsymbol{\Omega},E)+S(s^{\prime\prime},\boldsymbol{\Omega},E)\right)\,\exp\left(-\int_{s^{\prime\prime}}^{s}\sigma_{t}(s^{\prime},E)ds^{\prime}\right)ds^{\prime\prime}\;.$
(12)

Now the transport problem can be solved by sampling particles with phase space
coordinates $(\boldsymbol{r},\boldsymbol{\Omega},E)$ from the source function
$S(\boldsymbol{r},\boldsymbol{\Omega},E)$ and simulating their trajectories
and random interactions according to the respective probability distributions.
This corresponds to interpreting section 2.3 in the sense of equation 4 while
choosing the probability density function $p$ to represent the respective
particle source and physics. For the most accurate solution, this simulation
is done event by event for each individual particle, always sampling the
distance to the next collision, the type of collision, energy loss and
scattering angle [55]. However, in order to increase computation speed,
different variations of this method exist, such as condensed history
algorithms which only explicitly simulate so called catastrophic collision
events, while accounting for soft collisions collectively after a given
traveled distance. We will not elaborate this further here, the interested
reader might take a look at [55] or [47] for a more concise derivation of the
probability measures, the resulting solution algorithm and different Monte
Carlo variations.

When considering non-intrusive UQ methods, the Monte Carlo dose calculation
algorithm can be treated as a black box
$BB(\boldsymbol{r},\boldsymbol{\Omega},E)$, returning the dose
$D(\boldsymbol{r})$ at position $\boldsymbol{r}$ for the input vector
$\boldsymbol{X}=(\boldsymbol{r},\boldsymbol{\Omega},E)\sim S_{0}$ on
$\mathcal{D}_{S_{0}}$.

$D(\boldsymbol{r})=E_{S_{0}}[BB(\boldsymbol{r},\boldsymbol{\Omega},E)]=\int_{\mathcal{D}_{S_{0}}}BB(\boldsymbol{X})S_{0}(\boldsymbol{X})d\boldsymbol{X}\approx\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\;,$
(13)

where the realizations $\boldsymbol{x_{1}},...,\boldsymbol{x_{H}}$ are sampled
from $S_{0}(\boldsymbol{X})$.

Note, that the black box for fixed realizations of the input parameters still
describes a stochastic process, due to subsequent sampling of the path and
interactions of the particles. The input sampled from the source distribution
merely fixes the initial properties of released particles. Each response
$BB(\boldsymbol{r},\boldsymbol{\Omega},E)$ however also implicitly includes
realizations of the probability distributions governing the interactions of
these particles. Thus, the Monte Carlo method is applied to the complete
random particle trajectory, we however omit this, since we are only interested
in uncertainties affecting the phase space at a superficial level in the
source distribution.

##  3 Introducing uncertainties

We now introduce a random error vector $\boldsymbol{Z}$ with density $p_{Z}$
to the problem equation 1. The corresponding uncertain linear Boltzmann
equation takes the following form.

$\displaystyle\Omega\cdot\nabla\Phi(\boldsymbol{r},\boldsymbol{\Omega},E,\boldsymbol{Z})+\sigma_{t}(\boldsymbol{r},E,\boldsymbol{Z})\Phi(\boldsymbol{r},\boldsymbol{\Omega},E,\boldsymbol{Z})$
$\displaystyle=Q^{Sca}(\boldsymbol{r},\boldsymbol{\Omega},E,\boldsymbol{Z})+S(\boldsymbol{r},\boldsymbol{\Omega},E,\boldsymbol{Z})\;,$
(14)

We are interested in the expected value and variance of the dose estimate with
respect to uncertainties, which occur during radiation treatments.

###  3.1 Uncertainties in radiation therapy

Among the most important sources of uncertainties in proton therapy are such
that manifest in the patient set-up and particle range [41, 30, 42, 32, 33].
Set-up uncertainties directly affect the particle source $S_{0}$, by shifting
the initial particle positions $\boldsymbol{r}$ by a random amount
$\boldsymbol{Z_{r}}$. Range uncertainties are caused by a variety of factors
concerning the patient density, ranging from CT conversion errors to imaging
artifacts or changes in the patient geometry [54, 40, 32, 36]. Therefore, they
manifest in the material dependent quantities, such as the scattering cross
sections $\sigma_{s}$ and $\sigma_{t}$, as well as the density $\rho$ and the
stopping power $L$.

In the following, the method will be introduced for a general random error
vector $\boldsymbol{Z}$ which affects the initial particle parameters
$\boldsymbol{X}(\boldsymbol{Z})$ and the corresponding phase space
distribution. We assume that the conditional distribution
$S_{0,\boldsymbol{X}|\boldsymbol{Z}}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z})$
of $\boldsymbol{X}$ for a fixed error realization is known. An application of
the proposed method to set-up uncertainties will be discussed in section 6.

###  3.2 Solution of the uncertain problem using Monte Carlo

Most non-intrusive strategies involve computing the dose at different points
in the parameter space. In case of the Monte Carlo method, these points are
realizations randomly drawn from the probability distribution $p_{Z}$ of the
uncertain variable $\boldsymbol{Z}$.

Note that the introduction of an uncertain parameter adds a second stochastic
process around the black box transport solver. On the one hand, the dose
distribution depends on random input parameters following a distribution
defined by the physical radiation source. On the other hand, this distribution
now depends on the random error vector $\boldsymbol{Z}$. Thus, if we want to
compute the expected value $\mathcal{D}$ of the dose with respect to both
processes, we are looking for

$\mathcal{D}(\boldsymbol{r})=E_{p_{Z}}\left[E_{S_{0,\boldsymbol{X}|\boldsymbol{Z}}}\left[BB(\boldsymbol{X}(\boldsymbol{Z}))|\boldsymbol{Z}=\boldsymbol{z}\right]\right]=\int_{\mathcal{D}_{p_{Z}}}\int_{\mathcal{D}_{S_{0}}}BB(\boldsymbol{X}(\boldsymbol{Z}))\cdot
S_{0,\boldsymbol{X}|\boldsymbol{Z}}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z})\,d\boldsymbol{X}\cdot
p_{Z}(\boldsymbol{Z})\,d\boldsymbol{Z}$ (15)

The Monte Carlo method can now be applied to both random processes. Starting
with the error variable, for a fixed realization, the black box MC solver can
be applied analogously to the original problem. Thus, the dose estimate
$D_{i}(\boldsymbol{r})$ for realization $\boldsymbol{z_{i}}$ and the
corresponding conditional distribution
$S_{i}=S_{0,\boldsymbol{X}|\boldsymbol{Z}}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})$
is given by

$D_{i}(\boldsymbol{r})=E_{S_{i}}[BB(\boldsymbol{X}(\boldsymbol{Z}))|\boldsymbol{Z}=\boldsymbol{z_{i}}]=\int_{\mathcal{D}_{S_{i}}}BB(\boldsymbol{X}(\boldsymbol{z_{i}}))S_{i}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})d\boldsymbol{X}\approx\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}^{i}})$
(16)

with $\boldsymbol{x^{i}_{1}},...,\boldsymbol{x^{i}_{H}}$ sampled from $S_{i}$.
The Monte Carlo estimate of the expected value with respect to both processes
is then computed as the mean over
$D_{1}(\boldsymbol{r}),...,D_{N}(\boldsymbol{r})$ for
$\boldsymbol{z_{1}},...,\boldsymbol{z_{N}}$ sampled from $p_{Z}$

$\displaystyle\mathcal{D}(\boldsymbol{r})\approx\frac{1}{N}\sum_{i=1}^{N}D_{i}(\boldsymbol{r})\approx\frac{1}{N}\sum_{i=1}^{N}\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}^{i}})\;,$
(17)

and the variance as the sample variance, defined as

$S_{N}^{2}=\frac{1}{N-1}\sum_{i=1}^{N}(\boldsymbol{x_{i}}-\bar{\boldsymbol{x}})^{2}\,\text{
where }\bar{\boldsymbol{x}}=\frac{1}{N}\sum_{i=1}^{N}\boldsymbol{x_{i}}\,,$
(18)

for $\boldsymbol{x_{i}},\;i=1,...,N$ that are independent and identically
distributed. This is an unbiased estimator of the variance of the random
variable $\boldsymbol{X}$. Applied to the problem at hand, this results in the
following formula for the dose variance:

$\displaystyle Var_{p_{\boldsymbol{Z}}}(D(\boldsymbol{r},\boldsymbol{Z}))$
$\displaystyle\approx\frac{1}{N-1}\sum_{i=1}^{N}\left(D_{i}(\boldsymbol{r})-\mathcal{D}(\boldsymbol{r})\right)^{2}$
$\displaystyle\approx\frac{1}{N-1}\sum_{i=1}^{N}\left(\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}^{i}})-\frac{1}{N}\sum_{j=1}^{N}\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}^{j}})\right)^{2}\;.$
(19)

###  3.3 Direct computation of the expected value

When not only the conditional or marginal distributions of $\boldsymbol{X}$
and $\boldsymbol{Z}$, but also their joint distribution
$\mathcal{S}(\boldsymbol{X}(\boldsymbol{Z}))$ is known, it is possible to
compute the expected dose directly

$E_{p_{\boldsymbol{Z}}}[D(\boldsymbol{r},\boldsymbol{Z})]=\mathbb{E}_{\mathcal{S}}[BB(\boldsymbol{X}(\boldsymbol{Z}))]\approx\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\;,$
(20)

with $\boldsymbol{x}_{p},\;p=1,..,H$ sampled from $\mathcal{S}$.

##  4 (Re-)weighting Monte Carlo samples for uncertainty quantification

In the following section, we introduce the proposed (re-)weighting approach.
We exploit the fact, that the input vector $\boldsymbol{X}$ is non-
deterministic, i.e., has the probability density $S_{0}(\boldsymbol{X})$
regardless of uncertainties. Uncertainties usually manifest in a way that
merely changes this distribution. The concept of importance sampling can be
used to gain insight into dose realizations for different input distributions
$S_{i}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})$ corresponding to
error realizations $\boldsymbol{z}_{i},\,i=1,...,N$ from $p_{Z}$, from one
initial simulation.

We therefore now recall the basic principles of importance sampling.

###  4.1 Standard importance sampling

Importance sampling is a method originally introduced by [24] and later made
popular in a more general context by [18].

Let’s assume we are interested in $g(\boldsymbol{X})$. Let $\boldsymbol{X}$
have the probability density function $p(\boldsymbol{X})$ with support
$\mathcal{D}_{p}$, i.e., $\boldsymbol{X}\sim p(\boldsymbol{X})$. Then

$\int_{\mathcal{D}_{p}}g(\boldsymbol{X})p(\boldsymbol{X})d\boldsymbol{X}=\int_{\mathcal{D}_{p}}g(\boldsymbol{X})\frac{p(\boldsymbol{X})}{q(\boldsymbol{X})}q(\boldsymbol{X})d\boldsymbol{X}\;.$
(21)

So the integral can be approximated as

$\int_{\mathcal{Q}}g(\boldsymbol{X})\frac{p(\boldsymbol{X})}{q(\boldsymbol{X})}q(\boldsymbol{X})d\boldsymbol{X}\approx\sum_{i=1}^{N}g(\boldsymbol{x_{i}})\frac{p(\boldsymbol{x_{i}})}{q(\boldsymbol{x_{i}})}\;.$
(22)

Where $\boldsymbol{x_{i}}$, $i=1,...,N$ are realizations of
$\boldsymbol{X}\sim q(\boldsymbol{X})$ and $q$ is a distribution with support
$\mathcal{Q}\supseteq\mathcal{D}_{p}$.

Thus, the computation of $g(\boldsymbol{X})$ for realizations of
$p(\boldsymbol{X})$, can be replaced by computations of $g(\boldsymbol{X})$
with realizations $\boldsymbol{x_{i}}$, $i=1,...,N$ from $q(\boldsymbol{X})$
using additional weights $\frac{p(\boldsymbol{x_{i}})}{q(\boldsymbol{x_{i}})}$
in the aggregation of the results. This can be convenient if it is easier to
sample from $q(\boldsymbol{X})$ than $p(\boldsymbol{X})$ or because choosing
$q(\boldsymbol{X})$ cleverly can reduce the variance of the estimate. It can
also be helpful, when evaluating $g$ is much more expensive than evaluating
$\frac{p(\boldsymbol{x_{i}})}{q(\boldsymbol{x_{i}})}$, for example because $g$
is an expensive simulation. In this case, if we are interested in
$g(\boldsymbol{X})$ for different distributions of $\boldsymbol{X}$, we can
reuse the realizations $g(\boldsymbol{x_{1}}),...,g(\boldsymbol{x_{N}})$.

###  4.2 Importance (re-)weighting

For our purpose, this can be applied to mimic dose computations at different
points in the parameter space from a previously computed sample of particle
trajectories. In the following, we will assume random samples, the method can
however be analogously applied to quadratures or quasi-random numbers.

Instead of solving the transport problem $N$ times with slightly different
input distributions $S_{1},...,S_{N}$ (compare equation 16), it is solved once
and the results for other error realizations are determined by (re-)weighting
the deposited dose of this one sample of particles as explained in section
4.1. For an error realization $\boldsymbol{z_{i}}$ from $p_{Z}$, initial
parameter realizations $\boldsymbol{x_{1}},...,\boldsymbol{x_{H}}$ from the
sampling distribution $q(\boldsymbol{X})$ and target distributions
$S_{i}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})$, the dose is given
by

$\displaystyle\mathcal{D}_{i}(\boldsymbol{r})$
$\displaystyle=\int_{\mathcal{D}_{S_{i}}}BB(\boldsymbol{X}(\boldsymbol{z_{i}}))S_{i}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})d\boldsymbol{X}=\int_{\mathcal{D}_{S_{i}}}\frac{BB(\boldsymbol{X}(\boldsymbol{z_{i}}))\;S_{i}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{X})}\;\cdot
q(\boldsymbol{X})d\boldsymbol{X}$
$\displaystyle\approx\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{x_{p}})}\;.$
(23)

This yields the following procedure:

Algorithm 1 Calculate doses $\mathcal{D}_{S_{i}}$

(a) On-the-fly

for particles p=1:H do

Sample $\boldsymbol{x_{p}}\leftarrow q(\boldsymbol{X})$

Compute $BB(\boldsymbol{x_{p}})$

Sample $\boldsymbol{z_{i}}\leftarrow p_{Z}(\boldsymbol{Z})$

for target distributions
$S_{i}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})\;,\newline i=1,...,N$
do

$\mathcal{D}_{S_{i}}\mathrel{+}=\frac{1}{H}\frac{BB(\boldsymbol{x_{p}})\;S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{x_{p}})}$

end for

end for

(b) In post-processing

for particles p=1:H do

Sample $\boldsymbol{x_{p}}\leftarrow q(\boldsymbol{X})$

Compute $BB(\boldsymbol{x_{p}})$

end for

Sample $\boldsymbol{z_{i}}\leftarrow p_{Z}(\boldsymbol{Z})$

for target distributions
$S_{i}(\boldsymbol{X}|\boldsymbol{Z}=\boldsymbol{z_{i}})\;,\newline i=1,...,N$
do

$\mathcal{D}_{S_{i}}\mathrel{=}\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{x_{p}})}$

end for

The importance and target distributions have to be chosen according to the
quantities of interest. Here we show the derivation for the nominal dose
without uncertainty, the expected dose and the dose variance with respect to
uncertainties.

###  4.3 Nominal dose

The nominal dose is the solution to the problem introduced in section 2
without uncertainty, i.e., for the input distribution $S_{0}(\boldsymbol{X})$.
Let $q(\boldsymbol{X})$ be the simulated particle distribution, then the
nominal dose can be computed by replacing $p(\boldsymbol{X})$ with
$S_{0}(\boldsymbol{X})$ in section 4.1.

$\displaystyle D(\boldsymbol{r})$
$\displaystyle=E_{S_{0}}[BB(\boldsymbol{X})]=\int_{\mathcal{D}_{S_{0}}}BB(\boldsymbol{X})S_{0}(\boldsymbol{X})d\boldsymbol{X}=\int_{\mathcal{D}_{S_{0}}}\frac{BB(\boldsymbol{X})S_{0}(\boldsymbol{X})}{q(\boldsymbol{X})}\;\cdot
q(\boldsymbol{X})d\boldsymbol{X}$ (24)
$\displaystyle\approx\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;S_{0}(\boldsymbol{x_{p}})}{q(\boldsymbol{x_{p}})}$

This reduces to the regular dose computation with Monte Carlo described in
section 2, when $q$ is chosen to be $S_{0}$.

###  4.4 Expected value

Using section 4.2 to compute the dose for individual error realizations, we
can mimic the Monte Carlo method for uncertainties and compute the expected
dose as

$\displaystyle\mathcal{D}(\boldsymbol{r})\approx\frac{1}{N}\sum_{i=1}^{N}D_{i}(\boldsymbol{r})\approx\frac{1}{N}\sum_{i=1}^{N}\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{x_{p}})}\;,$
(25)

with $\boldsymbol{z_{i}},\,i=1,...,N$ sampled from $S_{i}$ and
$BB(\boldsymbol{x_{p}})$ again the pre-computed particle trajectories with
input parameters $\boldsymbol{x_{p}},\,p=1,...,H$ sampled from $q$.

For a known joint distribution of the uncertain input, it is possible to
directly compute the expected value, without sampling from $p_{Z}$, by
adapting equation 20

$\mathcal{D}(\boldsymbol{r})=\mathbb{E}_{\mathcal{S}}[BB(\boldsymbol{X}(\boldsymbol{Z}))]\approx\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;\mathcal{S}(\boldsymbol{x_{p}})}{q(\boldsymbol{x_{p}})}\;,$
(26)

where $\boldsymbol{x_{p}},\;p=1,...,H$ sampled from $\mathcal{S}$.

###  4.5 Variance

The variance can also be computed for error realizations
$\boldsymbol{z_{1}},...,\boldsymbol{z_{N}}$ from $p_{Z}$, using dose
realizations from equation 25 with the sample variance formula

$\sigma^{2}=Var_{p_{\boldsymbol{Z}}}(D(\boldsymbol{r},\boldsymbol{Z}))\approx\frac{1}{N-1}\sum_{i=1}^{N}\left(D_{i}(\boldsymbol{r})-\mathcal{D}(\boldsymbol{r})\right)^{2}\;.$
(27)

###  4.6 Choosing the importance distribution

The optimal choice of $q(\boldsymbol{X})$ depends largely on the target
quantities. It follows from sections 4.3, 4.4 and 4.5 that the nominal dose,
expected dose and the error scenarios required for the variance computation
correspond to $N+2$ different distributions of the input parameters. Thus,
when all of these quantities are computed from the same set of realizations,
$q(\boldsymbol{X})$ needs to be suitable for not just one but $N+2$ different
target distributions. Therefore it will most likely not be possible to choose
a function which is optimal for each individual computation and the choice
depends on how the accuracy of the different target quantities is prioritized.

For a given target distribution $p(\boldsymbol{X})$ there are two criteria, by
which $q(\boldsymbol{X})$ should be chosen (see e.g. [39]):

  * 1.

Let $Q=\\{\boldsymbol{X}|q(\boldsymbol{X})>0\\},\text{ if
}\;p(\boldsymbol{X})\cdot BB(\boldsymbol{X})\neq 0\text{ then
}\boldsymbol{X}\in Q$

  * 2.

To minimize the standard error of the approximation
$\sigma_{q}^{2}=\frac{1}{H}\sum_{p=1}^{H}\left[\frac{p(\boldsymbol{x_{p}})\cdot
BB(\boldsymbol{x_{p}})}{q(\boldsymbol{x_{p}})}-\mu\right]^{2}$ with
$\mu=E_{S_{0}}[BB(\boldsymbol{X})]$, $q(\boldsymbol{X})$ should be chosen to
be proportional to $p(\boldsymbol{X})\cdot BB(\boldsymbol{X})$

The first point is jointly attainable for all target values, for example if
the support of $q(\boldsymbol{X})$ includes the supports of
$S_{1}(\boldsymbol{X}),...,S_{N}(\boldsymbol{X})$ as well as those of $S_{0}$
and $\mathcal{S}$.

The second point creates a trade-off between the different targets. Choosing
$q(\boldsymbol{X})=\mathcal{S}$ or $q(\boldsymbol{X})=S_{0}$, leads to
accurate estimates of the expected dose and nominal dose respectively. It can
be expected, that $q(\boldsymbol{X})=\mathcal{S}$ is a better choice with
respect to the dose variance, since it has more density in the outer regions,
which is necessary for an acceptable estimate of the shifted densities $S_{i}$
for larger errors $\boldsymbol{z_{i}}$.

When prioritizing the computation of the variance estimate, a more fat-tailed
choice of distribution for $q(\boldsymbol{X})$ is beneficial. This ensures a
larger number of realizations and therefore higher accuracy in the outer
regions of the densities, where the variance is typically high. In general it
is not advisable to choose a $q$ which is more light-tailed than the target
distribution, as the variance can explode if $q$ is too small in regions where
the target distribution is not zero [39]. One option is using a mixture
distribution of the individual shifted distributions we want to reconstruct
[19]. Another common rule of thumb in case of a normal target distribution is
to use a Student’s t distribution, as they are similar but the latter has
heavier tails [39]. The accuracy of estimates for some exemplary choices of
$q$ will be investigated numerically in section 8.3.2.

##  5 Mathematical properties

###  5.1 Accuracy & Bias

An advantage of the proposed method is that it inherits some of the
mathematical properties of importance sampling. Namely, for the nominal dose,
expected dose and dose for the single error realizations, the estimates are
unbiased [39].

######  Proposition 1.

For $q(\boldsymbol{X}):\mathcal{D}_{X}\rightarrow\mathbb{R}_{+}$ with
$Q=\\{\boldsymbol{X}|q(\boldsymbol{X})>0\\},\text{ if
}\;p(\boldsymbol{X})\cdot BB(\boldsymbol{X})\neq 0\text{ let
}\boldsymbol{X}\in Q$ . Then the importance reweighting estimator
$\hat{D}(\boldsymbol{r})$ is unbiased, i.e.
$E_{q}[\hat{D}(\boldsymbol{r})]=D(\boldsymbol{r})$.

Proof: According to section 4.3, the importance sampling estimate is defined
as

$\hat{D}(\boldsymbol{r})=\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;p(\boldsymbol{x_{p}})}{q(\boldsymbol{x_{p}})}\;,$
(28)

where $\boldsymbol{x_{p}},\;p=1,...,H$ are realizations of $\boldsymbol{X}\sim
q(\boldsymbol{X})$. Then

$\displaystyle E_{q}[\hat{D}(\boldsymbol{r})]$
$\displaystyle=E_{q}\left[\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;p(\boldsymbol{x_{p}})}{q(\boldsymbol{x_{p}})}\right]=\frac{1}{H}\sum_{p=1}^{H}E_{q}\left[\frac{BB(\boldsymbol{x_{p}})\;p(\boldsymbol{x_{p}})}{q(\boldsymbol{x_{p}})}\right]$
(29)
$\displaystyle=\int_{\mathcal{Q}}\frac{BB(\boldsymbol{X})\;p(\boldsymbol{X})}{q(\boldsymbol{X})}\cdot
q(\boldsymbol{X})d\boldsymbol{X}=\int_{\mathcal{Q}}BB(\boldsymbol{X})\cdot
p(\boldsymbol{X})d\boldsymbol{X}\;,$
$\displaystyle=\int_{\mathcal{D}_{p}}BB(\boldsymbol{X})\cdot
p(\boldsymbol{X})d\boldsymbol{X}+\int_{\mathcal{Q}\cap\mathcal{D}_{p}^{c}}\underbrace{BB(\boldsymbol{X})}_{=0}\cdot
p(\boldsymbol{X})d\boldsymbol{X}+\int_{\mathcal{Q}^{c}\cap\mathcal{D}_{p}}BB(\boldsymbol{X})\cdot\underbrace{p(\boldsymbol{X})}_{=0}d\boldsymbol{X}\;,$
$\displaystyle=E_{S_{0}}[BB(\boldsymbol{X})]=D(\boldsymbol{r})\;$

∎.

Also, the standard error $\epsilon_{q}$ of these estimators can be computed
similarly to that of Monte Carlo estimates [39]. It depends on the number of
realizations as well as the similarity of the importance distribution and the
distribution of interest:

$\epsilon_{q}=\frac{\sigma_{q}^{2}}{H}\,\text{, where
}\sigma_{q}^{2}\approx\hat{\sigma}_{q}^{2}=\frac{1}{H-1}\sum_{p=1}^{H}\left[\frac{BB(\boldsymbol{x_{p}})S_{0}(\boldsymbol{x_{p}})}{q(\boldsymbol{x_{p}})}-\hat{D}(\boldsymbol{r})\right]^{2}\;.$
(30)

To lower the standard error, the importance distribution should resemble the
target distribution. In the variance computation, the estimates for more
extreme error realizations will thus be less accurate. For the type and extent
of errors in radiation therapy, the distributions are however mostly similar
enough, as we will demonstrate further in section 8.

In contrast to the dose estimates, the variance computed as in section 4.5 is
not generally unbiased. This is due to its computation from estimates of the
dose for shifted beams, which rely on the same single sample of particles.
What we consider as the realizations in the sample variance computation, while
individually unbiased, are thus no longer independent.

######  Proposition 2.

For $\boldsymbol{x_{p}}$ from
$S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})$, let
$w_{i}(\boldsymbol{x_{p}})=\frac{S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{x_{p}})}$
be independent of $BB(\boldsymbol{x_{p}})$. Then the importance reweighting
estimator for the dose variance is not generally unbiased and the bias is
given by

$b_{S_{N}^{2}}=\left(\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\cdot
E[w_{i}(\boldsymbol{x_{p}})]\right)^{2}-\frac{1}{H^{2}}\sum_{p=1}^{H}\sum_{h=1}^{H}BB(\boldsymbol{x_{p}})BB(\boldsymbol{x_{h}})\cdot
E[w_{i}(\boldsymbol{x_{p}})w_{j}(\boldsymbol{x_{h}})]$

Proof: In [4], general formulas for the moments of the sample variance
$S_{n}^{2}$ are derived. The first moment $E[S_{n}^{2}]$ for general non
independent or identically distributed (iid) samples $X_{1},...,X_{N}$ is
given by:

$E[S_{n}^{2}]=\frac{\sum_{i=1}^{n}\mathbb{E}[X_{i}^{2}]}{n}-\frac{\sum_{i\neq
j}\mathbb{E}[X_{i}X_{j}]}{n(n-1)}\;X_{1},...,X_{n}\text{not iid}\;.$

Applied to our context, the realizations are
$D_{i}=\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})w_{i}(\boldsymbol{x_{p}})$
for $i=1,...,N$, with
$w_{i}(\boldsymbol{x_{p}})=\frac{S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{x_{p}})}$.
Since all realizations use the same sample of computed particle trajectories
$BB(\boldsymbol{x_{p}})$ and initial particle properties
$\boldsymbol{x_{p}},\;p=1,...,H$, they cannot be assumed to be independent, in
particular
$D_{i}=D_{j}+\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\left(w_{i}(\boldsymbol{x_{p}}-w_{j}(\boldsymbol{x_{p}})\right)$.
Thus,

$\displaystyle E[S_{n}^{2}]$
$\displaystyle=\frac{\sum_{i=1}^{N}E[D_{i}^{2}]}{N}-\frac{\sum_{i\neq
j}E[D_{i}D_{j}]}{N(N-1)}\;.$ (31)

Since we assume that the $D_{i}$ are identically distributed, this reduces to

$\displaystyle E[S_{n}^{2}]$ $\displaystyle=E[D_{i}^{2}]-E[D_{i}D_{j}]$
$\displaystyle=\underbrace{E[D_{i}^{2}]-E[D_{i}]E[D_{j}]}_{=Var(D)}-\underbrace{\left(E[D_{i}]E[D_{j}]-E[D_{i}D_{j}]\right)}_{=b_{S_{N}^{2}}\neq
0\;,\text{ for}\;i\neq j}\;.$

Thus the estimator is not necessarily unbiased and the bias is

$\displaystyle b_{S_{N}^{2}}$
$\displaystyle=\underbrace{E[D_{i}]E[D_{j}]}_{=E[D_{i}]E[D_{i}]}-E[D_{i}D_{j}]$
$\displaystyle=E\left[\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})w_{i}(\boldsymbol{x_{p}})\right]^{2}-E\left[\left(\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})w_{i}(\boldsymbol{x_{p}})\right)\left(\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})w_{j}(\boldsymbol{x_{p}})\right)\right]$
$\displaystyle=\left(\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})E\left[w_{i}(\boldsymbol{x_{p}})\right]\right)^{2}-E\left[\frac{1}{H^{2}}\sum_{p=1}^{H}\sum_{h=1}^{H}BB(\boldsymbol{x_{p}})BB(\boldsymbol{x_{h}})w_{i}(\boldsymbol{x_{p}})w_{j}(\boldsymbol{x_{h}})\right]$
$\displaystyle=\left(\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})E\left[w_{i}(\boldsymbol{x_{p}})\right]\right)^{2}-\frac{1}{H^{2}}\sum_{p=1}^{H}\sum_{h=1}^{H}BB(\boldsymbol{x_{p}})BB(\boldsymbol{x_{h}})E\left[w_{i}(\boldsymbol{x_{p}})w_{j}(\boldsymbol{x_{h}})\right]$
(33)

∎.

Since this bias term is costly to compute, it is generally not feasible to
score this in addition to the quantities of interest. Computing the bias once
for a given problem class could however be of interest to judge its order of
magnitude. Its potential impact on our results will be assessed in greater
detail in section 8.

###  5.2 Bounds on the variance

While the discussed estimators can be expected to significantly lower the
expense of scenario computations, the variance still requires $N$ scoring
operations, either accumulated on-the-fly after each simulated particle
history, or in post processing as (re-)weighting of stored histories. For a
large number of simulated particles, run-times can therefore exceed the
required time-frame, especially in the context of robust optimization, where
the variance has to be computed in each of numerous iterations [46, 30, 7].

######  Proposition 3.

An upper bound to the variance estimator $\hat{\sigma}^{2}$ is given by

$\hat{\sigma}^{2}\leq\left(\;\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})^{2}\;\right)\left(\;\frac{1}{H}\sum_{p=1}^{H}S_{N}^{2}(w(\boldsymbol{x_{p}}))\;\right)$

Proof: Using the definition of the importance reweighing estimator for the
dose variance, the dose for a single error realization and the expected dose
(see equations 27, 4.2 and 26), we get

$\displaystyle\hat{\sigma}^{2}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}\;\left[\;\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;S_{i}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{i}})}{q(\boldsymbol{x_{p}})}-\frac{1}{N}\sum_{j=1}^{N}\frac{1}{H}\sum_{p=1}^{H}\frac{BB(\boldsymbol{x_{p}})\;S_{j}(\boldsymbol{x_{p}}|\boldsymbol{Z}=\boldsymbol{z_{j}})}{q(\boldsymbol{x_{p}})}\;\right]^{2}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}\;\left[\;\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\;w_{i}(\boldsymbol{x_{p}})-\frac{1}{N}\sum_{j=1}^{N}\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\;w_{j}(\boldsymbol{x_{p}})\;\right]^{2}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}\;\left[\;\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\;w_{i}(\boldsymbol{x_{p}})-\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\frac{1}{N}\sum_{j=1}^{N}\;w_{j}(\boldsymbol{x_{p}})\;\right]^{2}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}\;\left[\;\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})\left(w_{i}(\boldsymbol{x_{p}})-\frac{1}{N}\sum_{j=1}^{N}\;w_{j}(\boldsymbol{x_{p}})\right)\;\right]^{2}$
$\displaystyle\leq\left(\;\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})^{2}\right)\;\Biggl{(}\frac{1}{H}\sum_{p=1}^{H}\;\underbrace{\frac{1}{N-1}\sum_{i=1}^{N}\left[w_{i}(\boldsymbol{x_{p}})-\sum_{j=1}^{N}w_{j}(\boldsymbol{x_{p}})\right]^{2}}_{=\,S_{N}^{2}(w(\boldsymbol{x_{p}}))\,\approx\,Var(w(\boldsymbol{x_{p}}))}\Biggr{)}$
$\displaystyle=\left(\;\frac{1}{H}\sum_{p=1}^{H}BB(\boldsymbol{x_{p}})^{2}\;\right)\left(\;\frac{1}{H}\sum_{p=1}^{H}S_{N}^{2}(w(\boldsymbol{x_{p}}))\;\right)\;$

∎.

The advantage of this bound is, that the squared simulation response
$BB(\boldsymbol{x_{p}})^{2}$ is aggregated independently from the error
realizations. This includes the implicit assumption, that the variance in the
likelihood ratio between the simulated and target distribution, is enough to
describe the variance of the output. While this might hold for relatively
homogeneous background media, it is unclear whether it is sufficient in
heterogeneous media as present in many patient anatomies.

##  6 Modeling set-up uncertainties

We now consider an application of the proposed approach to the case of
uncertainties in the patient set-up or position. These can be modeled by
random shifts in the particle positions $\boldsymbol{r}$. The input
$\boldsymbol{X}$ is thus affected by the error additively, i.e., the input
$\boldsymbol{X(Z)}$ under uncertainty is

$\boldsymbol{X(Z)}=\boldsymbol{X}+\boldsymbol{Z}\;,\text{where
}\boldsymbol{X}=(\boldsymbol{r},\boldsymbol{\Omega},E)\text{ and
}\boldsymbol{Z}=(\boldsymbol{Z_{r}},\boldsymbol{0},0).$ (34)

We follow the common assumption, that errors as well as input parameters are
normally distributed [60, 54, 42, 1, 15, 57]

$p_{Z}=\mathcal{N}(\boldsymbol{\mu}_{Z},\boldsymbol{\Sigma})\;,\;S_{0}=\mathcal{N}(\boldsymbol{\mu}_{X},\boldsymbol{\Lambda})\;.$
(35)

In this case we can define the density $\mathcal{S}$ of the input including
the uncertain factors through convolution of $p_{Z}$ and $S_{0}$:

$\mathcal{S}:=\mathcal{N}(\boldsymbol{\mu}_{Z}+\boldsymbol{\mu}_{X},\boldsymbol{\Sigma}+\boldsymbol{\Lambda})\;.$
(36)

Thus, the expected value can be directly computed using section 4.4 and the
nominal dose as well as the dose variance are computed as described in
sections 4.3 and 4.5.

###  6.1 Beam correlation models

The most frequently used uncertainty model for set-up errors assumes one
global error
$\boldsymbol{Z_{r}}=(Z_{r_{x}},Z_{r_{y}},Z_{r_{z}})\sim\mathcal{N}(\boldsymbol{\mu}_{Z_{r}},\boldsymbol{\Sigma}_{r}).$
Therefore, it corresponds to movements or mispositioning of the patient which
affect the complete body or at least area of interest around the tumor.

However, in radiation therapy, often numerous beamlets with different energy
levels, positions and incoming angles are used in order to achieve better
coverage of the tumor. Here, the term beamlet or pencil beam refers to a
single instance of energy, angular and lateral position. An ensemble of
beamlets coming from the same geometrical set-up of nozzle and patient is then
referred to as beam, a group of beamlets with the same lateral positioning
within a beam as ray and a group of beamlets with the same energy within a
beam as energy level. In this case, it may be more reasonable to assume
beamlet-specific correlation patterns, e.g., depending on the differences in
time and position of application of the individual beamlets. Additionally,
treatments are usually administered in several sittings or fractions, which
can also affect how errors are correlated. A number of simple correlation
models, assuming full correlation between the beamlets belonging to the same
beam or ray, have been investigated by different authors [1, 57, 46] and also
discussed in our previous work [51]. However, depending on the dose
calculation and uncertainty quantification method, the incorporation of more
complex correlation models can be tedious to both set up and compute.

Fortunately, the introduced framework principally allows arbitrary correlation
models for uncertainties in the set-up of different beamlets. These can be
implemented using the covariance matrix of a multivariate error distribution.
We define the errors in each beamlet by
$\boldsymbol{Z^{b}_{r}}=(Z^{b}_{r_{x}},Z^{b}_{r_{y}},Z^{b}_{r_{z}})\sim\mathcal{N}(\boldsymbol{\mu}^{b}_{Z_{r}},\boldsymbol{\Sigma}_{r}^{b})$
and their correlations by

$C=\begin{pmatrix}\Sigma_{r}^{1}&\begin{matrix}\rho_{xx}^{12}&\rho_{xy}^{12}&\rho_{xz}^{12}\\\
\rho_{yx}^{12}&\rho_{yy}^{12}&\rho_{yz}^{12}\\\
\rho_{zx}^{12}&\rho_{zy}^{12}&\rho_{zz}^{12}\\\
\end{matrix}&\cdots&\begin{matrix}\rho_{xx}^{1B}&\rho_{xy}^{1B}&\rho_{xz}^{1B}\\\
\rho_{yx}^{1B}&\rho_{yy}^{1B}&\rho_{yz}^{1B}\\\
\rho_{zx}^{1B}&\rho_{zy}^{1B}&\rho_{zz}^{1B}\\\ \end{matrix}\\\
\begin{matrix}\rho_{xx}^{21}&\rho_{xy}^{21}&\rho_{xz}^{21}\\\
\rho_{yx}^{21}&\rho_{yy}^{21}&\rho_{yz}^{21}\\\
\rho_{zx}^{21}&\rho_{zy}^{21}&\rho_{zz}^{21}\\\
\end{matrix}&\Sigma_{r}^{2}&&\\\ \vdots&&\ddots&\vdots\\\
\begin{matrix}\rho_{xx}^{B1}&\rho_{xy}^{B1}&\rho_{xz}^{B1}\\\
\rho_{yx}^{B1}&\rho_{yy}^{B1}&\rho_{yz}^{B1}\\\
\rho_{zx}^{B1}&\rho_{zy}^{B1}&\rho_{zz}^{B1}\\\
\end{matrix}&\cdots&&\Sigma_{r}^{B}\end{pmatrix}\;.$ (37)

Here $\rho_{xy}^{ab}$ is the correlation between the set-up errors of beamlet
$a$ in dimension $x$ and beamlet $b$ in dimension $y$. It is then possible to
obtain error realizations for different choices of the correlations, using a
known probability distribution or, alternatively, a Copula with mean
$\boldsymbol{\mu}_{r}=(\boldsymbol{\mu}^{1}_{Z_{r}},...,\boldsymbol{\mu}^{B}_{Z_{r}})^{T}$
and covariance matrix $\boldsymbol{C}$.

Besides the discussed simple correlation models, this framework allows us to
model more complex correlation patterns. This can, for example, be used to
incorporate knowledge about the radiation process and even model certain time
dependencies. In the following, we demonstrate the derivation of the
covariance matrix $\boldsymbol{C}$ for two examples of autocorrelated set-up
errors.

####  6.1.1 Autocorrelated set-up errors

To show how the covariance matrix $\boldsymbol{C}$ can be obtained for a given
correlation model, we first consider the case of autocorrelated set-up errors,
which follow an AR(1) process. While irradiation itself usually only takes a
few milliseconds, beamlets are applied sequentially according to their energy
level and up to 1-2 seconds can pass between them. Thus, assuming the patient
is exhibiting random movement with some probability at all points in time, the
patient position corresponding to beamlets in adjacent energy groups should be
correlated to a higher degree than those of beamlets with a larger difference
in energies. Since first order autoregressive models are sufficient to capture
interfractional set-up errors [21], we choose such a model for the random
patient movements during treatment.

We thus describe movements by an autoregressive process of order one, meaning
that the overall deviation at time $t$ depends on the deviation at time $t-1$
as well as a random variable
$\boldsymbol{U}_{t}\mathrel{\stackon[4pt]{\sim}{\scalebox{0.7}{\mbox{$\displaystyle
i.i.d.$}}}}\mathcal{N}(\boldsymbol{\mu}_{Z_{r}},\boldsymbol{\Sigma}_{r})$,
which represents the movement error introduced between $t-1$ and $t$:

$\boldsymbol{Z}_{r;t}=\alpha\boldsymbol{Z}_{r;t-1}+\boldsymbol{U}_{t}\;.$ (38)

Here the factor $|\alpha|<1$ controls how strongly the position at the
previous time step affects the current position. It can for example be chosen
depending on how much time passes between the time steps. In our case each
beam $b$ corresponds to a point in time $t(b)$ according to its energy and we
assume that there is no additional correlation between the
(x,y,z)-coordinates, i.e.,

$\boldsymbol{\Sigma}=\begin{pmatrix}\boldsymbol{\sigma}_{x}^{2}&0&0\\\
0&\boldsymbol{\sigma}_{y}^{2}&0\\\ 0&0&\boldsymbol{\sigma}_{z}^{2}\\\
\end{pmatrix}\;.$ (39)

From this description of the process it is now possible to obtain the
covariance matrix $\boldsymbol{C}$ using the Autocovariance Function (ACVF):

$\rho_{ab}^{ij}=Cov(Z_{r_{a};t(i)},Z_{r_{b};t(j)})=\begin{cases}\frac{\alpha^{|t(i)-t(j)|}\sigma_{x}^{2}}{1-\alpha^{2}}\;\;,a=b=x\\\
\\\ \frac{\alpha^{|t(i)-t(j)|}\sigma_{y}^{2}}{1-\alpha^{2}}\;\;,a=b=y\\\ \\\
\frac{\alpha^{|t(i)-t(j)|}\sigma_{z}^{2}}{1-\alpha^{2}}\;\;,a=b=z\\\ \\\
0\;\;\;\;\;\;\;\;\;\;\;\;,a\neq b\end{cases}\;.$ (40)

Note that the variance is constant in time and the covariances only depend on
the time difference. To obtain a joint probability function we use a Gaussian
copula, which for Gaussian marginals reduces to the multivariate normal
distribution $\mathcal{N}(\boldsymbol{\mu}_{Z_{r}},\boldsymbol{C})$ with
$\boldsymbol{\mu}_{Z_{r}}=(\boldsymbol{\mu}_{Z_{r}}^{1},\boldsymbol{\mu}_{Z_{r}}^{2},...,\boldsymbol{\mu}_{Z_{r}}^{B})$
and

$\boldsymbol{C}=\begin{pmatrix}\frac{\sigma_{x}^{2}}{1-\alpha^{2}}&0&0&&&&\frac{\alpha^{|t(B)-t(1)|}\sigma_{x}^{2}}{1-\alpha^{2}}&0&0\\\
0&\frac{\sigma_{y}^{2}}{1-\alpha^{2}}&0&&\cdots&&0&\frac{\alpha^{|t(B)-t(1)|}\sigma_{y}^{2}}{1-\alpha^{2}}&0\\\
0&0&\frac{\sigma_{z}^{2}}{1-\alpha^{2}}&&&&0&0&\frac{\alpha^{|t(B)-t(1)|}\sigma_{z}^{2}}{1-\alpha^{2}}\\\
&&&&&&&&\\\ &\vdots&&&\ddots&&&\vdots&\\\ &&&&&&&&\\\
\frac{\alpha^{|t(1)-t(B)|}\sigma_{x}^{2}}{1-\alpha^{2}}&0&0&&&&\frac{\sigma_{x}^{2}}{1-\alpha^{2}}&0&0\\\
0&\frac{\alpha^{|t(1)-t(B)|}\sigma_{y}^{2}}{1-\alpha^{2}}&0&&\cdots&&0&\frac{\sigma_{y}^{2}}{1-\alpha^{2}}&0\\\
0&0&\frac{\alpha^{|t(1)-t(B)|}\sigma_{z}^{2}}{1-\alpha^{2}}&&&&0&0&\frac{\sigma_{z}^{2}}{1-\alpha^{2}}\\\
&\end{pmatrix}\;.$ (41)

This AR(1)-correlation model gives an intuitive description for random
movements of a patient over the time of irradiation. However, there are also
errors caused by periodic movements such as respiratory motion or heart beats.
Respiratory motion is often modeled using Gaussian processes with an
appropriate covariance kernel [11], however, there is also literature
exploring the use of a periodic autoregressive moving average (ARMA) model
[20]. Both could be incorporated into our uncertainty framework, we however
confine ourselves to showing an example for a Gaussian process with a (local)
periodic kernel here.

The Gaussian process (GP) implements the concept of data in a time series
following a joint multivariate Gaussian distribution, where the covariance
kernel defines the dependencies between different data points. As the type of
movement we are interested in has a dominant periodic component, which may
however vary over time, we choose the following kernel
$k(Z_{r;t(i)},Z_{r;t(j)})$, proposed similarly by [12],

$k(Z_{r;t(i)},Z_{r;t(j)})=\sigma^{2}\cdot\exp\left(-\frac{2\sin^{2}(\pi(t(i)-t(j))/p)}{l_{1}^{2}}\right)\cdot\exp\left(-\frac{(t(i)-t(j))^{2}}{2l_{2}^{2}}\right)\;.$
(42)

Here $\sigma^{2}$ is the variance of the stationary process, $p$ is the period
and $l_{1}$ and $l_{2}$ are scale parameters. The first exponential term
describes an exact periodic repetition, whereas the second term is a local
kernel which introduces the variations over time. Figure 1 illustrates the
covariance kernel as well as five exemplary realizations. The process starts
at a random point in the cycle and periodically repeats itself with some
random variations. From equation 42 the derivation of the covariance matrix
$\boldsymbol{C}$ is straightforward with the entries

(a)

(b)

Figure 1: Local periodic covariance kernel with $\sigma=5$, $p=3$,$l_{1}=1$
and $l_{2}=5$ (a) and five realizations of the associated GP (b)

$\rho_{ab}^{ij}=Cov(Z_{r_{a};t(i)},Z_{r_{b};t(j)})=\begin{cases}\sigma_{x}^{2}\cdot\exp(-\frac{2\sin^{2}(\pi(t(i)-t(j))/p)}{l_{1}^{2}})\cdot\exp(-\frac{(t(i)-t(j))^{2})}{2l_{2}^{2}})\;\;,a=b=x\\\
\\\
\sigma_{y}^{2}\cdot\exp(-\frac{2\sin^{2}(\pi(t(i)-t(j))/p)}{l_{1}^{2}})\cdot\exp(-\frac{(t(i)-t(j))^{2})}{2l_{2}^{2}})\;\;,a=b=y\\\
\\\
\sigma_{z}^{2}\cdot\exp(-\frac{2\sin^{2}(\pi(t(i)-t(j))/p)}{l_{1}^{2}})\cdot\exp(-\frac{(t(i)-t(j))^{2})}{2l_{2}^{2}})\;\;,a=b=z\\\
\\\ 0\;\;\;\;\;\;\;\;\;\;\;\;,a\neq b\end{cases}\;.$ (43)

The errors for all pencil beams can now be sampled from the multivariate
Gaussian $\mathcal{N}(\boldsymbol{\mu}_{\delta},\boldsymbol{C})$ and will
depend on the time at which their corresponding energy level is applied.

The two examples above demonstrate, that it is possible to implement complex
movement patterns and even to some degree model time dependencies, by changing
the way error realizations are sampled.

##  7 Implementation and Test cases

The described method was implemented as post-processing in Matlab. Radiation
plans were generated for proton dose calculations using the treatment planning
software matRad [59] and dose calculations were performed with the Monte Carlo
engine TOPAS [43]. For quicker convergence we use randomized quasi-Monte Carlo
with scrambled Sobol numbers when sampling the error scenarios for the dose
variance estimate (see section 2.2).

In the following, we present the results for a 3D water phantom as well as a
prostate and liver patient obtained from the open CORT (common optimization
for radiaion therapy) data set [8]. Overall, we consider three different
treatment plans, combined with the three error models discussed in section 6.
For the Gaussian set-up error we assume a standard deviation of
$3\text{\,}\mathrm{mm}$ and mean $\mu_{\delta}=0$, i.e., the errors do not
have a systematic component [comp. 42]. For two of the treatment plans,
intensity modulated proton therapy (IMPT) [25, 31] was employed. In these
cases, each source consists of a large number of narrow pencil beams to enable
a better and more flexible coverage of the target volume. More details
concerning the treatment plans and the respective error models applied to each
of them can be found in table 1.

Table 1: Overview of simulated plans and error models per test case. Test case | Water phantom | Prostate | Liver
---|---|---|---
Size | 60x25x25 voxels | 183x183x90 voxels | 217x217x168 voxels
Irradiation angles | ($None$,$None$) | ($None$,$None$)/($None$,$None$) | ($None$,$None$)
Pencil beams | 1 | 175 | 1 375/1 383 | 1408
Simulated particles | 100 000 | 2 566 453 | 16 992 193/16 748 034 | 13 528 430
Error model | global | global/AR(1)/GP | global | global/AR(1)/GP
Sample size | 500 | 500 | 100 | 100

###  7.1 Evaluation criteria

In order to evaluate the quality of our estimates, we compute a reference for
each test case using the regular randomized quasi-Monte Carlo method for the
uncertainties. The number of realizations was chosen to be equal to that used
for the importance (re-)weighting estimates, which is 500 for the water
phantom and 100 for the prostate patient.

We compare the results according to the $\gamma$ criterion [34], as well as a
difference map. For the difference maps, we compute

$\text{diff}_{i}(d^{\text{ref}}_{i},d_{i}^{\text{est}})=d_{i}^{\text{ref}}-d_{i}^{\text{est}}\;,$
(44)

for each voxel $i$ in the reference result $d^{\text{ref}}$ and (re-)weighting
estimate $d^{\text{est}}$.

The $\gamma$ criterion is frequently used for comparisons of three-dimensional
dose matrices in clinical applications. To compare doses $d$ and $d_{c}$, one
first computes

$\displaystyle\gamma(x)=\min_{x_{c}}\sqrt{\frac{|x_{c}-x|^{2}}{\Delta\tilde{x}^{2}}+\frac{|d_{c}(x_{c})-d(x)|^{2}}{\Delta\tilde{d}^{2}}}\;,$
(45)

where $x$ and $d(x)$ are a location and the dose at that location
respectively, $x_{c}$ and $d_{c}(x_{c})$ are locations and dose in the cube to
be compared and $\Delta\tilde{x}$ and $\Delta\tilde{d}$ are the user defined
distance and dose tolerances. The overall agreement is then measured according
to the $\gamma$-passrate over all points $x$. For each point

$\begin{cases}\gamma(x)\leq 1,\;x\text{ passes}\\\ \gamma(x)>1,\;x\text{
fails}\;.\end{cases}$ (46)

It implements the idea of agreement within a predefined neighborhood and
tolerance around a single coordinate or data point. In contrast to criteria
such as the difference or (root) mean square error, which compare the values
at an exact location, this does not take minor fluctuations in the spacial
dose distribution into account. Since the reference itself is also subject to
statistical errors, the dose distribution is not necessarily incorrect if it
deviates slightly from the reference. However, the $\gamma$ criterion is not a
metric in the mathematical sense if the tolerances are larger than zero. Thus,
two dose cubes could have $100\%$ agreement according to the $\gamma$-passrate
without actually being the same. Therefore the $\gamma$-passrates always have
to be interpreted together with the applied tolerances. In the following, we
choose values frequently used in literature: a distance tolerance of
$3\text{\,}\mathrm{mm}$ and a dose tolerance of $3\text{\,}\mathrm{\char
37\relax}$. Furthermore, the $\gamma$-passrate can depend on the grid
resolution as well as the implementation [22] . We use the matRad
implementation [59], which reduces the dependency on the resolution by
interpolating the comparison dose cube and comparing it at a higher
resolution.

##  8 Results

###  8.1 Fully correlated pencil beams

Figure 2 presents the expected dose and standard deviation for a single pencil
beam in a water phantom. In order to enable a direct comparison to a more
complex application case, we also include a recomputation of the respective
quantities for the prostate patient investigated in [51] with fully correlated
pencil beams, i.e., one global error, in figure 3. The importance distribution
was chosen to be the nominal distribution $S_{0}$. For this reason, the
nominal dose estimate is just a regular Monte Carlo estimate and is therefore
omitted in the following figures.

$E[\boldsymbol{d}]$

Estimate

Reference

Difference

$\gamma$-map

$\boldsymbol{\sigma}$

Figure 2: Expected dose and standard deviation w.r.t. set-up uncertainties
with $3\text{\,}\mathrm{mm}$ standard deviation for one pencil beam in a water
box. The left column shows the estimate computed with the proposed
(re-)weighting approach, the middle column a reference computed with Monte
Carlo and the right column a plot of the
$\gamma_{$3\text{\,}\mathrm{mm}$/$3\text{\,}\mathrm{\char
37\relax}$}$-indices.

The results using the importance (re-)weighting estimate are in good agreement
with the reference, with a
$\gamma_{$3\text{\,}\mathrm{mm}$/$3\text{\,}\mathrm{\char
37\relax}$}$-passrate of $100\text{\,}\mathrm{\char 37\relax}$ for the
expected value and over $99.9\text{\,}\mathrm{\char 37\relax}$ with respect to
the standard deviation in both waterbox and prostate patient.

$E[\boldsymbol{d}]$

Estimate

Reference

$\gamma$-map

Difference

$\boldsymbol{\sigma}$

Figure 3: Expected dose and standard deviation w.r.t. set-up uncertainties
with $3\text{\,}\mathrm{mm}$ standard deviation for two beams (gantry angles
90° and 270°, couch angle 0°) in a prostate patient. The left column shows the
estimate computed with the proposed (re-)weighting approach, the middle column
a reference computed with Monte Carlo and the right column the results of a
$\gamma_{$3\text{\,}\mathrm{mm}$/$3\text{\,}\mathrm{\char
37\relax}$}$-analysis.

###  8.2 Other correlation models

In section 6, we derived correlation matrices and error distributions for
different uncertainty models. Our implementation was validated against Monte
Carlo results for individual realizations of shifts over time. For a single
ray with $6$ energy levels in a water box, shifted by
$\pm$2\text{\,}\mathrm{m}\mathrm{m}$$ every $10\text{\,}\mathrm{m}\mathrm{s}$,
the dose results in each time step reached $100\text{\,}\mathrm{\char
37\relax}$ agreement in a $\gamma$-analysis with tolerances of
$$2\text{\,}\mathrm{mm}$/$2\text{\,}\mathrm{\char 37\relax}$$ in dose and
location, respectively. Figure 4 compares the dose standard deviation for
three such models in a water phantom with a beam consisting of several
superimposed pencil beams of seven different energy groups. The results
illustrate that already in this simple case the uncertainty assumption has a
significant impact on both the amount and shape of the dose standard
deviation. However, due to the low amount of energy levels and homogeneous
tissue, the characteristics of the different input models are not fully
visible.

This becomes apparent in figure 5, which compares the standard deviation for
the three error models in a realistic liver patient, with 21 energy levels and
166 rays. While the region of high standard deviation is still at the
gradients at the edge of the beam in the fully correlated and periodic model,
several local peaks can additionally be observed in the latter. In this case,
considering the accumulated dose over all time steps/energy levels results in
a lower magnitude of standard deviation, due to averaging of the errors. In
the AR(1) model, which mimics a random movement process over time, the
standard deviation is less concentrated to the beam edges. Here, pencil beams
are frequently shifted against instead with each other, due to the larger
random component of the error in each time step/energy level.

(a) Full correlation

(b) AR(1)

(c) Local periodic

Figure 4: Standard deviation of dose from a beam made up of 175 pencil beams
in a 3D waterbox for different error models: (a) one global error (b) AR(1)
model for movement between energy levels (c) seasonal AR(1) model for movement
between energy levels

(a) Full correlation

(b) AR(1)

(c) Local periodic

Figure 5: Standard deviation of dose from a beam with irradiation angle
(gantry angle $None$, couch angle $None$) made up of 1378 pencil beams in a
liver patient for different error models: (a) one global error (b) seasonal
AR(1) model for movement between energy levels

###  8.3 Practical investigation of mathematical properties

####  8.3.1 Bias

The above findings indicate that the bias in the dose variance estimate,
introduced in proposition 2, does not significantly affect the accuracy of the
results, as the agreement with the reference solutions is comparable for the
unbiased expected value estimator and the biased variance estimator. Further,
figure 6 shows that the bias reduces for an increasing number of particle
histories.

(a)

(b)

Figure 6: Convergence of the variance estimator’s bias in (a) number of error
realizations used to compute the bias and (b) the number of particle histories
used for the dose reconstructions for a single pencil beam in a water box.

####  8.3.2 Standard error

So far, the dose for all estimates was reconstructed from samples of the
distribution $S_{0}$. However, we have seen in equation 30, that the standard
error of the estimates depends on the relation of distribution from which the
original sample was drawn to the target distribution. When computing the dose
standard deviation, each error scenario corresponds to a Gaussian input
distribution with a shifted mean value (see section 6). Figure 7 shows that
the accuracy of the estimates decreases with growing distance of the mean of
the target distributions from that of the sampling distribution. This
development can be observed not only for the theoretical standard error, but
also for the mean square error and $\gamma$-passrates. This indicates that the
more extreme realizations of the uncertain parameters introduce relatively
high errors to the overall standard deviation estimate. We want to see,
whether the use of a better suited sampling distribution, such as a wider
Gaussian or a Gaussian mixture of some of the target distributions, can
improve the accuracy of the estimate.

For this, we sample input parameters once from the joint distribution
$\mathcal{S}$ and once from a mixture of the narrower Gaussian $S_{0}$,
shifted by $\pm\sigma$ and $\pm 2\sigma$ in all coordinate directions. In
order to better detect differences in the results, we compare the estimates
constructed from these samples with the results from $S_{0}$ using the
stricter criteria $2\text{\,}\mathrm{mm}$/$2\text{\,}\mathrm{\char 37\relax}$
for the $\gamma$-analysis.

(a)

(b)

(c)

Figure 7: Accuracy of dose estimates for increasing distance between
importance distribution $q$ and target distribution $\Phi_{i}$. (a) Measured
by mean square error (MSE), (b) by $\gamma$-passrate and (c) the mean of the
theoretical error (see equation 30) in voxels with non-zero dose.

Figure 8 illustrates that the $\gamma$-passrate increases when using both the
wider Gaussian and the Gaussian mixture. Thus, especially when one has an
uncertainty model involving errors with a high variance or if the expected
value and standard deviation of the dose are more important than the nominal
dose, it can be worthwhile to not use the regular parameter distribution to
generate the initial sample. For medical purposes, the high accuracy of the
nominal dose is however an integral part of the quality assurance of a method
and could often be prioritized over the quantification of uncertainties.

(a) $S_{0}$, $\gamma_{$2\text{\,}\mathrm{mm}$/$2\text{\,}\mathrm{\char
37\relax}$}$-passrate = $94.64\text{\,}\mathrm{\char 37\relax}$

(b) $\mathcal{S}$, $\gamma_{$2\text{\,}\mathrm{mm}$/$2\text{\,}\mathrm{\char
37\relax}$}$-passrate = $99.86\text{\,}\mathrm{\char 37\relax}$

(c) Gaussian mixture,
$\gamma_{$2\text{\,}\mathrm{mm}$/$2\text{\,}\mathrm{\char
37\relax}$}$-passrate = $99.98\text{\,}\mathrm{\char 37\relax}$

Figure 8: Estimates of the dose standard deviation in a water box using
different sampling distributions: (a) the nominal distribution $S_{0}$, (b)
the distribution $\mathcal{S}$ corresponding to the expected value and (c) a
Gaussian mixture distribution

####  8.3.3 Variance upper bound

Lastly, in section 5.2 an upper bound for the dose variance estimate was
derived, which is significantly less computationally expensive than the
variance estimate itself. Figure 9 shows a comparison of this upper bound
estimate to the dose variance. While the variance values are expectedly higher
in the upper bound estimate, both exhibit a similar structure with two
distinct variance peaks at positions with high dose gradients. This indicates
a potential use for optimization purposes. However, further analyses in more
complex test cases, in particular such with several irradiataion angles and in
heterogeneous materials, are necessary.

Upper bound of dose variance

Dose variance

Difference

Figure 9: Comparison of the upper bound to the dose variance estimate for a
beam made up of 175 pencil beams in a 3D water box.

###  8.4 Convergence

The proposed approach mimics a sampling-based uncertainty quantification
method. The solution of the modified problem for each error realization using
the black box solver is replaced by (re-)weighting the dose corresponding to
particle trajectories from a previous simulation. Thus, the convergence of the
method is equivalent to that of the underlying UQ method, in our case
randomized quasi-Monte Carlo. Figure 10 illustrates this for the dose standard
deviation of the single beam in a water box. While the convergence of
reference and estimate is very similar per computed realization, the speed-up
is achieved by the quicker computation of each iteration, that is, only
performing MC scoring operations compared to physical simulations. Therefore,
the amount of time which can be saved also depends on the specific patient
case, treatment plan and software, which determine how much overhead is
produced for example by the initialization and simulation of physical
processes.

(a)

(b)

Figure 10: Mean squared error of different methods compared (a) per iteration
and (b) per time.

For run-time comparisons, reference computations and the (re-)weighting
approach were run on the same virtual machine222Virtual machine including 64
CPUs with 1.995 GHz and 200GB RAM and timed during the computation of the dose
standard deviation using the global error model. As table 2 shows, we observe
reduced CPU times by a factor of approximately 183 and 80 for the two
treatment plans using the water phantom and factors of 32 and 23 for the liver
and prostate patient. Note, that the large difference in initialization times
between the IMPT plans and the plan involving just one beam is due to the size
of the sparse matrix containing the doses for all particle histories. For the
larger plans, this exceeded the available RAM such that the use of less
efficient tall arrays was necessary. Note, that these conclusions refer solely
to a post-processing implementation of the method. When implemented as on-the-
fly scoring no storage of histories is required, instead, one sparse dose cube
per considered scenario has to be stored. Also, the possibility to use
existing efficient structures of the Monte Carlo code might enable performance
improvements compared to the post-processing routine.

Table 2: CPU time comparison for the reference vs. (re-)weighting approach applied to different test cases and computed on the same machine. All values are given in seconds. Note that the times for 100 realizations include the initialization times, while the time for a single realization only refers to the dose computation time. |  | reference | (re-)weighting
---|---|---|---
Water phantom (1 beam) | Initialization | 2.27 | 1.34
One realization | 72.21 | 0.38
100 realizations | 7221.64 | 39.34
Water phantom (IMPT) | Initialization | 2.35 | 61.53
One realization | 2331.30 | 28.51
100 realizations | 233126.93 | 2912.53
Liver | Initialization | 2.44 | 2038.75
One realization | 39066.44 | 1198.74
100 realizations | 3906650.90 | 121912.75
Prostate | Initialization | 4.26 | 4867.75
One realization | 58762.40 | 2479.07
100 realizations | 5876253.86 | 252774.75

##  9 Discussion

In this paper we demonstrate how importance sampling can be adapted to lower
the computational costs of sampling-based uncertainty quantification in Monte
Carlo dose calculations using the example of proton irradiation. We derive an
unbiased estimator and the standard error for the expected value, as well as
an estimator, bias and upper bound for the standard deviation. We apply the
proposed method to Gaussian set-up uncertainties with $3\text{\,}\mathrm{mm}$
standard deviation in a homogeneous water phantom, a prostate and a liver
patient. For a simple, global uncertainty model we observe high agreements of
at least $99.9\text{\,}\mathrm{\char 37\relax}$ for both test cases (figures 2
and 3). We show that these results can be further improved by adapting the
distribution function used for the initial Monte Carlo simulation to the
quantity which is of highest interest. Here the user can choose between the
exactness of the nominal dose estimate and that of the expected dose or
standard deviation.

The sampling of error realizations is completely independent from the solution
of the transport problem, which allows us to incorporate complex correlation
models including movement patterns and time dependence. This is a clear
advantage compared to classic sampling based methods, which would require
elaborate set up and a large number of realizations to account for the degrees
of freedom in a multivariate model. Therefore, more complex correlation models
are currently often neglected in uncertainty studies [1, 46, 30], although it
is apparent in figures 4 and 5, that different uncertainty models can have a
significant impact on the dose standard deviation, already in a homogeneous
material. As an extension to this work, the time-dependent error modeling
introduced in section 6 might be extended to simulate spatial motion for 4D
treatment planning or even deformations in geometry components [1]. The
derived correlation models are of interest in their own right and can be
transferred to any dose calculation and uncertainty quantification framework
which supports multivariate error correlations [e.g. 1]. Also, the resulting
uncertainty estimates may be used for the purpose of robust optimization
during treatment planning.

While an extension to other types of uncertainties in the input parameters is
straightforward, the method is limited to errors which have a non-trivial
density function. Further work could explore the adaptation and application of
the model to uncertainties beyond errors in the patient set-up.

The main computational advantage of the proposed method is that it minimizes
the overhead associated with complex simulations. Therefore we cannot claim a
faster convergence rate and the speed-up is dependent on the specific
application and implementation. A combination with more sophisticated
scenario-based UQ methods, such as quadratures or a non-intrusive polynomial
chaos expansion, could however lead to further run-time improvements.

Further, we argue that computational performance can be improved through a
more efficient implementation. In particular, importance weighting can be
implemented on-the-fly using multiple sub-scoring routines, where each routine
accumulates the delivered dose for one scenario (see algorithm 1). Here,
performance gains should be achievable compared to the current post-processing
implementation, since the existing efficient structures of the Monte Carlo
code could be directly used. However, in this case the number of error
scenarios has to be predetermined, causing a loss in flexibility, and is also
limited with regard to memory constraints.

Lastly, the proposed (re-)weighting approach is not limited to proton therapy.
An application to uncertainties in radiation therapy using e.g. carbon ions or
photons is feasible and the method can in principal be relevant to any case of
input uncertainties in complex Monte Carlo simulations.

## Acknowledgments

The present contribution is supported by the Helmholtz Association under the
joint research school HIDSS4Health — Helmholtz Information and Data Science
School for Health.

## References

  * [1] M. Bangert, P. Hennig, and U. Oelfke.  Analytical probabilistic modeling for radiation therapy treatment planning.  Physics in Medicine and Biology, 58(16):5401–5419, July 2013.
  * [2] R. J. Beckman and M. D. McKay.  Monte carlo estimation under different distributions using the same simulation.  Technometrics, 29(2):153–160, 1987.
  * [3] J. L. Bedford.  Calculation of absorbed dose in radiotherapy by solution of the linear Boltzmann transport equations.  Physics in Medicine & Biology, 64(2):02TR01, Jan. 2019.
  * [4] E. Benhamou.  A few properties of sample variance.  working paper or preprint, Feb. 2019.
  * [5] R. E. Caflisch.  Monte Carlo and quasi-Monte Carlo methods.  Acta Numerica, 7:1–49, Jan. 1998.
  * [6] M. Casiraghi, F. Albertini, and A. Lomax.  Advantages and limitations of the ’worst case scenario’ approach in IMPT treatment planning.  Physics in medicine and biology, 58:1323–1339, Feb. 2013.
  * [7] M. Chu, Y. Zinchenko, S. G. Henderson, and M. B. Sharpe.  Robust optimization for intensity modulated radiation therapy treatment planning under uncertainty.  Physics in Medicine and Biology, 50(23):5463–5477, Dec. 2005.
  * [8] D. Craft, M. Bangert, T. Long, D. Papp, and J. Unkelbach.  Shared data for intensity modulated radiation therapy (IMRT) optimization research: The CORT dataset.  GigaScience, 3(1), Dec. 2014.
  * [9] B. Davison and J. B. Sykes.  Neutron Transport Theory.  Clarendon Press/Oxford University Press, Oxford/London, 1957.
  * [10] J. J. Duderstadt and W. R. Martin.  Transport Theory.  John Wiley & Sons, New York, 1979.
  * [11] R. Dürichen, M. A. F. Pimentel, L. Clifton, A. Schweikard, and D. A. Clifton.  Multitask Gaussian processes for multivariate physiological time-series analysis.  IEEE transactions on bio-medical engineering, 62(1):314–322, Jan. 2015.
  * [12] D. Duvenaud.  Automatic model construction with Gaussian processes.  PhD thesis, University of Cambridge, 2014.
  * [13] M. Frank.  Approximate models for radiative transfer.  Bulletin of the Institute of Mathematics. Academia Sinica. New Series, 2, Jan. 2007.
  * [14] A. Fredriksson and R. Bokrantz.  The scenario-based generalization of radiation therapy margins.  Physics in Medicine and Biology, 61, Oct. 2015.
  * [15] A. Fredriksson, A. Forsgren, and B. Hårdemark.  Minimax optimization for handling range and setup uncertainties in proton therapy.  Medical Physics, 38(3):1672–1684, Mar. 2011.
  * [16] K. A. Gifford, J. L. Horton, T. A. Wareing, G. Failla, and F. Mourtada.  Comparison of a finite-element multigroup discrete-ordinates code with Monte Carlo for radiotherapy calculations.  Physics in Medicine and Biology, 51(9):2253–2265, Apr. 2006.
  * [17] H. Hachiya, T. Akiyama, M. Sugiayma, and J. Peters.  Adaptive importance sampling for value function approximation in off-policy reinforcement learning.  Neural Networks, 22(10):1399–1410, 2009.
  * [18] W. K. Hastings.  Monte Carlo Sampling Methods Using Markov Chains and Their Applications.  Biometrika, 57(1):97–109, 1970.
  * [19] T. Hesterberg.  Weighted average importance sampling and defensive mixture distributions.  Technometrics, 37(2):185–194, 1995.
  * [20] N. Homma, M. Sakai, H. Endo, M. Mitsuya, Y. Takai, and M. Yoshizawa.  A new motion management method for lung tumor tracking radiation therapy.  WSEAS Transactions on Systems, 8:471–480, Apr. 2009.
  * [21] A. Hu and S. Qi.  Time Series Analysis of Interfraction Patient Setup in Image Guided Radiation Therapy.  International Journal of Radiation Oncology, Biology, Physics, 84(3):S736–S737, Nov. 2012.
  * [22] M. Hussein, C. Clark, and A. Nisbet.  Challenges in calculation of the gamma index in radiotherapy – Towards good practice.  Physica Medica, 36:1–11, Apr. 2017.
  * [23] K. Jabbari.  Review of Fast Monte Carlo Codes for Dose Calculation in Radiation Therapy Treatment Planning.  Journal of medical signals and sensors, 1:73–86, Jan. 2011.
  * [24] H. Kahn.  Random sampling (Monte Carlo) techniques in neutron attenuation problems–I.  Nucleonics, 6(5):27; passim, May 1950.
  * [25] T. Kanai, K. Kawachi, Y. Kumamoto, H. Ogawa, T. Yamada, H. Matsuzawa, and T. Inada.  Spot scanning system for proton radiotherapy.  Medical Physics, 7(4):365–369, 1980.
  * [26] T. Kollig and A. Keller.  Efficient multidimensional sampling.  Computer Graphics Forum, 21(3):557–563, 2002.
  * [27] A. C. Kraan, S. van de Water, D. N. Teguh, A. Al-Mamgani, T. Madden, H. M. Kooy, B. J. M. Heijmen, and M. S. Hoogeman.  Dose uncertainties in IMPT for oropharyngeal cancer in the presence of anatomical, range, and setup errors.  International Journal of Radiation Oncology, Biology, Physics, 87(5):888–896, Dec. 2013.
  * [28] S. Laine and T. Karras.  Stratified Sampling for Stochastic Transparency.  Computer Graphics Forum, 2011.
  * [29] P. L’Ecuyer.  Randomized quasi-monte carlo: An introduction for practitioners.  In A. B. Owen and P. W. Glynn, editors, Monte Carlo and Quasi-Monte Carlo Methods, pages 29–52, Cham, 2018. Springer International Publishing.
  * [30] W. Liu, X. Zhang, Y. Li, and R. Mohan.  Robust optimization of intensity modulated proton therapy.  Medical Physics, 39(2):1079–1091, Feb. 2012.
  * [31] A. Lomax.  Intensity modulation methods for proton radiotherapy.  Physics in Medicine and Biology, 44(1):185–205, Jan. 1999.
  * [32] A. J. Lomax.  Intensity modulated proton therapy and its sensitivity to treatment uncertainties 1: The potential effects of calculational uncertainties.  Physics in medicine and biology, 53(4):1027–1042, Feb. 2008.
  * [33] A. J. Lomax.  Intensity modulated proton therapy and its sensitivity to treatment uncertainties 2: The potential effects of inter-fraction and inter-field motions.  Physics in medicine and biology, 53(4):1043–1056, Feb. 2008.
  * [34] D. A. Low, W. B. Harms, S. Mutic, and J. A. Purdy.  A technique for the quantitative evaluation of dose distributions.  Medical Physics, 25(5):656–661, May 1998.
  * [35] J. Matoušek.  On the l2-discrepancy for anchored boxes.  Journal of Complexity, 14(4):527–556, 1998.
  * [36] S. E. McGowan, N. G. Burnet, and A. J. Lomax.  Treatment planning optimisation in proton therapy.  The British Journal of Radiology, 86(1021):20120288–20120288, Jan. 2013.
  * [37] A. B. Owen.  Randomly permuted (t, m, s)-nets and (t, s)-sequences.  In Monte Carlo and quasi-Monte Carlo methods in scientific computing, pages 299–317. Springer, 1995.
  * [38] A. B. Owen.  Monte carlo variance of scrambled net quadrature.  SIAM Journal on Numerical Analysis, 34(5):1884–1910, 1997.
  * [39] A. B. Owen.  Monte Carlo Theory, Methods and Examples.  2013\.
  * [40] H. Paganetti.  Range uncertainties in proton therapy and the role of Monte Carlo simulations.  Physics in Medicine and Biology, 57(11):R99–R117, June 2012.
  * [41] P. C. Park, J. P. Cheung, X. R. Zhu, A. K. Lee, N. Sahoo, S. L. Tucker, W. Liu, H. Li, R. Mohan, L. E. Court, and L. Dong.  Statistical assessment of proton treatment plans under setup and range uncertainties.  International journal of radiation oncology, biology, physics, 86(5):1007–1013, Aug. 2013.
  * [42] Z. Perkó, S. R. van der Voort, S. van de Water, C. M. H. Hartman, M. Hoogeman, and D. Lathouwers.  Fast and accurate sensitivity analysis of IMPT treatment plans using Polynomial Chaos Expansion.  Physics in Medicine and Biology, 61(12):4646–4664, May 2016.
  * [43] J. Perl, J. Shin, J. Schumann, B. Faddegon, and H. Paganetti.  TOPAS: An innovative proton Monte Carlo platform for research and clinical applications.  Medical Physics, 39(11):6818–6837, Nov. 2012.
  * [44] B. Perthame.  Transport Equations in Biology.  Frontiers in Mathematics. Birkhäuser Basel, 2007.
  * [45] L. Peshkin and C. R. Shelton.  Learning from scarce experience.  In Proceedings of the Nineteenth International Conference on Machine Learning, pages 498–505, 2002.
  * [46] D. Pflugfelder, J. J. Wilkens, and U. Oelfke.  Worst case optimization: A method to account for uncertainties in the optimization of intensity modulated proton therapy.  Physics in medicine and biology, 53(6):1689–1700, Mar. 2008.
  * [47] G. Poëtte.  A gPC-intrusive Monte-Carlo scheme for the resolution of the uncertain linear Boltzmann equation.  Journal of Computational Physics, May 2018.
  * [48] S. Shirakawa, Y. Akimoto, K. Ouchi, and K. Ohara.  Sample reuse via importance sampling in information geometric optimization.  arXiv preprint arXiv:1805.12388, 2018.
  * [49] I. Sobol’.  On the distribution of points in a cube and the approximate evaluation of integrals.  USSR Computational Mathematics and Mathematical Physics, 7(4):86–112, 1967.
  * [50] J. Spanier and E. M. Gelbard.  Monte Carlo Principles and Neutron Transport Problems.  Addison-Wesley Pub. Co., Reading, Mass., 1969.
  * [51] P. Stammer, L. Burigo, O. Jäkel, M. Frank, and N. Wahl.  Efficient uncertainty quantification for Monte Carlo dose calculations using importance (re-)weighting.  Physics in Medicine & Biology, 66(20):205003, Oct. 2021.
  * [52] A. D. Supanitsky and G. Medina-Tanco.  Effect of multiple reusing of simulated air showers in detector simulations.  Astroparticle Physics, 30(5):264–269, 2008.
  * [53] J. W. Tukey.  Configural polysampling.  Siam Review, 29(1):1–20, 1987.
  * [54] J. Unkelbach, T. C. Chan, and T. Bortfeld.  Accounting for range uncertainties in the optimization of intensity modulated proton therapy.  Physics in Medicine & Biology, 52(10):2755, 2007.
  * [55] O. N. Vassiliev.  The Boltzmann Equation.  In O. N. Vassiliev, editor, Monte Carlo Methods for Radiation Transport: Fundamentals and Advanced Topics, Biological and Medical Physics, Biomedical Engineering, pages 49–104. Springer International Publishing, Cham, 2017.
  * [56] M. T. Vilhena, C. F. Segatto, and L. B. Barichello.  A particular solution for the SN radiative transfer problems.  Journal of Quantitative Spectroscopy and Radiative Transfer, 53(4):467–469, Apr. 1995.
  * [57] N. Wahl, P. Hennig, H.-P. Wieser, and M. Bangert.  Efficiency of analytical and sampling-based uncertainty propagation in intensity-modulated proton therapy.  Physics in Medicine and Biology, 62(14):5790–5807, June 2017.
  * [58] X. Weng, Y. Yan, H. Shu, J. Wang, S. Jiang, and L. Luo.  A vectorized Monte Carlo code for radiotherapy treatment planning dose calculation.  Physics in medicine and biology, 48:N111–20, May 2003.
  * [59] H.-P. Wieser, E. Cisternas, N. Wahl, S. Ulrich, A. Stadler, H. Mescher, L.-R. Müller, T. Klinge, H. Gabrys, L. Burigo, A. Mairani, S. Ecker, B. Ackermann, M. Ellerbrock, K. Parodi, O. Jäkel, and M. Bangert.  Development of the open-source dose calculation and optimization toolkit matRad.  Medical physics, 44(6):2556–2568, June 2017.
  * [60] H.-P. Wieser, C. P. Karger, N. Wahl, and M. Bangert.  Impact of Gaussian uncertainty assumptions on probabilistic optimization in particle therapy.  Physics in Medicine & Biology, 65(14):145007, July 2020.
  * [61] D. Xiu and J. S. Hesthaven.  High-Order Collocation Methods for Differential Equations with Random Inputs.  SIAM Journal on Scientific Computing, 27(3):1118–1139, Jan. 2005\.