# The metallicity measurement of early-type galaxies

Yu-Zhong Wu

###### Abstract

We use data for 6048 early-type galaxies (ETGs) from Galaxy Zoo 1 that have
been cross-matched with the catalog of the MPA-JHU emission-line measurements
for the Sloan Digital Sky Survey Data Release 7. We measure the metallicity of
these ETGs by excluding various ionization sources, and study other properties
as well. We use the optimal division line of W2–W3=2.5 as a diagnostic tool,
and for the first time derive metallicity measurements for 2218 ETGs. We find
that these ETGs actually are closer to H II regions as defined by Kauffmann et
al. in the Baldwin–Philips–Terevich diagram, and they display younger stellar
populations. We present a full mass– metallicity relation and find that most
ETGs have lower metallicities than star-forming galaxies (SFGs) at a given
galaxy stellar mass. We use five metallicity calibrators to check our results.
We find that these metallicity indicators (R23, O32, and O3S2) give consistent
results. We suggest that the remaining two metallicity calibrators, which
increase metallicity by N-enrichment, can be used to calibrate metallicities
for SFGs, but not to estimate the metallicities of ETGs.

###### Subject headings:

galaxies: abundances — galaxies: evolution — galkaxies: accretion

††affiliationtext: Key Laboratory of Optical Astronomy, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China

##  1\. INTRODUCTION

Early-type galaxies (ETGs) generally appear as red, passive/retired objects,
almost without gas and dust. lenticular and elliptical galaxies are ETGs;
recent studies have demonstrated that ETGs usually contain a multiphase
interstellar medium (ISM; Herpich et al. 2018) comprising neutral hydrogen
(e.g., Krumm & Salpeter 1979; Goudfrooij et al. 1994; Oosterloo et al. 2010;
Serra et al. 2012; Lagos et al. 2014; Woods & Gilfanov 2014), molecule gas
(Combes, Young & Bureau 2007; Davis et al. 2015), halos of hot gas (Sarzi et
al. 2013), and dust (Goudfrooij & de Jong 1995). Gas in the ISM tend to be
accompanied by low efficientcy star formation (SF).

Mounting evidence reveals recent or ongoing SF in some ETGs. Yi et al. (2005)
demonstrated that recent SF activities could be identified in approximately
$15\%$ of their sample of 160 ETGs. Using a sample of roughly 2100 Sloan
Digital Sky Surcey (SDSS) ETGs, Kaviraj et al. (2007) have shown that $\sim
30\%$ of their sample experience low-level SF, deriving $\sim 1-3\%$ of the
total galaxy stellar mass. This implies that low-level SF is common in ETGs.

Much attention has been focused on explaining SF in ETGs. One possible cause
is gas accretion introduced by minor mergers (Kaviraj et al. 2009; Thilker et
al. 2010). Geréb et al. (2016) posited that an extend ($\sim$ 60 kpc) stellar
stream provides direct evidence for gas accretion, showing that a merger event
in GASS 3505 happened in the recent past. Another possible cause of SF is
accretion from the intergalatic medium (IGM). With regard to the formation of
Hoag’s Object, Finkelman et al. (2011) suggested that the core and H I disks
of the Hoag’s Object were formed at different evolutianary phases, and the
forming of the disks was prolonged by ‘cold’ accretion of pristine gas from
the IGM.

Gas-phase metallicity is crucial for studying various aspects of the evolution
of galaxies. Several attempts to obtain metallicity calibrations for H II
regions of ETGs have been made. Zhang et al. (2017) studied the impact of the
diffuse ionized gas (DIG) on metallicity measurements at kiloparsec scales
with the data from Mapping Nearby Galaxies at APO survey. Kumari et al. (2019)
estimated the effect of DIG contamination on metallicity measurements, and
obtained metallicity calibrators for the DIG and low-ionization emission
region.

To date, metallicity measurements of ETGs have only been provided by Athey &
Bregman (2009), Annibali et al. (2010), Bresolin (2013), and Griffith, Martini
& Conroy (2019); they have obtained dozens of oxygen abandances from optical
emission lines. Here, we investigate the metallicity of ETGs with a large
sample, and study the properties of these ETGs. In section 2, we provide a
brief description of the ETG sample and our data. Excluding several ionization
sources, we derive the metallicity of ETGs and explore the properties of these
ETGs in Section 3. Finally, our results are summarized in Section 4.

##  2\. THE DATA

In this work, we utilize the data of the SDSS Seventh Data Release (DR7;
Abazajian et al. 2009). The SDSS DR7 provides the spectra of more than one
billion galaxies, covering a wavelength range from 3600 to 9200 Å, with mean
spectral resolution $R\sim 1800$ and assembled from $3^{\prime\prime}$
diameter fibres. The measurements of stellar masses, emission line fluxes, and
star formation rates (SFRs) can be obtained from the SDSS DR7 catalog of Max
Planck Institute for Astrophysics$-$John Hopkins University (MPA-JHU), which
is publically available 111https://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/; this
catalog supplies the spectra measurements of about 900,000 galaxies.

Figure 1.— BPT diagnostic diagram. The black and red dots show initial
composite ETGs selected; the red dots are final ETGs with metallicity
measurements. The dotted blue curve represents the Kauffmann et al. (2003)
semi-empirical lower limit for the SFGs; the green dashed curve on this
diagram is the theoretical “extreme starburst line” obtained by Kewley et al.
(2001) as an upper boundary for SFGs.

First, we choose galaxies at $0.04<z<0.12$ to avoid the bias of mass-
metallicity (MZ) relations from the aperture (Kewley et al. 2005). Also, the
aperture-covering fractions are required to be $>20\%$ for all galaxies, and
the paremeters are computed from the r-band Petrosian and fiber magnitudes.
For our sample we select galaxies with a signal-to-noise ratio (S/N)$>3$ for
H$\alpha$ and H$\beta$, and with S/N $>2$ for [O II]$\lambda\lambda$ 3227,
3229, and [N II]$\lambda 6584$. Moreover, because the SFR FLAG keyword
represents the measuring status of SFRs in the catalog, the keyword must be
zero. As a result, we have an initial sample of 140,589 galaxies.

We first take the galaxies with $n_{\rm Sersic}>2.5$ to be ETGs, and then base
our galaxy morphologies on Galaxy Zoo 1 (Lintott et al. 2008, 2011) to match
our sample with that found in Table 2 of Lintott et al. (2011). The Sérsic
index comes from the New York University Value-Added Galaxy Catalog
222http://sdss.physics.nyu.edu/vagc-dr7/vagc2/sersic/ (NYU-VAGC; Blanton et
al. 2005). We cross-match the NYU-VAGC with our initial sample within
$2^{\prime\prime}$, and obtain 133,101 galaxies. Then we choose the galaxies
with $n_{\rm Sersic}>2.5$, showing a sample of 42,916 galaxies. Next, we
utilize Galaxy Zoo 1 to select those galaxies with an elliptical probability
that is higher than 0.5 (Herpich et al. 2018), and match them with those found
in Table 2 of Lintott et al. (2011). This produces a sample of 16,623 ETGs.

Figure 2.— Diagnostic diagram of [O I]/H$\alpha$ and [O III]/H$\beta$, which
identifies ionization from the stellar population with SD progenitor
contribution proposed by Woods & Gilfanov (2014). Galaxies with [O
I]/H$\alpha>0.5$ are dominated by SD ionization, and exist outside of the
plot.

We also need to assess the star-formating properties of our ETG sample and
their probability of their contribution of an active galactic nucleus (AGN) on
the Baldwin-Philips-Terevich (BPT) diagram (Baldwin et al. 1981). Based on
several models, Kewley et al. (2001) and Kewley et al. (2006) constructed
extreme starburst lines, which are an upper limit on the emission-line
strengths in star-forming galaxies (SFGs). Galaxies that lie below the two
curves are dominated by star formation (Griffith et al. 2019). Therefore, we
consider only those galaxies that are located in the composite region on the
BPT diagram in this work, which means that our sample contains 6048 composite
ETGs (see the BPT diagram in Figure 1).

To further study the properties of ETGs, we utilize the Wide-field Infrared
Survey Explorer (WISE; Wright et al. 2010) catalog to explore their
properties. This survey provides coverage of the whole sky in four bands: W1,
W2, W3, and W4, with the central wavelengths of 3.4, 4.6, 12, and 22 $\mu m$.
The accurate position, four-band fluxes, four-band instrumental profile-fit
phometry magnitudes, and their instrucmental profile-fit photometry S/N ratios
are publicly available 333https://irsa.ipac.caltech.edu/cgi-bin/Gator/nph-dd.
We cross-match the ETG sample using the ALLWISE source catalog within
$2^{\prime\prime}$ and S/N$>3$ for W2 and W3, and obtain 4177 galaxy sample.

In the MPA-JHU catalog for the SDSS DR7, $M_{\star}$ and SFRs measurements
assumed a Kroupa (2001) initial mass function (IMF) and are corrected by a
Chabrier (2003) IMF. In SFGs, gas-phase oxygen abundances are estimated by
using abundance measurements of extragalactic H II. Because ETGs include many
other ionization sources, such as AGN activities, cosmic rays, shocks, and
old, hot stars, Griffith et al (2019) obtained the metallicities of three ETGs
by considering the influence of these ionization sources on the abundance
measurements. In addition, Brown et al. (2016) used a sample of about 200,000
SFG sample to estimate the performance of some abundance indicators, and
suggested that the O3N2 method of Pettini & Pagel (2004, PP04) is the most
perfect calibrator. In the higher stellar masses, the MZ relation calibrated
by the PP04-O3N2 indicator can display a consistency with the direct method
($T_{e}$ method) MZ relation (Andrews & Martini 2013). In addition, the O3N2
indicator adopts flux ratios of more neighboring lines than other indicators,
avoiding more substantial effects of atmospheric dispersion (Griffith et al.
2019). In this Letter, we utilize PP04-O3N2 as our abundance estimator.

Figure 3.— W2-W3 vs. W1-W2 color-color diagram for composite ETGs. The blue
horizontal dotted-dashed line represents the mid-IR criterion to select AGNs
proposed by Stern et al. (2012). The green solid lines display the “AGN” wedge
as suggested by Mateos et al. (2012). The red vertical dashed line describes
the best boundary between retired galaxies and ongoing SF galaxies. The black
dots and purple triangles are composite ETGs, and the former ones have no SF
or AGNs, while the latter ones exhibit AGN activities. The cyan diamonds are
our study sample of ETGs with SF. The red “$*$” signs are those ETGs with
metallicity measurements from Bregman (2009), Annibali et al. (2010), Bresolin
2013, and Griffith et al. (2019).

In addition to the metallicity calibrator of PP04-O3N2, in this Letter we also
use another five metallicity calibrators to calculate oxygen abundance of
ETGs; Tremonti et al. 2004 (T04), Jones et al. (2015, Jon15), Curti et al.
(2017, Curti17), Sànchez-Almeida et al. (2018, Sánch18), and Sanders et al.
(2018, Sander18).

##  3\. Results

###  3.1. The Metallicity Measurement of ETGs

The gas-phase oxygen abundance of the warm ISM can be estimated by metallicity
measurements of extragalatic H II regions and photoionization models. In ETGs,
many ionization sources, such as AGN activities, shocks, post-asymptotic giant
brance (PAGB) stars, cosmic rays. In this work, the AGN activities can be
excluded due to the fact that our sample consists of composite ETGs, dominated
by SF (Griffith et al. 2019), and we do not need to consider the AGN
photoionization. Regarding the shock excitation mechanism, the observationed
line fluxes can be influenced by shocks. Compared with the typical densities
and velocities suggested by Sparks et al. (1989), the shock energy is too low
by two orders of magnitude (Athey & Bregman 2009). In Griffith et al. (2019),
we can see that [N II]/H$\alpha$ ratio excited by the shock from a Small
Magellanic Cloud-like galaxy (Allen et al. 2008) is far lower (about one order
of magnitude) than the ratio of galaxies lying in composite or AGN regions.
From Figure 4 of Griffith et al. (2019), we can also see that the line flux
ratio excited by the cosmic ray or extra heat is far lower than that excited
by galaxies located in composite or AGN region on the BPT diagram, and
excludes the two ionization sources. In addition, because PAGB stars adopting
the excitation source can provide energy on the same order of magnitude as
weak AGNs, Belfiore et al. (2016) suggested that many galaxies, dominated by
SF, are misclassified as active ones (Griffith et al. 2019).

Figure 4.— Comparison of the MZ relations of composite ETGs for different
metallicity calibrators. The green dotted-dashed lines represent a polynomial
fit of data from 93,089 SFGs (see the text).

Because the probability of the [Ne III]$\lambda 3869$ photoionization
transition is far lower than one of the collisional cross-section, Athey &
Bregman (2009) suggested that the line is an important index of collisinal
excitation; they excluded NGC 4125 and NGC 2768 from their sample due to the
fact that they have clear [Ne III]$\lambda 3869$ fluxes. This demonstrates
that the two objects have more a complicated excitation mechanism than
photoionization. We exclude two objects having significant [Ne III]$\lambda
3869$ line, and obtain 4175 ETGs.

Griffith et al. (2019) introduced a dianostic diagram of [O I]$\lambda
6300$/H$\alpha$ versus [O III]$\lambda 5007$/H$\beta$ to exclude an ionization
mechanism, which is excited by single-degenerate (SD) SNe Ia progenitors
(Woods & Gilfanov 2014). In this model, Woods & Gilfanov (2014) employed
diagnostic tools, utilizing [O I], [N I] emission lines increasingly ionized
by high-temperature, to test if WD progenitors could be a dominant ionization
sourece. In Figure 2, we show our sample on the diagnostic diagram of [O
I]$\lambda 6300$/H$\alpha$ versus [O III]$\lambda 5007$/H$\beta$, obtaining
4097 ETGs.

Figure 3 presents the relation between W2-W3 and W1-W2 colors for composite
ETGs. The optimal division line (red dashed line), W2-W3=2.5, is used as the
best demarcation between galaxies with and without SF (Herpich et al. 2016).
Stern et al. (2012) proposed that W1-W2$>0.8$ of Figure 2 be the mid-IR
standard by which to choose AGNs, and these measurements reveal that nuclear
activity provides almost all the IR emission (Caccianiga et al. 2015). The
‘AGN wedge’ proposed by Mateos et al. (2012) is displayed by the green lines
in Figure 2.

In Figure 3, the black dots at the left of the red dashed line W2-W3=2.5 show
1845 composite ETGs without SF. The red ‘$*$’ signs in Figure 3 represent
those ETGs with metallicity measurements from Athey & Bregman (2009), Annibali
et al. (2010), Bresolin 2013, and Griffith et al. (2019). One of the galaxies
at the right of the red dashed line of Figure 3 is NGC 4694, which appears in
the SF region of the BPT diagram of Figure 4 of Griffith et al. (2019). This
galaxy should be an SF ETG. The purple triangles locate in the ‘AGN’ wedge
(green solid lines) proposed by Mateos et al. (2012) are used to display ETGs
with nuclear activities; here there are 24 such ETGs. The cyan diamonds, which
lie to the right of this red dashed line of Figure 3, are thought to be
composite ETGs with SF, because Herpich et al. (2016) found that galaxies with
mid-IR color W2-W3$>2.5$ denote objects with SF. Less than half of all
composite ETGs have SF, and we suggest that SF is a dominant excitation source
in these ETGs. Finally, we derive our sample of 2218 ETGs.

Figure 5.— Distribution of $D_{n}4000$ for our ETG sample with metallicity
measurements.

###  3.2. Properties of these ETGs

In Figure 1, we use the BPT diagnostic diagram to show the ETG samples. The
blue dotted curve represents the Kauffmann et al. (2003) semiempirical lower
boundary for SFGs, while the green dashed curve on this diagram show the
theoretical “extreme starburst line” obtained by Kewley et al. (2001) as an
upper boundary for SFGs. The black and red dots present initial composite
ETGs; the red dots display our final ETG sample with metallicity measurements.
From Figure 1, we can see that most of the red dots are close to the blue
dotted curve, and a majority of the black dots are close to the green dashed
curve. This indicates that the ETGs represented by red dots may be dominated
by SF.

Due to the limited sample size of ETGs, we were previously unable to
completely present their MZ relation. Here, we have a chance to present it
using a large ETG sample. Figure 4 shows the MZ relation for our composite
ETGs with the six metallicity calibrators. We find that the range of most
galaxy stellar masses is from log$(M_{\star}/M_{\sun})\sim 10.0$ to
log$(M_{\star}/M_{\sun})\sim 11.0$. We find that these MZ relations calibrated
by the six metallicity indicators do not show clearly positive correlations,
and that most ETGs lie below the green fit curves in Figure 4. Next we discuss
the properties of the MZ relation in the following two cases: the MZ relation
of ETGs that lie below and above the median MZ relation of SFGs.

From Figures 4(a), 4(c), 4(d), and 4(f), the MZ relations for ETGs lie below
the median MZ relations for SFGs. They employ, respectively, the metallicity
indicators of R23, O32, O3N2, and O3S2 to calibrate the ETG metallicities. In
Figure 4(a), the green dotted-dashed line describes a polynomial fit of data
(deriving 93089 SFGs from the MPA-JHU catalog for the SDSS DR7 with the same
sample selection method of Wu et al. (2019), except S/N $>2$ for [O
II]$\lambda\lambda$ 3227, 3229, and [N II]$\lambda 6584$.), which have median
values of 30 bins, 0.1 dex in mass, and include more than 100 galaxies. The MZ
relation does not clearly show a positive correlation, with the Spearman
coefficient r=0.17. Also, the MZ relations calibrated by metallicity
indicators of Jon15, PP04, and Curti17 do not present a positive correlation,
with the Spearman coefficient r=0.1, r=0.26, and 0.21, respectively. We can
see that most of the composite ETGs displayed in Figure 4(a) lie below the
green fit line at a given stellar mass. Figures 4(c), 4(d), and 4(f) also show
the same result, with a majority of ETGs lying below the median MZ relation of
SFGs. The middle panel of Figure 9 of Griffith et al. (2019) shows the MZ
relation of ETGs using the O3N2 metallicity calibration of Pettini & Pagel
(2004). We can see that most ETGs may be below the MZ relation of SFGs.
Combined the results of our MZ relation and that of Griffith et al. (2019), we
find that the MZ relation of these ETGs may be different from that of late-
type galaxies; this indicates that the metallicities of ETGs are likely to be
lower than thoss of late-type galaxies at a fixed stellar mass.

Figures 4(b) and 4(e) do not show a significantly positive correlation, with
the Spearman coefficient r=0.28 and r=0.31, respectively. In these two Figures
we find that the metallicity calibrators of Sánch18 and Sander18 do not show
the result of Figures 4(a), (c), (d), and (f); in addition, almost all fo
theETGs lying above the green fit line in adopting N2 metallicity indicator.
In Figure 4(b), The Sánch18 calibrator employs the N2 indicator ([N
II]$\lambda 6584$/$\rm H\alpha$) and Equation (3) of Sánchez-Almeida et al.
(2018). Their metallicity increases with increasing N2. Similarly, in Figure
4(e), we find that about $52\%$ of ETGs lie below the green fit curve. We
suggest that this may originate from the metallicity calibrator of Sander18,
and the metallicity increases with increasing indicator N2O2 ([N II]$\lambda
6584$/([O II]$\lambda\lambda 3727,3729$)). In the top panel of Figure 9 of
Griffith et al. (2019), we find that most ETGs from Athey & Bregman (2009),
Annibali et al. (2010), and Griffith et al. (2019) have higher, about 0.2 dex,
than SFGs in O3N2 metallicity calibrated from the PP04-N2 metallicity
indicator. In the bottom panel, we also find that these ETGs have higher about
0.1 dex in metallicity than SFGs using PP04-O3N2 indicator calibrated from the
KD02-N2O2 estimator. We find that their results are almost consistent with our
results shown in Figures 4(b) and 4(d). According to Coziol et al. (1999)
model, the nitrogen production occures almost ceaselessly in SFGs, but the
oxygen enrichment is only produced in massive star ecolution, and the two
productions proceed alternately (Wu & Zhang 2013). It is difficult to observe
oxygen production in ETGs, while nitrogen production is easier to observe
because ETGs tend to experience low-level SF (Yi et al. 2005; Kaviraj et al.
2007). Therefore, we suggest that these metallicity calibrators, wherein
metallicity increases with N-enrichment, for example, the Sánch18 and Sander18
metallicity indicators, can be used to successfully calibrate metallicities
for SFGs, but may not to estimate the metallicities of ETGs.

In Figure 5, we show the distribution of $D_{n}4000$ for our final ETG sample.
We find that the majority of ETGs with metallicity measurements have younger
stellar populations, and about $85\%$ of these ETGs display $D_{n}4000<1.6$.
These ETGs with W2-W3$>2.5$ clearly show SF (Herpich et al. 2016), and
approach the Kauffmann et al. (2003) semiempirical line; therefore, most of
them should have younger stellar populations.

##  4\. Summary

In this Letter we derive the observational data of 6048 ETGs, and cross-match
Galaxy Zoo 1 data with SDSS DR7 MPA-JHU emission-line measurements. We exclude
various ionization sources to explore the metallicity measurement of ETGs, and
investigate the properties of these ETGs. We summarize our main results as
follows:

1\. We utilize the optimal division line of W2-W3=2.5 as the diagnose tool,
which selects ETGs (W2-W3$>2.5$) with SF as our composite ETG sample. We
obtain the final sample of 2 218 ETGs, and can calculate their metallicities.

2\. We find that our ETGs tend to be located close to the semiempirical lower
limit for SFGs proposed by Kauffmann et al. (2003) in the BPT diagram.
Moreover, we find that these ETGs have younger stellar populations, and $85\%$
of the ETGs have $D_{n}4000<1.6$.

3\. The MZ relation for ETGs is shown, and we find that the metallicity of
these ETGs tends to be lower than SFGs at the same galaxy stellar mass. We use
five metallicity calibrators to check the result. We find that three of these
metallicity indicators (R23, O32, and O3S2) can give consistent results, while
the N2O2 and N2 indicators do not obtain the same result. We suggest that
these metallicity calibrators, wherein metallicity increases with N-rerichment
(for instance, Sánch18 and Sander18 metallicity indicators), can be used to
successfully calibrate metallicities for SFGs, but may not be able to estimate
the metallicities of ETGs.

Y.Z. Wu thanks the anonymous referee for valuable suggestions and comments
that improved the quality of this Letter. This work was supported by the
Natural Science Foundation of China (NSFC; No. 11703044).

## References

  * Abazajian et al. 2009 (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A. et al. 2009, ApJS, 182, 543
  * Allen et al. 2008 (2008) Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20
  * Andrews & Martini 2013 (2013) Andrews, B. H & Martini, P. 2013, ApJ, 765, 140
  * Annibali et al. 2010 (2010) Annibali, A., Bressan, A., Rampazzo, R., Zeilinger, W. W., Vega, O & Panuzzo, P. 2010, A&A, 519, 40
  * Athey & Bregman 2009 (2009) Athey, A. E & Bregman, J. N. 2009, ApJ, 696, 681
  * Baldwin et al. 1981 (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Belfiore et al. 2016 (2016) Belfiore, F., et al. 2016, MNRAS, 461, 3111
  * Blanton 2005 (2005) Blanton, M. R., et al. 2005, AJ, 129, 2562
  * Bresolin 2013 (2013) Bresolin, F. 2013, ApJL, 772, 23
  * Brown 2016 (2016) Brown, J. S., Martini, P., Andrews, B. H. 2016, MNRAS, 458, 1529
  * Caccianiga et al. 2015 (2015) Caccianiga, A., Antón, S., Ballo, L., Foschini, L., Maccacaro, T., et al. 2015, MNRAS, 451, 1795
  * Chabrier 2003 (2003) Chabrier, G. 2003, PASP, 115, 763
  * Combes et al. 2007 (2007) Combes, F., Young, L. M & Bureau, M. 2007, MNRAS, 377, 1795
  * Coziol et al. 1999 (1999) Coziol, R., Reyes, R. E. C., considère, S., Davoust, E., & Contini, T. 1999, A&A, 345, 733
  * Curti et al. 2017 (2017) Curti, M., Cresci, G., Mannucci, F., Marconi, A., Maiolino, R., & Esposito, S. 20107, MNRAS, 465, 1384
  * Davis et al. 2015 (2015) Davis, T. A., et al. 2015, MNRAS, 449, 3503
  * Finkelman et al. 2011 (2011) Finkelman, I., Moiseev, A., Brosch, N & Katkov, I. 2011, MNRAS, 418, 1834
  * Geréb et al. 2016 (2016) Geréb, K., Catinella, B., Cortese, L., Bekki, K., Moran, S. M & Schiminovich, D. 2016, MNRAS, 462, 382
  * Goudfrooij & de Jong 1995 (1995) Goudfrooij, P & de Jong, T. 1995, A&A, 298, 784
  * Goudfrooij et al. 1994 (1994) Goudfrooij, P., de Jong, T., Hansen, L., Norgaard-Nielsen, H. U. 1994, MNRAS, 271, 833
  * Griffith et al. 2019 (2019) Griffith, E., Martini, P & Conroy, C. 2019, MNRAS, 484, 562
  * Herpich et al. 2016 (2016) Herpich, F., Mateus, A., Stasińska, G., Cid Fernandes, R., & Vale Asari, N. 2016, MNRAS, 462, 1826
  * Herpich et al. 2018 (2018) Herpich, F., Stasińska, G., Mateus, A., Vale Asari, N., & Cid Fernandes, R. 2018, MNRAS, 481, 1774
  * Jones et al. 2015 (2015) Jones, T., Martin,C., & Cooper, M. C. 2015, ApJ, 813, 126
  * Kauffmann et al. 2003 (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
  * Kaviraj et al. 2009 (2009) Kaviraj, S., Sebastien, P., Khochfar, S., Silk, J & Kay, S. 2009, MNRAS, 394, 1713
  * Kaviraj et al. 2007 (2007) Kaviraj, S., Schawinski, K., Devriendt, J. E. G., Ferreras, I., Khochfar, S., Yoon, S.-J., et al. 2007, ApJS, 173, 619
  * Kewley et al. 2001 (2001) Kewley, L. J., Dopita, M. A., Sutherland, R., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
  * Kewley et al. 2006 (2006) Kewley, L. J., Grovers, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961
  * Kewley et al. 2005 (2005) Kewley, L. J., Jansen, R. A., & Geller, M. J. 2005, PASP, 117, 227
  * Kroupa 2001 (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Krumm & Salpeter 1979 (1979) Krumm, N & Salpeter, E. E. 1979, ApJ, 228, 64
  * Kumari et al. 2019 (2019) Kumari, N., Maiolino, R., Belfiore, F & Curti, M. 2019, MNRAS, 485, 367
  * Lagos et al. 2014 (2014) Lagos, C. d. P., Davis, T. A., Lacey, C. G., Zwaan, M. A., Baugh, C. M., Gonzalez-Perez, V., Padilla, N. D. 2014, MNRAS, 443, 1002
  * Lintott et al. 2011 (2011) Lintott, C. J., Schawinski, K., Bamford, S., Slosar, A., Land, K., et al. 2011, MNRAS, 410, 166
  * Lintott et al. 2011 (2011) Lintott, C. J., Schawinski, K., Slosar, A., Land, K., Bamford, S., et al. 2008, MNRAS, 389, 1179
  * Mateos et al. 2012 (2012) Mateos, S., Alonso-Herreor, A., Carrera, F. J., Blain, A., Watson, M. G. et al. 2012, MNRAS, 426, 3271
  * Oosterloo 2010 (2010) Oosterloo, T., et al. 2010, MNRAS, 409, 500
  * Pettini & Pagel. 2004 (2004) Pettini, M., & Pagel,B. E. J. 2004, MNRAS, 348, 59
  * Sánchez-Almeida et al. 2018 (2018) Sánchez-Almeida, J., et al. 2018, MNRAS, 476, 4765
  * Sanders et al. 2018 (2018) Sanders, R. L., Shapley, A. E; Mariska, K., et al. 2018, ApJ, 858, 99
  * Sarzi et al. 2013 (2013) Sarzi, M., Alatalo, K., Blitz, L., Bois, M., Bournaud, F., et al. 2013, MNRAS, 432, 1845
  * Serra et al. 2012 (2012) Serra, P., Oosterloo, T., Morganti, R., Alatalo, K., Blitz, L., et al. 2012, MNRAS, 422, 1835
  * Sparks et al. 1989 (1989) Sparks, W. B., Macchetto, F., & Golombek, D. 1989, ApJ, 345, 153
  * Stern et al. 2012 (2012) Stern, D., Assef, R. J., Benford, D. J., Blain, A., Cutri, R. et al. 2012, ApJ, 753, 30
  * Thilker et al. 2010 (2010) Thilker, D. A., Bianchi, L., Schiminovich, D., Gil De Paz, A., Seibert, M., et al. 2010, ApJL, 714, 171
  * Tremonti et al. 2004 (2004) Tremonti, C. A., et al. 2004, ApJ, 613, 898
  * Woods et al. 2014 (2014) Woods, T. E & Gilfanov, M. 2014, MNRAS, 439, 2351
  * Wright et al. 2010 (2010) Wright, Edward L., Eisenhardt, Peter R M., Mainzer, Amy K., Ressler, Michael E., Cutri, Roc M. et al. 2010, AJ, 140, 1868
  * Wu et al. 2013 (2013) Wu, Y.-Z., & Zhang, S.-N. 2013, MNRAS, 436, 934
  * Wu et al. 2019 (2019) Wu, Y.-Z., Zhang, W., & Zhao, Y.-H. 2019, MNRAS, 486, 5310
  * Yi et al. 2005 (2005) Yi, S. K., Yoon, S.-J., Kaviraj, J.-M., Deharveng, R. M., Rich, R. M., et al. 2005, ApJ, 619, 111
  * Zhang et al. 2017 (2017) Zhang, K., et al. MNRAS, 466, 3217


# Convolution of Scale Invariant Continuous Ranked Probability Scores for
Testing Experts’ Statistical Accuracy

Tina Nane  Department of Applied Mathematics, Delft University of Technology
Roger Cooke  Department of Applied Mathematics, Delft University of Technology
Resources for the Future (ret)

###### Abstract

Computable solutions for expectations of Continuous Ranked Probability Scores
are presented. After deriving a scale invariant version of these scores, a
closed form for the convolutions of scores is presented. This closed form
enables the testing experts’ statistical accuracy. Results are compared with
tests using a familiar Chi-square goodness of fit test using a recent data set
of 6,761 expert probabilistic forecasts for which true values are known.

Keywords expert judgment, scoring rules, continuous ranked probability score,
probability interval score, geometric probability, classical model,
overconfidence, location bias.

##  1 Introduction

Scoring rules were introduced by de Finetti in 1937 as tools for encouraging
honesty in eliciting subjective probabilities (De Finetti, 1937) and have been
further developed by many authors including Shuford et al. (1966); Savage
(1971); Murphy (1977); Brown (1974); DeGroot and Fienberg (1983); Gneiting and
Raftery (2007). The latter reference gives an extensive overview. An expert
receives a score as a function of his/her probability assessment and the
realization. The score is strictly proper if the expert maximizes (for
negatively sensed rules, minimizes) his/her expected score per item by, and
only by, stating his/her true belief. Using a result of Murphy (1977), DeGroot
and Fienberg (1983) gave an additive decomposition of strictly proper rules
into ‘calibration’ and ‘refinement’ terms, thereby replacing Murphy’s
‘resolution’ (refinement applies only to well-calibrated experts). In the case
of the logarithmic rule, refinement becomes Kullback-Leibler divergence of the
sample distribution of realizations. Some authors (Hersbach, 2000) adopt a
framework in which nature picks a distribution for an unknown quantity and
forecasters attempt to predict this distribution. Scoring rules for sets of
variables play a key role in the classical model (CM, Cooke (1991)) for
combining expert judgments.

Scoring rules for individual variables were not designed for evaluating or
combining experts and are not generally fit for that purpose. Indeed,
rewarding honesty is not the same as rewarding quality. A simple example
illustrates this difference: Consider $100$ fair coin tosses. An expert
assesses the probability of heads on each toss as $1/2$. With the standard
scoring rules, the score for the outcome heads is the same as their score for
tails on each toss. If the score for all $100$ assessments is a function of
their $100$ scores for the individual tosses, then their score for $100$
tosses is independent of the outcome sequence; the outcome of $100$ heads
receives the same score as $50$ heads and $50$ tails. Equal scores do not
imply equal quality.

Another example concerns the quadratic rule for ‘rain / no rain’ events. This
rule is positively sensed on $[-1,1]$ and assigns the quadratic score
$2r–r^{2}–(1–r)^{2}$ if rain occurs, where $r$ is the expert’s probability of
rain. Interchange $r$ and $(1–r)$ in case it does not rain. Consider 1000 next
day forecasts of rain by two experts. Suppose the experts bin their forecasts
as shown below (Cooke, 2014).

Probability of rain next day | $5\%$ | $15\%$ | $25\%$ | $35\%$ | $45\%$ | $55\%$ | $65\%$ | $75\%$ | $85\%$ | $95\%$ | Totals
---|---|---|---|---|---|---|---|---|---|---|---
expert 1 | assessed | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 1000
realized | 5 | 15 | 25 | 35 | 45 | 55 | 65 | 75 | 85 | 95 | 500
expert 2 | assessed | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 1000
realized | 0 | 0 | 0 | 0 | 0 | 100 | 100 | 100 | 100 | 100 | 500
Table 1: Binned subjective probability of rain assessments for the 1000 next
days by two experts.

Ten probability bins are considered, each associated with a forecast
probability of rain. The experts’ assessments are equally informative in the
sense that they assign the same probabilities to the same number of days.
Expert 1 is statistically perfectly accurate, that is, the empirical frequency
of actual rainy days from the assessed 100 days is identical with the
probability associated with each bin. Expert 2 is massively inaccurate
statistically. The sample distributions bear little resemblance to his/her
assessed probabilities $(5\%,\ldots,95\%)$. Expert 1 has an average quadratic
score of 0.67 and Expert 2 an average quadratic score of 0.84. Expert 2 gets a
better quadratic score because the higher resolution in the sample out-weighs
statistical inaccuracy. Such examples make it difficult to explain intuitively
what the scores 0.67 and 0.84 mean. For more discussion, see (Cooke, 1991,
2014). In the context of expert judgment we would like to reward both honesty
and quality with scoring rules that are intuitive and easily explained. This
requires numerical insight into the rules’ behavior.

The negatively sensed Probability Interval Score ($PIS$) and its related
Continuous Ranked Probability Scores ($CRPS$) have recently been applied to
COVID-19 probabilistic predictions (Ray et al., 2020). Gneiting and Raftery
(2007) note: “Applications of the $CRPS$ have been hampered by a lack of
readily computable solutions to the integral (1)” (see below).

This article presents computable solutions to this integral which then allow
us to study its trade offs between statistical accuracy and ”sharpness”. The
$CRPS$ can be thrown into a scale invariant form which offers significant
advantages. After reviewing the $PIS$ and $CRPS$, we introduce a re-
parametrization of $CRPS$ by transforming the realizations according to the
probability integral transformation of an expert’s assessed cumulative
distribution function ($CDF$). An expert with $CDF\,\,F_{X}$ for continuous
variable $X$ is scored not with respect to the realization $y$ but with
$F_{X}(y)$, the quantile of the distribution of $X$ realized by $y$. The
proposed CRPS transformation has several advantages:

  1. (i)

The transformed $CRPS$ becomes scale invariant.

  2. (ii)

The expert’s sampling distribution of transformed $CRPS$ can be expressed in
closed form

  3. (iii)

The density of convolutions of transformed $CRPS$ scores for independent
variables is available in closed form.

  4. (iv)

Transformed $CRPS$ can then be used to test the expert’s statistical accuracy
without recourse to an asymptotic distribution.

On the downside, $CRPS$ is insensitive to location bias. If the experts assess
only certain quantiles, a second downside is that continuous $CDFs$ must be
interpolated before applying $CDFs$.

After introducing the Probability Interval Score and the Continuous Ranked
Probability Score, computable examples of the latter are given. This motivates
a scale invariant version of CRPS to be used in testing experts’ ”statistical
accuracy”, a term denoting goodness-of-fit tests adapted to expert judgments.
The closed form of the convolution of scale invariant $CRPS$ scores is
introduced. Using a recently compiled database of expert judgments with
realizations (Cooke et al., 2021), the results of this test are compared with
the statistical accuracy score of $CM$. The scores are also compared with
regard to rewarding proximity of the medians (considered as point forecasts)
to the realizations. We conclude that the scale invariant $CRPS$ better
rewards proximity to the median while failing to punish location bias. The
closed form convolution also confers advantages with respect to $CM$.

##  2 Probability Interval Scores

Numerical insight into the behavior of these scores requires a bit of effort.
For the $(1-\alpha)$ uncertainty interval $[L,U]$, with upper (lower) bound
$U(L)$, the $PIS$ (negatively sensed) (Aitchison and Dunsmore, 1968) for
realization $y$ is

$(U-L)+\frac{2}{\alpha}\times[(L-y)_{+}+(y-U)_{+}]$

where $X_{+}=X$ if $X>0$ and $X_{+}=0$ otherwise. Note that $s=2/\alpha$ is
the slope of the overconfidence penalty for $y\notin[L,U]$. The length
$\|U-L\|$ is called the “sharpness”; small values reward concentrated
probability mass.

To better understand the characteristics of $PIS$, consider $Y\sim U[0,1]$ and
the $(1-\alpha)$ uncertainty interval $[L,U]$. Then

$E_{Y}[PIS(y)]=U-L+\frac{2}{\alpha}\int_{0}^{L}(L-x)\,\mathrm{d}x+\frac{2}{\alpha}\int_{U}^{1}(x-U)\,\mathrm{d}x=U-L+\frac{1}{\alpha}\left[L^{2}+(U-1)^{2}\right].$

For the central $0.9$ interval $[0.05,0.95]$, the expected PIS is $0.95$.
Suppose an expert prefers to give an $80\%$ interval$,[0.1,0.9]$, then the
expected score is $0.9$. This is better than $0.95$ because the prediction
interval is sharper. An expert seeking to optimize (i.e., minimize) his/her
expected score might take a central $2\%$ prediction interval $[0.49,0.51]$
with expected score of $0.51$. The way in which the $PIS$ trades
overconfidence for sharpness may strike some as counter-intuitive. For
example, an expert claiming that the degenerate interval $[0.5,0.5]$ has
$40\%$ probability of catching the realization would achieve an expected score
of $0.833$, better than the score of the $90\%$ central interval. The
sharpness of an interval of zero length outweighs the overconfidence of
claiming $40\%$ mass at the point $0.5$.

##  3 Continuous Ranked Probability Score

Consider $y$ an unknown scalar quantity of interest. Suppose $y$ has a true
forecast cumulative distribution function ($CDF$) $F_{Y}$, characterizing the
distribution of a random variable $Y$, which is not known. An expert provides
$CDF\,\,F_{X}$ which (s)he believes to be the distribution of $Y$. We assume
both $F_{Y}$ and $F_{X}$ are continuous and strictly increasing on their
support. The continuous ranked probability score (CRPS) is defined as (Brown,
1974)

$CRPS(F_{X},y)=\int_{-\infty}^{\infty}[F_{X}(x)-1_{\\{x\geq y\\}}]^{2}dx.$ (1)

Lower values indicate better performance. CRPS is known to be strictly proper
relative to a class of Borel probability measures with finite first moment
(Gneiting and Raftery, 2007). As mentioned in the introduction, lack of
readily computable solutions for (1) have restricted the use of CRPS score. To
understand the behavior of the CRPS score, let us consider $X\sim U[L,H]$,
with $0<L<H<1$. As we will show later, this particular choice of distribution
is relevant for the development of our proposed score. For $y\in[0,1]$:

$CRPS(F_{X},y)=\left\\{\begin{array}[]{c}\begin{array}[]{ccc}L-y+\frac{H-L}{3}&for\
0\leq y<L\vspace{0.5mm}\\\
\frac{(y-L)^{3}}{3(H-L)^{2}}-\frac{(y-H)^{3}}{3(H-L)^{2}}&for\
y\in[L,H]\vspace{0.5mm}\\\ y-H+\frac{H-L}{3}&for\ H<y\leq
1.\end{array}\end{array}\right.$ (2)

Figure 1 shows the CRPS score as a function of $y$, for different cases of $L$
and $H$. Cases when $y$ falls within and outside the $F_{X}$ support are
highlighted.

(a)  L=0.3, H=0.4

(b) L=0.6, H=0.7

(c)  L=0.3, H=0.7

Figure 1: CRPS score for $X\sim U[L,H]$ and $y\in[0,1]$. Differences for $y<L$
(red), for $y\in[L,U]$(blue), and $y>U$(green) are highlighted.

The expectation of the CRPS score, which may be infinite, is given by

$E_{Y}[CRPS(F_{X},Y)]=\int_{y}\int_{x}[F_{X}(x)–1_{\\{x\geq
y\\}}]^{2}\mathrm{d}x\,\mathrm{d}F_{Y}(y).$ (3)

We discuss some computable solutions for this expectation.

###  3.1 Computable solutions

Consider $Y\sim U[0,1]$ and an assessment of $Y$’s distribution by an expert
as that of random variable $X$, $X\sim U[0,H],\,0<H<1$. The expert thinks
values greater than $H$ are impossible, although these can in fact arise. The
expected CRPS is computed based on the distribution of $Y$. The $CDF$ of
$X,F(x)=x/H$, for $x\in[0,H]$ and $F_{X}(x)=1,\,$for $x\geq H$, along with the
survivor function of $X,S(x)=1-F(x)$ are shown in Figure 2 (see also Candille
and Talagrand (2005)).

Figure 2: Cumulative distribution and survivor function of an uniformly
distributed random variable on [0,H].

Then $E_{Y}[CRPS(F_{X},Y)]=\int_{0}^{1}\int_{0}^{1}[F(x)-1_{x\geq
y}]^{2}\mathrm{d}x\,\mathrm{d}y$ is computed in 2 steps:

A) For $y<H:$

$\int_{0}^{H}\left[\int_{0}^{y}\frac{x^{2}}{H^{2}}dx+\int_{y}^{H}\left(\frac{H-x}{H}\right)^{2}dx\right]dy=\int_{0}^{H}\left[\frac{y^{3}}{3H^{2}}+\int_{H-y}^{0}\frac{z^{2}}{H^{2}}(-dz)\right]dy=\frac{H^{2}}{6}.$

B) For $y>H:$

$\displaystyle\int_{H}^{1}\left(\int_{0}^{H}\frac{x^{2}}{H^{2}}dx+\int_{H}^{y}\,dx+\int_{y}^{1}0\,dx\right)dy$
$\displaystyle=\int_{H}^{1}\left(\frac{H}{3}+y-H\right)\,dy=\frac{H(1-H)}{3}+\int_{0}^{1-H}\,zdz$
$\displaystyle=\frac{H(1–H)}{3}+\frac{(1–H)^{2}}{2}.$

Therefore:

$E_{Y}[CRPS(F_{X},Yy)]=\frac{H^{2}}{6}+\frac{H(1–H)}{3}+\frac{(1–H)^{2}}{2}.$
(4)

As noted by Hersbach (2000), these results acquire a physical dimension. The
result for (A) is the score an expert with $X\sim U[0,H]$ expects, namely
$H^{2}/6$, which has the physical dimension of $H^{2}$. If $X$ is in meters
and changes to centimeters, the expected score increases by a factor $10^{4}$.

If $X\sim U[L,H]$ , for $0<L<H<1$, then the same method of calculation applies
mutatis mutandis. If $L=1-H$, with $H\geq 0.5$, then the contributions from
$x<y$ and $y<x$ are equal and we need only to double the contribution from
$x<y$. If $y<L$, the contribution from $x<y$ is zero. We compute

$\displaystyle\int_{L}^{H}\int_{L}^{y}F(x)^{2}\,dx\,dy+\int_{H}^{1}\int_{L}^{y}F(x)^{2}\,dx\,dy$
$\displaystyle=$
$\displaystyle\int_{L}^{H}\int_{L}^{y}\frac{(x–L)^{2}}{(H–L)^{2}}\,dx\,dy+\int_{H}^{1}\left[\int_{L}^{H}\frac{(x-L)^{2}}{(H-L)^{2}}\,dx+\int_{H}^{y}\,dx\right]\,dy$
$\displaystyle=$
$\displaystyle\int_{L}^{H}\frac{y^{3}}{3(H-L)^{2}}\,dx+\int_{H}^{1}\left[\frac{(H-L)}{3}+(y-H)\right]\,dy$
$\displaystyle=$
$\displaystyle\frac{(H-L)^{2}}{12}+\frac{(H-L)(1-H)}{3}+\frac{(1-H)^{2}}{2}.$

Adding the identical contribution from $y<x$ gives:

$E_{Y}[CRPS(F,Y)]=\frac{(H-L)^{2}}{6}+\frac{2(H-L)(1-H)}{3}+(1-H)^{2}.$ (5)

Again, the score inherits a physical dimension from $X$. Substituting $L=1-H$
in (5) we find that $E_{Y}[CRPS(F,Y)]$ for $X\sim U[1-H,H]$ is equal to
$E_{Y}[CRPS(F,Y)]$ for $X\sim U[0,H]$ from (4). This holds for any $0<L<H<1$,
with $L=1-H$. Hence, for $X\sim
U[0,0.7],\,\,E_{Y}[CRPS(F,Y)]=0.1966=E_{Y}[CRPS(\tilde{F},Y)]$, for
$\tilde{X}\sim U[0.3,0.7]$. By the same token, $X\sim U[0,\,0.5]$ yields the
same expected score of $1/4$ as $\tilde{X}$ with distribution $\delta(0.5)$
concentrated at $0.5$.

If $F_{i}\rightarrow F$, then $E_{Y}[CRPS(F_{i},y)]\rightarrow
E_{Y}[CRPS(F,y)]$, by the Helly-Bray theorem (Billingsley, 2013). It follows
that, if $F_{n}\rightarrow U[0,\,0.5]$ and
$\tilde{F}_{n}\rightarrow\delta(0.5)$, then for all $\varepsilon>0$ and for
all sufficiently large $n$,
$|E_{Y}[CPRS(F_{n},Y)]-E_{Y}[CPRS(\tilde{F}_{n},Y)]|<\varepsilon$. This
illustrates how the $CRPS$ compensates loss of statistical accuracy by a gain
in “sharpness”, and again illustrates that equal scores do not entail equal
quality.

Note that if the probabilistic forecast is the distribution of $Y$ (uniform
$[0,1]$), then the expected $CRPS$ score is $1/6$, by (4). For $0<H<1$,
$E_{Y}[CRPS(F,Y)]>\frac{1}{6}$, by (4) and for $0.5<H<1,\,L=1-H$,
$E_{Y}[CRPS(F,Y)]>\frac{1}{6}$, by (5). An expert would receive a better
(lower) expected score if their probabilistic forecast were equal to the
distribution of $Y$.

##  4 Test for statistical accuracy

We would like to test the hypothesis that $Y$ follows the expert’s assessed
distribution $F_{X}$. Applying the probability integral transformation, let
$U=F_{X}(X)$ and define $V=F_{X}(Y)$. Then $F_{U}(u)=u$. The hypothesis that
$F_{X}=F_{Y}$ is equivalent to the hypothesis

$H_{0}:\,\,\,V\sim\,U[0,\,1].$

In this case $CRPS$ can be written, for the realization $v$ and for $U$
uniformly distributed on $[0,1]$, as

$\displaystyle CRPS(F_{U},v)$
$\displaystyle=\int_{-\infty}^{\infty}[u-1_{\\{u\geq v\\}}]^{2}\mathrm{d}{u}$
$\displaystyle=\int_{-\infty}^{0}(0-0)^{2}\mathrm{d}{u}+\int_{0}^{v}u^{2}\mathrm{d}{u}+\int_{v}^{1}(u-1)^{2}\mathrm{d}{u}+\int_{1}^{\infty}(1-1)^{2}\mathrm{d}{u}$
$\displaystyle=\frac{v^{3}}{3}-\frac{(v-1)^{3}}{3}.$ (6)

The range of the $CRPS(F_{U},v)$ is $\left[\frac{1}{12},\frac{1}{3}\right]$,
for $v\in[0,1]$ and $U\sim U[0,1]$. The distribution of $CRPS$ is the
distribution of the random variable

$\frac{1}{3}\left[V^{3}-(V-1)^{3}\right]=\frac{1}{3}-V+V^{2},$

taking values in $[\frac{1}{12},\frac{1}{3}]\,\,$(lower values are better).
Under the null hypothesis, $V$ is uniform $[0,\,1]$. For fixed
$Q\in[\frac{1}{12},\frac{1}{3}]$, to find the probability that
$CRPS(F_{U},v)\leq Q$, under the null hypothesis, we find the roots of
$V^{2}-V+(\frac{1}{3}-Q)=0:$

$v_{1,2}=\frac{1\pm\sqrt{1-4(\frac{1}{3}-Q)}}{2}=\frac{1\pm\sqrt{4Q–\frac{1}{3}}}{2}.$

Collecting the mass between the two roots, we obtain the CDF

$P(CRPS(F_{U},v)\leq x)=\sqrt{4x-\frac{1}{3}},$

with density

$\frac{2}{\sqrt{4x-\frac{1}{3}}},\,\,\,x\in\left(\frac{1}{12},\,\frac{1}{3}\right).$
(7)

Figure 3 shows $CRPS(F_{U},v)$, for $v\in[0,1]$ and $U\sim U[0,1]$, together
with its $CDF$ and density under the null hypothesis $H_{0}$.

Figure 3: CRPS function, for $U\sim U[0,1]$ and $v\in[0,1]$, together with its
cumulative and density functions, under the null hypothesis $H_{0}$.

From Figure 3 it is evident that the score $CRPS(F_{U},v)$ is symmetric around
the value $v=0.5$. This is different from the behavior of $CRPS(F_{U},v)$
exhibited in Figure 1a or 1b, and it illustrates a feature of the scale
invariant version of $CRSP$.

From equation (7), we can easily compute

$\displaystyle\centering E_{V}\left[CRPS(F_{U},V)\right]\@add@centering$
$\displaystyle=$ $\displaystyle\frac{1}{6}$ $\displaystyle
E_{V}\left[CRPS^{2}(F_{U},V)\right]$ $\displaystyle=$
$\displaystyle\frac{1}{16}\left(\frac{1}{5}+\frac{1}{3}\right)$ $\displaystyle
Var_{V}\left(CRPS(F_{U},V)\right)$ $\displaystyle=$ $\displaystyle
0.005555\underline{5}$

It is handier to consider the following transformation

$Z(U,V)=4CRPS(F_{U},V)-\frac{1}{3}.$ (8)

Then $Z$ has $CDF$ and density

$F_{Z}(x)=\sqrt{x},\,\,\,\text{
and}\,\,\,\,\,f_{Z}(x)=\frac{1}{2\sqrt{x}},\,\,x\in[0,1].$

Note that $f_{Z}$ is the density of $U^{2}$, where $U\sim U[0,1]$. So far,
only one unknown scalar quantity of interest, and expert’s resulting $CRPS$
score, have been considered.

Suppose an expert provides uncertainty assessments for $n$ random variables.
The emerging question is how to aggregate the $CRPS$ scores of each of the $n$
variables? The transformation (8), and the observation that the density
$f_{Z}$ is the density of a squared uniform random variable are again handy.

If we assume the $n$ variables to be independent, then we need to consider
$Z_{1},\ldots,Z_{n}$ independent variables, each with density $f_{Z}$. For
these, we need to find the density of $Z^{(n)}=Z_{1}+\cdots+Z_{n}$. Or, in
terms of the squared uniform random variables, we need to find the density of
$S_{n}=U_{1}^{2}+\cdots+U_{n}^{2}$.

Weissman (2017) provides closed form distributions for $S_{n}$, for
$n=3,4,5,6,8,10,12$ and their graphical representations. A connection is also
made with a topic of geometrical probability, that is, finding the cumulative
distribution function of $S_{n}$, $P(S_{n}\leq s)$ is equivalent to finding
the volume of the intersection between the unit n-cube and the ball of radius
$\sqrt{s}$, in $\mathbb{R}^{n}$, when both are centered at the origin. In his
comment to Weissman (2017), Forrester (2018) observes that the more generic
volume problem posed by Xu (1996), of finding the volume of the intersection
of a cube and a ball in n-space has already been solved by B. Tibken and D.
Constales (Rousseau and Ruehr (1997)). Weissman (2017) reports that Constales’
solution to the volume problem involves a method based on Fourier series and
implies that, for general $n$,

$F_{n}(s)=P(S_{n}\leq
s)=\frac{1}{6}+\frac{s}{n}+\frac{1}{\pi}\textbf{Im}\sum_{k=1}^{\infty}\left[\left(\frac{C(2\sqrt{k/n})-iS(2\sqrt{k/n})}{2\sqrt{k/n}}\right)^{n}\frac{e^{2\pi
iks/n}}{k}\right],$ (9)

where $S(x)=\int_{0}^{x}sin(t^{2})\,\mathrm{d}t$ and
$C(x)=\int_{0}^{x}cos(t^{2})\,\mathrm{d}t$ denote the Fresnel integrals and Im
is the imaginary part of a complex number.

(a)

(b)

Figure 4: The exact distribution (red) and empirical distribution function of
the sum of squared simulated uniform observations (blue), for n=2 (left) and
n=10 (right).

Figure 4 graphically compares the above cumulative distribution for n=2 (left)
and n=10 (right) with the empirical distribution function of the corresponding
sum of squared uniform observations. 100 observations were sampled for both
empirical distribution functions. The cumulative distribution function was
implemented in `R`, by making use of functions implementing the Fresnel
integrals in the `pracma` package (Borchers and Borchers, 2022).

##  5 Expert data

In this paper, we use data from $49$ studies involving 526 experts assessing
in total $580$ calibration variables from their fields for which realizations
are known (four experts from the original data were dropped because they did
not assess all calibration variables in their respective panels). In total
there are 6,761 expert probabilistic forecasts of variables from their fields
for which true values are known. The data is described and referenced in
(Cooke et al., 2021). The supplementary information for that reference gives a
description of the Classical Model, whose relevant aspects are briefly
reviewed here.

Figure 5: Number of experts per number of assessed calibration variables.

The number of assessed calibration variables differ per study and Figure(5)
provides information about this. Experts assessed at least 7 and at most 21
calibration variables during all studies, and 10 calibration were used in 21
of the 49 studies.

##  6 Results of Statistical Tests

This section compares statistical accuracy ($SA$) tests based on the classical
model $(CM,\,$Cooke (1991)) versus based on $CPRS$. In the classical model,
$SA$ is measured as the probability of falsely rejecting the hypothesis that a
probabilistic assessor is statistically accurate. It is, in other words, the
$P$-value of rejection for this hypothesis. We hasten to add that $CM$ does
not test and reject expert hypotheses but, in compliance with proper scoring
rule theory for sets of assessments, uses this $P$-value to measure the degree
of correspondence between assessments and data in forming weighted
combinations of expert distributions.

In the data used for this analysis, expert assessments take the form of fixed
quantiles, $5\%,\,50\%,$, $95\%$, from the assessor’s subjective distribution
for a continuous unknown quantity. When $n$ true values for a number of such
quantities are observed, we compute the sample distribution $s$ of inter-
quantile relative frequencies and compare this with the theoretical inter-
quantile mass function $p=(0.05,0.45,0.45,0.05)$. The test statistic is
$2nI(s,p)$, where $I$ is the Shannon relative information (log likelihood
ratio) and $n$ is the number of calibration variables. Assuming that the
realizations are independently sampled from the assessor’s distributions, this
statistic is asymptotically chi-square distributed with degrees of freedom
equal the number of assessed quantiles. Thus, $SA$ is measured as
$1-F_{\chi^{2}}(2nI(s,p))$ where $F_{\chi^{2}}$ is the $CDF$ of a $\chi^{2}$
distribution with three degrees of freedom. High scores (near $1$) are good,
low scores (near $0$) mean it is unlikely that the divergence between $s$ and
$p$ should arise by chance. Note that $CM$ relies on an asymptotic
approximation which for a small number of observations is not very good
(Cooke, 2014). Simulations for ten calibration variables are provided in
(Hanea and Nane, 2021). It is deemed capable of detecting only large
differences in experts’ performances. It uses only the assessed quantiles and
does not rely on an interpolated $CDF$.

When applied to this expert judgment data, a test based on $CRPS$ must
interpolate an expert’s $CDF$. For this purpose we follow CM and adopt a
minimally informative distribution relative to a uniform support (chosen by
the analyst111In principle, any background measure supporting the experts’
quantiles and realization may be used, this data imposes the uniform
background for convenience) which complies with the expert’s quantile
specification. The resulting CDF can be found in Hanea and Nane (2021). $CM$
uses this interpolation for computing an expert’s informativeness, but not for
computing SA.

Consider $n$ observations of continuous variables assessed by an expert $e$.
The following procedure calculates $SA$ for the $CRPS$ statistic

  1. 1.

For each realization $y_{i}=1,\ldots,n$, compute $q_{e}^{i}=F_{i,e}(y_{i})$,
the quantile of $y_{i}$ in expert $e$’s $CDF$ $F_{i,e}$.

  2. 2.

The $SA$ hypothesis $H_{0}$ entails that these quantiles are independent
samples from a uniform distribution. Under this hypothesis,
$CRPS_{i}(e)=CRPS(F_{U},q_{e}^{i})$ can be computed from (6).

  3. 3.

For each $i=1,\ldots,n$, compute $z_{i}(e)=4CRPS_{i}(e)-1/3$, which is given
in (8)

  4. 4.

Compute $s(e)=F_{n}(\sum_{i=1}^{n}z_{i}(e))$, where $F_{n}$ is the exact
distribution of the sum of $n$ independent squared uniform variables, given in
(9).

Note that the procedure can be applied for continuous and invertible CDFs.
$CRPS$ uses an exact instead of an asymptotic distribution for the convolution
of these $CDF$s. From Figure 3 it is evident that the score $CRPS(v)$ for
value $v$ is symmetric around the value $0.5$. The distribution of the sum of
such variables is insensitive to location bias in the following sense: the
score for $2n$ observations of $0.4$ is the same as for $n$ observations of
$0.4$ and $n$ observations of $0.6$. Figure 3 also shows that $CRPS$ is also
insensitive to under-confidence: An expert whose probability transformed
realizations are all $0.5$ scores better than one for whom the hypothesis of
Section 4 that $F_{X}=F_{Y}$ holds. Under-confidence is rare with expert
judgment. Over-confidence, on the other hand, is not rare and the $CRPS$ score
is sensitive to over-confidence (see Figure 6).

Figure 6: Frequency of percentiles of realizations for $6761$ expert
probabilistic forecasts from 49 studies.

Figure 7 plots the SA scores for $526$ experts based on $CM$ and on $CRPS$.
Although the drift of the two scores is similar, there is substantial scatter.
$CRPS$’s log geomean $SA$ score is $-2.76$ while that of $CM$ is $-3.47$; in
this sense $CRPS$ is less severe.

Figure 7: Statistical accuracy of 526 experts with respect to CRPS (red) and
CM (blue). SA scores are ordered by CRPS SA scores.

We define an expert’s location bias as the absolute difference between the
percent of realizations above the medians and $50\%$. Location bias of 50%
means that all realizations are above or all realizations are below the
medians. Figure 8 black circles the $CM$ scores of those experts for whom the
location bias is greater or equal to $20\%$. $CM\,SA$ scores of experts for
whom the location bias is $50\%$ (all realizations were either below or all
above the medians) are red circled. The location bias of these circled experts
is missed by $CRPS$ and helps explain some of the downward scatter.

Figure 8: Statistical accuracy of 526 experts with CM and CRPS with location
biases circled. Location bias is the absolute difference between the
percentage of medians above the realizations and $50\%$. Black circles denote
$CM$ scores of those experts for whom the location bias is greater or equal to
$20\%$. Red circles denote location bias of $50\%$.

For 70 experts, the location bias is 0%. These are termed experts without
location bias (though of course there could be location bias in the lowest and
highest inter-quantile intervals). Figure 9 plots these $70$ $CM\,SA$ scores
against all the $CRPS\,SA$ scores. On this subset, $CRPS$’s log geomean $SA$
score is $-2.50$ while that of $CM$ is $-2.68$.

Figure 9: Statistical accuracy of 526 experts with CRPS and 70 experts with CM
without location bias.

The number of calibration variables assessed in the 49 panels ranges from $7$
to $21$, see Figure 5. This influences $CM\,SA$ in two ways, a larger number
(i) increases the accuracy of the $\chi^{2}$ approximation and (ii) tends to
lower SA scores of poorly calibrated experts as the test of statistical
accuracy has greater power. The first should decrease the differences between
$CM\,SA$ and $CRPS\,SA$, at least for location unbiased experts, whereas the
second enables a greater range of $CM\,SA$ scores and might therefore increase
the differences. A multiple regression of $\log(\frac{CRPS\,SA}{CM\,SA})$
against location bias and number of calibration variables explains $30\%$ of
the variance (adjusted $R^{2}$) and both explanatory variables have a
significant positive effect on the dependent variable. Thus, the influence of
(ii) exceeds that of (i). The Pearson correlations of the dependent variable
with location bias and with number of calibration variables are $0.41$ and
$0.29$ respectively. The correlation of the two explanatory variables is
$-0.17$. A more detailed analysis might better explain the differences in the
two $SA$ scores but at this point it appears that location bias is the major
factor.

##  7 MAPE

Statistical accuracy is not the only scoring variable of interest; the
proximity of the median to the realizations is also important. There are many
measures for such proximity (Morley et al., 2018; Gneiting et al., 2007), each
with benefits and drawbacks. Perhaps the most popular is the Mean Absolute
Percentage Error (MAPE), defined for forecasts $x_{i}$ and realization $r_{i}$
$\,i=1\dots n$ as

$\frac{1}{n}\sum_{i=1}^{n}|\frac{x_{i}-r_{i}}{r_{i}}|.$

This is evidently unstable for very small $r_{i}$. Instability arises on this
data set, as the largest $MAPE$ is over one million. The $MAPE$s of $326$
experts were less than $2$ and we focus on this subset. Figure 10 plots these
$MAPE$ scores (left axis) and also plots the corresponding values of $CM\,SA$
and $CRSP\,SA$ (right axis). Although not overwhelmingly clear in the figure,
the $CRSP\,SA$ scores tend to be higher than those of $CM\,SA$, especially for
very low $MAPE$s (see trend lines). The Spearman correlation of $CM\,SA$ and
$MAPE$ on this data subset is $-0.15$ while that of $CRSP\,SA$ and $MAPE$ is
$-0.26$. This results from the fact that $CRPS$ uses the (interpolated) $CDF$
whereas $CM$ is based on inter-quantile hit-rates. It is reasonable to expect
that weighing experts according to $CRPS$ scores will produce better $MAPE$
values for the combination of experts than $CM$. Other researchers (Flandoli
et al., 2011) have used likelihood scores based on interpolated $CDF$s and
achieved better $MAPE$ performance than with $CM$, but such scores are
notoriously improper. The great advantage of $CRPS$ in this regard is that it
is based on a strictly proper scoring rule.

Figure 10: 326 MAPE scores $<2$ (left axis) and SA for CM and CRPS (right
axis), with trend lines.

##  8 Conclusion

A scale invariant version of the Continuous Ranked Probability Score (CRPS)
has been developed and applied to an expert judgment data base involving $49$
studies with $526$ experts assessing in total $580$ calibration variables from
their fields for which realizations are known. The transformed CRPSc yields a
test for experts’ statistical accuracy which has the advantage of a closed
form solution without appeal to an asymptotic distribution. Compared to the
statistical accuracy test used in the classical model it has the advantage of
better rewarding proximity of a median point forecast to the realization.
Nonetheless, the CRPS is insensitive to location and under-confidence bias. A
future study will focus on combinations of experts’ judgments, comparing the
performance of $CRPS$ with other tests based on the Chi Square, the Kolmogorov
Smirnov and the Cramer Von Mises statistics.

## References

  * Aitchison and Dunsmore (1968) J. Aitchison and I. Dunsmore.  Linear-loss interval estimation of location and scale parameters.  _Biometrica_ , 55(1):141–148, 1968.
  * Billingsley (2013) P. Billingsley.  _Convergence of probability measures_.  John Wiley & Sons, 2013.
  * Borchers and Borchers (2022) H. W. Borchers and M. H. W. Borchers.  Package ‘pracma’.  _accessed on_ , 4, 2022.
  * Brown (1974) T. A. Brown.  Admissible scoring systems for continuous distributions.  _RAND Corporation_ , 1974.
  * Candille and Talagrand (2005) G. Candille and O. Talagrand.  Evaluation of probabilistic prediction systems for a scalar variable.  _Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography_ , 131(609):2131–2150, 2005.
  * Cooke (1991) R. M. Cooke.  _Experts in uncertainty: opinion and subjective probability in science_.  Oxford University Press on Demand, 1991.
  * Cooke (2014) R. M. Cooke.  Validating expert judgment with the classical model.  _Experts and consensus in social science_ , pages 191–212, 2014.
  * Cooke et al. (2021) R. M. Cooke, D. Marti, and T. Mazzuchi.  Expert forecasting with and without uncertainty quantification and weighting: What do the data say?  _International Journal of Forecasting_ , 37(1):378–387, 2021.
  * De Finetti (1937) B. De Finetti.  La prévision: ses lois logiques, ses sources subjectives.  In _Annales de l’institut Henri Poincaré_ , volume 7, pages 1–68, 1937.
  * DeGroot and Fienberg (1983) M. H. DeGroot and S. E. Fienberg.  The comparison and evaluation of forecasters.  _Journal of the Royal Statistical Society: Series D (The Statistician)_ , 32(1-2):12–22, 1983.
  * Flandoli et al. (2011) F. Flandoli, E. Giorgi, W. P. Aspinall, and A. Neri.  Comparison of a new expert elicitation model with the classical model, equal weights and single experts, using a cross-validation technique.  _Reliability Engineering & System Safety_, 96(10):1292–1310, 2011.
  * Forrester (2018) P. J. Forrester.  Comment on “sum of squares of uniform random variables” by i. weissman.  _Statistics & Probability Letters_, 142:118–122, 2018\.
  * Gneiting and Raftery (2007) T. Gneiting and A. E. Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477):359–378, 2007.
  * Gneiting et al. (2007) T. Gneiting, F. Balabdaoui, and A. E. Raftery.  Probabilistic forecasts, calibration and sharpness.  _J. R. Statist. Soc. B Part 2_ , 69:243–368, 2007.
  * Hanea and Nane (2021) A. M. Hanea and G. F. Nane.  An in-depth perspective on the classical model.  _Expert Judgement in Risk and Decision Analysis_ , pages 225–256, 2021.
  * Hersbach (2000) H. Hersbach.  Decomposition of the continuous ranked probability score for ensemble prediction systems.  _Weather and Forecasting_ , 15:559–570, 2000.
  * Morley et al. (2018) S. K. Morley, T. V. Brito, and D. T. Welling.  Measures of model performance based on the log accuracy ratio.  _Space Weather_ , 16(1):69–88, 2018.
  * Murphy (1977) A. H. Murphy.  The value of climatological, categorical and probabilistic forecasts in the cost-loss ratio situation.  _Monthly Weather Review_ , 105(7):803–816, 1977\.
  * Ray et al. (2020) E. L. Ray, N. Wattanachit, J. Niemi, A. H. Kanji, K. House, E. Y. Cramer, J. Bracher, A. Zheng, T. K. Yamana, X. Xiong, et al.  Ensemble forecasts of coronavirus disease 2019 (covid-19) in the us.  _MedRXiv_ , 2020.
  * Rousseau and Ruehr (1997) C. Rousseau and O. Ruehr.  Problems and solutions. subsection: The volume of the intersection of a cube and a ball in n-space. two solutions by bernd tibken and denis constales.  _SIAM Review_ , 39(4):779–786, 1997.
  * Savage (1971) L. J. Savage.  Elicitation of personal probabilities and expectations.  _Journal of the American Statistical Association_ , 66(336):783–801, 1971.
  * Shuford et al. (1966) E. H. Shuford, A. Albert, and H. Edward Massengill.  Admissible probability measurement procedures.  _Psychometrika_ , 31(2):125–145, 1966.
  * Weissman (2017) I. Weissman.  Sum of squares of uniform random variables.  _Statistics & Probability Letters_, 129:147–154, 2017\.
  * Xu (1996) L. Xu.  The volume of the intersection of a cube and a ball in n-space.  _SIAM Review_ , 38(4):669–669, 1996.