- Research
- Open Access

# Global recording rates for large eruptions

- Jonathan Rougier
^{1}Email author, - Stephen R. Sparks
^{2}and - Katharine V. Cashman
^{2}

**5**:11

https://doi.org/10.1186/s13617-016-0051-4

© The Author(s) 2016

**Received:**22 May 2015**Accepted:**18 September 2016**Published:**10 October 2016

## Abstract

A non-parametric statistical approach is used to assess the global recording rate for large (M4+) stratovolcano eruptions in a modern database, LaMEVE (v3.1). This approach imposes minimal structure on the shape of the recording rate through time. We find that recording rates have declined rapidly, going backwards in time. Prior to the year 1600 they are below 50 *%*, and prior to 1100 they are below 20 *%*. Even in the recent past, e.g. the 1800s, they are likely to be less than 100 *%*. The assessment for very large (M5+) eruptions is more uncertain, due to the scarcity of events.

## Keywords

- Under-recording
- Non-parametric
- Confidence set

## Background

Under-recording is an unavoidable issue in empirical studies of large volcanic eruptions, as has been discussed by, e.g., Simkin (1993), Siebert et al. (2010, pp. 31–34) and Brown et al. (2014). As Simkin (1993, p. 436) puts it, the rise in the number of recorded eruptions follows from “more observers, in a wider geographical distribution, with better communication, and broader publication.” Statisticians refer to this type of absence as ‘missing not at random’, and it must be accounted for when making inferences from observations, in order to remove, as far as possible, systematic biases. But this under-recording poses a major challenge, because we have few quantitative beliefs about the way in which the timing, location, and magnitude of an eruption combine together to affect its missingness from the modern record.

For an initial analysis, we can make one drastic simplification at the outset, and consider global volcanism above a specified large magnitude threshold. This reduces three dimensions of missingness to just one—timing. It also simplifies the statistical modelling: first, because we avoid complicated issues regarding the modelling and accuracy of magnitude measurements; second, because the process of global volcanism is likely to be more stable than its component parts. These issues are discussed further below (Section “Our statistical non-parametric approach”).

This paper provides an empirical analysis of the global recording rate of large eruptions of stratovolcanoes; the definition of the global recording rate is discussed and clarified in Section “Definition of the recording rate”. ‘Large’ is taken to be at least 100 million tonnes of ejected matter, or M4 according to the magnitude scale of Mason et al. (2004), for which *M* := log10*m*−7, where *m* is the mass erupted in kilograms. *M* thus provides a continuous measure of eruption magnitude based on estimates of erupted mass; more common is the use of the Volcanic Explosivity Index (VEI), which is discrete and assigned on the basis of either erupted volume or eruption intensity (average mass eruption rate).

Our key innovation is to use a non-parametric approach (Section “Our statistical non-parametric approach”), and to present our results visually in terms of a 95 % confidence set (Section “Results and discussion”). We are able to quantify recording rates through time in terms of upper bounds of 95 % confidence intervals, and show, for example, that this upper bound is below 50 % prior to the year 1600.

## Methods

### Definition of the recording rate

The notion of a ‘global recording rate’ is imprecise. We provide a definition which highlights the difference between the global recording rate and the mean recording rate across volcanoes.

Imagine a process, extensive in space and time, which indicates with a 1 or a 0 whether there is an information link from a large eruption at a given location and date, through to today and an appearance in the database. Index the active stratovolcanoes with *i*=1,…,*n*, and write *c*
_{
i
}(*t*) for this linkage function at the location of volcano *i* and time *t*. If we observed this function for a given *i* going backwards through time, we would see it at ‘on’ (i.e. 1) for today, but then at some point in the past it would start to flicker between ‘on’ and ‘off’, indicating that a link to today is contingent on there being a suitable person present at the eruption, or traces (e.g. a near-field tephra deposit) that have survived and been discovered and interpreted. As we go back further, the function spends longer and longer amounts of time ‘off’, until finally it is permanently ‘off’. This sequence, from currently ‘on’, through flickering, to ‘off’, varies from volcano to volcano. A global recording function reflects the amalgamation of such functions for the global collection of volcanoes.

*t*as

where *λ*
_{
i
} is the large-eruption rate of volcano *i*, and *λ* :=*λ*
_{1}+⋯+*λ*
_{
n
}, the global large-eruption rate; these rates are taken to be constant in time—see the next paragraph for the generalisation. The numerator of (1) is the database’s global large-eruption rate for time *t*, which can be written as *π*(*t*)·*λ*, where *λ* is today’s global large-eruption rate. Note that *π*(*t*)→1 as *t* tends to today, because each *c*
_{
i
}(*t*)→1. In (1), we are not assuming that large eruptions follow a Poisson process for each volcano, although we are assuming that we can combine rates in the same way that we would combine rates were the volcanoes to be following independent Poisson processes. In Section “Our statistical non-parametric approach”, we will treat the global large-eruption process as a time-homogeneous Poisson process, but we will not need to assert that the individual processes are either Poisson or stationary. This issue has been discussed by De la Cruz-Reyna (1991).

*c*

_{ i }(·)’s. Let

*λ*

_{ i }be today’s large-eruption rate for volcano

*i*, and then write

*c*

*i*′(

*t*) :=

*c*

_{ i }(

*t*)·

*λ*

_{ i }(

*t*)/

*λ*

_{ i }. We expect that

*c*

_{ i }(

*t*)≈

*c*

*i*′(

*t*), because this is exactly true when

*c*

_{ i }(

*t*)=0 or

*t*tends to today, and approximately true elsewhere, if the fluctuations in

*λ*

_{ i }(

*t*) are not too large, in proportionate terms.

*j*of width

*Δ*centered at time

*t*

_{ j }. Starting with (1), take the time-average over the period to find the global large-eruption rate for period

*j*:

The value *π*
_{
ij
} is the recording rate of volcano *i* in period *j*. It is not a probability, but it can be interpreted as one, under additional conditions. If a large eruption of volcano *i* is randomly and uniformly distributed in period *j*, then *π*
_{
ij
} is the probability of that eruption being recorded.

*j*is approximately the mean recording rate across all volcanoes. In general, though, (2a) shows that the global recording rate is biased away from the mean recording rate towards the recording rates of the most active volcanoes. So, for example, if recordings of a very active volcano start at the end of period

*j*, then

*π*

_{ j }to

*π*

_{ j+1}will show a large increase, larger than the increase that would be expected for the addition of one extra volcano. Just monitoring

*π*

_{ j }alone, it would not be possible to distinguish between, say, three newly-recorded volcanoes with modest eruption rates, and one newly-recorded volcano with a high eruption rate.

This distinction between the global recording rate and the mean recording rate is an important part of the definition and interpretation of the global recording rate, but whether it really matters is moot. Large eruptions are scarce. We do not expect to constrain our uncertainty about *π*
_{
j
} to within the scale of the distinction, and from this point of view our uncertainty about *π*
_{
j
} is more-or-less the same as our uncertainty about the mean recording rate, which is how we will interpret our results below.

### Review of previous approaches

*t*

_{ j }, of width

*Δ*

_{ j }, which has

*x*

_{ j }recorded large eruptions. Equating the expectation with the observation suggests that

*Δ*

_{ j }·

*π*

_{ j }·

*λ*≈

*x*

_{ j }, where

*λ*is the global large-eruption rate, defined after (1), and taken to be effectively constant over the time-frame of the analysis;

*π*

_{ j }is defined in (2a). Let the most recent period, period

*k*say, have a recording rate of

*π*

_{ k }=1. Then

*Δ*

_{ k }·

*λ*≈

*x*

_{ k }, and substituting for

*λ*gives

Naïve though this is, it will turn out to be the Maximum Likelihood (ML) estimator in the model proposed in Section “Our statistical non-parametric approach”. It is, of course, a woeful estimator, being very over-fitted: it has one parameter per observation, and no constraints.

There are several ways to proceed. One approach is to smooth the *x*
_{
j
}’s, or equivalently the *π*
_{
j
}’s. This approach was implemented in Guttorp and Thompson (1991), who discussed the critical issue of choosing a smoother and a bandwidth (their p. 579). There was no assessment of uncertainty in their estimated *π*
_{
j
}’s, although this was not their main concern. We present a similar estimate below, at the start of Section “Results and discussion”.

A second approach is to propose a specific functional form for *π*(·), which is controlled by a small number of unknown parameters; this imposes smoothness on *π*(·) and reduces over-fitting. This approach was originally suggested in a more general context by Solow (1993), who proposed a monotonic four-parameter functional form based on the distribution function of the Beta distribution. The parameters of the proposed function are fitted to the *x*
_{
j
}’s, which typically requires a statistical model for the number of recorded large-eruptions, if uncertainty is to be assessed. This is the parametric approach that has been used in the volcanological literature to date.

*π*(·), used by Coles and Sparks (2006) and Deligne et al. (2010), is

where *t*
_{0} and *b* are the parameters, *t*
_{0}≤*t*
_{1}, and *b*>0. This functional form encodes two basic beliefs: that the recording rate is currently 1, and that it decreases smoothly backwards in time. A very similar approach is used by Kiyosugi et al. (2015), who propose an exponential curve.

where *π*
_{0} and *t*
_{
s
} are the parameters, with 0≤*π*
_{0}≤1 and *t*
_{1}≤*t*
_{
s
}≤*t*
_{
k
}. Equation (5) defines a step function that moves abruptly from *π*
_{0} to 1 at time *t*
_{
s
}. This functional form encodes the same type of belief, except replacing ‘smoothly backwards in time’ with ‘dramatically at some time previous to today’. Furlan (2010) uses a single global *π*
_{0} and *t*
_{
s
}, while Mead and Magill (2014) allow *π*
_{0} and *t*
_{
s
} to vary by region.

Mead and Magill’s (2014) model is consistent with smooth functions such as (4) at the global level, as our analysis in Section “Definition of the recording rate” shows. In our terms, their model asserts that \(\pi _{ij} = \pi _{0}^{(i)}\) for \(t_{j} \leq t_{s}^{(i)}\), and 1 thereafter, the superscript *i*’s indicating that these parameter values differ by volcano or region. If different regions have different \(t_{s}^{(i)}\) values, then the total effect on the global large-eruption rate may look like a smooth curve, or a smooth curve with kinks, where a volcano or region with a high *λ*
_{
i
} starts to be fully recorded. But a smooth curve would also arise in the more general case where, instead of jumping directly from \(\pi _{0}^{(i)}\) to 1, the *π*
_{
ij
}’s increased more gradually.

The existence of these two different functions for *π*(·) indicates a difficulty with this parametric approach: there is little guidance about the shape of the functional form. In fact there is more disagreement than we have presented here, because both functions also allowed for the effect of magnitude, but in different ways.

A third approach, adopted in this paper and new to volcanology, is to use a non-parametric representation of *π*(·), subject to some simple conditions. Because the conditions do not prescribe a tight relationship between the *π*
_{
j
}’s, this option is sometimes known as a ‘soft’ constraint, compared to the second (parametric) approach, which would be a ‘hard’ constraint. A hard constraint is more effective if the parametric form of *π*(·) is known, but a soft constraint is more widely applicable, and more suitable for dealing with observations affected by numerous complex factors related to recording of historical and geological events. As we show in Section “Our statistical non-parametric approach”, a soft constraint can generalize both of the models given above. And as shown in Section “Results and discussion”, although the uncertainty is large, the soft constraint still entails useful statements about temporal changes in the global recording rate for large eruptions.

### Our statistical non-parametric approach

The notation is the same as the previous sections. Let there be *k* non-overlapping periods, with known lengths *Δ*
_{1},…,*Δ*
_{
k
}, uncertain recording rates π :=(*π*
_{1},…,*π*
_{
k
}), and for which the number of recorded large eruptions is x :=(*x*
_{1},…,*x*
_{
k
}). A statistical assessment of π, including a quantification of uncertainty, requires a probabilistic model for x.

We are considering global sums of large eruptions, and therefore we treat x as arising from a homogeneous Poisson process with unknown rate denoted *λ* (units of /yr). This simple model is suggested by De la Cruz-Reyna (1991), who also noted that it would be less applicable for individual volcanoes. It is not possible to test the time-homogeneity of this model if the recording rate varies by period. Instead, homogeneity becomes the condition under which the recording rates are identified (Solow 1993, makes a similar point).

Most of the papers discussed in Section “Review of previous approaches” (Coles and Sparks 2006; Deligne et al. 2010; Furlan 2010; Mead and Magill 2014) used a 2D Poisson process in time and magnitude, which is homogeneous in time. Because this implies a 1D Poisson process which is homogeneous in time, our model is more general. In lumping all eruptions of at least M4 together, we have ignored information in the database which these other papers have used. This is a deliberate choice on our part, due to the difficulty of incorporating the magnitude information, which we consider more than offsets the benefits of the extra information.

First, a 2D Poisson process over time and magnitude has many more possible implementations than a 1D process over time alone (which has only one). The common choice of all four papers cited above is to parameterise this process according to the peaks over threshold (POT) model for extremes, based on asymptotic theory (Pickands 1971). This requires three parameters rather than just one, and there is also the delicate issue of whether the asymptotic conditions hold, or whether some other justification might be found for the POT model.

*a fortiori*, this accuracy has not changed through time. But magnitude is always imperfectly inferred. Mead and Magill (2014) and Kiyosugi et al. (2015), for example, implicitly assume that VEIs have been accurately categorized for several thousands of years, despite the very different ways in which VEIs have been assessed for eruptions occurring at different times over the Holocene and the Quaternary, respectively. Evidence of misrecording of magnitudes over the last millenium is given in Fig. 1, which shows the magnitudes ‘piling up’ on the integer values 4 and 5. Furlan (2010, sec. 5.4) shows that inferences are sensitive to misrecording of magnitudes.

So, in suppressing the magnitude information, our results are robust to the specification of the 2D Poisson process over time and magnitude, and robust to magnitude measurement errors. This is not the case for the results of the other papers. The smooth-in-time model used by Coles and Sparks (2006) and Deligne et al. (2010) is also sensitive to timing errors, while the abrupt-in-time model used by Furlan (2010) and Mead and Magill (2014) is insensitive to timing errors, as is our model (because we bin the times into intervals).

*λ*) is

*λ*is a nuisance parameter which can be eliminated, in order to focus on π. Let

*λ*have a Gamma distribution with shape

*a*and rate

*b*(units of yr), where both these values are specified. Then integrating out

*λ*gives the likelihood function for π:

where *Γ* is the Gamma function, and the final line discards all multiplicative constants not involving π see, e.g., Lunn et al. (2013, ch. 3). A similar approach was used by Solow (2001).

We need to specify *a* and *b* in the Gamma distribution for *λ*. The functional form of (7) suggests that *a* and *b* have a natural interpretation, in terms of dataset augmentation. The marginal distribution for *λ* carries the same weight as an additional period of length *b*yr which has *a* large eruptions. Informed by previous studies (e.g. Deligne et al. 2010), we assess the global rate of large eruptions at about one every two years, although acknowledge that this value is ‘contaminated’ by our exposure to the dataset x. So we set *a*=1 and *b*=2yr, with this small value for *b* indicating that, putting aside the information in the dataset, we have only vague beliefs about *λ*. Very similar results would follow in passing to the (improper) limit of *a*=0 and *b*=0yr, which is completely vague. However, it is not possible to simulate *λ* from this limit, and simulation is required for code verification and for finding confidence sets.

which can be verified by substituting into the first-order conditions. This expression is intuitive, because *a*/*b* is the expectation of *λ*, and so \({\Delta _{j} \cdot \hat {\pi }_{j} \cdot E(\lambda) = x_{j}}\), similar to the starting-point in Section “Review of previous approaches”. Replacing *E*(*λ*) with the estimate *x*
_{
k
}/*Δ*
_{
k
} then leads to (3), although, as already discussed, this is a very poor estimator. As a Reviewer has noted, the ML estimator in (8) is insensitive to the strength of the information in the marginal distribution for *λ*, because it depends only on *a*/*b*, rather than on the value of *b* directly.

- 1.
The global large-eruption recording rate since 1980 is 1.

- 2.
For sufficiently long periods, the global large- eruption recording rate of an earlier period is never higher than that of a later period.

taking the final period to be 1980–2014. This parameter space generalizes the parametric models given in (4) and (5), since it is possible for the *π*
_{
j
}’s to decline smoothly, or abruptly, or both. The constraints in the parameter space allow observations from later periods to influence the estimates of *π*
_{
j
} in earlier periods; in this way these early *π*
_{
j
}’s can ‘borrow strength’ from the later periods, to make up for the paucity of observations.

Although *Ω* is hard to explore in a sequential fashion, not having a simple interior, it is easy to enumerate a large set of elements of *Ω*, from which \(\hat {\boldsymbol {\pi }}\) can be approximated using the maximum value from the set. This algorithm is described in the Appendix, along with a method for computing a 95 *%* confidence set for π∈*Ω*, which implies a 95 *%* confidence interval for each *π*
_{
j
}.

As the Appendix demonstrates, it is quite hard to assess uncertainty about π. A Bayesian approach is more straightforward, but requires a prior distribution on *Ω*. However, it is not possible to specify such a prior distribution without introducing additional beliefs to the two basic beliefs given above. ‘Obvious’ priors such as uniform on *Ω* are actually highly constricting, since they impose a linear relationship on the prior expectations of the *π*
_{
j
}’s, namely *E*(*π*
_{
j
})=*j*/*k*. Thus, for a large *k*, the prior distribution for *π*
_{1} is forced towards zero. We do not want to incorporate any additional beliefs about π, and therefore we have eschewed a Bayesian approach.

## Results and discussion

The dataset is drawn from the LaMEVE database, see Crosweller et al. (2012) and Brown et al. (2014), version 3.1, downloaded Oct 2015. Our code, including the dataset, is available from the first author, implemented in the statistical computing environment R (R Core Team 2013).

*%*confidence set for π for large (

*M*≥4) eruptions. Technically this is a Parallel Coordinates Plot (see, e.g., Venables and Ripley 2010, sec. 11.1); each grey line represents a point in the 95

*%*confidence set. The initial value of 1 and non-positive slope (going back in time) are imposed by

*Ω*, but the broadly exponential shape reflects the dataset, as does the low recording rate c1100. These features are consistent with the Coles and Sparks (2006) model, given in (4) above, with a value of

*b*>1. As anticipated, there seems to be a discontinuity around 1550, which is not a feature of the Coles and Sparks (2006) model. This discontinuity was picked up in the parametric step-change model fitted by Furlan (2010), given in (5) above, which put the step-change just prior to 1600. Furlan’s (2010) dataset is similar but not identical to ours, as she uses all large eruptions, not just those for stratovolcanoes, and we have extended the dataset to include recent eruptions. But her claim that “the under-recording process largely disappears in the most recent 400 years” (Abstract, p. 113) is refuted in our analysis.

The confidence intervals for individual *π*
_{
j
}’s are large, but the main message from Fig. 4 is clear: the global recording rates for large eruptions decline rapidly going back in time. Prior to 1600 they are below 50 *%*, and prior to 1100 they are below 20 *%*. Even in the recent past, e.g. the 1800s, they are likely to be significantly less than 100 *%*. As a technical point, the error bars in Fig. 4 are margins of a 95 *%* confidence set, and as such their marginal coverage is *at least* 95 *%*. So the actual 95 *%* confidence intervals for each *π*
_{
j
} are likely to be narrower than the error bars, and these statements are made with at least 95 *%* confidence.

*M*=4 (and also at

*M*=5). So possibly some of the eruptions recorded as

*M*=4 eruptions are smaller than this, and should be excluded. To check the sensitivity of our results to this source of mis-recording, we repeated the analysis for eruptions with recorded

*M*≥4.1; the results are shown in Fig. 5. The uncertainty is larger, compared to Fig. 4, because 57 observations out of 184 have been dropped. But the qualitative conclusions are unchanged, and we are satisfied that our simple quantitative assessment of recording rates is robust to this source of mis-recording.

*M*≥5.0), and redid the calculation. Figure 6 shows the result: the number of M5+ eruptions in the database is low (42), and our uncertainty about the recording rates is large. We also found that using eruptions

*M*≥5.1 (24 eruptions) gave a similar confidence set, although the point estimates were different. According to Fig. 6, the recording rate is likely to be higher for M5+ than for M4+, but we do not think the uncertainty assessment in Fig. 6 is reliable.

The confidence set figures include an additional curve, the centroid of the 50 *%* confidence set, as an alternative point estimate to the ML estimate. This centroid will tend to be smoother than the ML estimate, which can get closer to the edges of *Ω*. But the difference between the two point estimates is mostly less than 10 percentage points for M4+, which is small in relation to the width of the confidence set. Overall, we hesitate to provide a point estimate for the recording rates from this dataset, given the volume of the 95 *%* confidence set.

## Conclusions

There is empirical and narrative evidence that the recording rate for volcanic eruptions was lower in the past that it is today. We would like to estimate historical under-recording rates, to better interpret the historical record, and, by extension, to better predict the future. Yet we have few beliefs about the quantitative features of the under- recording process. Hence in this paper we have adopted a non-parametric approach to estimate the recording rate as a function of time, in which the global large-eruption rate is treated as constant in time, and the only two constraints on the recording rates are that the global recording rate is now one (since 1980), and that the global recording rates are non-increasing, going backwards in time. Compared to the parametric approach used previously in the volcanological literature (Sections “Review of previous approaches” and “Our statistical non-parametric approach”), our approach has the advantage of being far more general, but may not reduce our uncertainty much below its current levels. This possibility is especially acute given our focus on large (M4+) stratovolcano eruptions, which occur at a global rate of only about one every two years.

Fortunately, the dataset from the LaMEVE database is sufficiently rich that we can derive a clear message about historical under-recording (Section “Results and discussion”), including that the global recording rate drops below 100 *%* even in the recent past (e.g. in the 1800 s), is below 50 % before 1600, and is below 20 % before 1100 (all with 95 *%* confidence). These statements, which are formally about the global recording rate, apply, approximately, to the average recording rate across volcanoes (Section “Definition of the recording rate”). We also find that the sequence of global recording rates appears to combine aspects of both of the proposed parametric models, being generally exponential in shape, but also having more abrupt changes, notably around 1550. We extended our approach to very large (M5+) eruptions, but we judged our results to be unreliable.

Our results are consistent with the narrative evidence in Siebert et al. (2010, pp. 31–34) and Brown et al. (2014). Siebert et al. (2010) attribute the step-change around 1550 to the Spanish and Portuguese explorations at the end of the 15th century, which “opened Latin America and much of the western Pacific to European record-keeping” (p. 31), and also to the development of printing.

Finally the methods we have used here may also be useful in assessing under-recording of other hazards, e.g. earthquakes and floods. Like volcanoes, these infrequent hazards suffer from a sparse record of larger events, that is further thinned by low recording rates in the past. Our methods and computer code should be easily adaptable to other applications, and we would be happy to discuss this further.

## Appendix: Computing a 95 % confidence set

Putting aside *π*
_{
k
}, which is equal to 1, the order statistics of *k*−1 independent uniform random quantities are uniformly distributed on *Ω*, see Cox and Hinkley (1974, Appendix 2). This result provides a (stochastic) method to generate a point uniformly in *Ω*: generate *k*−1 independent uniform random quantities, order them from smallest to largest, and then append 1.

^{ k−1}, in this case a Sobol sequence. This sequence, which does not suffer from random variability, will have better coverage of any large connected subset of

*Ω*, because it decreases the clumpiness of the points, relative to a random sample. Figure 7 illustrates this property in 2D.

- 1.
Generate an

*m*-point Sobol sequence in*k*−1 dimensions. Order the values in each point from smallest to largest, and then append 1. This is now a deterministic*m*-spoint space-filling design on*Ω*. - 2.
Compute the log-likelihood for each point using (7), denoted

*ℓ*_{ i }for point π^{(i)}. Let \(\hat {i} = \arg \max _{i}\ell _{i}\). - 3.
For some value

*c*, keep all points π^{(i)}for which \(\ell _{i} \geq \ell _{\hat {\imath }} - c\). That is, all points for which the log-likelihood is within*c*of the maximum log-likelihood. These points form a space-filling design in some subset of*Ω*, denoted \(\mathcal {C}\).

A Sobol sequence can be generated using the sobol function from the randtoolbox package in the statistical computing evironment R (R Core Team 2013).

*d*degrees of freedom. According to the standard asymptotics, setting

will imply that \(\mathcal {C}\) is a 95 *%* confidence set for π; see Cox (2006, ch. 6) or van der Vaart (1998, ch. 16). But the standard asymptotics do not hold here, because the sample size is small, and the space *Ω* is irregular. In a nutshell, the 95 *%* level set of the likelihood function is not elliptical, because it is clipped by the edges of *Ω*. Therefore *c* must be adjusted to correct for ‘level error’, which is where the actual coverage is not equal to the nominal coverage.

*c*

^{∗}, is found empirically, in a simulation experiment using synthetic datasets; see DiCiccio and Efron (1996, notably sec. 7) for bootstrap methods for approximating confidence sets. Set \(\hat {\boldsymbol {\pi }} := \boldsymbol {\pi }^{(\hat {i})}\), and generate

*n*independent samples x

^{(1)},…,x

^{(n)}from the model with \(\boldsymbol {\pi } \leftarrow \hat {\boldsymbol {\pi }}\). This involves, for each sample:

- 1.
Simulate

*λ*∼Gamma(*a, b*), - 2.
For

*i*=1,…,*k*, simulate \(x_{i} \sim \text {Poisson}(\Delta _{i} \cdot \hat {\pi }_{i} \cdot \lambda)\).

For each sample, evaluate the log-likelihood for all points π
^{(1)},…,π
^{(m)}, denote this *ℓ*
_{
ij
} := log*p*(x
^{(j)}|π
^{(i)}). Also identify the maximum log-likelihood value, \(\hat {\ell }_{j} \mathrel {:=}\nolinebreak \max _{i} \ell _{ij}\). The value *c*
^{∗} is chosen so that 95 *%* of the sample contain the ‘true’ value \(\hat {\boldsymbol {\pi }}\). In other words, *c*
^{∗} solves

where
is the indicator function. A simple rearrangement shows that *c*
^{∗} is the 95th percentile of the histogram of \(\{ \hat {\ell }_{j} - \ell _{\hat {i} j} \}_{j=1}^{m}\).

Denote the resulting level set based on *c*
^{∗} as \(\mathcal {C}^{*}\). This empirical approach ensures that the coverage of \(\mathcal {C}^{*}\) is very close to 95 *%* at \(\boldsymbol {\pi } = \hat {\boldsymbol {\pi }}\), and, one hopes, approximately 95 *%* in the region of *Ω* around \(\hat {\boldsymbol {\pi }}\). Thus the resulting \(\mathcal {C}^{*}\) is approximately an exact 95 *%* confidence set.

In our calculation we used *m*←10^{4} and *n*←1001, and the 95 *%* confidence set cut-off at \(\hat {\boldsymbol {\pi }}\) was found to be *c*
^{∗}=4.61. By way of contrast, the asymptotic cut-off from (11) is *c*=9.15. In order to check the coverage of this cut-off at plausible values other than \(\hat {\boldsymbol {\pi }}\), we chose ten other points at random from the 95 *%* confidence set, and evaluated the coverage of the confidence set at these π’s using the above values of *m*, *n*, and *c*
^{∗}; the coverages were 94,92,94,92,94,93,95,93,95,92 *%*, which is reassuring.

## Declarations

### Acknowledgements

This work was supported by the Natural Environment Research Council (NERC) funded Consortium on Risk in the Environment: Diagnostics, Integration, Benchmarking, Learning and Elicitation (CREDIBLE); grant number NE/J017450/1. Sparks acknowledges funding from the European Research Council in the VOLDIES project which led to the development of the LaMEVE database. Cashman acknowledges support from the AXA Research Fund. We would like to thank two anomymous Reviewers who helped improve the description of our approach and how it differs from past parametric approaches.

### Authors’ contributions

JR carried out the statistical analysis, SS and KC provided volcanological expertise. All authors read and approved the final manuscript.

### Competing interests

The authors declare they have no competing interests.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## Authors’ Affiliations

## References

- Brown, S, et al. Characterisation of the Quaternary eruption record: Analysis of the Large Magnitude Explosive Volcanic Eruptions (LaMEVE) database. J Appl Volcanol. 2014; 3(5). http://www.appliedvolc.com/content/3/1/5.
- Coles, S, Sparks R. In: Mader H, Coles S, Connor C, Connor L, (eds).Extreme value methods for modelling historical series of large volcanic magnitudes. Geological Society of London: Special Publication of IAVCEI; 2006, pp. 47–56.Google Scholar
- Cox, D. Principles of Statistical Inference. Oxford University Press. Oxford: UK; 2006.View ArticleGoogle Scholar
- Cox, D, Hinkley D. Statistics Theoretical. Chapman and Hall. London: UK; 1974.Google Scholar
- Crosweller, H, et al. Global database on Large Magnitude Explosive Volcanic Eruptions (LaMEVE). J Appl Volcanol. 2012; 1(4). http://www.appliedvolc.com/content/1/1/4.
- De la Cruz-Reyna, S. Poisson-distributed patterns of exposive eruptive activity. Bull Volcanol. 1991; 54:57–67.View ArticleGoogle Scholar
- Deligne, N, Coles S, Sparks R. Recurrence rates of large explosive volcanic eruptions. J Geophys Res. 2010; 115(B06):203.Google Scholar
- DiCiccio, T, Efron B. Bootstrap confidence intervals. Stat Sci. 1996; 11(3):189–212. with discussion and rejoinder, 212–228.View ArticleGoogle Scholar
- Furlan, C. Extreme value methods for modelling historical series of large volcanic magnitudes. Stat Modell. 2010; 10(2):113–32.View ArticleGoogle Scholar
- Guttorp, P, Thompson M. Estimating second-order properties of volcanicity from historical data. J Am Stat Assoc. 1991; 86:578–83.View ArticleGoogle Scholar
- Kiyosugi, K, Connor C, Sparks R, Crosweller H, Brown S, Siebert L, Wang T, Takarada S. How many explosive eruptions are missing from the geologic record? Analysis of the quaternary record of large magnitude explosive eruptions in Japan. J Appl Volcanol. 2015; 4(1). doi:10.1186/s13617-015-0035-9.
- Lunn, D, Jackson C, Best N, Thomas A, Spiegelhalter D. The BUGS Book: A Practical introduction to Bayesian Analysis. Boca Raton FL, USA: CRC Press; 2013.Google Scholar
- Mason, B, Pyle D, Oppenheimer C. The size and frequency of the largest explosive eruptions on earth. Bull Volcanol. 2004; 66(8):735–48.View ArticleGoogle Scholar
- Mead, S, Magill C. Determining change points in data completeness for the Holocene eruption record. Bull Volcanol. 2014; 76(14).Google Scholar
- Pickands, J. The two-dimensional Poisson process and extremal processes. J Appl Probab. 1971; 8:745–56.View ArticleGoogle Scholar
- R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria, ; 2013. http://www.R-project.org/.
- Siebert, L, Simkin T, Kimberly P. Volcanoes of the World. Berkeley and Los Angeles CA, USA: University of California Press; 2010.Google Scholar
- Simkin, T. Terrestrial volcanism in space and time. Ann Rev Earth Planetary Sci. 1993; 21:427–52.View ArticleGoogle Scholar
- Solow, A. Estimating record inclusion probability. Am Stat. 1993; 47(3):206–8.Google Scholar
- Solow, A. An Empirical Bayes analysis of volcanic eruptions. Math Geology. 2001; 33(1):95–102.View ArticleGoogle Scholar
- van der Vaart, A. Asymptotic Statistics. Cambridge, UK : Cambridge University Press; 1998.View ArticleGoogle Scholar
- Venables, W, Ripley B. Modern Applied Statistics with S, 4th edn. USA: Springer, New York NY; 2010.Google Scholar