Importance Sampling for Cost-Optimized Estimation of Burn-Probability Maps in Wildfire Monte Carlo Simulations

1. Introduction

Monte Carlo simulations [1] are an approach to studying random processes by drawing samples from them. This approach has often been applied to studying the potential impacts of wildfire [2]: burn probability maps can be produced from Monte Carlo simulations, by averaging the burned areas of many randomly ignited fires simulated by a wildfire behavior model [3].

Simulating a large number of fires is required to ensure that the Monte Carlo estimator converges to an asymptotic value; otherwise, running the simulations again may produce a significantly different result. This random variability is known as the variance of the Monte Carlo estimator. Reducing the variance is required for making the estimation precise.

Importance Sampling is a general strategy for reducing the computational cost of Monte Carlo simulations. Importance Sampling redistributes the sampling probability of simulation outcomes, while inversely scaling their contribution to simulation results. This preserves the expected value of the Monte Carlo estimator while altering its variance, enabling faster convergence to the asymptotic estimand. In the case of wildfire simulations, importance sampling boosts or depresses the probability that some fires will be simulated, and accordingly decreases or increases their contribution to the aggregated simulation results. In [1] importance sampling is described as “the most fundamental variance reduction technique”. Therefore, it appears important for fire modelers to consider using Importance Sampling in their simulations; this article aims to help them understand its underlying principles, assess its potential benefits, and implement it when they are significant.

Importance Sampling has been used for sub-aspects of wildfire simulation. [4] have studied the use of Importance Sampling for data assimilation. [5] use Importance Sampling as part of the simulation of spotting behavior. We know of no academic study of Importance Sampling for accelerating the estimation of risk maps.

This paper is structured as follows. Section 2 will develop the conceptual methods introduced by this paper: this will mostly consist of framing the problem of computational efficiency in wildfire simulations, through concepts like cost-to-precision ratio and efficiency gain, introducing a generalized notion of variance (based on the $L_2$ norm) suited to estimating burn-probability maps (Section 2.1, Equation (2)), then recalling the definition of Importance Sampling (Section 2.4) as a technique for accelerating estimator convergence in Monte-Carlo simulations, adapting it to the context of wildfire simulations with $L_2$ variance, calculating its performance characteristics, considering various generalizations, and suggesting approaches to implementing the approach (Section 2.12). Section 3 will illustrate these methods will various case studies. Section 4 discusses the limitations of this approach. Table 1 summarizes the symbols most commonly encountered in formulas; their meaning will be defined in Section 2.

2. Materials and Methods

This analysis is conducted in the context of estimating burn probability maps over a large Area of Interest $\mathcal{A}$.

2.1. Background: Mean, Variance and Precision of Random Maps

For measuring the difference between multidimensional outputs like burn-probability maps, we use the $L_2$ norm $\parallel z\parallel$, defined by the spatial integral:

$${\parallel z\parallel }^2:={\int }_{s\in \mathcal{A}}{\left|z\left(s\right)\right|}^{2}ds$$

(1)

The distance between two maps $z_1$ and $z_2$ is then $\parallel z_1-z_2\parallel$.

We generalize the notion of variance from random numbers to random maps as follows. If $\widehat{z}$ is a random map, then we define the variance of $\widehat{z}$ to be:

$$\mathrm{Var}\left[\widehat{z}\right]:=\mathbb{E}\left[{∥\widehat{z}-\mathbb{E}\left[\widehat{z}\right]∥}^2\right]=\mathbb{E}\left[\parallel \widehat{z}{\parallel }^2\right]-{\parallel \mathbb{E}\left[\widehat{z}\right]\parallel }^2$$

(2)

Note that this can also be expressed in terms of pointwise variances - for each spatial location s, $\widehat{z}\left(s\right)$ is a random number of variance $\mathrm{Var}\left[\widehat{z}\left(s\right)\right]$, and we have

$$\mathrm{Var}\left[\widehat{z}\right]={\int }_{s\in \mathcal{A}}\mathrm{Var}\left[\widehat{z}\left(s\right)\right]ds$$

(3)

Many properties of scalar variance are retained: variance is additive on independent variables, and a random variable is constant if and only if its variance is zero.

A useful related concept is the precision, defined as the inverse variance:

$$\mathrm{Prec}\left[\widehat{z}\right]:={\displaystyle \frac{1}{\mathrm{Var}\left[\widehat{z}\right]}}$$

(4)

Aiming for low variance is equivalent to aiming for high precision.

2.2. Mean Estimation from Monte Carlo Simulations

Assume that we want to estimate the unknown expectation $\mu :=\mathbb{E}\left[\widehat{y}\right]$ of a distribution of interest, represented by a random map $\widehat{y}$. (Note that $\mu$ itself is a map, not a number). When all we know about $\widehat{y}$ is how to draw samples from it, a Monte Carlo simulation lets us estimate $\mu$ by drawing many i.i.d samples ${\widehat{y}}_k\sim \widehat{y}$, and computing their empirical mean ${\widehat{m}}_y$:

$${\widehat{m}}_y:={\displaystyle \frac{1}{N}}\sum _{k=1}^N{\widehat{y}}_k$$

(5)

Expectation being linear, it is easily seen that $\mathbb{E}\left[{\widehat{m}}_y\right]=\mathbb{E}\left[\widehat{y}\right]=\mu$. The point of drawing a large number of samples N is to increase the precision of that estimation, so that we can rely on ${\widehat{m}}_y\approx \mu$. Indeed, from the properties of variance, the variance and precision of ${\widehat{m}}_y$ are related to the sample size N as follows:

$$\begin{array}{cc}\hfill \mathrm{Var}\left[{\widehat{m}}_y\right]& ={\displaystyle \frac{1}{N^2}}\sum _{k=1}^N\mathrm{Var}\left[{\widehat{y}}_k\right]={\displaystyle \frac{1}{N}}\mathrm{Var}\left[\widehat{y}\right]\hfill \\ \hfill \mathrm{Prec}\left[{\widehat{m}}_y\right]& =N\mathrm{Prec}\left[\widehat{y}\right]\hfill \end{array}$$

(6)

The precision formula is insightful: we can think of each sample contributing $\mathrm{Prec}\left[\widehat{y}\right]$ to estimator precision. Thus, $\mathrm{Prec}\left[\widehat{y}\right]$ can be viewed as a precision per sample performance characteristic of our sampling distribution.

Assume now that with each sample ${\widehat{y}}_k$ is associated a computational cost${\widehat{c}}_k^{\left(y\right)}$ of drawing it - typically ${\widehat{c}}_k^{\left(y\right)}$ is expressed in computer time. For example, in wildfire simulation, large fires are often more costly to simulate than small fires. Then the cost of drawing a large number N of i.i.d samples will be approximately $N\mathbb{E}\left[{\widehat{c}}^{\left(y\right)}\right]$. This motivates the following concept of cost-to-precision ratio, which quantifies the overhead of our estimator:

$$r_{C/P}\left[\widehat{y}\right]={\displaystyle \frac{\mathbb{E}\left[{\widehat{c}}^{\left(y\right)}\right]}{\mathrm{Prec}\left[\widehat{y}\right]}}=\mathbb{E}\left[{\widehat{c}}^{\left(y\right)}\right]\mathrm{Var}\left[\widehat{y}\right]$$

(7)

2.3. Burn-Probability Maps from Monte Carlo Simulations

We assume the following structure for how simulations run. Each simulation k draws initial conditions ${\widehat{I}}_k$ (typically, ${\widehat{I}}_k$ encompasses the space-time location of ignition, and the surrounding weather conditions), from which a fire ${\widehat{F}}_k$ is simulated to yield a binary burn map $b_{{\widehat{F}}_k}$, in which

$$b_{{\widehat{F}}_k}\left(s\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{\widehat{F}}_k\mathrm{burned}\mathrm{location}s\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

Thus, each $b_{{\widehat{F}}_k}\left(s\right)$ is a Bernoulli random variable, and its expectation $\mu \left(s\right):=\mathbb{E}\left[b_{{\widehat{F}}_k}\left(s\right)\right]$ is the probability that one simulation burns location s.

We assume that the ${\widehat{F}}_k$ are sampled i.i.d., and so will write $\widehat{F}$ to denote a random variable representing their common distribution. Likewise for the ${\widehat{I}}_k$.

The goal of the simulations is to estimate the expected burn probability, which is defined as the expectation of the burn map, therefore yielding a map of burn probabilities:

$$\begin{array}{cc}\hfill \mu & :=\mathbb{E}\left[b_{\widehat{F}}\right]\hfill \\ \hfill \mu \left(s\right)& =\mathbb{P}\left[b_{\widehat{F}}\left(s\right)=1\right]\hfill \end{array}$$

(8)

That estimation is typically done by using the Crude Monte Carlo (CMC) estimator ${\widehat{m}}_{\mathrm{CMC}}$, which draws a large number N of fires and computes the average of the $b_{{\widehat{F}}_k}$:

$${\widehat{m}}_{\mathrm{CMC}}:={\displaystyle \frac{1}{N}}\sum _{k=1}^N{b}_{{\widehat{F}}_k}$$

(9)

Because each $b_{\widehat{F}}\left(s\right)$ is a Bernoulli random variable (taking only 0 and 1 as values), $\mathrm{Var}\left[b_{\widehat{F}}\right]$ has an interesting expression in terms of expected fire size:

$$\mathrm{Var}\left[b_{\widehat{F}}\right]=\mathbb{E}\left[\parallel b_{\widehat{F}}{\parallel }^2\right]-{\parallel \mu \parallel }^2=\mathbb{E}\left[a\left(\widehat{F}\right)\right]-{\parallel \mu \parallel }^2$$

(10)

in which $a\left(\widehat{F}\right)\ge 0$ is the (random) fire size (say, in ${\mathrm{m}}^2$ or acres, or whatever unit was chosen for $ds$). More precisely, $a\left(\widehat{F}\right)$ is the area burned by the random fire within the AOI $\mathcal{A}$.

The above derivation relies on $\parallel b_{\widehat{F}}{\parallel }^2=a\left(\widehat{F}\right)$, which we now prove. Denoting $\widehat{\mathcal{B}}\subseteq \mathcal{A}$ the (random) set of locations burned by a simulated fire, and observing that $|b_{\widehat{F}}{\left(s\right)|}^2=b_{\widehat{F}}\left(s\right)$, we have:

$$\parallel b_{\widehat{F}}{\parallel }^2={\int }_{s\in \mathcal{A}}|b_{\widehat{F}}{\left(s\right)|}^{2}ds={\int }_{s\in \mathcal{A}}{b}_{\widehat{F}}\left(s\right)ds={\int }_{s\in \widehat{\mathcal{B}}}1ds+{\int }_{s\notin \widehat{\mathcal{B}}}0ds=a\left(\widehat{F}\right)+0$$

(11)

Another insightful relationship is that the expected fire size is the spatial integral of the expected burn probability:

$$\mathbb{E}\left[a\left(\widehat{F}\right)\right]={\int }_{s\in \mathcal{A}}\mu \left(s\right)ds$$

(12)

Finally, we note a useful approximation. In the common settings where the burn probabilities $\mu \left(s\right)$ are low (each location s is burned only by a small minority of simulated fires), then $\mu \left(s\right)\ll 1$ implies that $\mathrm{Var}\left[b_{\widehat{F}}\left(s\right)\right]=\mu \left(s\right)(1-\mu \left(s\right))\approx \mu \left(s\right)\gg \mu {\left(s\right)}^2$, and since $\mu \left(s\right)=\mathbb{E}\left[|b_{\widehat{F}}{\left(s\right)|}^2\right]$ it follows that ${\parallel \mu \parallel }^2\ll \mathbb{E}\left[a\left(\widehat{F}\right)\right]$, and

$$\mathrm{Var}\left[b_{\widehat{F}}\right]\approx \mathbb{E}\left[a\left(\widehat{F}\right)\right]\gg {\parallel \mu \parallel }^2$$

(13)

In other words, the variance of the binary burn map is well approximated by the expected fire size.

2.4. Using Importance Sampling for Reducing the Cost-to-Precision Ratio

We now borrow from the concept of Importance Sampling from the generic Monte Carlo literature [1]. This will enable us to design other Monte Carlo estimators for the burn-probability $\mu$, with lower cost-to-precision ratio than the CMC estimator, thus requiring less computation to achieve convergence.

Denote $p\left(F\right)$ the “natural” probability distribution from which simulated fires are sampled (i.e. $\widehat{F}\sim p$), which is hopefully representative of real-world fire occurrence. The CMC estimator consists of drawing samples ${\widehat{F}}_k$ from p then averaging the $b_{{\widehat{F}}_k}$ maps. Importance Sampling modifies this procedure by drawing samples ${\widehat{F}}_k^{\left(q\right)}$ from another distribution $q\left(F\right)$ (called the proposal distribution) and reweighting the samples to remain representative of p, yielding the Likelihood Ratio estimator1${\widehat{m}}_q$:

$${\widehat{m}}_q:={\displaystyle \frac{1}{N}}\sum _{k=1}^N{w}_q\left({\widehat{F}}_k^{\left(q\right)}\right)b_{{\widehat{F}}_k^{\left(q\right)}}$$

(14)

in which we defined the importance weight as the ratio of probabilities:

$$w_q\left(F\right)={\displaystyle \frac{p\left(F\right)}{q\left(F\right)}}$$

(15)

It is easily verified that this change of estimator preserves the expected value:

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{m}}_q\right]& =\mathbb{E}\left[{\displaystyle \frac{p\left(b_{{\widehat{F}}^{\left(q\right)}}\right)}{q\left(b_{{\widehat{F}}^{\left(q\right)}}\right)}}{b}_{{\widehat{F}}^{\left(q\right)}}\right]\hfill \\ \hfill & ={\int }_{F}dFq\left(F\right){\displaystyle \frac{p\left(F\right)}{q\left(F\right)}}{b}_F\hfill \\ \hfill & ={\int }_{F}dFp\left(F\right)b_F\hfill \\ \hfill & =\mathbb{E}\left[b_{\widehat{F}}\right]=\mu =\mathbb{E}\left[{\widehat{m}}_{\mathrm{CMC}}\right]\hfill \end{array}$$

(16)

What can make this new estimator valuable is its variance:

$$\begin{array}{cc}\hfill \mathrm{Var}\left[{\widehat{m}}_q\right]& ={\displaystyle \frac{1}{N}}\left(\mathbb{E}\left[w_q{\left({\widehat{F}}^{\left(q\right)}\right)}^2{\parallel b_{{\widehat{F}}^{\left(q\right)}}\parallel }^2\right]-{\parallel \mu \parallel }^2\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{N}}\left(\mathbb{E}\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]-{\parallel \mu \parallel }^2\right)\hfill \end{array}$$

(17)

The approach of Importance Sampling consists of seeking a proposal distribution q which minimizes the cost-to-precision ratio $r_{C/P}^{\left(q\right)}:=r_{C/P}\left[{\displaystyle \frac{p\left({\widehat{F}}^{\left(q\right)}\right)}{q\left({\widehat{F}}^{\left(q\right)}\right)}}{b}_{{\widehat{F}}^{\left(q\right)}}\right]$. When the cost $\widehat{c}$ is constant, this amounts to minimizing variance, in which case a classical result of the Monte Carlo literature says that the optimal proposal distribution $q_{\mathrm{MinVar}}^{*}$ and the corresponding variance are given by:

$$\begin{array}{cc}\hfill q_{\mathrm{MinVar}}^{*}\left(F\right)& \propto \parallel b_F\parallel p\left(F\right)\hfill \\ \hfill & =\mathbb{E}{\left[\parallel b_{\widehat{F}}\parallel \right]}^{-1}\parallel b_F\parallel p\left(F\right)\hfill \\ \hfill \mathrm{Var}\left[{\widehat{m}}_{q_{\mathrm{MinVar}}^{*}}\right]& ={\displaystyle \frac{1}{N}}\left(\mathbb{E}{\left[\parallel b_{\widehat{F}}\parallel \right]}^2-{\parallel \mu \parallel }^2\right)\hfill \end{array}$$

(18)

This is an improvement over $\mathrm{Var}\left[{\widehat{m}}_{\mathrm{CMC}}\right]$ because $\mathbb{E}{\left[\parallel b_{\widehat{F}}\parallel \right]}^2\le \mathbb{E}\left[\parallel b_{\widehat{F}}{\parallel }^2\right]$, as can be seen for example by applying Jensen’s inequality [6].

We place ourselves in a more general context, corresponding to a non-constant cost $\widehat{c}$. More precisely, we will assume:

• the Poisson process approximation (justified in Section 2.14):

  $${\parallel \mu \parallel }^2\ll \mathbb{E}\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]$$

  (19)
• the cost $\widehat{c}$ being a function of the simulated fire:

  $$\widehat{c}=c\left(\widehat{F}\right)$$

  (20)

The cost-to-precision ratio achieved by a proposal distribution q is then:

$$\begin{array}{cc}\hfill r_{C/P}^{\left(q\right)}& =\mathbb{E}\left[c\left({\widehat{F}}^{\left(q\right)}\right)\right]\mathbb{E}\left[{\displaystyle \frac{p\left(b_{\widehat{F}}\right)}{q\left(b_{\widehat{F}}\right)}}{\parallel b_{\widehat{F}}\parallel }^2\right]=\mathbb{E}\left[{\displaystyle \frac{q\left(b_{\widehat{F}}\right)}{p\left(b_{\widehat{F}}\right)}}c\left(\widehat{F}\right)\right]\mathbb{E}\left[{\displaystyle \frac{p\left(b_{\widehat{F}}\right)}{q\left(b_{\widehat{F}}\right)}}{\parallel b_{\widehat{F}}\parallel }^2\right]\hfill \\ \hfill r_{C/P}^{\left(q\right)}& =\mathbb{E}\left[w_q{\left(\widehat{F}\right)}^{-1}c\left(\widehat{F}\right)\right]\mathbb{E}\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]\hfill \end{array}$$

(21)

For practical purposes, we often need to know the weight function $w_q(·)$ only up to a proportionality constant. For example, we can compute the efficiency gain $g_q$ of proposal distribution q as:

$$g_q:={\displaystyle \frac{r_{C/P}^{\left(p\right)}}{r_{C/P}^{\left(q\right)}}}={\displaystyle \frac{\mathbb{E}\left[c\left(\widehat{F}\right)\right]\mathbb{E}\left[\parallel b_{\widehat{F}}{\parallel }^2\right]}{\mathbb{E}\left[w_q{\left(\widehat{F}\right)}^{-1}c\left(\widehat{F}\right)\right]\mathbb{E}\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]}}$$

(22)

The efficiency gain $g_q$ quantifies how the cost-to-precision ratio improves from the Crude Monte Carlo estimator (${\widehat{m}}_{\mathrm{CMC}}$) to Importance Sampling (${\widehat{m}}_q$). We see that replacing w by $\lambda w$ in Equation (22) leaves the result unchanged, because the constant $\lambda$ gets simplified away in the denominator. Likewise for equation .

2.5. Computing Sampling Frequency Multipliers

To better understand the effect of Importance Sampling with some reweighting function $w:=p/q$, we introduce the notion of frequency multipliers, which describe how the sampling frequencies of simulated fires change when we use Importance Sampling.

The same-cost frequency multipliers correspond to the following procedure. Assume that we initially planned to simulate N fires with the natural distribution p, which means an expected computational cost of $N\mathbb{E}\left[c\left(\widehat{F}\right)\right]$. We then choose to adopt Importance Sampling with reweighting function w (or equivalently, with proposal density q), and consequently adjust the number of simulated fires to $N^{\prime }$ to keep the same expected computational cost, i.e. $N^{\prime }\mathbb{E}\left[c\left({\widehat{F}}^{\left(q\right)}\right)\right]=N\mathbb{E}\left[c\left(\widehat{F}\right)\right]$. For any potential simulated fire b, we can then ask the question: how much more frequently does b get drawn by Importance Sampling than by our original sampling design? This number is the same-cost frequency multiplier of b, and can be computed as:

$${\displaystyle \frac{N^{\prime }q\left(F\right)}{Np\left(F\right)}}={\displaystyle \frac{\mathbb{E}\left[c\left(\widehat{F}\right)\right]}{\mathbb{E}\left[w_q{\left(\widehat{F}\right)}^{-1}c\left(\widehat{F}\right)\right]}}{w}_q{\left(F\right)}^{-1}$$

(23)

Similarly, we can conceive of the same-precision frequency multiplier. Then the number of simulated fires is adjusted from N to $N"$ so that precision remains the same, i.e. $N^{\prime }\mathbb{E}{\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]}^{-1}=N\mathbb{E}{\left[{\parallel b_{\widehat{F}}\parallel }^2\right]}^{-1}$. Then the frequency multiplier gets computed as:

$${\displaystyle \frac{N^{\prime }q\left(F\right)}{Np\left(F\right)}}={\displaystyle \frac{\mathbb{E}\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]}{\mathbb{E}\left[{\parallel b_{\widehat{F}}\parallel }^2\right]}}{w}_q{\left(F\right)}^{-1}$$

(24)

Confusion pitfall: note that the frequency multipliers are inversely proportional to the reweighting function w.

2.6. The Optimal Proposal Distribution

This section shows that there is an optimal proposal distribution, and that it can be described quantitatively. Applying the probabilistic Cauchy-Schwarz inequality (Appendix A.2) to Equation (2.4) yields an expression for the optimal proposal distribution $q_c^{*}$ and its cost-to-precision ratio:

$$w_{q_c^{*}}{\left(F\right)}^2\propto {\displaystyle \frac{c\left(F\right)}{\parallel b_F{\parallel }^2}}$$

$$\begin{array}{c}\hfill q_c^{*}\left(F\right)=\mathbb{E}{\left[{\displaystyle \frac{\parallel b_{\widehat{F}}\parallel }{c{\left(\widehat{F}\right)}^{\frac{1}{2}}}}\right]}^{-1}{\displaystyle \frac{\parallel b_F\parallel }{c{\left(F\right)}^{\frac{1}{2}}}}p\left(F\right)\end{array}$$

(25)

$$r_{C/P}^{\left(q_c^{*}\right)}=\mathbb{E}{\left[c{\left(\widehat{F}\right)}^{{\displaystyle \frac{1}{2}}}\parallel b_{\widehat{F}}\parallel \right]}^2$$

(26)

We define the maximum potential efficiency gain $g^{*}$ as the efficiency gain (Equation (22)) achieved by the optimal proposal distribution $q_c^{*}$. By combining equations 22 and 26, we obtain the following expression for $g^{*}$:

$$g^{*}={\displaystyle \frac{\mathbb{E}\left[c\left(\widehat{F}\right)\right]\mathbb{E}\left[\parallel b_{\widehat{F}}{\parallel }^2\right]}{\mathbb{E}{\left[c{\left(\widehat{F}\right)}^{\frac{1}{2}}\parallel b_{\widehat{F}}\parallel \right]}^2}}$$

(27)

2.7. Application to Burn-Probability Maps

In the case of burn maps, recalling that $\parallel b_{\widehat{F}}{\parallel }^2=a\left(\widehat{F}\right)$, it follows that

$$q_c^{*}\left(F\right)=\mathbb{E}{\left[{\left({\displaystyle \frac{a\left(\widehat{F}\right)}{c\left(\widehat{F}\right)}}\right)}^{{\displaystyle \frac{1}{2}}}\right]}^{-1}{\left({\displaystyle \frac{a\left(F\right)}{c\left(F\right)}}\right)}^{{\displaystyle \frac{1}{2}}}p\left(F\right)$$

(28)

In other words, the best proposal distribution reweights the probability of each fire by the square root of its size-to-cost ratio. Two special cases of cost functions $c\left(F\right)$ stand out:

• when the cost of simulating a fire is proportional to its size, then Importance Sampling makes no difference: the initial distribution $p\left(F\right)$ is already the optimal proposal distribution;
• when the cost of simulating a fire is constant, then the optimal proposal distribution reweights the probability of each fire by the square root of its size.

It is also interesting to consider a mixture of these cases. The cost function can be modeled as the sum of a constant term and a term proportional to fire size: $c\left(F\right)=c_0(1+a\left(F\right)/A_0)$. For very large fires ($A\gg A_0$), the size-to-cost ratio will be very close to 1, such that no reweighting will occur: large fires will be sampled according to their natural distribution. For very small fires ($A\approx 0$), the size-to-cost ratio will be very small, such that these fires will have virtually no chance of being sampled; and since they do not contribute much to the mean, the weight will not be high enough to compensate this sampling scarcity. This provides a principled justification for the common heuristic of dispensing with simulating small fires: provided that $A_0$ is not too large, that heuristic is very close to the optimal proposal distribution.

To further illustrate this phenomenon, Figure 1 displays the maximum potential efficiency gain (Equation (27)) for an imaginary model in which the cost c is determined by the fire size A as $c\propto (1+{(a\left(F\right)/A_0)}^{ϵ})$, based on the distribution of fire sizes recorded in the Fire Occurrence Database [7]. When the constant term is high (high $A_0$), a lot of efficiency can be gained by reducing the sampling frequency of small fires. When the cost is a superlinear function of size ($ϵ>1$), efficiency can be gained by reducing the sampling frequency of large fires.

2.8. Generalization: Optimal Proposal Distribution for Ignitions

We now assume that we do not have enough mastery over $p\left(F\right)$ to find the optimal proposal distribution $q_c^{*}$. Instead, we restrict ourselves to seeking a proposal distribution for the ignition $\widehat{I}$, such that $q(F,I)$ is of the form:

$$q(F,I)=p\left(F\right|I)q^I\left(I\right)$$

(29)

This decomposition is consistent with the idea that we’re going to change the distribution of ignitions without changing the fire spread algorithm itself. For example, Section 3.4 will change the sampling distribution of ignitions based on fire duration. Under this constraint, the cost-to-precision ratio of the proposal distribution can then be expressed as:

$$r_{C/P}^{\left(q\right)}=\mathbb{E}\left[{\displaystyle \frac{q^I\left(\widehat{I}\right)}{p\left(\widehat{I}\right)}}{c}_I\left(\widehat{I}\right)\right]\mathbb{E}\left[{\displaystyle \frac{p\left(\widehat{I}\right)}{q^I\left(\widehat{I}\right)}}{a}_I\left(\widehat{I}\right)\right]$$

(30)

in which the functions $a_I\left(I\right),c_I\left(I\right)$ compute the conditional expectation of fire size and computational cost:

$$\begin{array}{cc}\hfill a_I\left(I\right)& :=\mathbb{E}\left[\parallel b_{\widehat{F}}{\parallel }^2|\widehat{I}=I\right]=\mathbb{E}\left[a\left(\widehat{F}\right)|\widehat{I}=I\right]\hfill \\ \hfill c_I\left(I\right)& :=\mathbb{E}\left[c\left(\widehat{F}\right)|\widehat{I}=I\right]\hfill \end{array}$$

(31)

The proof to Equation (25) can then be generalized to show that the optimal proposal distribution $q_c^{I*}$ is given by:

$$q_c^{I*}\left(\widehat{I}\right)\propto {\left({\displaystyle \frac{a_I\left(\widehat{I}\right)}{c_I\left(\widehat{I}\right)}}\right)}^{{\displaystyle \frac{1}{2}}}p\left(\widehat{I}\right)$$

(32)

It can then be seen that $q_c^{I*}$ yields the following density and cost-to-precision ratio:

$$\begin{array}{cc}\hfill q_c^{I*}(F,I)& =\mathbb{E}{\left[{\left({\displaystyle \frac{a_I\left(\widehat{I}\right)}{c_I\left(\widehat{I}\right)}}\right)}^{{\displaystyle \frac{1}{2}}}\right]}^{-1}{\left({\displaystyle \frac{a_I\left(I\right)}{c_I\left(I\right)}}\right)}^{{\displaystyle \frac{1}{2}}}p(F,I)\hfill \\ \hfill r_{C/P}^{\left(q_c^{I*}\right)}& =\mathbb{E}{\left[{\left(c_I\left(\widehat{I}\right)a_I\left(\widehat{I}\right)\right)}^{{\displaystyle \frac{1}{2}}}\right]}^2\hfill \end{array}$$

(33)

Note that any random variable $\widehat{I}$ can be chosen for these results to apply. Here we have described $\widehat{I}$ as the “ignition” or “initial conditions”, but it can be any predictor of fire size and simulation cost. Two extreme cases of $\widehat{I}$ are noteworthy:

• If $b_{\widehat{F}}$ is fully determined by $\widehat{I}$ (deterministic simulation), then we recover $q_c^{I*}=q_c^{*}$.
• If $\widehat{I}$ is uninformative (independent) w.r.t $b_{\widehat{F}}$, then we obtain a no-progress optimum of $q_c^{I*}=p$.

2.9. Generalization: Impact-Weighted $L_2$ Distance

An interesting generalization consists of using another measure than area for computing $L_2$ distances, by introducing a density function $u\left(s\right)$ which quantifies the “value” or “concern” per unit area:

$${\parallel z\parallel }_u^2:={\int }_s{\left|z\left(s\right)\right|}^{2}u\left(s\right)ds$$

(34)

Then all the results of the previous section apply (with $ds$ being replaced by $u\left(s\right)ds$ in integrals), except that the interpretation changes: quantities like $\parallel b_{\widehat{F}}{\parallel }_u^2$, $a\left(\widehat{F}\right)$ and $a_I(·)$ no longer measure fire size, but burned value. In particular, when using the optimal proposal distribution for importance sampling, each fire sees its sampling probability reweighted by the square root of its burned value/simulation cost ratio.

Weighting by value instead of area can dramatically change the optimal proposal distribution. It can boost the probability of sampling small fires near valuable assets, while depressing the probability of large fires in low-concern areas. From this perspective, “Area of Interest” takes on a very precise meaning: the set of locations s where $u\left(s\right)>0$.

The choice of value measure $u\left(s\right)ds$ depends on how the burn-probability map will be used: where do we want precision? For example, the burn-probability maps might be used to plot the risk to properties and infrastructure, and also be aggregated to estimate GHG and smoke emissions. To that end, a good value function could place strong measure on properties and infrastructure, and mix that with small measure on fuel loads (the latter do not need as high a measure because emissions are destined to be aggregated, which already reduces variance).

2.10. Generalization: Non-Geographic Maps and Non-Binary Outcomes

Another way to understand Section 2.9 is to pretend that the variable $s\in \mathcal{A}$ ranges not over geographic space, but over “assets space”. Attentive readers may have noticed that nothing in the above analysis is specific to the spatial nature of the set $\mathcal{A}$. We can conceive of other type of burn-probability maps: for example, $\mathcal{A}$ may be:

• a space-time region, like $\mathcal{A}:=\mathrm{California}\times [2024,2034)$ or $\mathcal{A}:=\mathrm{Alberta}\times \mathrm{Months}$, in which $\mathrm{Months}:=\{\mathrm{January},\mathrm{February},\cdots ,\mathrm{December}\}$;
• a combination of geographic space and fire severity buckets, like $\mathcal{A}:=\mathrm{Australia}\times \{\mathrm{Severity}\mathrm{Class}1,\mathrm{Severity}\mathrm{Class}2,\cdots \}$;
• a combination of geographic space and paired scenario outcomes, like:

  $$\begin{array}{cc}\hfill \mathcal{A}:=& \mathrm{Area}\mathrm{of}\mathrm{Interest}\hfill \\ \hfill & \times \{“\mathrm{burned}\mathrm{without}\mathrm{fuel}\mathrm{treatment}”,“\mathrm{did}\mathrm{not}\mathrm{burn}\mathrm{without}\mathrm{fuel}\mathrm{treatment}”\}\hfill \\ \hfill & \times \{“\mathrm{burned}\mathrm{with}\mathrm{fuel}\mathrm{treatment}”,“\mathrm{did}\mathrm{not}\mathrm{burn}\mathrm{with}\mathrm{fuel}\mathrm{treatment}”\}\hfill \end{array}$$

The above paragraph generalizes $\mathcal{A}$, that is the domain of function $s\mapsto b_{\widehat{F}}\left(s\right)$; we can also generalize its codomain, that is the simulation outcomes. In particular:

• $b_{\widehat{F}}\left(s\right)$need not be binary: for example, $b_{\widehat{F}}\left(s\right)$ might be a density of greenhouse-gas emissions;
• $b_{\widehat{F}}\left(s\right)$need not be about fire.

What matters for the analysis of this article to apply is that the $b_{\widehat{F}}$ outcomes be aggregated by averaging, and for the Poisson process approximation to be valid (Section 2.14). Note that when the outcomes are not binary, $\parallel b_F{\parallel }^2$ can no longer be interpreted as the impacted area, i.e. Equation (11) ($a\left(F\right)=\parallel b_F{\parallel }^2$) no longer holds.

2.11. Variant: Sampling from a Poisson Point Process

The above analysis assumes that the number of fires is fixed, and that the fires are sampled i.i.d. A variant of this sampling scheme is to sample from a Poisson Point Process (PPP) [8] instead: this is equivalent to saying that the number of fires is no longer fixed, but follows a Poisson distribution. The PPP can be appealing because it is algorithmically convenient (subregions of space-time can be sampled independently) and yields mathematical simplifications, as we will see below.

In the context of a PPP, a random number ${\widehat{N}}_q\sim \mathrm{Poisson}(\lambda =N)$ of fires are sampled, and the Likelihood Ratio Estimator, which we now denote ${\ddot{m}}_q$, is virtually unchanged:

$${\ddot{m}}_q:={\displaystyle \frac{1}{N}}\sum _{k=1}^{{\widehat{N}}_q}{w}_q\left({\widehat{F}}_k^{\left(q\right)}\right)b_{{\widehat{F}}_k^{\left(q\right)}}$$

(35)

By writing $\mathbb{E}\left[{\ddot{m}}_q\right]=\mathbb{E}\left[\mathbb{E}\left[{\ddot{m}}_q\right|{\widehat{N}}_q]\right]$, it is easily verified that the mean of this estimator is correct: $\mathbb{E}\left[{\ddot{m}}_q\right]=\mu$.

An interesting aspect of the PPP is that the variance approximation of Equation (19) becomes exact2, such that:

$$\mathrm{Var}\left[{\ddot{m}}_q\right]={\displaystyle \frac{1}{N}}\mathbb{E}\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]$$

(36)

It follows that the results of the previous section about the optimal proposal distribution and cost-to-precision ratio remain unchanged.

2.12. Practical Implementation

The previous sections described the optimal proposal distribution which can be achieved with perfect knowledge of the impact and cost functions $a_I(·)$ and $c_I(·)$. In practice, these functions will be unknown, so that Equation (32) is not applicable. Instead, we suggest using a calibration approach to estimate an approximately optimal proposal distribution, through the following steps:

• Run simulations to draw a sample of ignition and simulated fire $(\widehat{I},\widehat{F})$, yielding an empirical distribution which we will represent by the random variables $(\tilde{I},\tilde{F})$.
• Apply Equation (2.4) to $\tilde{F}$ instead3 of $\widehat{F}$ to estimate the potential gain in cost-to-precision ratio. If that gain is not significant, abandon Importance Sampling.
• Design a parametric family of weights functions $w(I;\theta )$.
• Find the optimal parameter vector ${\theta }^{*}$ by solving the following optimization problem, inspired by Equation (2.4):

  $${\theta }^{*}:=\mathrm{argmin}\left(\theta \mapsto \mathbb{E}\left[w{(\tilde{I};\theta )}^{-1}c\left(\tilde{F}\right)\right]\mathbb{E}\left[w(\tilde{I};\theta )a\left(\tilde{F}\right)\right]\right)$$

  (37)
  If ${\theta }^{*}$ does not yield a significant improvement compared to $w\left(I\right)=1$, try a new family of weights functions, or abandon Importance Sampling.
• Adopt $q^I\left(I\right)\propto w{(I,{\theta }^{*})}^{-1}p\left(I\right)$ as the proposal distribution for Importance Sampling.

Care should be taken to appropriately constrain the ${\left(w(·;\theta )\right)}_{\theta }$ parametric family to avoid any problem of identifiability, since the objective function in Equation (37) will yield exactly the same score for a reweighting function w and a scaled variant $\lambda w$: therefore, the family of reweighting functions must not be stable by rescaling. It might also be sensible to add a regularization term to the objective function.

In addition to using a parametric family, it can be sensible to use a curve-fitting approach: one can compute the weights which are optimal for the (discrete) empirical distribution, then use some interpolation algorithm to derive a continuous weight function, like a supervised learning algorithm or nearest-neighbor method [9].

2.13. Analysis of Geographic Inequalities in Precision

The above analysis describes how to optimize the global precision (in the sense of the $L_2$ norm); however, doing so might alter the relative geographic distribution of local precision across the map. In particular, we can expect that this optimization will favor absolute precision over relative precision, spending the computational cost mostly in areas of high burn probability. Section 3.3 provides an example that illustrates this phenomenon.

2.14. Justifying the Poisson Process Approximation

The proof to Equation (25) relies on the Poisson process approximation (Equation (19)), which we now justify. It follows from Equation (18) that $\mathbb{E}{\left[\parallel b_{\widehat{F}}\parallel \right]}^2$ is the smallest possible value for $\mathbb{E}\left[w_q\left(\widehat{F}\right){\parallel b_{\widehat{F}}\parallel }^2\right]$. Therefore, it is enough to show that $\mathbb{E}{\left[\parallel b_{\widehat{F}}\parallel \right]}^2\gg {\parallel \mu \parallel }^2$. We suggest that this generally holds, by considering the case where $\mu$ is spatially uniform. Denoting $\left|\mathcal{A}\right|$ the total area, we then have $\mu \left(s\right)={\left|\mathcal{A}\right|}^{-1}\mathbb{E}\left[a\left(\widehat{F}\right)\right]$, and ${\parallel \mu \parallel }^2\ll \mathbb{E}{\left[a{\left(\widehat{F}\right)}^{{\displaystyle \frac{1}{2}}}\right]}^2$ is equivalent to the following condition on the burned fraction$a\left(\widehat{F}\right)/\left|\mathcal{A}\right|$:

$$\mathbb{E}\left[{\displaystyle \frac{a\left(\widehat{F}\right)}{\left|\mathcal{A}\right|}}\right]\ll \mathbb{E}\left[{\left({\displaystyle \frac{a\left(\widehat{F}\right)}{\left|\mathcal{A}\right|}}\right)}^{{\displaystyle \frac{1}{2}}}\right]$$

(38)

... which will be the case in the typical situation where the burned fraction is almost certainly much smaller than 1. The general phenomenon here is that ${\parallel \mu \parallel }^2$ becomes an insignificant component of variance as the AOI is made much larger than the fires; this corresponds to the Binomial-to-Poisson approximation we showed in a previous paper [10].

Another perspective to shed light on this phenomenon is that making the AOI very large means that the sampling procedure behaves like a Poisson Point Process (Section 2.11) in any minority subregion of the AOI; then Equation (36) holds approximately. This is the same phenomenon as the Binomial-to-Poisson approximation: the number of points drawn in a small subregion follows a Binomial distribution with a small probability parameter, which is well approximated by a Poisson distribution.

2.15. Importance Sampling versus Stratified or Neglected Ignitions

Stratification [11] is another strategy for reducing the variance of Monte Carlo simulations. In the case of burn-probability estimation, stratification might consist of partitioning the set of possible ignitions into several subsets (for example, low-risk and high-risk), and simulating more fire seasons in some subsets than others. The burn probability maps from each subset are then combined by a weighted sum, the weights being inversely proportional to the number of fire seasons. A similar strategy consists of dividing the space of ignitions between a low-risk and high-risk subsets, and simulating only the high-risk subset, thus neglecting the other.

Observe that both strategies are highly similar to a special case of Importance Sampling, in which the weight function $w_q(·)$ is piecewise constant. Therefore, Importance Sampling is useful for reasoning about these other strategies, and provides an upper bound to their potential efficiency gains.

3. Results

We illustrate the above methods by applying them to an existing dataset of results from a wildfire Monte-Carlo simulation. We will implement various forms of Importance Sampling, estimate the efficiency gains, and discuss the trade-offs. Due to limited access to computed resources, we did not run new simulations based on the reweighted sampling scheme; instead, we use simulation results that were drawn without Importance Sampling, and quantity the efficiency gains that would have been obtained by implementing Importance Sampling.

3.1. Simulation Design and Results

The dataset of simulated fires studied here was produced as part of the 3rd version of the FireFactor project [12]. We will not try to discuss the validity and implementation details of these simulations here - that would require its own publication - and we will instead focus solely on aspect relevant to our study of Importance Sampling. 800 repetitions of the 13 years from 2011 to 2023 were simulated with historical weather. The ignition locations, times, and fire durations were randomly drawn as a Poisson Point Process [8], following a density function that was influenced by predictor variables such as fuel type, weather, distance to built areas, and a spatial density of ignitions derived from historical data. Compared to previous versions of the FireFactor project, a salient feature of this version is a heavy tail of fire durations, in an attempt to represent the occurrence of very large fires more faithfully. As already noted, we did not use Importance Sampling in these simulations.

This simulation procedure produced a dataset of approximately 31.7 million simulated fires. Because this size of data was impractical for the experiments described below, we thinned this dataset into a random sub-sample of approximately 2.7 million fires: for each fire, we independently drew a random binary variable to determine whether it would be selected in the sub-sample. To improve precision, the larger fires had a higher probability of being selected in the sub-sample4. Based on a rough examination of a fire size histogram and some judgement-based trial and error, the selection probability $p_s$ was chosen to be the following function of the fire size A:

$$p_s\left(A\right):=min\left(1,{\left({\displaystyle \frac{2000\mathrm{ac}+A}{2000\mathrm{ac}+30000\mathrm{ac}}}\right)}^{1.2}\right)$$

(39)

Thus, the selection probability was almost uniform when $A\ll 2000\mathrm{ac}$ (809.3 ha), approximately proportional to $A^{1.2}$ when $A\gg 2000\mathrm{ac}$, and saturating to 1 when $A\ge 30000\mathrm{ac}$ (12140.6 ha).

Denote ${\left(F_i\right)}_i$ the original sample of 31.7 million fires, and ${\left(F_k^{\prime }\right)}_k$ the thinned sub-sample. Then, for any function $u\left(F\right)$, we will be able to approximate the sum of the $u\left(F_i\right)$ as follows:

$$\sum _{i}u\left(F_i\right)\stackrel{\mathbb{E}}{\approx }\sum _k{s}_{k}u\left(F_k^{\prime }\right)$$

(40)

in which the symbol $\stackrel{\mathbb{E}}{\approx }$ denotes equality in expectation, and $s_k$ is the appropriate subsampling weight derived from the selection probability:

$$s_k:={\displaystyle \frac{1}{p_s\left(A_k^{\prime }\right)}}$$

(41)

3.2. Maximum Potential Efficiency Gain

By applying Equation (27) to the dataset described in Section 3.1, we obtained $g^{*}=1.56$ (standard error $\sigma =0.005$). The standard error $\sigma$ was estimated by bootstrap-resampling [13] using 100 bootstrap replications.

The above result is an upper bound on the efficiency gain that can be obtained by Importance Sampling; however, it does not tell us what reweighting function to use in order to approach that upper bound. The next section will explore various reweighting strategies.

3.3. Reweighting Pyromes Uniformly

We start with a reweighting function that is constant within each pyrome. The approach is therefore similar to stratifying (Section 2.15) by pyrome. Pyromes [14] are a spatial partition of the Conterminous United States into 128 areas of relatively homogeneous fire occurrence characteristics. This reweighting strategy means that instead of simulating the same number of fire seasons in all pyromes, we simulate more fire seasons in some pyromes than others; however, the relative frequencies of simulated fires ignited within a given pyrome remain the same.

By applying Equation (33) with variable I defined as the pyrome of ignition, we know the importance weight for a pyrome is proportional to $\sqrt{c\left(I\right)/a\left(I\right)}$, the square root of expected cost divided by expected burned area in that pyrome. Thus, we estimated $a\left(I\right)$ and $c\left(I\right)$ for each pyrome I, yielding the optimal weighting function, up to some proportionality constant. Using Equation (22), we computed the efficiency gain $g_q$, yielding a value of $g_q=1.09$.

We then assessed the geographic inequality in precision implied by this reweighting scheme (Section 2.13). To do so, we computed the same-cost frequency multipliers (Equation (23)), telling us how the number of simulated seasons in each pyrome has increased or decreased. All else being equal, precision is proportional to the number of simulated years, therefore the gain in precision within each pyrome is equal to this multiplier. Table 2 shows the multipliers for each pyrome, along with some other descriptive statistics; these numbers are plotted in Figure 2. We see that most pyromes have a multiplier between 0.7 and 1.3, with a few as low as 0.4 and none higher than 2.

In our opinion from looking at these results, this particular reweighting scheme is more trouble than it is worth: the efficiency gain is rather small, and the reduce precision in some pyromes is concerning.

3.4. Duration-Based Parametric Reweighting in Pyrome 33

As another example, we restrict our attention to only one pyrome and use a parametric weight function (as suggested by Section 2.12) based on simulated fire duration. The assumption here is that the simulated fire duration is randomly drawn at the same time as the ignition location (before the fire is simulated), and that we can oversample or undersample fires based on their duration.

The pyrome under consideration is Pyrome 33 (South Central California Foothills and Coastal Mountains). For this pyrome, we computed (Equation (27)) a maximum potential efficiency gain of $g^{*}=1.38$.

We will use a parametric reweighting function, of the form:

$$ln{w}_1(t;\alpha )=\mathrm{const}+\alpha exp\left(-{\displaystyle \frac{t}{150\mathrm{hr}}}\right)$$

(42)

in which $\alpha$ is a parameter controlling how short fires will get over-sampled ($\alpha <0$) or under-sampled ($\alpha >0$) by Importance Sampling.

To estimate the optimal value ${\alpha }^{*}$, we used the scipy.optimize.minimize function of the SciPy library [15] to solve the following optimization problem, corresponding to Equation (37):

$${\alpha }^{*}=\mathrm{argmin}\left(\alpha \mapsto \left(\sum _k{s}_k{w}_1{(t_k;\alpha )}^{-1}c\left(F_k\right)\right)\left(\sum _k{s}_k{w}_1(t_k;\alpha )a\left(F_k\right)\right)\right)$$

(43)

in which k indexes the simulated fires of the calibration dataset, $t_k$ is the duration of the k-th simulated fire, and $s_k$ is the subsampling weight (see Equation (41)).

The best-fitting parameter value is ${\alpha }^{*}=-1.95$, achieving an efficiency gain of 1.25 (as computed by Equation (22)). Shorter fires get a smaller weight (${\alpha }^{*}<0$) which means that Importance Sampling will improve efficiency by over-representing them. The curve in Figure 3 shows how Importance Sampling alters the sampling frequency of fires under this reweighting function, by displaying the same-cost frequency multipliers (Equation (23)) corresponding to $w(·,{\alpha }^{*})$, along with a point cloud of the frequency multipliers yielded by the best possible reweighting scheme (given by Equation (25)).

A side effect of the reweighted function fitted above is that it increases the sampling frequency of very short fires. There is cause to be unsatisfied with this, in particular because simulating a very large number of tiny fires can make the results unwieldy in ways not captured by our cost function. To address this limitation, we now fit a more refined reweighting function, expressed as:

$$ln{w}_2(t;{\alpha }_1,{\tau }_1,{\alpha }_2,{\tau }_2)=\mathrm{const}+{\alpha }_{1}exp\left(-{\displaystyle \frac{t}{{\tau }_1}}\right)+{\alpha }_{2}exp\left(-{\displaystyle \frac{t}{{\tau }_2}}\right)$$

(44)

The above function is more flexible, having 4 parameters instead of 1; in particular, the duration lengthscales are no longer fixed, thanks to the ${\tau }_i$ parameters. We fitted this function with initial parameters $({\alpha }_1=0,{\tau }_1=5\mathrm{hr},{\alpha }_2=0,{\tau }_1=150\mathrm{hr})$: in practice, this means that $({\alpha }_1,{\tau }_1)$ influence short fires while $({\alpha }_2,{\tau }_2)$ influence long fires. The best-fit values were $({\alpha }_1^{*}=2.53,{\tau }_1^{*}=3.63\mathrm{hr},{\alpha }_2^{*}=-2.15,{\tau }_2^{*}=81.2\mathrm{hr})$, achieving an efficiency gain of 1.27 (confirmed with K-fold cross-validation). Note that ${\alpha }_1^{*}>-{\alpha }_2>0$, thus reducing the sampling frequency of very short fires. Figure 4 displays the fitted frequency multipliers: in summary, Importance Sampling in this case will increase the sampling frequency of medium-duration fires, at the expense of low-duration and high-duration fires.

3.5. Fitting a Machine-Learning Model to Empirical Optimal Weights

We now illustrate the last suggestion of Section 2.12: fitting a supervised learning to the optimal weights derived from the empirical sample. Using the dataset of 2.7 million simulated fires described in Section 3.1, we trained a gradient-boosted regression tree algorithm [9] to predict the optimal weight $w^{*}$ (computed from the empirical sample using Equation (25)). More precisely, we used the HistGradientBoostingRegressor class of the Scikit-Learn library [16], which is inspired by LightGBM [17]. We used $ln{w}^{*}$ as the target variable, and we weighted samples by fire size. As predictor variables, we used the fire duration, pyrome number, ignition-point fuel model number, ignition-time ERC percentile, some fire behavior metrics (rate of spread, reaction intensity in ${\mathrm{kW}/\mathrm{m}}^2$, fire-line intensity in kW/m) derived from the fuel model in reference weather conditions5 using the Rothermel model [18], and spatial coordinates of ignition point $(x,y)$.

We fitted two models, differing only in which subset of predictor variables was used:

• Model A: used all predictor variables listed above.
• Model B: used all predictor variables except for the spatial coordinates $(x,y)$.

The models were tested by a method similar to K-fold cross-validation: training on a random subset of 70% of the fires, computing the efficiency gain on the remaining 30%, repeating this procedure 30 times and computing the mean and standard deviation of the efficiency gain. The test performance of the model is summarized by Table 3. We see that Model A achieves a much higher efficiency gain than Model B, revealing that it learns fine-grained spatial features from the spatial coordinates of ignition. While its superior efficiency gain is appealing, we are wary of using Model A, because of the subtle potential side effects of having estimator precision that varies at fine spatial scales.

4. Discussion

Perhaps the most questionable aspect of these methods is the use of $L_2$ as the metric of deviation. While this choice is mathematically convenient, and has various theoretical justifications, it may create geographic disparities in precision, and excessively favor absolute precision over relative precision (as noted in Section 2.13).

When the computational cost of a simulated fire is asymptotically a super-linear function of its size, then Importance Sampling will reduce the sampling frequency of large fires (Section 2.7). This is unintuitive, and can be problematic: the largest fires are typically those we are most interested in. Reducing the sampling frequency of large fires may improve precision in the sense of $L_2$, but it might at the same time conflict with other precision objectives. This issue can be mitigated in several ways:

• by weighting the $L_2$ metric by value (Section 2.9), e.g. by giving more importance of Highly Valued Resources and Assets;
• by including non-geographic outcomes like burn severity (Section 2.10) and weighting them unevenly;
• by constraining the family of reweighting functions under consideration to be uniform in the tail of large fires.

Consequently, the main recommendation of this article is to consider using some form of importance sampling, not necessarily the one which is optimal in the sense of $L_2$ convergence. We find the approach most promising for reducing the sampling frequency of many small inconsequential fires in situations where the simulation pipeline makes them disproportionately costly.

To make such approaches practical, we recommend that wildland fire simulation software enable users to provide weights alongside ignitions, these weights being taken into account for aggregating simulation results into outputs. It should be immaterial whether such weights originate from importance sampling or similar techniques like stratified sampling.

5. Conclusions

This article has studied the application of Importance Sampling to Monte Carlo wildfire simulations. This technique can accelerate the convergence of Monte Carlo simulations to their asymptotic burn-probability map; Importance Sampling does not introduce bias, and therefore cannot improve or degrade the realism of that asymptotic limit. Simulating an infinite number of fires would yield the same burn-probability map with or without Importance Sampling.

Our analysis starts from the classical generic result of the Monte Carlo literature about importance sampling for a scalar estimator, which shows that the variance-minimizing proposal distribution reweights the probability of each sample by its absolute value (see Equation (18)). We have made this result applicable to the estimation of burn-probability maps through a series of generalizations:

• estimating a map rather than a scalar (Section 2.2), using the $L_2$ distance as the metric of deviation underlying variance (Section 2.1);
• accounting for the variability of computational costs across simulated fires, adopting cost-to-precision as the performance metric (Equation (7)) instead of simply variance;
• deriving the optimal cost-to-precision and the corresponding proposal distribution (Section 2.4), relying on the Poisson process-approximation (Section 2.14) that each fire affects only a small fraction of the area of interest;
• generalizing to proposal distributions constrained to an upstream variable (typically the ignition), see Section 2.8;
• generalizing to impact-weighted $L_2$ distance (Section 2.9);
• generalizing from burn-probability estimation to other types of maps (Section 2.10).

We illustrated these methods in Section 3 with various case studies. Here is a summary of the key findings:

• Section 2.1 observes that the relevant metric for quantifying the efficiency of Monte Carlo simulation is the cost-to-precision ratio $r_{C/P}$, that is the product of variance to expected computational cost (Equation (7)).
• Without Importance Sampling, the variance of the burn-probability estimator is simply the expected burned area divided by the number of simulated fires (Equations (13) and (6)).
• The efficiency characteristics of the Monte Carlo sampling, and in particular the proposal distribution which maximizes convergence speed, are fully determined by the joint distribution of fire size and computational cost (Equation (2.4)). Beyond these two variables, the diversity of shapes and behaviors of simulated fires is irrelevant.
• The optimal proposal distribution (lowest cost-to-precision ratio) reweights the natural probability of each fire by a factor proportional to the square root of the ratio of burned area to computational cost (Equation (25)). When only the distribution of ignitions can be reweighted, the burned area and computational cost must be replaced by their expectation conditional on the ignition (Equation (33)).
• In particular, when computational cost is strictly proportional to burned area, importance sampling can achieve no progress: the natural distribution is already the optimal proposal distribution. This linearity assumption is most dubious at the “small fires” tail of the natural distribution.
• Finding a good proposal distribution is not trivial, as it requires predictive power. Section 2.12 suggests a machine learning approach based on calibration runs.
• However, the best achievable cost-to-precision ratio is easily estimated based on a calibration sample of fire sizes and computational costs, using Equation (26). This allows for quickly assessing whether Importance Sampling is worth pursuing.
• Stratified or neglected ignitions can be seen as approximations to special cases of Importance Sampling, in which the reweighting is piecewise constant (Section 2.15).

We do not necessarily recommend applying the methods and formula of this paper to the letter. In our opinion, the most useful contributions of this work are to provide a relevant framing of the computational efficiency problem, to raise awareness of Importance Sampling in general to the fire modeling community, and to showcase some calculation techniques associated with it. The particular flavor of Importance Sampling studied here, which uses $L_2$ distance as the metric of convergence, is more questionable and should be assessed on a case-by-case basis (see Section 4). Ultimately, the optimal proposal distribution depends on intended downstream use of simulations results.

Appendix A.1. Probability Basics

Notation: we write $\mathbb{E}\left[Y\right]$ for the expected value of random variable Y, and $\mathrm{Var}\left[Y\right]$ for its variance. $Y\sim Y^{\prime }$ means that Y and $Y^{\prime }$ follow the same probability distribution; if $p\left(y\right)$ is a probability density function, $Y\sim p$ means that Y follows the corresponding probability distribution. $Y\stackrel{\mathbb{E}}{\approx }{Y}^{\prime }$ means that Y and $Y^{\prime }$ have the same expected value, i.e. $\mathbb{E}\left[Y\right]=\mathbb{E}{\left[Y\right]}^{\prime }$.

If Y is a random variable with probability density $p\left(y\right)$, and $u\left(y\right)$ is a function, then the expected value $\mathbb{E}\left[u\right(Y\left)\right]$ can be written as an integral:

$$\mathbb{E}\left[u\left(Y\right)\right]={\int }_{y}dyp\left(y\right)u\left(y\right)$$

(A1)

Appendix A.2. The Probabilistic Cauchy-Schwarz Inequality

The probabilistic Cauchy-Schwarz inequality is the following theorem. For any two real-valued random variables U and V, we have:

$$\mathbb{E}\left[UV\right]\le \mathbb{E}{\left[U^2\right]}^{{\displaystyle \frac{1}{2}}}\mathbb{E}{\left[V^2\right]}^{{\displaystyle \frac{1}{2}}}$$

(A2)

This inequality is an equality if and only if U and V are proportional to one another with a non-negative coefficient, i.e. if $U=\lambda V$ or $V=\lambda U$ for some $\lambda \ge 0$ (more rigorously, if there exists $\lambda \ge 0$ such that either $\mathbb{P}[U=\lambda V]=1$ or $\mathbb{P}[V=\lambda U]=1$).