ABSTRACT

We study the random injury outcome caused by multiple flash bang submunitions on a crowd. We are particularly interested in the fluctuations in injury outcome among individual realizations. Previously, to simulate the distribution of the actual number of injured, we developed a comprehensive Monte Carlo model. While the full computational model is important for thorough theoretical investigations, in practical operations, it is desirable to characterize the phenomenological behavior of injury outcome using a concise model with only one or two parameters. Conventionally, the injury outcome is indicated by the average fraction of injured, which is called the risk of significant injury (RSI). The single metric RSI description fails to capture fluctuations in the injury outcome. The number of injured in the crowd is influenced by many random factors: the aiming error of flash bang mortar, the dispersion of submunitions after mortar burst, the amount of acoustic dose reaching individual subjects, and the biovariability of individual subjects’ reactions to a given acoustic dose. We aim to include these random factors properly and concisely. In this study, we represent the random injury outcome as a compound binomial model, in which the hidden injury probability is drawn from a two-parameter model distribution. We formulate and examine six model distributions for the injury probability. The best performer is a mixture of uniform and triangle distributions, parameterized by (RSI, dp) where dp is the standard deviation of the hidden injury probability. This mixture model predicts the behavior of injury outcome with uncertainty, based solely on the two parameters (RSI, dp) in the flash bang description. For example, we can predict the probability of the injury outcome not exceeding a prescribed tolerance. We advocate the adoption of this two-parameter characterization for flash bangs to replace the single-parameter RSI description. Whenever we need to give a high level coarse description of a flash bang situation, we state that the injury risk is represented by (RSI, dp).

Keywords:

Flash Bangs, Actual Number of Injured in a Crowd, Fluctuations in the Hidden Injury Probability, Simplified Phenomenological Models

1. Background

Flash bangs have been widely used for crowd control, riot control and area denial [1] [2] [3] [4] . In order to use flash bangs effectively, it is important to assess the associated injury risk. For a crowd of n_c subjects uniformly distributed in a circle of diameter d_c, we study the injury outcome on the crowd caused by a flash bang mortar with multiple submunitions. Previously, we constructed a comprehensive Monte Carlo model for the crowd injury outcome, simulating the detailed effects of many model parameters [5] . In the current study, we consider the problem of how to summarize a flash bang situation using a properly designed two-parameter model so that the distribution of the random injury outcome can be accurately reconstructed solely based on the two parameters supplied. We first review briefly the problem setup in the full computational model.

To establish the coordinate system, we put the origin at the center of the crowd circle. We select the range direction of a flash bang mortar as the x-axis and the deflection direction as the y-axis. In this coordinate system, the crowd is uniformly distributed in circle $\text{Cir}\left(\mathrm{0,}{d}_c/2\right)$. The intended burst position of flash bang mortar is in the air directly above the crowd center. The actually realized burst position is random, affected by the aiming error. There is more randomness in the problem: the dispersion of flash bang submunitions upon the mortar burst is random [5] . In addition, even at a fixed acoustic dose (SELA) reaching the surroundings of a subject, the injury outcome is still random, affected by the uncertainty of dose propagating into the injury site [6] , and affected by the biovariability of subjects’ individual responses to a given dose at the injury site [7] . Because of these uncertainties, the actual number of injured N out of n_c subjects will be a random variable fluctuating in individual realizations. Previously we developed a comprehensive computational framework simulating random realizations of crowd injury outcome in a situation specified by the crowd and flash bang parameters [5] . Some of the parameters are summarized below:

· $d_c=20\text{m}$, diameter of the crowd circle;

· $n_c=100$, number of subjects in the crowd;

· $N=\text{randomvariable}$, actual number of injured;

· $n_s=20$, number of submunitions in a flash bang mortar;

· $\sigma =\text{variable}\in \left(0,30\text{m}\right)$, standard deviation of aiming error;

· $d_{\text{rng}}=25\text{m}$, dimension in the range direction, of the submunitions dispersion area;

· $d_{\text{defl}}=15\text{m}$, dimension in the deflection direction, of the submunitions dispersion area;

· ${F|}_{r=1\text{m}}=150\text{dB}$, acoustic intensity of a single flash bang submunition at 1 meter away;

· ${\text{ID}}_{50}=163\text{dB}$, median acoustic dose for a hearing loss of 30 dB or more in noise-induced permanent threshold shift (PTS);

· $c_{\alpha }=0.1{\left(\text{dB}\right)}^{-1}$, steepness coefficient in the dose-injury function for a hearing loss of 30 dB or more in noise-induced permanent threshold shift (PTS).

Mathematically, the distribution of random variable N is fully described by the Monte Carlo model. For practical operations, it is convenient to have a high-level summary of the flash bang situation using only one or two parameters. Conventionally, the overall injury effect of a flash bang mortar on a crowd is indicated by the average fraction of injured, which is called the risk of significant injury (RSI). The average of N is then given by $n_c\times \text{RSI}$. Unfortunately, the single metric RSI value does not provide information about fluctuations in N among individual realizations. In particular, the RSI description does not reflect the fluctuations of aiming error, which swings the injury probability up and down simultaneously for all individual subjects in a crowd.

In this study, we seek to construct a two-parameter probabilistic description of random crowd injury outcome caused by multiple flash bang submunitions. In this high level coarse description, each flash bang situation is summarized in a two-parameter characterization, and the distribution of N is predicted solely from the two parameters in the given characterization, without resorting to the full Monte Carlo simulations on all model parameters. With the predicted distribution of N, we can make meaningful predictions regarding the random crowd injury outcome. For example, we can estimate the probability of having more than 10 injured in a crowd of 100 subjects, $\mathrm{Pr}\left({N>10|}_{n_c=100}\right)$. These probabilities are very informative and are sought after when making real operation decisions. We point out that probabilities of this type cannot be estimated solely from the single metric RSI value. For example, an RSI value of 6% may yield $\mathrm{Pr}\left({N>10|}_{n_c=100}\right)=0$ when N is narrowly distributed within interval [3, 9]; the same 6% RSI value may lead to $\mathrm{Pr}\left({N>10|}_{n_c=100}\right)>50\%$ if N alternates between 0 and 11. To achieve a reasonably accurate estimate of these probabilities in a concise model, we consider two-parameter descriptions of N.

2. Mathematical Formulation

We start by introducing relevant variables and functions.

· $\omega$ = a random realization of submunitions distribution relative to the crowd center. Here symbol $\omega$ represents both the randomness in the burst position of a mortar round (aiming error) and the randomness in submunitions dispersion after mortar burst.

· N = the actual number of injured out of n_c subjects. N is a random variable, with outcome influenced by three factors:

i) distribution of submunitions relative to the crowd center ( $\omega$ ),

ii) distribution of n_c subjects in the crowd, and

iii) fluctuations in injury outcome for a subject exposed to a given acoustic dose.

Note that the dose-injury model predicts only the injury probability corresponding to a given dose. It does not completely determine the binary injury outcome. As a result, even when the distribution of submunitions ( $\omega$ ) is given and fixed, the conditional random variable $\left(N\mathrm{<mo>\; |\; </mo>}\omega \right)$ is still random and far from being deterministic.

· $P\left(x\mathrm{,}y\mathrm{,}\omega \right)$ = injury probability of a subject at location $\left(x\mathrm{,}y\right)$. $P\left(x\mathrm{,}y\mathrm{,}\omega \right)$ is a random variable, varying with submunitions distribution $\omega$, and with subject location $\left(x\mathrm{,}y\right)$.

· $P\left(\omega \right)$ = injury probability of a random subject, averaged over subject location. $P\left(\omega \right)$ is a random variable, varying with submunitions distribution $\omega$. Mathematically, $P\left(\omega \right)$ has the expression:

$P\left(\omega \right)\equiv E_{\left(x\mathrm{,}y\right)}\left[P\left(x\mathrm{,}y\mathrm{,}\omega \right)\right]=\frac{1}{\text{π}{\left(d_c/2\right)}^2}{\displaystyle {\int }_{\text{Cir}\left(\mathrm{0,}{d}_c/2\right)}}P\left(x\mathrm{,}y\mathrm{,}\omega \right)\text{d}x\text{d}y$

· $I^{\left(k\right)}$ = indicator function that subject k is injured. $I^{\left(k\right)}$ is the random injury outcome of subject k, which is affected by factors i), ii) and iii) listed above.

We model the locations of $n_c$ subjects as independently and identically distributed (i.i.d.) with uniform distribution over circle $\text{Cir}\left(\mathrm{0,}{d}_c/2\right)$. Furthermore, we assume the crowd subjects distribution is independent of the submunitions distribution ( $\omega$ ). It follows that, at any given submunitions distribution ( $\omega$ ), the conditional random injury outcomes of individual subjects satisfy several properties.

Property 1:

Given $\omega$, the conditional individual injury outcomes of $n_c$ subjects, $\left\{\left(I^{\left(k\right)}|\omega \right),k=1,2,\cdots ,n_c\right\}$, are independently and identically distributed (i.i.d.).

Property 2:

$\left(I^{\left(k\right)}|\omega \right)$, the conditional injury outcome for subject k given submunitions distribution $\omega$, is a Bernoulli random variable with probability $P\left(\omega \right)$. Since the n_c subjects are i.i.d., the distribution of $\left(I^{\left(k\right)}|\omega \right)$ is independent of k.

Property 3:

$\left(N|\omega \right)={\displaystyle {\sum }_{k=1}^{n_c}}\left(I^{\left(k\right)}|\omega \right)$, the conditional crowd injury outcome for n_c subjects given $\omega$, is a binomial random variable with parameter n_c and $P\left(\omega \right)$.

$\left(N\mathrm{<mo>\; |\; </mo>}\omega \right)\sim {\text{Bin}}_{(n_c\mathrm{,}P(\omega ))}$

Note that these statements are valid conditional on a given submunitions distribution $\omega$. Over random $\omega$, injury outcomes of the n_c subjects $\left\{I^{\left(k\right)}\left(\omega \right),k=1,2,\cdots ,n_c\right\}$ are no longer independent of each other, and probability $P\left(\omega \right)$ itself is a random variable. For example, when the position of mortar burst is completely off the crowd area due to a large aiming error, the binary outcomes for the n_c subjects are all very likely to be “not injured”.

We introduce quantities $p_a$ and $\delta p$, the two parameters that we will use to specify phenomenological models for the crowd injury outcome.

$p_a\equiv E\left[P\right]\equiv E_{\omega }\left[P\left(\omega \right)\right]$

$\delta p\equiv \text{std}\left[P\right]\equiv {\text{std}}_{\omega }\left[P\left(\omega \right)\right]$

where $E_{\omega }\left[\cdot \right]$ denotes the average over the distribution of $\omega$. Parameter $p_a$ is the same as RSI (the average fraction of injured), and $\delta p$ measures fluctuations in the hidden injury probability $P\left(\omega \right)$, caused by randomness in submunitions distribution $\omega$. We write the second moment of $P\left(\omega \right)$ in terms of $p_a$ and $\delta p$ :

$E\left[P^2\right]={\left(\delta p\right)}^2+p_a^2$

Note that the injury probability $P\left(\omega \right)$ is a hidden variable, not directly measurable. As a result, parameters $p_a$ and $\delta p$ cannot be directly estimated as the mean and standard deviation of $P\left(\omega \right)$. On the other hand, the crowd injury outcome N is directly observable. We like to estimate $p_a$ and $\delta p$ based on observed samples of N.

We first calculate the mean and variance of injury outcome $I^{\left(k\right)}$ of subject k, and then we calculate those of the crowd injury outcome $N={\displaystyle {\sum }_k}{I}^{\left(k\right)}$.

$E\left[I^{\left(k\right)}\mathrm{<mo>\; |\; </mo>}\omega \right]=P(\omega )$

$\text{var}\left[I^{\left(k\right)}\mathrm{<mo>\; |\; </mo>}\omega \right]=P\left(\omega \right)(1-P(\omega ))$

$E\left[N\mathrm{<mo>\; |\; </mo>}\omega \right]=n_{c}E\left[I^{\left(k\right)}\mathrm{<mo>\; |\; </mo>}\omega \right]=n_{c}P(\omega )$

$\text{var}\left[N\mathrm{<mo>\; |\; </mo>}\omega \right]=n_c\text{var}\left[I^{\left(k\right)}\mathrm{<mo>\; |\; </mo>}\omega \right]=n_{c}P\left(\omega \right)(1-P(\omega ))$

$E\left[N\right]=E_{\omega }\left[E\left[N\mathrm{<mo>\; |\; </mo>}\omega \right]\right]=n_c{E}_{\omega }\left[P\left(\omega \right)\right]=n_c{p}_a$ (1)

$\begin{array}{c}E\left[{\left(N\mathrm{<mo>\; |\; </mo>}\omega \right)}^2\right]=\text{var}\left[N\mathrm{<mo>\; |\; </mo>}\omega \right]+{\left(E\left[N\mathrm{<mo>\; |\; </mo>}\omega \right]\right)}^2\\ =n_{c}P\left(\omega \right)\left(1-P\left(\omega \right)\right)+n_c^{2}P{\left(\omega \right)}^2\\ =n_{c}P\left(\omega \right)+n_c\left(n_c-1\right)P{\left(\omega \right)}^2\end{array}$

$\begin{array}{c}E\left[N^2\right]=E_{\omega }\left[E\left[{\left(N\mathrm{<mo>\; |\; </mo>}\omega \right)}^2\right]\right]\\ =n_c{E}_{\omega }\left[P\left(\omega \right)\right]+n_c\left(n_c-1\right)E_{\omega }\left[P{\left(\omega \right)}^2\right]\\ =n_c{p}_a+n_c\left(n_c-1\right)\left({\left(\delta p\right)}^2+p_a^2\right)\end{array}$

$\text{var}\left[N\right]=E\left[N^2\right]-{\left(E\left[N\right]\right)}^2=n_c{p}_a\left(1-p_a\right)+n_c\left(n_c-1\right){\left(\delta p\right)}^2$ (2)

Since N is directly observable in field operations and in Monte Carlo simulations, quantities $E\left[N\right]$ and $\text{var}\left[N\right]$ can be estimated from measured samples of N. The results derived above offer a way of expressing $p_a\equiv E\left[P\right]$ and $\delta p\equiv \text{std}\left[P\right]$ in terms of the measurable quantities $E\left[N\right]$ and $\text{var}\left[N\right]$. Using results (1) and (2), we write $p_a$ and $\delta p$ as

$p_a=\frac{E\left[N\right]}{n_c}$ (3)

${\left(\delta p\right)}^2=\frac{1}{n_c-1}\left[\frac{\mathrm{var}\left[N\right]}{n_c}-\frac{E\left[N\right]}{n_c}\left(1-\frac{E\left[N\right]}{n_c}\right)\right]$ (4)

Note that the injury probability $P\left(\omega \right)$ for a random subject is not affected by n_c, the number of subjects in the crowd. Mathematically, both $p_a\equiv E\left[P\right]$ and $\delta p\equiv \text{std}\left[P\right]$ are independent of n_c. In Monte Carlo simulations and in field operations, however, all observable quantities vary with n_c. In expressions (3) and (4), the n_c-independent quantities $p_a$ and $\delta p$ are estimated based on parameter n_c and n_c-dependent observable quantities $E\left[N\right]$ and $\text{var}\left[N\right]$. The key significance of (3) and (4) is that they express non-observable quantities $p_a$ and $\delta p$ in terms of observable quantities, and thus, make $p_a$ and $\delta p$ effectively measurable. The estimation of $p_a$ and $\delta p$ opens the possibility of parametrically representing the non-observable distribution of injury probability $P\left(\omega \right)$ using measurable quantities. If the distribution of random variable $P\left(\omega \right)$ is adequately captured in the parametric representation, it will produce a phenomenological model for the distribution of N, which, in turn, provides a practical method for making meaningful predictions regarding the random crowd injury outcome N. For example, we can predict the probability of the crowd injury outcome not exceeding a prescribed tolerance.

Our goal in this study is to construct a probabilistic description of N based solely on the measurable parameters $p_a$ and $\delta p$. The key step in the construction is to formulate a two-parameter model distribution for the non-observable injury probability $P\left(\omega \right)$. In next section, we explore six models for $P\left(\omega \right)$.

3. Six Models for the Distribution of $P(\omega )$

We consider six models for representing the distribution of non-observable injury probability $P\left(\omega \right)$. The parametric representation in each model is based on measurable quantities $p_a\equiv E\left[P\right]$ and $\delta p\equiv \text{std}\left[P\right]$. Two factors motivate us to adopt $p_a$ and $\delta p$ as the key quantities for characterizing the non-observable $P\left(\omega \right)$.

· $p_a$ and $\delta p$ are effectively measurable, which is important for any practical model.

· $p_a$ and $\delta p$ capture the holistic effects of multiple flash bang submunitions in causing injury on a crowd; $p_a$ is the RSI (average fraction of injured); the addition of $\delta p$ in the formulation incorporates the effects of aiming error and submunitions distribution uncertainty.

Here, we emphasize that $P\left(\omega \right)$ is a random variable, varying with submunitions distribution $\omega$. Our objective in this section is to phenomenologically represent the distribution of $P\left(\omega \right)$ as accurately as possible, using only parameters $p_a$ and $\delta p$. Below we discuss and examine six models, one in each subsection.

3.1. Model 1: Delta Function Distribution for $P(\omega )$

Injury probability $P\left(\omega \right)$ is a random variable, influenced by submunitions distribution $\omega$. Let ${\rho }_P\left(p\right)$ be the probability density function (PDF) of random variable $P\left(\omega \right)$. We model ${\rho }_P\left(p\right)$ as a delta function, namely,

${\rho }_P\left(p\right)=\delta \left(p-p_a\right)$ (5)

The corresponding model for N, the actual number of injured out of n_c subjects, is a binomial distribution of parameters n_c and p_a. The probability mass function (PMF) of N has the expression,

$P\left(\omega \right)\sim \delta \left(p-p_a\right)\mathrm{,}N\sim {\text{Bin}}_{(n_c\mathrm{,}P(\omega ))}$

${\text{Pr}}_{\left[N=k\right]}={\text{Bin}}_{\left(n_c\mathrm{,}{p}_a\right)}\left(k\right)\equiv C\left(n_c\mathrm{,}k\right)p_a^k{\left(1-p_a\right)}^{\left(n_c-k\right)}$ (6)

When the distribution of $P\left(\omega \right)$ is not concentrated in a narrow region near p_a, the delta function representation of $P\left(\omega \right)$ is not expected to yield a good approximation distribution for N. The advantage of the delta function model is its simplicity: it uses only one parameter p_a; it does not require $\delta p$.

3.2. Model 2: Beta Distribution for $P(\omega )$

We model ${\rho }_P\left(p\right)$ as a beta distribution with shape parameters $\alpha$ and $\beta$.

${\rho }_P\left(p\right)={\text{Beta}}_{\left(\alpha \mathrm{,}\beta \right)}\left(p\right)\equiv \frac{p^{\alpha -1}{\left(1-p\right)}^{\beta -1}}{B\left(\alpha \mathrm{,}\beta \right)}$ (7)

where $B\left(u\mathrm{,}v\right)$ is the beta function defined as

$B\left(u\mathrm{,}v\right)\equiv {\displaystyle {\int }_0^1}{t}^{u-1}{\left(1-t\right)}^{v-1}\text{d}t$ (8)

$B\left(u\mathrm{,}v\right)$ is related to the gamma function, $\Gamma \left(z\right)\equiv {\displaystyle {\int }_0^{\infty }}{s}^{z-1}{\text{e}}^{-s}\text{d}s$, by

$B\left(u\mathrm{,}v\right)=\frac{\Gamma \left(u\right)\Gamma \left(v\right)}{\Gamma \left(u+v\right)}$ (9)

We determine parameters $\alpha$ and $\beta$ in model (7) by matching the mean and variance of the beta distribution to those of $P\left(\omega \right)$.

$\frac{\alpha }{\alpha +\beta }=p_a$ (10)

$\frac{\alpha \beta }{{\left(\alpha +\beta \right)}^2\left(\alpha +\beta +1\right)}={\left(\delta p\right)}^2$ (11)

Solving for $\alpha$ and $\beta$ from (10) and (11), we obtain

$\alpha =p_a\left(\frac{p_a\left(1-p_a\right)}{{\left(\delta p\right)}^2}-1\right)$ (12)

$\beta =\left(1-p_a\right)\left(\frac{p_a\left(1-p_a\right)}{{\left(\delta p\right)}^2}-1\right)$ (13)

The actual number of injured out of n_c subjects, N, has a beta-binomial distribution, which is the result of compounding the binomial distribution with a beta distribution for the probability parameter.

$P\left(\omega \right)\sim {\text{Beta}}_{\left(\alpha \mathrm{,}\beta \right)}\mathrm{,}N\sim {\text{Bin}}_{(n_c\mathrm{,}P(\omega ))}$

${\text{Pr}}_{\left[N=k\right]}=C\left(n_c\mathrm{,}k\right)\frac{B\left(k+\alpha \mathrm{,}{n}_c-k+\beta \right)}{B\left(\alpha \mathrm{,}\beta \right)}$ (14)

The probability mass function (PMF) of N is directly described by $\left(n_c\mathrm{,}\alpha \mathrm{,}\beta \right)$. Via Equations (12) and (13), ${\text{Pr}}_{\left[N=k\right]}$ is indirectly specified by the two flash bang parameters $p_a$ and $\delta p$, along with the number of subjects n_c.

3.3. Model 3: Mixture of Uniform Distributions for $P\left(\omega \right)$, Version A

Note that ${\rho }_P\left(p\right)$, the density of random variable $P\left(\omega \right)$, has a domain restricted to interval $\left[\mathrm{0,1}\right]$ since $P\left(\omega \right)$ has the meaning of injury probability and thus is constrained to interval $\left[\mathrm{0,1}\right]$. We model ${\rho }_P\left(p\right)$ as a uniform distribution and adjust it, if necessary, to make the support of model distribution fall inside interval $\left[\mathrm{0,1}\right]$. We start with a uniform distribution over $\left(a\mathrm{,}b\right)$. Parameters a and b are determined by matching the mean and variance of the uniform distribution to those of $P\left(\omega \right)$ :

$\frac{1}{2}\left(a+b\right)=p_a$ (15)

$\frac{1}{12}{\left(b-a\right)}^2={\left(\delta p\right)}^2$ (16)

Solving for a and b from (15) and (16), we obtain

$a=p_a-\sqrt{3}\left(\delta p\right)$ (17)

$b=p_a+\sqrt{3}\left(\delta p\right)$ (18)

The values of $p_a$ and $\delta p$ for flash bangs are typically well below 20%. As a result, we don’t need to worry about the situation of $b>1$. The situation of $a<0$, however, occurs when $\sqrt{3}\left(\delta p\right)>p_a$. For flash bang mortars with random aiming errors of significant magnitude, $\sqrt{3}\left(\delta p\right)>p_a$ is actually very likely. The aiming errors increase $\delta p$ and reduce $p_a$ at the same time, resulting in $\sqrt{3}\left(\delta p\right)>p_a$, which in turn leads to $a<0$. With a negative value of a, the proposed uniform distribution is an invalid representation for the density of $P\left(\omega \right)$. When we have $\sqrt{3}\left(\delta p\right)>p_a$, we need to change the model distribution for $P\left(\omega \right)$. We discuss the cases of $\sqrt{3}\left(\delta p\right)\le p_a$ and $\sqrt{3}\left(\delta p\right)>p_a$ separately.

Case 1: $\sqrt{3}\left(\delta p\right)\le p_a$

For $\sqrt{3}\left(\delta p\right)\le p_a$, we have $a\ge 0$, and the uniform distribution over $\left(a\mathrm{,}b\right)$, denoted by ${\text{unif}}_{\left(a\mathrm{,}b\right)}\left(p\right)$, is a valid model for the density of $P\left(\omega \right)$. In this case, we use ${\text{unif}}_{\left(a\mathrm{,}b\right)}\left(p\right)$ to represent ${\rho }_P\left(p\right)$, without modification.

${\rho }_P\left(p\right)={\text{unif}}_{\left(a\mathrm{,}b\right)}\left(p\right)\equiv \frac{1}{b-a}\{\begin{array}{ll}\mathrm{1,}\hfill & \text{for}p\in \left[a\mathrm{,}b\right]\hfill \\ \mathrm{0,}\hfill & \text{otherwise}\hfill \end{array}$ (19)

For the crowd injury outcome, N, the distribution corresponding to model (19) is a uniform-binomial distribution. The compound distribution draws $P\left(\omega \right)$ from uniform distribution (19) and then generates random variable N using the binomial distribution with parameters $\left(n_c\mathrm{,}P\left(\omega \right)\right)$. The PMF of the compound distribution is expressed as an integral,

$P\left(\omega \right)\sim {\text{unif}}_{\left(a\mathrm{,}b\right)}\mathrm{,}N\sim {\text{Bin}}_{(n_c\mathrm{,}P(\omega ))}$

$\begin{array}{c}{\text{Pr}}_{\left[N=k\right]}=C\left(n_c\mathrm{,}k\right)\frac{1}{b-a}{\displaystyle {\int }_a^b}{p}^k{\left(1-p\right)}^{n_c-k}\text{d}p\\ =C\left(n_c\mathrm{,}k\right)\frac{B\left(b\mathrm{<mo>\; ;\; </mo>}k+\mathrm{1,}{n}_c-k+1\right)-B\left(a\mathrm{<mo>\; ;\; </mo>}k+\mathrm{1,}{n}_c-k+1\right)}{b-a}\end{array}$ (20)

where $B\left(s\mathrm{<mo>\; ;\; </mo>}u\mathrm{,}v\right)$ is the incomplete beta function defined as

$B\left(s\mathrm{<mo>\; ;\; </mo>}u\mathrm{,}v\right)\equiv {\displaystyle {\int }_0^s}{t}^{u-1}{\left(1-t\right)}^{v-1}\text{d}t$ (21)

For moderately large $\left(u+v\right)$, function value of $B\left(s\mathrm{<mo>\; ;\; </mo>}u\mathrm{,}v\right)$ is extremely small. Numerically, it is inconvenient to work with $B\left(s\mathrm{<mo>\; ;\; </mo>}u\mathrm{,}v\right)$ directly in finite precision arithmetics. Most software packages (such as Matlab) offer a built-in function for calculating the regularized incomplete beta function, $I_s\left(u\mathrm{,}v\right)$, which is related to the incomplete beta function as

$I_{\left(s\right)}\left(u\mathrm{,}v\right)\equiv \frac{B\left(s\mathrm{<mo>\; ;\; </mo>}u\mathrm{,}v\right)}{B\left(u\mathrm{,}v\right)}=\frac{{\displaystyle {\int }_0^s}{t}^{u-1}{\left(1-t\right)}^{v-1}\text{d}t}{{\displaystyle {\int }_0^1}{t}^{u-1}{\left(1-t\right)}^{v-1}\text{d}t}$ (22)

The regularized incomplete beta function has the nice property that as input s increases, the function value increases monotonically from ${I_{\left(s\right)}\left(u,v\right)|}_{s=0}=0$ to ${I_{\left(s\right)}\left(u,v\right)|}_{s=1}=1$. It can be conveniently evaluated by calling built-in function “betainc” in Matlab. In (20), we use $B\left(s\mathrm{<mo>\; ;\; </mo>}u\mathrm{,}v\right)=B\left(u\mathrm{,}v\right)I_{\left(s\right)}\left(u\mathrm{,}v\right)$ and the identity

$\begin{array}{l}C\left(n_c\mathrm{,}k\right)\cdot B\left(k+\mathrm{1,}{n}_c-k+1\right)\\ =\frac{\Gamma \left(n_c+1\right)}{\Gamma \left(k+1\right)\Gamma \left(n_c-k+1\right)}\cdot \frac{\Gamma \left(k+1\right)\Gamma \left(n_c-k+1\right)}{\Gamma \left(n_c+2\right)}=\frac{1}{n_c+1}\end{array}$ (23)

With these identities, we write the PMF of N as

${\mathrm{Pr}}_{\left[N=k\right]}=\frac{I_{\left(b\right)}\left(k+1,n_c-k+1\right)-I_{\left(a\right)}\left(k+1,n_c-k+1\right)}{\left(n_c+1\right)\left(b-a\right)}$ (24)

This expression of PMF for uniform-binomial distribution serves as a building block in the discussion of case 2 below.

Case 2: $\sqrt{3}\left(\delta p\right)>p_a$

For $\sqrt{3}\left(\delta p\right)>p_a$, we have $a<0$, and the uniform distribution over $\left(a\mathrm{,}b\right)$ is not a valid model for the density of $P\left(\omega \right)$. To match the mean and variance of $P\left(\omega \right)$ with a valid model, we need to broaden the form of candidate distributions. We use a mixture of two uniform distributions, which generally has 5 degrees of freedom: two parameters for each uniform distribution mode, plus the weighing coefficient. Since a viable candidate model needs to be completely specified by the two given flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$, we need to reduce the degrees of freedom to two in the model distribution. We consider a distribution of the form

${\rho }_P\left(p\right)=\frac{1}{q+1}{\text{unif}}_{\left(\mathrm{0,}b\right)}\left(p\right)+\frac{q}{q+1}{\text{unif}}_{\left(\mathrm{0,}\mu b\right)}\left(p\right)$ (25)

Distribution (25) has three parameters: $\mu$, b and q. In this subsection, we fix coefficient $\mu =0.15$ and treat $\left(b\mathrm{,}q\right)$ as two tunable parameters. The resulting formulation is labeled version A. In next subsection, we consider version B, in which b is estimated without using the equations of first two moments, and $\left(\mu \mathrm{,}q\right)$ are treated as tunable parameters and they are determined from the equations of first two moments. For analytical clarity and convenience, when discussing version A, we still denote coefficient $\mu$ as a variable, instead of writing out explicitly $\mu =0.15$. This will allow a unified formulation in variables $\left(\mu \mathrm{,}b\mathrm{,}q\right)$ for both version A and version B.

We formulate the equations for $\left(\mu \mathrm{,}b\mathrm{,}q\right)$ by matching the first two moments. This set of equations is good both for solving version A where $\mu =0.15$, and for solving version B where b is estimated separately. We use the second moment instead of variance because it is simpler to work with the second moment in a mixture distribution. The first two moments of ${\text{unif}}_{\left(\mathrm{0,}b\right)}$ are:

$X\sim {\text{unif}}_{\left(\mathrm{0,}b\right)}\mathrm{<mo>\; :\; </mo>}E\left[X\right]=\frac{1}{2}b$

$E\left[X^2\right]=\frac{1}{3}{b}^2$

For a mixture of uniform distributions, we have

$X\sim \frac{1}{q+1}{\text{unif}}_{\left(\mathrm{0,}b\right)}+\frac{q}{q+1}{\text{unif}}_{\left(\mathrm{0,}\mu b\right)}$

$E\left[X\right]=\frac{q}{q+1}\cdot \frac{1}{2}\left(\mu b\right)+\frac{1}{q+1}\cdot \frac{1}{2}b=\frac{\mu q+1}{2\left(q+1\right)}b$

$E\left[X^2\right]=\frac{q}{q+1}\cdot \frac{1}{3}{\left(\mu b\right)}^2+\frac{1}{q+1}\cdot \frac{1}{3}{b}^2=\frac{{\mu }^{2}q+1}{3\left(q+1\right)}{b}^2$

Equating the first two moments of $P\left(\omega \right)$ to those of model distribution (25), we obtain two equations for $\left(\mu \mathrm{,}b\mathrm{,}q\right)$ :

$\frac{\mu q+1}{2\left(q+1\right)}b=p_a$ (26)

$\frac{{\mu }^{2}q+1}{3\left(q+1\right)}{b}^2={\left(\delta p\right)}^2+p_a^2$ (27)

In this subsection, we discuss version A of model distribution (25), which has $\mu =0.15$. We solve for $\left(b\mathrm{,}q\right)$ from (26) and (27), and express them in terms of flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$. The derivation is presented in Appendix 1A. The solution result is

$\mu =0.15$

$q=\frac{2c_A}{{\left(1-\mu \right)}^2\left(1-\frac{2c_A\mu }{{\left(1-\mu \right)}^2}+\sqrt{1-\frac{4c_A\mu }{{\left(1-\mu \right)}^2}}\right)}$ (28)

$b=\frac{2\left(q+1\right)}{\mu q+1}{p}_a$ (29)

where $c_A\equiv \frac{3{\left(\delta p\right)}^2-p_a^2}{4p_a^2}.$

Corresponding to model distribution (25) for injury probability $P\left(\omega \right)$, the crowd injury outcome N has a compound distribution. It draws $P\left(\omega \right)$ from the mixture of uniform distributions given in (25); then it generates the actual number of injured N using the binomial distribution with parameters $\left(n_c\mathrm{,}P\left(\omega \right)\right)$. The PMF of N is a mixture of two PMFs of the form (24) with $a=0$ :

$P\left(\omega \right)\sim \frac{1}{q+1}{\text{unif}}_{\left(\mathrm{0,}b\right)}+\frac{q}{q+1}{\text{unif}}_{\left(\mathrm{0,}\mu b\right)}\mathrm{,}N\sim {\text{Bin}}_{(n_c\mathrm{,}P(\omega ))}$

${\text{Pr}}_{\left[N=k\right]}=\frac{I_{\left(b\right)}\left(k+\mathrm{1,}{n}_c-k+1\right)+\frac{q}{\mu }{I}_{\left(\mu b\right)}\left(k+\mathrm{1,}{n}_c-k+1\right)}{\left(n_c+1\right)\left(q+1\right)b}$ (30)

This expression of PMF is good for both version A and version B of model distribution (25). In version A, $\mu =0.15$ is fixed, and $\left(q\mathrm{,}b\right)$ are expressed in terms of $\left(p_a\mathrm{,}\delta p\right)$ in (28) and (29). The PMF of N based on version A of (25) is completely specified by the two flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$, along with the number of subjects n_c. This statement is also true for version B, as we will see in next subsection.

3.4. Model 4: Mixture of Uniform Distributions for $P\left(\omega \right)$, Version B

In the previous subsection, a mixture of uniform distributions given in (25) is used as the model distribution. In the absence of additional condition, (25) has three parameters: $\mu$, b and q. In version A discussed in the previous subsection, we fixed $\mu =0.15$, and treated $\left(b\mathrm{,}q\right)$ as unknowns in Equations (26) and (27). From there we solved for $\left(b\mathrm{,}q\right)$ and express them in terms of the given flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$. In this subsection, we consider version B in which we first find a reasonable estimate of b without using Equations (26) and (27), and without using q and $\mu$. With the estimated value of b, we then solve for the two remaining unknowns $\left(q\mathrm{,}\mu \right)$ in Equations (26) and (27).

The single mode uniform distribution given in (19) has support interval $\left(p_a-\sqrt{3}\left(\delta p\right)\mathrm{,}{p}_a+\sqrt{3}\left(\delta p\right)\right)$. When $\sqrt{3}\left(\delta p\right)>p_a$, the support interval extends into the negative side, which make it invalid as a domain for probability values. As a remedy, when $\sqrt{3}\left(\delta p\right)>p_a$, the model distribution is changed to the mixture of uniform distributions given in (25), which has support interval $\left[\mathrm{0,}b\right]$. That motivates us to use the size of interval $\left(p_a-\sqrt{3}\left(\delta p\right)\mathrm{,}{p}_a+\sqrt{3}\left(\delta p\right)\right)$ as an approximation of b in the absence of knowing q and $\mu$. In version B, we set $b=2\sqrt{3}\left(\delta p\right)$, and solve for $\left(q\mathrm{,}\mu \right)$ from Equations (26) and (27). The derivation is presented in Appendix 1B. The solution result is

$b=2\sqrt{3}\left(\delta p\right)$ (31)

$q=\frac{4\left(1-c_B\right)}{1+c_B}$ (32)

$\mu =c_B+\frac{c_B-1}{q}$ (33)

where $c_B\equiv \frac{p_a}{\sqrt{3}\left(\delta p\right)}.$

The corresponding distribution for crowd injury outcome N has the same general expression as given in (30). For version B, the PMF of N is indirectly specified, via (31), (32) and (33), by the two flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$, along with the number of subjects n_c.

3.5. Model 5: Mixture of Triangle and Uniform Distributions for $P\left(\omega \right)$, Version A

As we have seen in previous two subsections, when $\sqrt{3}\left(\delta p\right)\le p_a$, there is a single mode uniform distribution that matches the first two moments of $P\left(\omega \right)$ and that has a support interval appropriate for modeling the distribution of non-negative $P\left(\omega \right)$. Thus, for $\sqrt{3}\left(\delta p\right)\le p_a$, it is reasonable to use this valid uniform distribution to represent $P\left(\omega \right)$. When $\sqrt{3}\left(\delta p\right)>p_a$, however, there is no valid uniform distribution that matches the first two moments of $P\left(\omega \right)$. For $\sqrt{3}\left(\delta p\right)>p_a$, we need to change the form of candidate distributions. The option used in previous two subsections is a mixture of uniform distributions given in (25). In this and next subsections, we study another model distribution.

We consider a mixture of triangle and uniform distributions. Since a general triangle distribution has three degrees of freedom, there are many triangle distributions even when the support interval $\left[a\mathrm{,}b\right]$ is specified. We first introduce the particular triangle distribution over interval $\left[\mathrm{0,}b\right]$ that we will use in the new model distribution.

${\text{tria}}_{\left(\mathrm{0,}b\right)}\left(p\right)\equiv \frac{2}{b}\{\begin{array}{l}\frac{b-p}{b}\text{for}p\in \left[\mathrm{0,}b\right]\\ 0\text{otherwise}\end{array}$ (34)

The first two moments of ${\text{tria}}_{\left(\mathrm{0,}b\right)}$ are:

$X\sim {\text{tria}}_{\left(\mathrm{0,}b\right)}\mathrm{<mo>\; :\; </mo>}E\left[X\right]=\frac{1}{3}b$

$E\left[X^2\right]=\frac{1}{6}{b}^2$

When $\sqrt{3}\left(\delta p\right)>p_a$, we use a mixture of triangle and uniform distributions to represent the distribution of $P\left(\omega \right)$. The mixture takes the form

${\rho }_P\left(p\right)=\frac{1}{q+1}{\text{unif}}_{\left(\mathrm{0,}b\right)}\left(p\right)+\frac{q}{q+1}{\text{tria}}_{\left(\mathrm{0,}\mu b\right)}\left(p\right)$ (35)

Similar to the situation with (25), model distribution (35) has three parameters: $\left(\mu \mathrm{,}b\mathrm{,}q\right)$. To make the model distribution completely specified by flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$, we consider two versions of (35), each having only two degrees of freedom. In version A (discussed in this subsection), coefficient $\mu =0.2$ is fixed, and $\left(b\mathrm{,}q\right)$ are treated as tunable parameters. In version B (discussed in next subsection), b is estimated separately without using the equations of first two moments, and $\left(\mu \mathrm{,}q\right)$ are treated as tunable parameters. Again, we work with the general formulation in three variables $\left(\mu \mathrm{,}b\mathrm{,}q\right)$, to accommodate both version A and version B. We construct the equations for $\left(\mu \mathrm{,}b\mathrm{,}q\right)$ by matching the first two moments of $P\left(\omega \right)$ and the candidate distribution. The first two moments of mixture (35) are:

$X\sim \frac{1}{q+1}{\text{unif}}_{\left(\mathrm{0,}b\right)}+\frac{q}{q+1}{\text{tria}}_{\left(\mathrm{0,}\mu b\right)}$

$E\left[X\right]=\frac{q}{q+1}\cdot \frac{1}{3}\left(\mu b\right)+\frac{1}{q+1}\cdot \frac{1}{2}b=\frac{2\mu q+3}{6\left(q+1\right)}b$

$E\left[X^2\right]=\frac{q}{q+1}\cdot \frac{1}{6}{\left(\mu b\right)}^2+\frac{1}{q+1}\cdot \frac{1}{3}{b}^2=\frac{{\mu }^{2}q+2}{6\left(q+1\right)}{b}^2$

Equating the first two moments of mixture (35) to those of $P\left(\omega \right)$ gives us two equations:

$\frac{2\mu q+3}{6\left(q+1\right)}b=p_a$ (36)

$\frac{{\mu }^{2}q+2}{6\left(q+1\right)}{b}^2={\left(\delta p\right)}^2+p_a^2$ (37)

In version A of mixture (35), we proceed with $\mu =0.2$, and solve for $\left(b\mathrm{,}q\right)$ from (36) and (37). The derivation is given in Appendix 2A. The solution result is

$\mu =0.2$

$q=\frac{2c_A}{\varphi \left(\mu \right)\left(1-\frac{4c_A\mu }{3\varphi \left(\mu \right)}+\sqrt{1-\frac{8c_A\mu }{3\varphi \left(\mu \right)}+\frac{2c_A{\mu }^2}{9{\varphi }^2\left(\mu \right)}}\right)}$ (38)

$b=\frac{6\left(q+1\right)}{2\mu q+3}{p}_a$ (39)

where $c_A\equiv \frac{3{\left(\delta p\right)}^2-p_a^2}{4p_a^2}$ and $\varphi \left(\mu \right)\equiv 1-\frac{4\mu }{3}+\frac{{\mu }^2}{2}.$

When the injury probability $P\left(\omega \right)$ is governed by (35), the corresponding crowd injury outcome N has a compound distribution. It draws $P\left(\omega \right)$ from the mixture of triangle and uniform distributions given in (35). Then it generates the crowd injury outcome N using the binomial distribution with parameters $\left(n_c\mathrm{,}P\left(\omega \right)\right)$. To obtain the PMF of binomial compounded with model distribution (35), we first derive the PMF of binomial compounded with ${\text{tria}}_{\left(\mathrm{0,}b\right)}$, which is expressed as an integral of binomial PMF weighted by density

${\text{tria}}_{\left(\mathrm{0,}b\right)}\left(p\right)=\frac{2}{b}-\frac{2}{b^2}{p}^2$.

$P\left(\omega \right)\sim {\text{tria}}_{\left(\mathrm{0,}b\right)}\mathrm{,}N\sim {\text{Bin}}_{(n_c\mathrm{,}P(\omega ))}$

$\begin{array}{c}{\mathrm{Pr}}_{\left[N=k\right]}=C\left(n_c,k\right){\displaystyle {\int }_0^b}\left(\frac{2}{b}-\frac{2}{b^2}p\right)p^k{\left(1-p\right)}^{\left(n_c-k\right)}\text{d}p\\ =C\left(n_c,k\right)\left(\frac{2}{b}B\left(b;k+1,n_c-k+1\right)-\frac{2}{b^2}B\left(b;k+2,n_c-k+1\right)\right)\end{array}$

We recall identity (23) and the identity

$B\left(b;u+1,v\right)=\frac{\Gamma \left(u+1\right)\Gamma \left(v\right)}{\Gamma \left(u+v+1\right)}=\frac{u}{u+v}B\left(b;u,v\right)$ (40)

Using these identities, we write out the PMF of binomial compounded with ${\text{tria}}_{(0,b)}$

${\mathrm{Pr}}_{\left[N=k\right]}=\frac{2I_{\left(b\right)}\left(k+1,n_c-k+1\right)-\frac{2\left(k+1\right)}{\left(n_c+2\right)b}{I}_{\left(b\right)}\left(k+2,n_c-k+1\right)}{\left(n_c+1\right)b}$ (41)

The PMF of binomial compounded with distribution (35) is a mixture of two PMFs: one of form (24) with $a=0$ and the other of form (41).

$P\left(\omega \right)\sim \frac{1}{q+1}{\text{unif}}_{\left(\mathrm{0,}b\right)}+\frac{q}{q+1}{\text{tria}}_{\left(\mathrm{0,}\mu b\right)}\mathrm{,}N\sim {\text{Bin}}_{(n_c\mathrm{,}P(\omega ))}$

$\begin{array}{l}{\text{Pr}}_{\left[N=k\right]}=\frac{1}{\left(n_c+1\right)\left(q+1\right)b}(I_{\left(b\right)}\left(k+\mathrm{1,}{n}_c-k+1\right)+\frac{2q}{\mu }{I}_{\left(\mu b\right)}\left(k+\mathrm{1,}{n}_c-k+1\right)\\ -\frac{2q\left(k+1\right)}{\left(n_c+2\right){\mu }^{2}b}{I}_{\left(\mu b\right)}\left(k+\mathrm{2,}{n}_c-k+1\right))\end{array}$ (42)

PMF expression (42) of crowd injury outcome N is good for both version A and version B of (35). In version A, $\mu =0.2$ is fixed, and $\left(q\mathrm{,}b\right)$ are expressed in terms of $\left(p_a\mathrm{,}\delta p\right)$ in (38) and (39). Thus, PMF (42) is indirectly specified by the two flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$, along with the number of subjects n_c. This statement is also true for version B of (35), as we will see in next subsection.

3.6. Model 6: Mixture of Triangle and Uniform Distributions for $P\left(\omega \right)$, Version B

The mixture of triangle and uniform distributions given in (35) for modeling $P\left(\omega \right)$ has three parameters: $\left(\mu \mathrm{,}b\mathrm{,}q\right)$. To determine $\left(\mu \mathrm{,}b\mathrm{,}q\right)$ from the two given flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$, we need to introduce additional condition. In version A (discussed in previous subsection), we solved for $\left(b\mathrm{,}q\right)$ at fixed $\mu =0.2$. In this subsection, we consider version B in which we set $b=2\sqrt{3}\left(\delta p\right)$ and solve for the two remaining unknowns $\left(q\mathrm{,}\mu \right)$ from Equations (36) and (37). The derivation is presented in Appendix 2B. The solution result is

$b=2\sqrt{3}\left(\delta p\right)$ (43)

$q=\frac{3\left(1-c_B\right)}{2c_B-1+\sqrt{c_B^2-4c_B+3}}$ (44)

$\mu =\frac{3}{2}\left(c_B+\frac{c_B-1}{q}\right)$ (45)

where $c_B\equiv \frac{p_a}{\sqrt{3}\left(\delta p\right)}$.

The corresponding distribution of crowd injury outcome N has the same general expression as given in (42). For version B of mixture (35), PMF (42) is indirectly specified by the two flash bang parameters $\left(p_a\mathrm{,}\delta p\right)$, via (43), (44) and (45), along with the crowd size n_c.

4. Comparison of the Six Models

We examine the six model distributions for $P\left(\omega \right)$ formulated in the previous section, each parameterized by $\left(p_a\mathrm{,}\delta p\right)$. We compare their performances in predicting the distribution of crowd injury outcome N. We judge the performance using the difference between the true distribution of N and the distribution of N predicted from $\left(p_a\mathrm{,}\delta p\right)$ using each of the six models. At each problem setup, we first run Monte Carlo simulations of the full computational model [5] to generate large number of independent realizations of crowd injury outcome N, and we use the collection of independent samples to approximate the true distribution of N. In a two-parameter phenomenological model for N, the injury effect of the flash bang situation is characterized by $\left(p_a\mathrm{,}\delta p\right)$, calculated from the Monte Carlo samples of N. The pair $\left(p_a\mathrm{,}\delta p\right)$ serves as a two-parameter summary of the flash bang situation. This is similar to RSI serving as a single-parameter characterization of the flash bang situation. A single-parameter or two-parameter characterization of the flash bang situation is a high level coarse description without the details of the full computational model. We study how well the behavior of crowd injury outcome can be predicted solely based on such a high level coarse description. We apply each of the six models to predict the distribution of N using only the values of $\left(p_a\mathrm{,}\delta p\right)$, and we calculate the prediction accuracy of each model accordingly.

We test the six models over various problem setups that yield a diverse range for characterization parameters $\left(p_a\mathrm{,}\delta p\right)$. Specifically, let $\sigma =\text{std}\left(\text{aimingerror}\right)$, the standard deviation of random aiming error. We carry out tests for $\sigma \in \left[0,30\text{m}\right]$. At each value of σ, Monte Carlo runs are repeated to generate a set of $m={10}^6$ independent realizations of N. This large sample set serves two roles in our comparison study: a) it provides the observed distribution of N (which we regard as the true distribution), and b) we use it to calculate the characterization parameters $\left(p_a\mathrm{,}\delta p\right)$. The predicted distribution of N uses only the information of $\left(p_a\mathrm{,}\delta p\right)$, not the Monte Carlo sample set. The performance of each model is evaluated by comparing the predicted distribution and the observed distribution. Table 1 displays values of $\sigma =\text{std}\left(\text{aimingerror}\right)$ vs parameters $\left(p_a\mathrm{,}\delta p\right)$ and coefficient $\sqrt{3}\left(\delta p\right)/p_a$. Recall that in models 3 - 6, condition $\sqrt{3}\left(\delta p\right)/p_a>1$ is used as the indicator when switching from a single mode uniform distribution for $P\left(\omega \right)$ to a mixture distribution. The results of the six candidate models formulated in previous section are displayed in Figures 1-6. For each candidate model, we calculate the errors in the predicted distribution of N relative to the observed distribution. We calculate both errors in the probability mass function (PMF) and those in the cumulative distribution function (CDF). In each figure, we compare the observed and the predicted distributions of N at aiming errors of various magnitudes, and we study the prediction errors vs. $\sigma =\text{std}\left(\text{aimingerror}\right)$, the standard deviation of random aiming error. Each of Figures 1-6 has six panels listed below.

Top left panel: comparison of the observed distribution of N at σ = 10 m and the distribution predicted from the values of $\left(p_a\mathrm{,}\delta p\right)$ using a two-parameter model distribution for injury probability $P\left(\omega \right)$. The observed distribution of N in Monte Carlo simulations of the full model is regarded as the true distribution.

Top right panel: error in PMF and error in CDF at σ = 10 m.

Middle left panel: comparison of the observed distribution of N at σ = 25 m and the distribution predicted from the values of $\left(p_a\mathrm{,}\delta p\right)$ using a two-parameter model.

Middle right panel: error in PMF and error in CDF at σ = 25 m.

Bottom left panel: maximum errors over all $n\ge 0$, respectively, in PMF and in CDF, as functions of σ.

Bottom right panel: Comparison of the observed and the predicted $p\left(N>10\right)$, as functions of σ. Here $p\left(N>10\right)\equiv \mathrm{Pr}\left({N>10|}_{n_c=100}\right)$ is the probability of having more than 10 injured in a crowd of 100 subjects.

Model 1 is actually a single-parameter description of flash bang situation specified

by RSI. As shown in Figure 1, Model 1 has large errors in predicted distribution and in predicted probabilities; in particular, the predicted probability $\mathrm{Pr}\left(N>10\right)$ is very different from the observed one. Models 2 - 6 are two-parameter descriptions of flash bang situation specified by $\left(p_a=\text{RSI}\mathrm{,}\delta p\right)$. In comparison with Model 1, models 2 - 6 all show significant improvements in predictions (Figures 2-6). Figure 7 compares the six models by plotting the maximum error in the predicted PMF of N as a function of $\sigma$ for each of the six models. Overall, the highest accuracy

is achieved by model 6. For aiming errors of magnitude $\sigma \le 28\text{m}$, the error in PMF $p\left(n\right)$ for model 6 is less than 2% over all $n\ge 0$. The maximum error in PMF $p\left(n\right)$ over $n>2$ is smaller than the unconstrained maximum. Figure 8 examines the accuracies of the six models in predicting probabilities $p\left(N>10\right)$ and $p\left(N>5\right)$ where $p\left(N>n_{\text{tolerance}}\right)$ is the probability of having more than $n_{\text{tolerance}}$ injured in a crowd of 100 subjects. Again, the highest accuracy is achieved by model 6, with error in $p\left(N>10\right)$ bounded by 0.0054 (0.54%) and

error in $p\left(N>5\right)$ bounded by 0.016 (1.6%), for aiming errors of magnitude $\sigma \le 30\text{m}$.

5. Summary

We studied the injury outcome caused by multiple flash bang submunitions. The injury outcome for a crowd of n_c subjects is measured by the actual number of injured N. The crowd injury outcome N is random, influenced by two factors: i)

the hidden injury probability $P\left(\omega \right)$ corresponding to a realized spatial distribution $\omega$ of n_s submunitions, and ii) the randomness other than the submunitions distribution. Factor i) contains the aiming error of flash bang mortar and the uncertainty in submunitions dispersion upon mortar burst. Factor ii) includes the uncertainty in the amount of acoustic dose propagating toward and reaching individual subjects and the biovariability of individual subjects’ reactions to a given acoustic dose received.

The distribution of crowd injury outcome N is fully described by the Monte Carlo computational model with all relevant parameters [5] . For practical field operations, it is desirable to avoid the complexity of the Monte Carlo model, and to do prediction based on a concise description of N using only one or two parameters. Conventionally, the overall injuring effect of the flash bang situation is simply characterized by the average fraction of injured, which is called the risk of significant injury (RSI). While the RSI description achieves the simplicity of only one parameter, it lacks information about fluctuations and uncertainties in the

actual number of injured. In particular, the RSI description excludes the effect of random aiming error, which tends to change the injury probability up and down simultaneously for all individual subjects in a crowd. To capture the uncertainty in the crowd injury outcome, we adopt a two-parameter framework in which an injury probability $P\left(\omega \right)$ is drawn from a two-parameter model distribution and then given the injury probability drawn, the conditional crowd injury outcome is modeled as a binomial random variable. The model distribution for injury probability $P\left(\omega \right)$ is specified by parameters $p_a\equiv E\left[P\left(\omega \right)\right]$ and $\delta p\equiv \text{std}\left[P\left(\omega \right)\right]$. The selection of these two parameters is motivated by their unique properties. First, $\left(p_a\mathrm{,}\delta p\right)$ as the mean and standard deviation, are a natural choice for characterizing the distribution of $P\left(\omega \right)$. Second, despite the fact that the hidden injury probability $P\left(\omega \right)$ is not observable, quantities $\left(p_a\mathrm{,}\delta p\right)$ can be calculated based on observed samples of crowd injury outcome N. The resulting two-parameter model for N is completely specified by $\left(p_a\mathrm{,}\delta p\right)$ and it serves as a high level coarse description for the flash bang situation.

To clarify, the two-parameter description of N is not meant to completely replace the comprehensive Monte Carlo model as a tool for investigating the effects of various parameters on injury outcome. Rather, it is meant to provide an improved alternative to the single-parameter RSI description as a high level coarse description. Theoretically, parameters $\left(p_a\mathrm{,}\delta p\right)$ are the mean and standard deviation of the non-observable injury probability $P\left(\omega \right)$. Operationally, parameters $\left(p_a\mathrm{,}\delta p\right)$ are estimated from observed samples of random injury outcome N, and the flash bang situation is specified by the estimated values of $\left(p_a\mathrm{,}\delta p\right)$ in a high level coarse description. We formulated and explored six candidate models for approximating the distribution of hidden injury probability $P\left(\omega \right)$, each parameterized by $\left(p_a\mathrm{,}\delta p\right)$. Once the distribution of $P\left(\omega \right)$ is constructed, the crowd injury outcome N is predicted as the compound distribution $N\sim {\text{Bin}}_{\left(n_c\mathrm{,}P\left(\omega \right)\right)}$. Thus, a two-parameter model for $P\left(\omega \right)$ leads to a two-parameter description for the flash bang situation. We tested six models in their performances of predicting the distribution of N. The best performer is model 6 described in subsection 3.6, which uses a mixture of uniform and triangle distributions for $P\left(\omega \right)$. The corresponding two-parameter description for N yields fairly accurate prediction, using only the values of $\left(p_a\mathrm{,}\delta p\right)$. This two-parameter description allows us to estimate the probability of crowd injury outcome exceeding a prescribed tolerance, for example, the probability of having more than 10 injured in a crowd of 100 subjects. We advocate the adoption of this two-parameter description to replace the conventional single-parameter RSI description. Whenever we need to give a high level coarse description of a flash bang situation, instead of stating that the average injury risk is RSI, we say that the injury risk is represented by the pair $\left(p_a=\text{RSI}\mathrm{,}\delta p\right)$. We will explore extending and revising this framework of two-parameter description to accommodate other situations with uncertain outcomes.

Acknowledgements and Disclaimer

The authors thank the Joint Non-Lethal Weapons Directorate of U.S. Department of Defense and the Naval Postgraduate School for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the U.S. Government.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Wang, H.Y., Simms, M., Labaria, G., Burgei, W.A. and Zhou, H. (2019) A Two-Parameter Description for the Injury Risk of Flash Bangs with Uncertainty. American Journal of Operations Research, 9, 79-108. https://doi.org/10.4236/ajor.2019.93005

References

Appendices

Appendix 1A

Solving for $\left(b\mathrm{,}q\right)$ from (26) and (27) at a given $\mu .$

We divide (27) by the square of (26) to eliminate b. The resulting equation for q is

$\frac{4\left(q+1\right)\left({\mu }^{2}q+1\right)}{3{\left(\mu q+1\right)}^2}=\frac{{\left(\delta p\right)}^2+p_a^2}{p_a^2}$

We rewrite it as

$\frac{\left(q+1\right)\left({\mu }^{2}q+1\right)}{{\left(\mu q+1\right)}^2}=1+c_A$

where $c_A\equiv \frac{3{\left(\delta p\right)}^2-p_a^2}{4p_a^2}.$

Multiplying the equation for q by its denominator and moving all terms to one side, we get a quadratic equation

$c_A{\mu }^2{q}^2-\left[{\left(1-\mu \right)}^2-2c_A\mu \right]q+c_A=0$

Solving the quadratic equation gives us two solutions for q

$q=\frac{2c_A}{\left[{\left(1-\mu \right)}^2-2c_A\mu \right]\pm \sqrt{{\left[{\left(1-\mu \right)}^2-2c_A\mu \right]}^2-4c_A^2{\mu }^2}}$

As functions of $c_A$, the two branches of solution satisfy respectively

$\left(\text{Branch}+\right)\mathrm{<mo>\; :\; </mo>}\underset{c_A\to 0+}{\lim }q\left(c_A\right)=0+$

$\left(\text{Branch}-\right)\mathrm{<mo>\; :\; </mo>}\underset{c_A\to 0+}{\lim }q\left(c_A\right)=+\infty$

In the mixture of uniform distributions, $\frac{1}{q+1}{\text{unif}}_{\left(\mathrm{0,}b\right)}\left(p\right)+\frac{q}{q+1}{\text{unif}}_{\left(\mathrm{0,}\mu b\right)}\left(p\right)$, coefficient q is the relative weight of the second uniform mode. Recall that this second uniform mode is used only in the case of $\sqrt{3}\left(\delta p\right)>p_a$, which corresponds to $c_A>0$. For $c_A\le 0$, only the primary uniform mode, ${\text{unif}}_{\left(\mathrm{0,}b\right)}\left(p\right)$, is used; the second uniform mode is absent. Mathematically, when $c_A\le 0$, the second uniform mode has zero weight. For $c_A$ slightly above zero, the second uniform mode is present but should have a small relative weight $q\left(c_A\right)$. Following this reasoning, we select (Branch+). Simplifying the expression of q in (Branch+) and then using equation (26) to express b in terms of q, we obtain the solution of $\left(q\mathrm{,}b\right)$

$q=\frac{2c_A}{{\left(1-\mu \right)}^2\left(1-\frac{2c_A\mu }{{\left(1-\mu \right)}^2}+\sqrt{1-\frac{4c_A\mu }{{\left(1-\mu \right)}^2}}\right)}$

$b=\frac{2\left(q+1\right)}{\mu q+1}{p}_a$

which are reported as equations (28) and (29) in the main text.

Appendix 1B

Solving for $\left(q\mathrm{,}\mu \right)$ from (26) and (27) at $b=2\sqrt{3}\left(\delta p\right).$

Using (26), we express $\mu$ in terms of q:

$\mu =c_B+\frac{c_B-1}{q}$

where

$c_B\equiv \frac{2p_a}{b}=\frac{p_a}{\sqrt{3}\left(\delta p\right)}$ (46)

Note that $c_B$ defined in (46) is always positive. Substituting $\mu$ into (27), we get

${\left(c_B+\frac{c_B-1}{q}\right)}^{2}q+1=\left(q+1\right)c_2$