# Generalized maximum likelihood estimation of the mean of parameters of
mixtures, with applications to sampling.

Eitan<EMAIL_ADDRESS>[ Israel Census Bureau
of Statistics;  Israel Census Bureau of Statistics  Ya’acov Ritov
<EMAIL_ADDRESS>[ University of Michigan;  University of
Michigan.

###### Abstract

Let $f(y|\theta),\;\theta\in\Omega$ be a parametric family, $\eta(\theta)$ a
given function, and $G$ an unknown mixing distribution. It is desired to
estimate $E_{G}(\eta(\theta))\equiv\eta_{G}$ based on independent observations
$Y_{1},...,Y_{n}$, where $Y_{i}\sim f(y|\theta_{i})$, and $\theta_{i}\sim G$
are iid.

We explore the Generalized Maximum Likelihood Estimators (GMLE) for this
problem. Some basic properties and representations of those estimators are
shown. In particular we suggest a new perspective, of the weak convergence
result by Kiefer and Wolfowitz (1956), with implications to a corresponding
setup in which $\theta_{1},...,\theta_{n}$ are fixed parameters. We also
relate the above problem, of estimating $\eta_{G}$, to non-parametric
empirical Bayes estimation under a squared loss.

Applications of GMLE to sampling problems are presented. The performance of
the GMLE is demonstrated both in simulations and through a real data example.

Keywords: GMLE, Mixing distribution, Non-parametric Empirical Bayes, Sampling.

##  1 Introduction

Let $f_{\theta}(y)\equiv f(y|\theta),\;\theta\in\Omega,$ be a parametric
family of densities with respect to some measure $\mu$. Let $G$ be an unknown
mixture distribution. We observe $Y_{1},...,Y_{n}$, which are realizations of
the following process. Let $\theta_{i}\sim G$, $i=1,...,n,$ be independent,
where conditional on $\theta_{1},...,\theta_{n}$, $Y_{i}\sim f(y|\theta_{i})$
are independent. Given a function $\eta(\theta)$, consider the functional
$\eta^{*}$, where $\eta^{*}(G)=E_{G}\eta(\theta)$. Denote,

$\eta^{*}(G)=\eta_{G}.$

Based on our observations, it is desired to estimate $\eta_{G}$.

Consider for example the case $Y_{i}\sim N(\theta_{i},1)$, $\theta_{i}\sim G$
are iid, $i=1,...,n$. In the case $\eta(\theta)=E_{\theta}(Y)$, the natural
estimator for $\eta_{G}$ is $\frac{1}{n}\sum Y_{i}$, which is consistent and
efficient. Assume the support of $G$ is bounded away from zero, and consider
$\eta(\theta)=\frac{1}{\theta}$. In the later case, deriving a consistent
estimator for $\eta_{G}=E_{G}\frac{1}{\theta}$ is less obvious. Yet, the
Generalized Maximum Likelihood Estimator (GMLE), defined in the following,
yields a consistent estimator for $\eta_{G}=E_{G}\frac{1}{\theta}$, and more
generally.

Let $\hat{G}$ be a GMLE estimator for $G$, defined through (2). A GMLE
estimator of $\eta_{G}\equiv E_{G}\eta(\theta)$, is defined as:

$\hat{\eta}_{G}\equiv E_{\hat{G}}\eta(\theta)$ (1)

GMLE for ${G}{G}$. Given a family $\\{f(y|\theta),\;\theta\in\Omega\\}$ and
$G$ as above, define

$f_{G}(y)=\int f(y|\vartheta)dG(\vartheta).$

Given observations $Y_{1},...,Y_{n}$, a GMLE $\hat{G}$ for $G$ (Kiefer and
Wolfowitz, 1956) is defined as:

$\hat{G}=\mathop{\mathrm{argmax}}_{{G}}\;\Pi f_{{G}}(Y_{i}).$ (2)

Traditionally, a GMLE estimator $\hat{G}$ for a mixture $G$ is computed via
EM-algorithm, see the seminal paper, Laird (1978). Koenker and Mizera (2014)
suggested exploitation of convex optimization techniques.

Much of the research on GMLE is about characterization of the GMLE $\hat{G}$
in terms of the size of its support, see, e.g., Lindsay (1995). Asymptotics
and Efficiency of the estimator $\hat{\eta}_{G}$ for $\eta_{G}$ is studied in
semi-parametric theory, see e.g., Bickel, et.a.l (1995). One goal of this
paper is to present some new problems and applications where GMLE estimators
are useful. We elaborate on examples of ‘sampling with no response’ and on
‘post-stratification’.

We also relate GMLE estimators $\hat{\eta}_{G}$ to non-parametric empirical
Bayes, and provide useful and insightful representations. In non-parametric
empirical-Bayes the goal is to estimate the specific values
$\eta_{i}\equiv\eta(\theta_{i})$, $i=1,...,n$, under a given loss, mainly, as
in our case, a squared loss. A mixing distribution $G$ induces a joint
distribution of $(Y,\theta)$. We denote the corresponding conditional
expectation of $\theta$ conditional on $Y$, by $E_{G}(\theta|Y)$. The Non
Parametric Empirical Bayes approach is to estimate $\eta_{i}$, $i=1,...,n$ by
$\hat{\eta}_{i}=E_{\hat{G}}(\eta(\theta)|Y_{i})$, see, e.g., Jiang and Zhang
(2009), Koenker and Mizera (2014), and $g$-modeling in Efron(2014). In
subsection 3.2 we show that $\hat{\eta}_{G}$ may be represented as
$\hat{\eta}_{G}=\frac{1}{n}\sum
E_{\hat{G}}(\eta(\theta)|Y_{i})=\frac{1}{n}\hat{\eta}_{i}$. In the case of
$Y_{i}\sim N(\theta_{i},1)$ with $\eta(\theta)=E_{\theta}Y$, $\theta_{i}\sim
G$ are iid (and more generally) we show, that $\hat{\eta}_{G}\equiv
E_{\hat{G}}\eta(\theta)=\frac{1}{n}\sum Y_{i}$.

Asymptotics of GMLE estimators $\hat{\eta}_{G}$, mainly in terms of
consistency, is explored both theoretically and through simulations. There are
cases where $\hat{G}$ is not unique, e.g., in cases where $G$ is non-
identifiable. When $\hat{G}$ is not unique, we consider as GMLE any
$\hat{\eta}_{G}=E_{\hat{G}}\eta(\theta)$ that corresponds to any GMLE
$\hat{G}$. In such cases, as $n\rightarrow\infty$, often, not every sequence
of GMLE is consistent; we argue that GMLE estimators are still plausible and
worthwhile, and suggest ”‘GMLE related”’ Confidence-Intervals for $\eta_{G}$.

Given realizations $Y_{1},...,Y_{n}$ as above, a related problem, studied in
Zhang (2005), is the prediction of the realized $\sum\eta(\theta_{i})$, where
$\theta_{i}\sim G$ are iid. His approach is also related to GMLE. Zhang
elaborates on the difference between the estimation and the prediction
problems in terms of efficiency. He also presents interesting applications
where it is desired to estimate the realized $\sum\eta(\theta_{i})$.

###  1.1  More on motivation and formulation. Fixed parameters
$\theta_{1},...,\theta_{n}$.

Our problem may be formalized without appealing to a mixture distribution $G$.
Given independent $Y_{i}\sim f(y|\theta_{i}),\;i=1,...,n$ with $\theta_{i}$
fixed, one may want to estimate $\sum\eta(\theta_{i})$. Let $G^{n}$ be the
empirical distribution of $\theta_{1},...,\theta_{n}$. Let $\hat{G}^{n}$ be
the GMLE based on $Y_{1},...,Y_{n}$, as defined in (2). Under suitable
tringular array formulation and conditions, the sequence of signed measures
$(\hat{G}^{n}-G^{n})$, converges weakly to the zero-measure. However, a
reference to such a general result is not known to us, and we give a proof in
Section 3.3. This weak convergence result motivates estimating
$\frac{1}{n}\sum\eta(\theta_{i})=E_{G^{n}}\eta(\theta)=\eta_{G^{n}}$ by
$\;E_{\hat{G}^{n}}\eta(\theta).$

As in most of the literature, we often appeal to a mixing distribution $G$ and
random $\theta_{i}\sim G$, which makes the formulation more convenient. The
two approaches of appealing to a mixing distribution $G$, versus avoiding it,
are analogous to Empirical-Bayes versus Compound-Decision approaches. See,
e.g., Zhang (1997), Brown and Greenshtein (2009). It is also related to the
distinction between incidental and random nuisance parameters, see, e.g.,
Pfanzagl (1993).

In Section 2 we present some motivating examples. In Section 3 we present some
theoretical results concerning representations of $\hat{\eta}_{G}$, and some
asymptotics. In Section 4 we present simulations. In Section 5 we sketch a
‘GMLE related’ Confidence Interval for $\eta_{G}$. In Section 6 we present a
real data example.

##  2 Examples.

A main motivation for our study of GMLE are the following sampling models and
problems in the context of stratification.

In stratified sampling and in poststratification, it is desired to have a very
fine stratification. Then, conditional on a fine stratum it appeals that
missing observations are missing completely at random. Similarly in
observational study (or, in convenience sampling), fine stratification makes
the assumption that conditional on the strata the success or failure of a
treatment are independent of the, often unknown, allocation mechanism to
treatment and control. However, under very fine stratification many strata are
likely to be empty, i.e., with no corresponding observations. The GMLE method
is a natural way of handling the resulting difficulties of empty strata. In
practice, as the number of observations is increased, the stratification
typically becomes finer.

We elaborate on the above. Suppose that it is desired to estimate the
proportion of unemployed in the population. A sampled subject will respond to
the survey with unknown probability. The proportion of unemployed among the
responders in the survey is strictly a biased estimator of their proportion in
the population, since there is a strong correlation between employment status
and response. It is known that in the relevant surveys, neighborhoods (or,
more generally, strata) with higher unemployment rates have lower response
rates. We may hope, however, that if we consider small enough and homogenous
stratum, the willingness to answer and the answer will be, conditionally on
the stratum, practically independent. Thus, the missing responses would be
practically missed completely at random (MCAR) within each stratum. See, e.g.,
Little and Rubin (2002).

Let $K_{i}$ be the number of observations in Stratum $i$ (i.e., $K_{i}$ were
sampled and responded). Let $X_{i}$ be the number of unemployed among the
observations from Stratum $i$. Formally, we have $n$ strata with equal
weights. Our observations are $Y_{i}=(X_{i},K_{i})$, $i=1,\dots,n$, where the
conditional distribution of $X_{i}$ conditional on $K_{i}$ is
$B(K_{i},p_{i})$. It is desired to estimate the population’s proportion
$p=n^{-1}\sum p_{i}$.

Similar considerations apply in observational studies. Suppose that
conditional on a fine stratification, within each stratum the probability of
success or failure, under treatment or control, is independent of the unknown
allocation mechanism to treatment or control. Given a treatment, we thinks of
$X_{i}$, the number of successes of the treatment in stratum $i$, and $K_{i}$,
the number of times the treatment was applied in Stratum $i$. We observe
$Y_{i}=(X_{i},K_{i})$. One approach of analyzing such observational data is to
pair treatments based on their propensity score. Using our method we do not
need pairing, and we may handle even the extreme case, where one treatment was
applied in one subset of strata, while the other treatment was applied in
another, disjoint, subset of strata.

###  2.1 Sampling Models

We consider the following two realistic scenarios, where the sample sizes
$K_{i}$ are random. In light of Section 3.2, the following motivating models
apply to both setups of random and fixed $\theta_{1},...,\theta_{n}$. We use
the notations and formulation of ramdom $\theta_{i}$, $i=1,...,n$, for
convenience.

Model (i) _Stratified sampling with non-response._ A random sample of
$\kappa_{i}$ subjects from stratum $i$ is sampled. The probability of a random
subject from stratum $i$ to respond is $\pi_{i}\leq 1$. Thus, the number of
actual responses, $K_{i}$, $K_{i}\sim B(\kappa_{i},\pi_{i})$ is a random
variable.

Model (ii): _Post Stratification_. We assume that $K_{i}$ has a
$Poisson(\lambda_{i})$ distribution, $i=1,2,\dots,n$. This scenario is
reasonable whenever we have a convenience sample, or, very low response rate.

In Model (i), $\theta_{i}=(\theta_{i1},\theta_{i2})\equiv(\pi_{i},p_{i})$,
where $\pi_{i}$ is the probability of response, while $p_{i}$ is the
proportion (of, say, unemployed) in stratum $i$. Conditional on $\theta_{i}$,
$X_{i}\mid K_{i}\sim B(K_{i},p_{i})$ and $K_{i}\sim B(\kappa_{i},\pi_{i})$. In
Model (ii), $\theta_{i}=(\theta_{i1},\theta_{i2})\equiv(\lambda_{i},p_{i})$.
Given $\theta_{i}$, $X_{i}\mid K_{i}\sim B(K_{i},p_{i})$ as before, and
$K_{i}\sim Poisson(\lambda_{i})$.

Suppose we are interested in estimating $\eta_{G}=E_{G}\eta(\theta)$, for
$\eta(\theta_{i})=p_{i}$. If there are no empty strata, i.e., $K_{i}>0$,
$i=1,...,n$ then the obvious and naive estimator:

$\hat{p}_{N}=\frac{1}{n}\sum_{i=1}^{n}\frac{X_{i}}{K_{i}},$

is applicable. The problem with the above estimator is that many strata may be
empty, i.e., with corresponding $K_{i}=0$. When there are some empty strata, a
common ad-hoc approach is of ‘collapsing strata’ where after the data is
observed some empty strata are unified with non-empty ones. An extreme
collapsing is to a single stratum, which yields the Extreme Collapsing
estimator:

$\hat{p}_{EC}=\frac{\sum X_{i}}{\sum K_{i}},$

which is, in fact, desirable when $p_{1}=...=p_{n}$.

Another way of handling the difficulty of empty strata, is to assume that
strata are ‘missing at random’, that is $K_{i}$ is uncorrelated with $p_{i}$,
and strata with $K_{i}=0$ may be simply ignored.

Our approach, of applying GMLE, handles cases with many empty strata in a
natural way, which is not ad-hoc and does not rely heavily on missing at
random type of assumptions.

Finding a GMLE for a two dimensional distribution $G$, outside of the normal
setup, only recently appears in the literature, see, e.g., Gu and Koenker
(2017), Feng and Dicker (2018). One reason might be its recent popularity due
to Koenker and Mizera‘s computational methods. Our Model (ii) is presented in
Feng and Dicker, however, their motivation and context is very different than
ours, and is not related to sampling. In addition, they study the estimation
of the the individual $p_{i}$, $i=1,...,n$ via non parametric empirical Bayes,
while we study the estimation of $\eta_{G}=E_{G}\eta(\theta)$, where
$\eta(\theta_{i})=p_{i}$. We also provide consistency results in the
estimation of $\eta$.

Further sampling models. In Greenshtein and Itskov (2018), related scenarios
and sampling models are studied.

We first elaborate on their “censored-case”. In the censored-case, $n$
individuals are attempted to be interviewed, where there are at most
$\kappa^{0}$ attempts for each individual $i,\;i=1,...,n$. The probability of
a response from individual $i$ in any single attempt is $\pi_{i}$, and
attempts are independent. Let $\kappa_{i}\leq\kappa^{0}$ be the number of
interviewing attempts of individual $i$, then $\kappa_{i}$ is a truncated
Geometric variable. Let $Z_{i}$ be a 0-1 valued r.v., e.g., indicator of
unemployment, $p_{i}=P(Z_{i}=1)$. Suppose $\theta_{i}\equiv(\pi_{i},p_{i})$
are i.i.d distributed $G$. It is desired to estimate
$E_{G}Z_{i}=E_{G}(\eta(\theta))$, where $\eta(\theta_{i})=p_{i}$, for a
randomly sampled item $i$. The possible outcomes, of the attempted interviews
of individual $i$, are $(Z_{i},\kappa_{i})=(1,\kappa_{i}),\;(0,\kappa_{i})$,
$\kappa_{i}\leq\kappa^{0}$, and $NULL$, where $NULL$ indicates non-response in
$\kappa^{0}$ attempts. This scenario is within our setup. One may think of
individuals as strata of size one.

The corresponding ”‘truncated-case”’ in Greenshtein and Itskov (2018) is
derived from the above scenario, with the twist that we do not know about the
existence of items that did not respond. Hence, inference is based only on
responded items. Denote by $G^{t}$ the conditional distribution of $\theta$,
conditional on response. Then: $dG^{t}(\theta)\propto
P_{\theta}(\kappa_{i}\leq\kappa^{0})\times dG(\theta)$, i.e., $G^{t}$ and $G$
determine each other when $\pi_{i}>0$ w.p.1, similarly $\hat{G}$ and
$\hat{G}^{t}$. Let $\hat{G}^{t}$ be the GMLE for $G^{t}$, which is based only
on the (non-truncated) observations $\\{i|\kappa_{i}\leq\kappa^{0}\\}$. Then,
when $\pi_{i}>0$ w.p.1 under $\hat{G}$, we suggest the GMLE estimator

$\hat{\eta}=\int\eta(\theta)d\hat{G}(\theta)\equiv\int\frac{\eta(\theta)[P_{\theta}(\kappa_{i}\leq\kappa^{0})^{-1}]d\hat{G}^{t}(\theta)}{\int
P_{\theta}(\kappa_{i}\leq\kappa^{0})^{-1}d\hat{G}^{t}(\theta)}.$

##  3 Asymptotic Results and equivalent representations of GMLE estimators.

###  3.1 Consistency and Asymptotics in Models (i) and (ii).

####  3.1.1  Consistency in Model (ii)

We consider the function $\eta(\theta_{i})=p_{i}$, where
$\theta_{i}=(\lambda_{i},p_{i})$. The results are under the setup of
$\theta_{i}\sim G$, $i=1,...n$ are iid random variables. In light of
subsection 3.3, analogous results may be derived for fixed
$\theta_{1},...,\theta_{n}$.

In Model (ii) we have consistency in the estimation of $\eta_{G}$, when
$P_{G}(\lambda=0)\equiv G(\\{\lambda=0\\})=0$, as proved in the following
theorem.

We first remind the notion of identifiability of mixtures.

Identifiability. Let $f_{G}(y)=\int f_{\vartheta}(y)dG(\vartheta)$. Then, $G$
is identifiable if $f_{G}=f_{\tilde{G}}$ implies $G=\tilde{G}$.

It is convenient to re-parametrize the problem, as follows. Given $X_{i}\sim
B(K_{i},p_{i})$ conditional on $K_{i}$, and $K_{i}\sim Poisson(\lambda_{i})$.
Denote $W_{i1}\equiv X_{i}$, and $W_{i2}=K_{i}-X_{i}$. Then $W_{i1}$ and
$W_{i2}$ are independent Poissons conditional on $(\lambda_{i},p_{i})$, with
corresponding parameters $\xi_{1i}\equiv p_{i}\lambda_{i}$ and
$\xi_{2i}\equiv(1-p_{i})\lambda_{i}$.

######  Theorem 1.

Consider $G_{\xi}$ the distribution of $(\xi_{1},\xi_{2})$ as above. Let
$\hat{G}_{\xi}$ be the GMLE based on iid $(W_{i1},W_{i2}),\;i=1,...,n$ from
the mixture $G_{\xi}$. Then $\hat{G}_{\xi}$ converges weakly to $G_{\xi}$.

###### Proof.

The proof follows from Kiefer and Wolfowitz (1956). Checking the conditions is
standard, the identifiability condition is verified, e.g., by Karlis and
Xekalak (2005).  $\square$

######  Corollary 1.

In Model (ii) for any function $\eta(\theta)$, such that
$\eta(\theta)=\psi(\xi_{1},\xi_{2})$, for $\psi$ which is continuous and
bounded on the support of $(\xi_{1},\xi_{2})$ under $G$,
$E_{\hat{G}}\eta(\theta)\rightarrow E_{G}\eta(\theta)$.

In our setup
$\eta_{\epsilon}(\theta)=\psi_{\epsilon}(\xi_{1},\xi_{2})\equiv\frac{\xi_{1}}{\xi_{1}+\xi_{2}+\varepsilon}1\\{\xi_{1}+\xi_{2}>0\\}$
satisfies the requirements for any $\epsilon>0$. If $G(\\{\lambda=0\\})=0$
then $\lim_{\epsilon\rightarrow 0}E_{G}\eta_{\epsilon}=E_{G}\eta(\theta)$, and
$E_{\hat{G}}\eta(\theta)\rightarrow E_{G}\eta(\theta)$.

When confining ourselves to the class of distributions in $\Gamma=\\{G\mid
G(\\{\lambda=0\\})=0\\}$, we have consistency of the estimator
$E_{\hat{G}}\eta(\theta)$ for $E_{G}\eta(\theta)$ for any $G\in\Gamma$.
However, the convergence may be arbitrarily slow.

####  3.1.2  (In)Consistency in Model (i)

Again, we consider the function $\eta(\theta_{i})=p_{i}$. Now
$\theta_{i}=(\pi_{i},p_{i})$.

In Model (i) there is no consistency in estimating $G$ for every $G$, neither
a consistency in estimating $E_{G}(\eta(\theta))$. This is a result of the
non-identifiability, and it is demonstrated in the following example for
$\kappa_{i}\equiv 1$.

######  Example 1.

When $\kappa_{i}\equiv 1$ there are $M=3$ possible outcomes of the
$i^{\prime}th$ observation, we list them as: $X_{i}=1$, $X_{i}=0$, and
$K_{i}=0$; the corresponding probabilities are:
$(\pi_{i}p_{i},\pi_{i}(1-p_{i}),(1-\pi_{i}))$. Suppose the outcomes of $n$
realization have $n_{1}$ occurrences of $X_{i}=1$, $n_{2}$ occurrences of
$X_{i}=0$, and $n_{3}$ occurrences of $K_{i}=0$. Note $(n_{1},n_{2},n_{3})$ is
multinomial. Suppose $\frac{1}{n}(n_{1},n_{2},n_{3})=(0.25,0.25,0.5)$,
obviously the following $\hat{G}_{1}$ and $\hat{G}_{2}$ are both GMLE. Let
$\hat{G}_{1}$ be degenerate at $(\pi,p)=(0.5,0.5)$. Let $\hat{G}_{2}$ be the
distribution whose support is $(0,1),(1,0),(1,1)$ with corresponding
probabilities 0.5, 0.25, 0.25. Then, obviously both $\hat{G}_{1}$ and
$\hat{G}_{2}$ are GMLE, while $E_{\hat{G}_{1}}\eta(\theta)=0.5\neq
E_{\hat{G}_{2}}\eta(\theta)=0.75.$

###  3.2  Equivalent representations.

Given a function $\eta(\theta)$, suppose it is desired to estimate
$\eta_{G}=E_{G}\eta(\theta)$. Again, for conveniemce we consider the notations
of random pramaters $\theta_{1},...,\theta_{n}$, but, in light of Section 3.2
the results apply also to the setup of fixed $\theta_{1},...,\theta_{n}$. By
the following theorem, the estimator $\hat{\eta}_{G}=E_{\hat{G}}\eta(\theta)$
equals to the average of $E_{\hat{G}}(\eta(\theta)\mid Y_{i})$, $i=1,...,n$.
The appeal of this fact is that if we estimate the values of the individual
$\eta(\theta_{i})$ via non-parametric empirical Bayes under squared loss,
specifically by $E_{\hat{G}}(\eta(\theta)\mid Y_{i})$, there is a consistency
and agreement between the estimates of the individual parameters and the
estimate of their total, or, of their average. In Zhang (2005), the problem of
estimating random sums involving a latent variable is explored. One approach
in Zhang (2005) is to estimate the random sum, by the sum of the estimated
conditional expectations of the summands. This is analogous to estimate
$\sum_{i}\eta(\theta_{i})$ by $\sum_{i}E_{\hat{G}}(\eta(\theta)\mid Y_{i})$,
the last term equals to $nE_{\hat{G}}\eta(\theta)$ by the following theorem.

######  Theorem 2.

Assume $\eta(\theta)$ is a bounded function, then:

$\hat{\eta}_{G}\equiv
E_{\hat{G}}\eta(\theta)=\frac{1}{n}\sum_{i}E_{\hat{G}}(\eta(\theta)\mid
Y_{i}).$

###### Proof.

Let

$d\hat{G}_{t}(\theta)=\bigl{(}1+t(\eta(\theta)-\hat{\eta}_{G})\bigr{)}d\hat{G}(\theta).$

Since $\eta(\theta)$ is bounded, it follows that for $\varepsilon>0$ small
enough, ${\mathcal{G}}=\\{G_{t},\;t\in(-\varepsilon,\varepsilon)\\}$ is a
curve of cdf’s (i.e., $G_{t}$ is a positive measure with total mass 1).

Since $\hat{G}$ is a GMLE, it maximizes the likelihood within ${\mathcal{G}}$.
Hence:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{1}{n}\frac{d}{dt}\sum_{i}\log(\int
f(Y_{i}\mid\theta)d\hat{G}_{t}(\theta))\mid_{t=0}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i}\frac{\int\bigl{(}\eta(\theta)-\hat{\eta}_{G}\bigr{)}f(Y_{i}\mid\theta)d\hat{G}(\theta)}{\int
f(Y_{i}\mid\theta)d\hat{G}(\theta)}=\frac{1}{n}\sum_{i}\int\eta(\theta)d\hat{G}(\theta\mid
Y_{i})-\hat{\eta}_{G}.$

This concludes the proof.

$\square$

####  3.2.1  GMLE for the mean of exponential mixtures

Let $Y_{i}\sim N(\theta_{i},1),\;i=1,...,n$ be independent observations, the
obvious estimator for $\sum\theta_{i}$ is $\sum Y_{i}$. One may wonder about a
comparison between the obvious estimator $\sum Y_{i}$ and the estimator
$nE_{\hat{G}}\theta=\sum E_{\hat{G}}(\theta|Y_{i})$. This is especially in
light of the ‘wonders’ done by non-parametric empirical Bayes, in the point-
wise estimation of $\theta_{i}$ by $E_{\hat{G}}(\theta|Y_{i})$, see, e.g.,
Brown and Greenshtein (2009), Jiang and Zhang (2009), Koenker and Mizera
(2014). Similarly one may wonder about the analogous situation where
$Y_{i}\sim Poisson(\lambda_{i})$, and it is desired to estimate
$\sum\lambda_{i}$. In the following we will show that the two estimators are
in fact equal. Thus, GMLE does not improve over the obvious estimator, yet, it
does not harm.

In the following theorem the notations are for a one dimensional exponential
family, but it applies for a general exponential family.

######  Theorem 3.

Let $Y_{1},...,Y_{n}$, be independent observations, $Y_{i}\sim
f_{\theta_{i}}(y)$. Suppose that the density of $Y$ is of the form
$f_{\theta}(y)\equiv f(y|\theta)=\exp(\theta y-\psi(\theta))$,
$\theta\in\Omega$, with respect to some dominating measure $\mu$. Let
$\hat{G}$ be a GMLE supported on the interior of the parameter set $\Omega$.
Let $\eta(\theta)=E_{\theta}Y$. Then:

$nE_{\hat{G}}\eta(\theta)=\sum Y_{i}.$

###### Proof.

Given the observations $Y_{1},...,Y_{n}$, let $\hat{G}$ be a GMLE. Define the
translation $\hat{G}_{\Delta}(\;.\;)=\hat{G}(\;.-\Delta)$.

Note:

$\int\exp(\theta
y-\psi(\theta))d\hat{G}_{\Delta}(\theta)=\int\exp((\theta-\Delta)y-\psi(\theta-\Delta))d\hat{G}(\theta).$

Since $\hat{G}$ is GMLE,

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{d}{d\Delta}\sum\log(\int\exp(\theta
Y_{i}-\psi(\theta))d\hat{G}_{\Delta}(\theta))|_{\Delta=0}$ $\displaystyle=$
$\displaystyle\frac{d}{d\Delta}\sum\log(\int\exp((\theta-\Delta)Y_{i}-\psi(\theta-\Delta))d\hat{G}(\theta))|_{\Delta=0}$
$\displaystyle=$
$\displaystyle\sum\frac{\int(-Y_{i}+\psi^{\prime}(\theta))\exp(\theta
Y_{i}-\psi(\theta))d\hat{G}(\theta)}{\int\exp(\theta
Y_{i}-\psi(\theta))d\hat{G}(\theta)}$ $\displaystyle=$ $\displaystyle-\sum
Y_{i}+\sum E_{\hat{G}}(\psi^{\prime}(\theta)|Y_{i})$ $\displaystyle=$
$\displaystyle-\sum Y_{i}+nE_{\hat{G}}\eta(\theta)$

The last equality follows by Theorem 2, and since for exponential family
$\eta(\theta)=\psi^{\prime}(\theta)=E_{\theta}Y$.

$\square$

The above theorem implies that under an exponential family setup,
$E_{\hat{G}}\eta(\theta)$ is unique for $\eta(\theta)=E_{\theta}Y$ even if
$\hat{G}$ is not. See the following example.

######  Example 2.

Let $Y_{i}\sim Bernoulli(p_{i})$ be independent where $p_{i}\sim G$ are
independent, $i=1,...,n$. Suppose $Y_{1}=...=Y_{\frac{n}{2}}=0$ and
$Y_{\frac{n}{2}+1}=...=Y_{n}=1$. Then, obviously both $\hat{G}_{1}$ and
$\hat{G}_{2}$ are GMLE, where $\hat{G}_{1}$ has half of its mass at $p=0.1$
and the other half at $p=0.9$, while $\hat{G}_{2}$ is degenerate at $p=0.5$.
Obviously $\hat{G}_{1}\neq\hat{G}_{2}$, however,
$E_{\hat{G}_{1}}p=E_{\hat{G}_{2}}p=0.5$.

###  3.3  Weak convergence

The seminal paper of Kiefer and Wolfowitz (1956), appeals to a completely
unknown distribution $G$, assuming that $\theta_{i}\sim G$, $i=1,2,...$ are
iid. They prove weak convergence of $\hat{G}$ to $G$. It may be seen in the
introduction of that paper, that the authors are somewhat uncomfortable with
this (very weak) assumption.

The weak convergence in KW (1956) is explored from a broader perspective in
Chen, J. H. (2017).

In this subsection we attempt to present a weak convergence result that does
not depend on the assumption that $\theta_{i},\;i=1,2,...$ are iid, and we
consider fixed $\theta_{1},...,\theta_{n}$. We consider a general triangular
array formulation with $(\theta_{1}^{n},...,\theta_{n}^{n})$, as also
elaborated in the sequel.

Recently there are results about rates of convergence of $\hat{G}$ and other
estimators to $G$. See, e.g., Saha and Guntuboyina (2020) for the multivariate
Gaussian case, Philippe and Kahn (2018) for finite mixtures. Those papers do
not cover the weak convergence for triangular arrays presented in this
subsection. In Philippe and Kahn (2018) a non-GMLE estimator, $\tilde{G}$ for
$G$, is studied and they provide minmax rates for the Wasserstein distance
between $\tilde{G}$ and $G$, for mixtures with $m$ components. However, those
rates do not imply our weak convergence result in the triangular array
formulation. Specifically, in our triangular array setup, as the number of
observations $n$ increases, the number of components $m$ in the mixture is
also increased since we allow $m=n$.

Assumptions

A1) The parameter space $\Omega$ is compact.

Let

$h(y)=\sup_{\theta}|\log(f(y|\theta))|,$

A2) $\sup_{\theta\in\Omega}E_{\theta}h(Y)^{4}<\infty$.

A3) The functions $f(y_{0}|\theta)$ are continuous in $\theta$ for every
$y_{0}$, and uniformly bounded.

A4) The class of densities $\\{f_{G}(y)\\}$, that correspond to the set of all
possible mixtures $\\{G\\}$, is identifiable.

We now define a truncation of $\log(f_{G}(Y))$, given a constant $M>0$.

$\displaystyle\log((f_{G}^{(M)}(Y))$ $\displaystyle=$
$\displaystyle\log(f_{G}(Y))\times I(|\log(f_{g}(Y))|<M)$ $\displaystyle+$
$\displaystyle M*(sign(\log(f_{G}(Y)))\times I(|\log(f_{G}(Y))|\geq M),$

here $I$ is an indicator function.

A5) For every $G_{0}\in\\{G\\}$, every $\epsilon>0$, and every $M>0$, the set

$\\{G\;|\;||\log(f_{G}^{(M)})-\log(f_{G_{0}}^{(M)})||_{\infty}<\epsilon\\},$

is an open set under the weak convergence topology.

The purpose of A5 is to ensure an open cover for $\\{G\\}$, from which a
finite covering subset may be extracted by compactness.

A5 is implied under A1-A4, if in addition it is assumed that, for every $M>0$,
the functions $\log((f_{\theta}^{(M)}(y))$, $\theta\in\Omega$ are uniformly
continuous.

The setup in this subsection is of a triangular array, as described in the
following. Given any triangular array sequence of sets
$\\{\theta_{1}^{n},...,\theta_{n}^{n}\\}$, $n=1,2,...$, denote by $G^{n}$ its
corresponding empirical distribution. Given independent observations
$Y_{1}^{n},...,Y_{n}^{n}$, $Y_{i}^{n}\sim f(y|\theta_{i}^{n})$, we denote the
corresponding GMLE by $\hat{G}^{n}$, as before. In the sequel we omit the
superscript $n$ and write $\theta_{i}$, $Y_{i}$.

Note that, in our current setup the joint likelihood
$f_{G}^{n}(Y_{1},...,Y_{n})$ does not satisfy
$f_{G}^{n}(Y_{1},...,Y_{n})=\Pi_{i}f_{G}(Y_{i})$ (since $Y_{i}$ are not
independent). We introduce new notations. Denote by $G^{n}_{1}$ the joint
distribution of $Y_{1},...,Y_{n}$ when the parameters are sampled
independently from $G^{n}$. The later is the setup in the previous sections;
we may think of it as sampling the parameters with replacement from the set
$\\{\theta_{1},...,\theta_{n}\\}$. We denote by $G^{n}_{2}$ the joint
distribution of $Y_{1},...,Y_{n}$ under the current setup, where sampling is
without replacement. Note that, for any function $\psi$,

$E_{G_{1}^{n}}\sum_{i}\psi(Y_{i})=E_{G_{2}^{n}}\sum_{i}\psi(Y_{i}).$ (3)

We keep the notations of previous sections and for $\theta\sim G$, $Y\sim
f(y|\theta)$, and a function $\psi$, we write $E_{G}\psi(Y)$.

######  Theorem 4.

Assume A1-A5, then for any bounded and continuous function $f$, the sequence
$a_{n}=(\int fd\hat{G}^{n}-\int fdG^{n})$, $n=1,2,...$ converges to 0 with
probability 1.

###### Proof.

By A1 we trivially have tightness, and hence any sequence $G^{n}$ of measures,
has a weakly converging sub-sequence. Our plan is to show that every sub-
sequence of $a_{n}$ has a further sub-sequence that converges to zero w.p.1.
The later implies the assertion of the theorem.

Given any sub-sequence of $a_{n}$, take a corresponding further sub-sequence
$n_{j}$, so that $G^{n_{j}}\Rightarrow_{w}G^{0}$ and
$\hat{G}^{n_{j}}\Rightarrow_{w}\hat{G}^{0}$, for some $\hat{G}^{0}$ and
$G^{0}$.

By Lemma 1 in the Appendix, under both sequences of joint measures
$G^{n_{j}}_{1}$ and $G^{n_{j}}_{2}$,

$\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\log(f_{\hat{G}^{n_{j}}}(Y_{i}))\rightarrow_{a.s.}E_{G^{0}}\log(f_{\hat{G}^{0}}(Y)).$
(4)

Note, the above is an almost sure convergence for a triangular array.

Similarly, under both sequences of joint measures $G^{n_{j}}_{2}$, and
$G^{n_{j}}_{1}$

$\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\log(f_{{G}^{n_{j}}}(Y_{i}))\rightarrow_{a.s.}E_{G^{0}}\log(f_{G^{0}}(Y)).$
(5)

By definition of $\hat{G}^{n_{j}}$,
$\frac{1}{n_{j}}\sum_{i}\log(f_{\hat{G}^{n_{j}}}(Y_{i}))\geq\frac{1}{n_{j}}\sum_{i}\log(f_{{G}^{n_{j}}}(Y_{i}))$,
thus under the sequence of measures $G^{n_{j}}_{2}$:

$E_{G^{0}}\log(f_{\hat{G}^{0}}(Y))=\lim\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\log(f_{\hat{G}^{n_{j}}}(Y_{i}))\geq\lim\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\log(f_{G^{n_{j}}}(Y_{i}))=E_{G^{0}}\log(f_{G^{0}}(Y)),$
(6)

Since $f_{G^{0}}=\mathop{\mathrm{argmax}}_{f_{G}}E_{G^{0}}\log(f_{G}(Y))$, we
obtain that:

$E_{G^{0}}\log(f_{\hat{G}^{0}}(Y))\leq E_{G^{0}}\log(f_{{G}^{0}}(Y)).$ (7)

From the last two inequalities we obtain:

$E_{G^{0}}\log(f_{G^{0}}(Y))=E_{G^{0}}\log(f_{\hat{G}^{0}}(Y)).$

By the identifiability assumption A4, and by concavity of the $\log$ function,
$\mathop{\mathrm{argmax}}_{\\{f_{G}\\}}E_{G^{0}}\log(f_{G}(Y))$ is unique, and
$\hat{G}^{0}=G^{0}$.

This concludes our proof.

$\square$

Remark. In compound decision, where $\theta_{1},...,\theta_{n}$ are fixed, it
is desired to estimate $\eta_{i}=\eta(\theta_{i})$. Then, the goal is to
approximate the optimal sperable estimator, some times also termed optimal
simple symmetric estimator. Under a squared loss, the later estimator for
$\eta_{i}$, given an observation $Y_{i}$, is
$E_{G^{n}}(\eta(\theta_{i})|Y_{i})$, see, e.g., Zhang (1997), Brown and
Greenshtein (2009). In our triangular array setup, our weak convergence
result, implies under suitable conditions
$[E_{G^{n}}(\eta(\theta_{i})|Y_{i})-E_{\hat{G}^{n}}(\eta(\theta_{i})|Y_{i})]\rightarrow
0$. A simple such possible additional condition in our triangular array setup,
is that $\liminf f_{G^{n}}(Y_{i}^{n})>0$.

##  4 Simulations.

In this section, we simulate the estimation of $\eta_{G}=E_{G}\eta(\theta)$,
in the case $\eta(\theta_{i})=p_{i}$, under models (i) and (ii). In all of the
following simulations, for convenience, $\eta_{G}$ equals 0.5. For any
parameter configuration, the number of repetitions is 50.

We compute the GMLE $\hat{G}$ via EM algorithm on a grid. As suggested by
Koenker and Mizera(2014), we search for a ”‘confined GMLE”’ where we confine
our search to distributions with a finite and specific support. Our choice of
a specific support is ad-hoc, see some rigorous treatment in Dicker and Zhao
(2014). The search for GMLE among distributions on the specific grid via EM-
algorithm, is also the approach in Feng and Dicker (2018), see further
elaboration there.

The grids for $(\theta_{i1},\theta_{i2})$ contain $40\times 40=1600$ equally
spaced grid points in a range that fits the relevant problem. The
parametrization in Model (ii) is via a two dimensional Poisson, as explained
in Section 3.1.1. The EM algorithm started with $G$ uniform on the grid and
ran for 1000 iterations.

###  4.1  Poisson sample sizes.

Table 1 presents the results of Model (ii) with 2-points support distribution
$G$. There are 500 strata corresponding to each of the two points in the
support of $G$. The table presents the mean and the sample standard deviation
of the naive and the GMLE estimators for $\eta_{G}$.

Table 1: Poisson Simulation discrete $G$. The estimated mean and standard deviations o two estimators.  Support points of $(\lambda,p))$ | The naive estimator | The GMLE
---|---|---
$(2,0.4),(1,0.6)$ |  0.486, (0.014) |  0.503, (0.020)
$(2,0.2),(1,0.8)$ |  0.453, (0.015) |  0.496, (0.018)
$(2,0.2),(0.5,0.8)$ |  0.385, (0.013) |  0.505, (0.022)

iIn Table 2 we report on one such set of simulations. In the three cases
presented in Table 2 there are again 500 strata of two types. The probability
$p^{I}$ of first 500 strata is fixed while for the rest of the strata the
probability is $p^{II}=(1-p^{I})$. The Poisson parameter $\lambda$ is
continuous in all the three cases and is chosen from a uniform distribution,
$U(0.5,1)$ in the first 500 strata and from $U(0.5,2)$ in each of the
remaining 500 strata.

Table 2: Poisson simulation with continuous $G$. The distribution of $\lambda$ is $U(0.5,1)$, for the first 500 strata, and $U(0.5,2)$ for the remaining strata. The binary probabilities $p^{I}$ corresponding to the first 500 strata are given the table. The binary probabilities of in the rest of the strata are $1-p^{I}$. The table gives estimates of the mean and standard deviation of the two estimators.  $p^{I}$ | Naive | GMLE
---|---|---
$0.4$ |  0.513, (0.015) |  0.500, (0.023)
$0.3$ |  0.529, (0.017) |  0.501, (0.027)
$0.2$ |  0.538, (0.016) |  0.491, (0.027)

###  4.2  Binomial sample sizes.

We next study Model (i), where $K_{i}$, the realized sample size from stratum
$i$, is distributed $B(\kappa_{i},\pi_{i})$. Again, our simulated populations
have two types of strata, 500 of each type. In the simulations reported in
Table 3, $\kappa_{i}=4$ for all of the 1000 strata, and $\eta_{G}=0.5$
throughout the three simulations summarized in the table. In 500 strata
$\pi_{i}=p_{i}=0.5-\delta$ while in the other 500 strata,
$\pi_{i}=p_{i}=0.5+\delta$.

Table 3: The mean and standard deviation of two estimators. Binomial Simulation. $\kappa\equiv 4$ and $\pi_{i}=p_{i}=0.5\pm\delta$.  $\delta$ | Naive | GMLE
---|---|---
$0.3$ |  0.559, (0.012) |  0.502, (0.014)
$0.2$ |  0.522, (0.011) |  0.504, (0.012)
$0.1$ |  0.504, (0.010) |  0.501, (0.010)

The final reported simulations are summarized in Table 4. We study Binomial
sampling with various values of $\kappa$, fixed at $\kappa=1,\dots,5$. Again,
there are two types of strata, 500 strata for each of the two types; for 500
strata $p_{i}$ and $\pi_{i}$ are sampled independently from $U(0.1,0.6)$,
while for the rest of the strata they are i.i.d. from $U(0.4,0.9)$.

Table 4: Binomial Simulations with continuous $G$. $\kappa$=1,2,3,4,5. $\kappa$ | Naive | GMLE
---|---|---
$1$ |  0.544, (0.019) |  0.530, (0.015)
$2$ |  0.528, (0.014) |  0.502, (0.021)
$3$ |  0.522, (0.014) |  0.498, (0.022)
$4$ |  0.517, (0.012) |  0.499, (0.020)
$5$ |  0.512, (0.009) |  0.501, (0.013)

It is surprising how well the GMLE is doing already for $\kappa=2,3$, in spite
of the non-identifiability of $G$ and the inconsistency of the GMLE.

##  5 Confidence Interval for $\eta_{G}$

Although the GMLE $\hat{\eta}_{G}$ is an appealing estimator, we do not have a
clue about how good is its performance, beside the consistency reults we
established; this is especially under non-identifiability where consistency is
not implied. In the following we suggest an asymptotically level-$(1-\alpha)$
conservative confidence interval for $\eta_{G}=E_{G}\eta(\theta)$. We only
elaborate on defining the corresponding convex optimization problem.
Implementing and solving the sketched convex optimization problem, could be
challenging.

Let $Z=Z(Y)$ be a random variable with $M$ possible outcomes,
$z_{j},\;j=1,...,M$. One may think of $M$ ”chosen cells” in a goodness of fit
test. The considerations for the choice and for the number of cells is beyond
the scope of this section, in particular $M(n)\equiv M$ is fixed. Denote
$p_{G}^{j}=P_{G}(Z=z_{j})$; denote
$\hat{p}^{j}=\frac{\\{\\#i|Y_{i}=z_{j}\\}}{n}\equiv\frac{n_{j}}{n}$, where
$n_{j}$ is implicitly defined. Recall that:

$2(\sum_{j}n_{j}\log(\hat{p}^{j})-\sum_{j}n_{j}\log(p_{G}^{j}))\Rightarrow_{G}\chi^{2}_{(M-1)}.$

Let $\chi^{2}_{(M-1),1-\alpha}$, be the $(1-\alpha)$ quantile of a
$\chi^{2}_{(M-1)}$ distribution. Define

$\Gamma_{\alpha}=\\{G\;:\;\sum_{j}n_{j}\log({p}_{G}^{j})\geq\sum_{j}n_{j}log(\hat{p}^{j})-\frac{1}{2}\chi^{2}_{(M-1),1-\alpha}\\}.$

Observe that $\Gamma_{\alpha}$ is a convex set of distributions.

Assume that the functional $\eta^{*}(G)=E_{G}\eta(\theta)$, is a linear
functional. Note, this is the case in all the examples presented in Section 2.

Let

$\eta^{L}=\min_{G\in\Gamma_{\alpha}}\eta^{*}(G);\;\eta^{U}=\max_{G\in\Gamma_{\alpha}}\eta^{*}(G).$

By the above, finding $(\eta^{L},\eta^{U})$ is a convex problem, and in
addition the later interval is a $(1-\alpha)$ conservative confidence interval
for $E_{G}\eta(\theta)$.

It may be shown that the above is also a conservative $(1-\alpha)$-level CI
for $\sum\eta(\theta_{i})$ for fixed $\theta_{1},...,\theta_{n}$, and also a
$(1-\alpha)$-level conservative credible set for $\sum\eta(\theta_{i})$ when
$\theta_{1},...,\theta_{n}$, are independent realizations $\theta_{i}\sim G$.

##  6  Real Data Example.

A rough description of the ‘Social-Survey’, conducted yearly by the Israeli
census bureau, is the following. We randomly sample a $1/1000$ fraction of the
age 20 or older individuals in the registry, then verify their home address
and interview them in person.

We study here the ‘social survey’ data, accumulated for Tel-Aviv, in the
surveys collected during 2015–2017. The total sample size in those three years
is 1256. There are 156 ‘statistical-areas’ in Tel-Aviv, very roughly of equal
size, about 3000 individuals in each. Statistical-areas are considered
homogeneous in many respects, and we take them as our strata.

Let $K_{i}$ be the sample size in stratum $i$, $i=1,...156$, then our data
satisfy $K_{i}>0$, $i=1,...,156$. ( Admittedly, we neglected a few small
statistical-areas that actually had zero sample sizes ). In our analysis,
$p_{i}$ is the proportion of individuals in stratum $i$ that own their living
place (or, it is owned by a member of their household). The goal is to
estimate $n^{-1}\sum p_{i}$, roughly the proportion of individuals that own
their living-place.

The naive estimator is applicable since $K_{i}>0$, $i=1,...,156$. The
estimated proportion obtained by the naive estimator is:

$\frac{1}{n}\sum_{i=1}^{n}\frac{X_{i}}{K_{i}}={\bf 0.434}.$

On the other hand, the ‘extreme collapse’ estimator, satisfy:

$\frac{\sum X_{i}}{\sum K_{i}}={\bf 0.488}.$

The significant difference between the two estimates has to do also with the
fact that strata are, in fact, not of equal size. But, more importantly for
us, the ‘extreme collapse’ seems to over estimate the proportion, since owners
are over represented in the sample. One reason is that their address in the
registry is more accurate and thus it is easier to find them. In other words,
individuals are not MCAR (missing completely at random). This phenomena is
partially corrected by the stratification, when MAR (missing at random)
conditional on the fine strata is approximately right.

In Table 5 we compare the GMLE and the naive estimators in semi-real simulated
scenarios in which only a randomly sub-sampled with probability $\gamma$, of
the described sample, is retained. The average estimates based on 25
simulations applied on the real data with $\gamma=0.1,0.2,0.25$, are presented
in Table 5. The (random) number of simulated strata with zero sample sizes,
corresponding to $\gamma=0.1,0.2,0.25$, are around 70, 40, and 30,
correspondingly.

It is reasonable to assume that the number of observations in the strata are
Poisson random variables. Thus, the sample sizes $K_{i}$ in the simulated sub-
sample are ${\mbox{P}oisson}(\lambda_{i})$ too, $i=1,...,156$. The naive
estimate based on the entire data equals 0.434, which is a reasonable
benchmark.

In the semi-real simulations the GMLE seems to perform much better than the
naive estimator.

Table 5:  Home ownership in Tel-Aviv. Performance of the two estimators as function of sampling rate.  $\gamma$ | Naive | GMLE
---|---|---
$0.1$ | 0.471 | 0.457
$0.2$ | 0.467 | 0.443
$0.25$ | 0.447 | 0.434

##  7 Appendix

######  Lemma 1.

Assume A1-A5, and the setup of Theorem 4. Then, under both $G^{n_{j}}_{1}$ and
$G^{n_{j}}_{2}$,

i)

$\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\log(f_{\hat{G}^{n_{j}}}(Y_{i}))\rightarrow_{a.s.}E_{G^{0}}\log(f_{\hat{G}^{0}}(Y)).$
(8)

ii)

$\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\log(f_{{G}^{n_{j}}}(Y_{i}))\rightarrow_{a.s.}E_{G^{0}}\log(f_{G^{0}}(Y)).$
(9)

###### Proof.

We prove part i).

First we show that:

$\lim
E_{{G}^{n_{j}}}\log(f_{\hat{G}^{n_{j}}}(Y))=E_{{G}^{0}}\log(f_{\hat{G}^{0}}(Y)).$
(10)

Note that by A2 for every $\epsilon>0$, there exists large enough $m>0$, such
that for every $G^{\prime},G^{\prime\prime}\in\\{G\\}$,

$|E_{G^{\prime}}\log f_{G^{\prime\prime}}(Y)\times I(h(Y)\leq
m)-E_{G^{\prime}}\log f_{G^{\prime\prime}}(Y)|<\epsilon,$

here $I$ is an indicator function.

Hence, (10) is implied by the following,

$\displaystyle\lim_{n_{j}}E_{{G}^{n_{j}}}\log(f_{\hat{G}^{n_{j}}}(Y))\times
I(h(Y)\leq m)$ $\displaystyle=$
$\displaystyle\lim_{n_{j}}E_{{G}^{0}}\log(f_{\hat{G}^{n_{j}}}(Y))\times
I(h(Y)\leq m)\times\frac{f_{{G}^{n_{j}}}(Y)}{f_{{G}^{0}}(Y)}$ $\displaystyle=$
$\displaystyle E_{G^{0}}\lim_{n_{j}}\log(f_{\hat{G}^{n_{j}}}(Y))\times
I(h(Y)\leq m)\times\frac{f_{{G}^{n_{j}}}(Y)}{f_{{G}^{0}}(Y)}$ $\displaystyle=$
$\displaystyle E_{{G}^{0}}\log(f_{\hat{G}^{0}}(Y))\times I(h(Y)\leq m).$

The above is by Lebesgue dominating convergence theorem, applying A2 . Note
that, the likelihood ratio $\frac{f_{{G}^{n_{j}}}(Y)}{f_{{G}^{0}}(Y)}$ is
bounded when $h(Y)\leq m$. Equation (10) follows by letting
$m\rightarrow\infty$ and utilizing A2.

Under $G^{n_{j}}_{1}$, the proof will follow, by Borel Cantelli, if we show
that for every $\epsilon>0$:

$\sum_{n_{j}}P_{{G}^{n_{j}}_{1}}(|\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\log(f_{\hat{G}^{n_{j}}}(Y_{i}))-E_{{G}^{n_{j}}}\log(f_{\hat{G}^{n_{j}}}(Y))|>\epsilon)<\infty.$
(11)

The complication in the above is that in the terms
$f_{\hat{G}^{n_{j}}}(Y_{i})$, the random function $f_{\hat{G}^{n_{j}}}$
depends on its argument $Y_{i}$. We will prove (11) through the following
steps.

In fact the convergence in equation (11) is of a sub-series $n_{j}$, but also
the original series converges. Indeed, in the sequel we consider the original
series with a general index $n$. Let
$\log(\Psi)\in\\{\log(f_{G}),\;G\in\\{G\\}\\}$, be a fixed function, in
particular, independent of the observations. Then for every $\epsilon>0$

$P_{G_{1}^{n}}(|\frac{1}{n}\sum_{i=1}^{n}\log(\Psi(Y_{i}))-E_{G^{n}}\log(\Psi(Y))|)>\epsilon)<\kappa/n^{2},$

for a suitable large enough $\kappa$. The above is by is by Markov inequality,
utilizing the bounded 4’th moment assumption A2, and the fact that under
${G_{1}^{n}}$ (sampling with replacement), $Y_{i},\;i=1,...,n,$ are
independent.

Note, that by tightness, $\\{G\\}$ is compact under the weak convergence
topology, since that every sequence $G^{k}$ has a converging subsequence. For
every $\epsilon>0$ and $M>0$, by A5 and compactness of $\\{G\\}$, we may find
a finite cover of $L$ open sets, and corresponding functions
$\Psi_{1},...,\Psi_{L},\;\Psi_{i}\in\\{f_{G}\\}$, such that for every
$f_{0}\in\\{f_{G}\\}$

$\inf_{\Psi\in\\{\Psi_{1},...,\Psi_{L}\\}}||\log(f_{0}^{(M)})-\log(\Psi^{(M)})||_{\infty}<\epsilon.$
(12)

Next, since $L$ is finite, the above considerations coupled with Bonferroni,
imply that:

$P_{{G}^{n}_{1}}(\max_{\Psi\in\\{\Psi_{1},...,\Psi_{L}\\}}(|\frac{1}{n}\sum_{i=1}^{n}\log(\Psi(Y_{i}))-E_{{G}^{n}}\log(\Psi(Y))|>\epsilon)<\kappa/n^{2},$
(13)

for a suitable large enough $\kappa$.

Now we show that in Equation (11), for large enough $M$, we may neglect the
terms where $|\log(f_{G^{n_{j}}}(Y_{i}))|>M$. For every $\epsilon>0$, for
large enough $M$, a.s.,

$\limsup\frac{1}{n}\sum_{\\{i\;|\;|\log(f_{\hat{G}^{n}}(Y_{i}))|\geq
M\\}}|\log(f_{\hat{G}^{n}}(Y_{i}))|<\epsilon.$

This follows since $|\log(f_{\hat{G}^{n}}(Y_{i}))|\leq h(Y_{i})$, and by A2,
coupled with Borel Cantelli. Hence, in the following, considering large $M$,
we may neglect

$\frac{1}{n}\sum_{\\{i\;|\;|\log(f_{\hat{G}^{n}}(Y_{i}))|\geq
M\\}}|\log(f_{\hat{G}^{n}}(Y_{i}))|.$

By (12) we may conclude, by taking large enough $M$, that for every
$\epsilon>0$,

$P_{{G}^{n}_{1}}(\sup_{\Psi\in\\{\
f_{G}\\}}(|\frac{1}{n}\sum_{i=1}^{n}\log(\Psi(Y_{i}))-E_{{G}^{n}}\log(\Psi(Y))|>\epsilon)<\kappa/n^{2},$
(14)

for a suitable large enough $\kappa$.

Equation (11) follows by the last equation and the above.

The Lemma is now implied by Borel Cantelli, for the case $G^{n_{j}}_{1}$, of
sampling with replacement.

The convergence under $G^{n_{j}}_{2}$ is obtained by comparing sampling with
replacement to sampling without replacement. Let
$S_{n_{j}}=\sum_{i=1}^{n_{j}}\log(\Psi(Y_{i}))$. Then,
$E_{G^{n_{j}}_{1}}S_{n_{j}}^{4}\geq E_{G^{n_{j}}_{2}}S_{n_{j}}^{4}$. It may be
seen by applying Rao-Blackwell on the convex function
$g(S_{n_{j}})=S_{n_{j}}^{4}$, when conditioning on $k_{i},\;i=1,...,n_{j}$,
where $k_{i}$ is the number of times $\theta_{i}$ was sampled in the sampling
with replacement process, that defines $G^{n_{j}}_{1}$. The implied bounded
4’th moment of $S_{n_{j}}$ under $G^{n_{j}}_{2}$, coupled with (3), implies
(11), along the lines of the above, under $G^{n_{j}}_{2}$. This concludes our
proof.  $\square$

References:

  * P.J. Bickel, C. A. Klassen, Y.Ritov and J.A. Wellner (1998). Efficient and Adaptive Estimation For Semiparametric Models. Springer-Verlag NY.

  * Brown, L.D. and Greenshtein, E. (2009). Non parametric empirical Bayes and compound decision approaches to estimation of high dimensional vector of normal means. Ann. Stat. 37, No 4, 1685-1704.

  * Chen, J. (2017). ”Consistency of the MLE under Mixture Models.” Statist. Sci. 32 (1) 47 - 63.

  * Dicker, L. and Zhao, S. (2014). Nonparametric empirical Bayes and maximum likelihood estimation for high-dimensional data analysis. ArXiv preprint arXiv:1

  * B. Efron. (2014). Two modeling strategies for empirical Bayes estimation. Statistical science: a review journal of the Institute of Mathematical Statistics, 29(2):285.

  * Greenshtein, E. and Itskov, T (2018), Application of Non-Parametric Empirical Bayes to Treatment of Non-Response. Statistica Sinica 28 (2018), 2189-2208.

  * Gu, J. and Koenker, R. (2017). Unobserved Heterogeneity in Income Dynamics: An Empirical Bayes Perspective. Journal of Business and Economics Statistics. Volume 35, 2017 - Issue 1

  * Heinrich, Philippe, and Jonas Kahn. ”Strong identifiability and optimal minimax rates for finite mixture estimation.” Annals of Statistics 46.6A (2018): 2844-2870.

  * Jiang, W. and Zhang, C.-H., (2009), General maximum likelihood empirical Bayes estimation of normal means. Ann.Stat.. 37 No. 4, 1647-1684.

  * Karlis, D. and Xekalak, E. (2005), Mixed Poisson Distributions. International Statistical Review, 73, 1, 35–58, Printed in Wales by Cambrian Printers.

  * Kiefer, J. and Wolfowitz, J. (1956). Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters. Ann.Math.Stat. 27 No. 4, 887-906.

  * Koenker, R. and Mizera, I. (2014). Convex optimization, shape constraints, compound decisions and empirical Bayes rules. JASA  109, 674-685.

  * Laird, N. (1978). Nonparametric maximum likelihood estimation of a mixing distribution. JASA 78, No 364, 805-811.

  * Lindsay, B. G. (1995). Mixture Models: Theory, Geometry and Applications. Hayward, CA, IMS.

  * Little, R.J.A and Rubin, D.B. (2002). Statistical Analysis with Missing Data. New York: Wiley

  * Long, Feng and Lee, H. Dicker (2018). Approximate nonparametric maximum likelihood for mixture models: A convex optimization approach to fitting arbitrary multivariate mixing distributions. Computational Statistics & Data Analysis 122: 80-91.

  * J.Pfanzagl (1993). Incidental Versus Random Nuisance Parameters. Ann Stat, 21, 1663-1691.

  * Sujayam Saha. Adityanand Guntuboyina. ”On the nonparametric maximum likelihood estimator for Gaussian location mixture densities with application to Gaussian denoising.” Ann. Statist. 48 (2) 738 - 762, April 2020. https://doi.org/10.1214/19-AOS1817

  * Zhang, C.-H. (1997). Empirical Bayes and compound estimation of a normal mean. Statist. Sinica 7 181–193.

  * Zhang, C-H. (2005). Estimation of sums of random variables: Examples and information bounds. Ann. Stat. 33, No.5. 2022-2041.


# Accelerated electrons observed down to $<$7 keV in a NuSTAR solar microflare

Lindsay Glesener11affiliation: University of Minnesota, Minneapolis, USA , Säm
Krucker22affiliation: University of California at Berkeley, Berkeley, USA
33affiliation: University of Applied Sciences and Arts Northwestern
Switzerland, Windisch, Switzerland , Jessie Duncan11affiliation: University of
Minnesota, Minneapolis, USA , Iain G. Hannah44affiliation: University of
Glasgow, Glasgow, UK , Brian W. Grefenstette55affiliation: California
Institute of Technology, Pasadena, USA , Bin Chen66affiliation: New Jersey
Institute of Technology, Newark, USA , David M. Smith77affiliation: University
of California at Santa Cruz, Santa Cruz, USA , Stephen M. White88affiliation:
Air Force Research Laboratory, Albuquerque, USA , Hugh Hudson22affiliation:
University of California at Berkeley, Berkeley, USA 44affiliation: University
of Glasgow, Glasgow, UK

###### Abstract

We report the detection of emission from a non-thermal electron distribution
in a small solar microflare (GOES class A5.7) observed by the Nuclear
Spectroscopic Telescope Array (NuSTAR), with supporting observation by the
Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI). The flaring
plasma is well accounted for by a thick-target model of accelerated electrons
collisionally thermalizing within the loop, akin to the “coronal thick target”
behavior occasionally observed in larger flares. This is the first positive
detection of non-thermal hard X-rays from the Sun using a direct imager (as
opposed to indirectly imaging instruments). The accelerated electron
distribution has a spectral index of 6.3$\pm$0.7, extends down to at least 6.5
keV, and deposits energy at a rate of $\sim 2\times$1027 erg s-1, heating the
flare loop to at least 10 MK. The existence of dominant non-thermal emission
in X-rays down to $<$5 keV means that RHESSI emission is almost entirely non-
thermal, contrary to what is usually assumed in RHESSI spectroscopy. The ratio
of non-thermal to thermal energies is similar to that of large flares, in
contrast to what has been found in previous studies of small RHESSI flares. We
suggest that a coronal thick target may be a common property of many small
microflares based on the average electron energy and collisional mean free
path. Future observations of this kind will enable understanding of how flare
particle acceleration changes across energy scales, and will aid the push
toward the observational regime of nanoflares, which are a possible source of
significant coronal heating.

Sun: corona — Sun: flares — Sun: X-rays, gamma rays

##  1 Introduction

While large solar flares of classes M and X garner most public and scientific
attention, small flares occur far more frequently. Microflares of GOES class A
and B occur exclusively in active regions (e.g. Hannah et al., 2008; Christe
et al., 2008). Small, sub-A-class brightenings have been observed in active
and in quiet regions (e.g. Kuhar et al., 2018), but it is still unclear as to
whether faint quiet-Sun brightenings signify the same physical processes as
those that occur in active regions. Microflares are generally observed to be
similar in nature to larger flares, with impulsive phases followed by gradual
cooling (e.g. Glesener et al., 2017).

There is a great deal of interest in studying the parameters of flare-
accelerated electrons across a wide range of flare magnitudes, including to
very small flare energies. The observed scaling properties can help to assess
theories of particle acceleration. In particular, since nanoflares are
possible coronal heating candidates, it is of great interest to examine
particle acceleration in small events to determine if they are similar to
large flares and how much energy they could deposit in the corona. One factor
that could affect the efficiency of particle acceleration across energy scales
is the magnitude of the guide field relative to the reconnecting field. Dahlin
et al. (2016) and Dahlin et al. (2017) have shown that the presence of a
large-order guide field (greater than the reconnecting field) suppresses
particle acceleration, which would lead one to expect low acceleration
efficiency for nanoflares. While the small flares currently observable in hard
X-rays (HXRs) are still far from the nanoflare scale, each advance in
sensitivity pushes further toward this observational regime.

HXR observations are ideal for characterizing the hot thermal plasma and any
non-thermal emissions generated in small energy releases. However, the most
advanced solar HXR spacecraft instrument, the Reuven Ramaty High Energy Solar
Spectroscopic Imager (RHESSI) spacecraft, was limited in its sensitivity to
small flares due to its indirect imaging method, although GOES class A
microflares could be observed (Hannah et al., 2008; Christe et al., 2008). The
recent advent of directly focusing HXR instruments in the form of the Nuclear
Spectroscopic Telescope Array (NuSTAR) spacecraft and the Focusing Optics
X-ray Solar Imager (FOXSI) sounding rocket have enabled the observation of
small microflares orders of magnitude fainter than those observed by RHESSI
(Krucker et al., 2014; Glesener et al., 2016, 2017; Wright et al., 2017; Kuhar
et al., 2018; Athiray et al., 2020; Hannah et al., 2019; Vievering et al.,
2020). However, while some of these observations implied a high-energy excess
that might arise from accelerated electron distributions, a clear, distinct
measurement of flare-accelerated electrons by focusing HXR instruments has
heretofore been prevented by the faintness of these flares and the rarity of
observational opportunity (since neither NuSTAR nor FOXSI observes the Sun
often).

Here we report the first direct detection of non-thermal emission from a solar
flare using a focusing HXR imager. We analyze the microflare’s thermal and
non-thermal properties and compare these to larger flares.

##  2 Observations

NuSTAR is a NASA Astrophysics Small Explorer launched in 2012 (Harrison et
al., 2013). Unlike all previous HXR-observing spacecraft, NuSTAR utilizes
directly focusing HXR optics to achieve far better sensitivity than any
previous HXR instrument. Although it is an astrophysics mission, NuSTAR can
measure faint solar emission during relatively quiet times, when best use is
made of the instrument’s limited throughput (Grefenstette et al., 2016).
Observations of the Sun are performed several times per year for one to
several hours at a time.

On 2017 August 21, NuSTAR observed the Sun for one orbit just before the solar
disk was partly occulted by the Moon (NuSTAR observation IDs 20312001001 and
20312002001). For a few minutes on this day, the Sun was totally occulted as
viewed from several locations in North America, an event commonly referred to
as the “Great American Eclipse.” Due to the high level of public excitement
and scientific interest generated by this event, most telescopes capable of
observing the Sun did so on that day.

NuSTAR’s target of interest on 2017 August 21 was an active region with NOAA
number 12671, which was observed from 18:49:58 to 19:50:03 UTC,111Summary
plots and information can be found at http://ianan.github.io/nsigh_all/. a
total of 3605 s, ($\sim$2940 s before the lunar occultation of the region
began at $\sim$19:39). NuSTAR’s livetime during the non-occulted observation
was, on average, 0.3% (with a minimum of 0.1% at microflare peak) for an
effective exposure of $\sim$9 s. This region produced a few C-class flares
before and after the NuSTAR observation and produced several evident
microflares during the observation. Here, we concentrate on a microflare
occurring in the west of the active region, at a location approximately
[360,45] arcsec west and north of the solar center. (Future work will analyze
in detail the entire set of NuSTAR microflares in this region.)

###  2.1 Microflare temporal and spatial observations

Figure 1:  Lightcurves from NuSTAR and GOES. Panels (a) and (b) show GOES soft
X-ray flux, with an exponentially falling background subtracted for Panel (b).
In Panel (a), the shorter-wavelength GOES channel flux has been multiplied by
a factor of 10. Panels (c)-(f) show livetime-corrected NuSTAR counts from both
telescopes in 10 s time bins and several energy bands, with 1-sigma
statistical uncertainties shown. Pale blue lines show data integrated over the
entire active region, while black lines show a 2 arcmin FoV around the
microflare site. Panel (g) shows that the high-energy NuSTAR emission mimics
the derivative of the GOES long wavelength flux in the first minutes of the
flare, indicative of the Neupert effect. The gray box indicates the 3 min
interval chosen for spectroscopy, and the 2-4 keV background-subtracted curve
is shown for comparison.

NuSTAR has two quasi-identical telescopes. These record X-rays on a single-
photon basis; events can then be arbitrarily binned in space, time, and
energy. Since NuSTAR pointing knowledge carries large uncertainties during
solar observation (see Grefenstette et al., 2016), we coaligned NuSTAR images
empirically to data from the Atmospheric Imaging Assembly (AIA) aboard the
Solar Dynamics Observatory (SDO). A linear combination of AIA data in the 94Å,
171Å, and 211Å filters was taken to isolate the Fe XVIII contribution, as in
Del Zanna (2013); this line has a formation temperature of log T $\approx$ 6.9
and is sensitive to a temperature range that overlaps that of NuSTAR. We
coaligned NuSTAR and AIA Fe XVIII data at the peak of the microflare
($\sim$18:55 UTC) and then cross-correlated the NuSTAR images to each other;
this method assumes slow (or no) source motion. Uncertainties in this
coalignment are estimated to be $\sim$10 arcsec.

Figure 1 shows (Panels (a)-(b)) soft X-ray lightcurves from the GOES X-ray
Spectrometer (XRS). After subtracting an exponentially falling background from
the long-wavelength flux, the GOES class is A5.7. Panels (c)-(f) show NuSTAR
lightcurves in several energy bands for a 2 arcmin region centered on the
flare, and also for a 12 arcmin region that corresponds to NuSTAR’s entire
field of view (FoV). The active region was contained well within this 12
arcmin FoV. At lower NuSTAR energies, emission from the entire active region
is evident, while at higher energies, emission emanates from the flare only.
(Residual differences between blue and black curves in the 8-10 keV energy
band are due to the wings of the instrumental point spread function, half
power diameter $\sim$1 arcmin.) High time variability is evident, especially
at higher energies. NuSTAR high-energy emission closely follows the derivative
of the flux in the GOES low-energy channel in the first few minutes of the
microflare (see Panel (g)). In this panel, GOES data have been smoothed over a
boxcar interval of 2 min before taking the derivative, and both NuSTAR and the
GOES derivative have been normalized to their maximum values over the plot
time range. The gray box in Panel (g) shows a 3 min interval at the beginning
of the flare (18:50–18:53 UT) on which we concentrate our efforts in this
paper; this interval was chosen because it covers the impulsive phase of the
flare and because the NuSTAR pointing was relatively steady over this
interval; the last $\sim$minute of the impulsive phase was excluded due to
pointing motion, which would complicate analysis.

Figure 2:  SDO/AIA images overlaid with NuSTAR, RHESSI, and EOVSA emission.
All AIA images are integrated from 18:50 to 18:53 UT. (Panel (a)) NuSTAR and
EOVSA emission overlaid on an AIA image of the entire active region, from
which the Fe XVIII component has been isolated. The NuSTAR image has been
deconvolved (50 iterations) and contour levels are 5, 10, 20, 40, 60, and 80%
of the maximum. EOVSA contours are 30, 50, 70, and 90% of the maximum. The
black box indicates the region shown in other panels. (Panel (b)) NuSTAR
emission deconvolved for 800 iterations (2-6 keV) and 100 iterations (6-10
keV) overlaid on an AIA image from the same time interval. In Panels (a) and
(b), NuSTAR images have been coaligned to AIA images, and NuSTAR data shown
are from Focal Plane Module B only. (Panels (c)-(e)) RHESSI vis_fwdfit images
overlaid on various AIA filter images. (No coalignment is necessary.) The HXR
emission from both RHESSI and NuSTAR shows good agreement with the AIA flaring
loop morphology.

Images of NuSTAR HXR emission are shown overlaid on AIA images in Figure 2.
NuSTAR’s two detector assemblies are termed Focal Plane Modules A and B (FPMA
and FPMB). For this event, FPMB recorded higher-quality data because some of
the bright emission fell in the small gap between the detector quadrants of
FPMA, so FPMB data are utilized for the images. In Panel (a), FPMB data from 2
to 10 keV have been integrated over the three-minute interval indicated in
Figure 1, have had the instrument point spread function deconvolved for 50
iterations using the IDL procedure max_likelihood.pro222Available within the
IDL Astronomy User’s Library at https://idlastro.gsfc.nasa.gov, and have been
coaligned to AIA data as previously described. This figure also includes data
from the Extended Owens Valley Solar Array (EOVSA), which is sensitive to
microwave emission from flare-accelerated electrons (e.g. Gary et al., 2018).
Panel (b) shows a zoomed-in image of FPMB emission in two energy bands after
800 iterations of deconvolution for the 2-6 keV band and 100 iterations for
the 6-10 keV band. (Different iteration numbers are chosen based on the
statistics available in each image.) Both NuSTAR sources are shown at 10, 30,
50, 70, and 90% level contours. The NuSTAR and AIA source shapes are similar,
and all HXR emission (in all available energy ranges) emanates primarily from
the flaring loop(s).

The microflare was observed by RHESSI, although it was too faint for inclusion
in the auto-generated RHESSI flare list. At this late stage in RHESSI’s
mission (15.5 years post-launch), only detectors 1, 3, 6, and 8 were
operating. Analysis of a microflare this faint and at such low energies is
challenging with RHESSI’s performance at the time, but the array of spatial
frequencies covered by this set of four subcollimators is sufficient to
produce an image of the microflare using the vis_fwdfit imaging method333See
https://hesperia.gsfc.nasa.gov/rhessi3/software/imaging-software/image-
algorithm-summary/index.html for a summary of RHESSI imaging algorithms., as
was used in Hannah et al. (2008). This method presupposes a source shape (in
this case a loop) and forward-fits the source parameters to the observed
visibilities. The result of this method for energies 6–9 keV is shown in
Panels (c)-(e) of Figure 2. RHESSI images produced using other imaging
algorithms (e.g. Clean; not shown) and at higher energies (e.g. 8-10 keV) all
show similar results; the HXR loop matches the position, loop shape, and
rotation angle of the AIA loop. RHESSI provides highly accurate source
locations, and so no co-alignment with AIA was necessary. In summary, all HXR
observations, from both NuSTAR and RHESSI, reveal that HXR emission emanates
primarily from the flaring coronal loop(s).

Figure 3:  Spectral fits to NuSTAR data in the 18:50–18:53 UTC time interval
for (a) an isothermal model, (b), a double-thermal model, and (c) an
isothermal + broken power law model. Count spectra shown are livetime- and
pileup-corrected, and a gain correction has been included in the fit. In each
spectrum, the top panel shows the count spectrum (and fit) for the two
telescopes (FPMA, red, and FPMB, black) and the bottom panel shows the error-
normalized residual distribution. The residuals indicate a more appropriate
fit for the isothermal + broken power law model than for either of the purely
thermal distributions, revealing a flare-accelerated electron distribution.
Fit parameter values are listed in Table 1.  Table 1: Fit parameters for all
spectral models in Section 2.2, along with best fit statistic values.

| XSPEC Models | OSPEX Models
---|---|---
| vapeca | vapec$+$ vapeca | vapec$+$ bknpower | vth$+$ thick2 | warm_thick
Temperature [MK] | 12.6${}^{+0.1}_{-0.3}$ | 11.8${}^{+0.4}_{-0.3}$ | 10.4${}^{+0.4}_{-0.2}$ | 10.3${}^{+0.7}_{-0.7}$ | 10.2${}^{+0.7}_{-0.7}$
EM [1045 cm${}^{-3}]$ | 4.0${}^{+0.1}_{-0.1}$ | 4.5${}^{+0.1}_{-0.1}$ | 4.2${}^{+0.5}_{-0.4}$ | 5.0${}^{+1.3}_{-1.3}$ |
Density [109 cm-3] |  |  |  |  | 6.0${}^{+2.0}_{-2.0}$
Temperatureb [MK] |  | 400${}^{+60}_{-180}$ |  |  |
EM [1045 cm${}^{-3}]$ |  | 0.015${}^{+0.003}_{-0.003}$ |  |  |
Break Energy [keV] |  |  | 5.0${}^{+0.1}_{-0.1}$ |  |
Photon Indexc $\gamma$ |  |  | 5.5${}^{+0.3}_{-0.2}$ |  |
Norm [phot keV-1 cm-2 s-1 at 1 keV] |  |  | 530${}^{+60}_{-50}$ |  |
Cutoff Energy [keV] |  |  |  | 6.2${}^{+0.9}_{-0.9}$ | 6.5${}^{+0.9}_{-0.9}$
Electron Index $\delta$ |  |  |  | 6.2${}^{+0.6}_{-0.6}$ | 6.3${}^{+0.7}_{-0.7}$
Electron Flux [1035 e- s-1] |  |  |  | 2.1${}^{+1.2}_{-1.2}$ | 1.8${}^{+0.8}_{-0.8}$
Loop Half Length [Mm] |  |  |  |  | 15 (fixed)
C-Statistic (reduced) | 2.5 | 1.8 | 1.2 |  |
$\chi^{2}$ value (reduced) | 2.0 | 2.0 | 1.1 | 0.8 | 0.7

  * a

The purely thermal fits are included for fit comparison purposes only; they
are poorer fits than the thermal$+$non-thermal fits and are not good
representations of this flare.

  * b

Parameter was allowed to vary only up to 40 keV (464 MK).

  * c

Index above the break. The index below the break was fixed at 2.

###  2.2 HXR spectral fitting

An examination of the NuSTAR spectral data over the first three minutes of the
flare indicated a small pileup component, and so a pileup correction was
performed on the count spectra (see Appendix C in Grefenstette et al., 2016).
Since each NuSTAR event is assigned a “grade,” comparison of events of various
grades gives an estimate of the pileup contribution, which can then be
subtracted from the spectrum. The necessary pileup correction is no more than
6.25% in any 0.64 keV energy band. Statistical uncertainties were widened to
account for uncertainties in the subtracted components. Additionally, a gain
slope correction was included to account for variations in the NuSTAR gain
that occur only at extremely high rates (e.g. livetime of a few percent)
encountered at the Sun (Duncan et al., 2019); this parameter was allowed to
vary for XSPEC fits and was fixed to 0.95 (the best fit value) for OSPEX fits.
The gain adjustment resulted in a 5% increase in temperature and a 20%
decrease in emission measure.

Following these pileup and gain corrections, fitting of various spectral
models was performed using the XSPEC software package (Arnaud, 1996). Fits
were performed simultaneously to FPMA and FPMB using the C-statistic to assess
goodness of fit. Figure 3 shows the results of fitting for three different
spectral models: an isothermal model (vapec), a double thermal model
(vapec+vapec), and an isothermal model plus a broken power law
(vapec+bknpower). The fit parameters are summarized in Table 1. For the
bknpower component, the spectral index below the break (where emission is
dominated by thermal plasma) was fixed at 2. The simplest model (isothermal
only) shown in Panel (a) exhibits a high-energy excess that is not well fit by
the model. The double thermal model shown in Panel (b) picks up this high-
energy component but retains systematic residuals that are not well
distributed across energy; the iron line intensity is poorly predicted.
Additionally, the temperature required by the hotter component (400 MK) is
unreasonably high. When the temperature was restricted to a more physical
range (e.g. $<$100 MK), the fit value was always driven to the highest allowed
temperature; the iron line was still poorly predicted; and the fit statistic
values were worse. The third model, which includes a broken power law in
photon space arising from a non-thermal electron distribution, exhibits well
distributed residuals and a significantly better fit statistic than the
thermal-only fits. This model is selected as the best fit to the data and
reveals the presence of an accelerated electron distribution.

As is often the case in fitting HXR spectra, the power-law spectral break
(which is related to the low-energy cutoff of the accelerated electron
distribution) is poorly constrained in the presence of bright thermal plasma.
The fit value is best viewed as an upper limit; see Section 3 of Holman et al.
(2011) for a thorough discussion. In fact, for this microflare, even a thermal
component plus an unbroken power law produced only a slight worsening of the
fit, bringing the temperature up to $T=11.9\pm(0.9,0.6)$ MK and emission
measure down to $EM=2.3\pm(0.3,0.6)\times 10^{45}$cm-3. Although we do not
consider a single power law likely, since it would require an accelerated
electron distribution extending down to extremely low energies, we use it to
set one bound on the flare thermal parameters. The resulting range of allowed
parameters (for the thermal+power-law models) is $T=10.2-12.8$ MK and
$EM=(1.7-4.7)\times 10^{45}$ cm-3. We consider this a range of allowed
parameters, but the most likely ones are those in the last column of Table 1.

AIA data were examined for consistency with the NuSTAR temperature. A simple
ratio was taken of data from AIA filters 131Å and 94Å (filters with
significant and relatively isolated response to hot flare-temperature plasma).
This ratio yields temperatures of 9.2–10.5 MK between 18:50 and 18:53 UTC and
an emission measure consistent with that obtained via the NuSTAR spectral fit.

Figure 4: Results of thick-target spectral fitting in OSPEX using models
(left) thick2 and (right) thick_warm, which model an accelerated electron
distribution propagating in a cold or warm plasma target, respectively. Fits
were performed to FPMB data only. The warm target model fits the data well
with no additional thermal component needed, indicating that the thermal
plasma arises from energetic electron thermalization within the loop. For the
warm-target model, the loop half length was fixed to 15 Mm from AIA images and
both temperature and density were allowed to vary.

A more physical model than a broken power law is to directly fit an electron
distribution to the X-ray data. To accomplish this, a thick-target model of
X-rays emitted by accelerated electrons (thick2) was fit to NuSTAR data along
with an isothermal component (vth) in the spectral fitting package OSPEX444See
https://hesperia.gsfc.nasa.gov/rhessi3/software/spectroscopy/spectral-
analysis-software/index.html., which is commonly used to fit solar HXR flares.
Since the emission is integrated over the spatial extent of the flare and over
a few minutes, we assume that the energetic electron distribution must
completely thermalize, so no thin-target fit was performed. Because OSPEX
performs only single-instrument fits and only uses a chi-squared value as a
fit parameter, data from FPMB (which has better coverage of this flare) were
selected and rebinned to ensure at least 10 counts in each energy bin. The
result of this fitting is shown in Panel (a) of Figure 4. The parameters
obtained from this fit are shown in Table 1 and are consistent with the
parameters obtained from the broken power law fit in Figure 3. (Note that in a
thick-target model, the electron and photon spectral indices $\delta$ and
$\gamma$ are related by $\delta=\gamma+1$.) An energy-integrated electron rate
of $(2.1\pm 1.2)\times 10^{35}$ electrons s-1 is obtained, along with a low-
energy cutoff of $6.2\pm 0.9$ keV and an electron power-law index of $6.2\pm
0.6$.

The thick2 model assumes a cold plasma target, requiring any thermal plasma to
be added as a separate component. On the other hand, the warm thick-target
model (WTTM) implemented in OSPEX as thick_warm provides a self-consistent fit
assuming the presence of a thermal component due to electrons thermalized out
of the accelerated population (see Kontar et al., 2015, 2019). We find that
the thick_warm model can well fit the entire NuSTAR spectrum for this flare,
with no additional thermal component required; see panel (b) of Figure 4. For
this fit, the warm plasma length was fixed to 15 Mm (the loop half-length from
AIA images) and both temperature and density were allowed to vary, as were the
spectral index, low-energy cutoff, and normalization of the accelerated
electron distribution555We note that our fitting approach differs from that
suggested by Kontar et al. (2019) because precision on the low energy cutoff
is not our goal.. The obtained results are almost identical to those obtained
via the cold thick-target model (CTTM), indicating that the CTTM is a valid
approximation for this case.

In summary, a multitude of spectral fitting tools and models were utilized to
fit the observed NuSTAR spectrum, including the X-ray fitting tools commonly
used by both the solar and non-solar astrophysics communities. All of these
tools and models point toward a thermal component plus a non-thermal, flare-
accelerated electron population, and the different methods produce consistent
quantitative results. The gain and pileup corrections mentioned at the start
of this section had small impacts on the fit parameters but did not have any
impact on which model provided the best fit. The next section will interpret
the spectral fitting results for their consequences on the dynamic evolution
of this flare.

##  3 Discussion

NuSTAR spectral fitting indicates that the non-thermal bremsstrahlung X-ray
flux intensity is greater than or equal to the thermal flux intensity above
$\sim$4.5 keV (see Figures 3 and 4). Further evidence for the non-thermal,
rather than thermal, nature of the emission comes from the close relationship
between the NuSTAR high-energy emission and the derivative of the GOES
emission (see Figure 1, Panel (g)). This relationship, known as the Neupert
effect, is commonly observed in larger flares and is often interpreted to
indicate that the energy collisionally deposited by non-thermal electrons is
responsible for the generation of the observed thermal plasma (Neupert, 1968;
Veronig et al., 2002).

NuSTAR does not have fine enough angular resolution to resolve much spatial
detail of the flare source shape at high energies, but the 6–9 keV RHESSI
source map in Figure 2 shows a loop matching the location, elongation, and
orientation of the loop observed by AIA. Therefore, the non-thermal emission
observed by NuSTAR and RHESSI emanates from the flare loop, not from its
footpoints, akin to the “coronal thick target” flares studied by Veronig &
Brown (2004); Veronig et al. (2005); Fleishman et al. (2016). This differs
from the standard thick-target flare model in which energy is transported from
the corona to the chromosphere primarily by accelerated electron beams, which
emit strong bremsstrahlung radiation at their thermalization locations in the
chromosphere and much fainter HXR emission from collisions in the tenuous
corona.

A look at the energetics supports the coronal thick target scenario. The non-
thermal fit to NuSTAR data utilizing the warm-target model (described in
Section 2.2) reveals that the accelerated electron distribution extends down
to a cutoff energy of $\sim$$6.5$ keV, or possibly even lower given that the
cutoff is not well constrained in the particular way we are utilizing
thick_warm, and has a spectral index of $\sim$$6.3$. This distribution has an
average electron energy of $\sim$8 keV. To investigate the plasma within which
this distribution propagates, we made geometric estimates from AIA. We
measured the area of the Fe XVIII AIA loop, using two thresholds on pixel
brightness (one liberal and one conservative). Taking the loop to be a bent
cylinder with (line-of-sight) depth equal to its (plane-of-sky) width, the
estimated volume is (1.2–2.2)$\times 10^{26}$ cm3. The full loop length,
measured directly from Fe XVIII AIA images, is 24–32 Mm. Assuming that volume
is filled with the thermal plasma indicated by the cold-target fit, and
assuming a filling factor of unity, this leads to a loop density of
(2.8–6.3)$\times 10^{9}$ cm-3. No projection effects have been corrected for
in these estimates, and if the filling factor is less than unity, then the
density is even higher. However, we note the calculated density is consistent
with that fit in the warm-target model ($6\times 10^{9}$ cm-3), showing good
correspondence between the two models.

If the flare-accelerated electrons are injected into the very top of the loop
and must travel one half loop length to reach the chromosphere, the column
density encountered is high enough to collisionally thermalize all electrons
below 5.4–7.1 keV. If the electrons are mirrored by the relatively stronger
magnetic fields encountered at lower altitudes, then they encounter even
higher column density. Most electrons below 9.4–12.4 keV would be
collisionally stopped by even 1.5 bounces in the loop. Therefore, it is
reasonable that the observed accelerated electron distribution deposits most
of its energy via collisions in the corona. This is different from the
expected behavior in larger-energy flares, in which electrons are observed up
to many tens or hundreds of keV, and low-energy electron cutoff energies are
often tens of keV; those distributions would easily penetrate loops of typical
coronal densities and produce the bright HXR footpoints typically observed in
large RHESSI flares. We note that even in larger flares, densities can
sometimes be high enough to produce this effect, as in the Veronig & Brown
(2004) flares.

From the spectral fits, the non-thermal electrons in this flare deposit energy
at a rate of $\sim$$2\times 10^{27}$ erg/s, for a total of $\sim$$4\times
10^{29}$ ergs when integrated over the three minutes chosen for spectroscopy.
The average thermal energy over this time is estimated from NuSTAR
spectroscopy to be $\sim$$3\times 10^{27}$ erg; these and other energetics
values are summarized in Table 2. We have compared this to thermal energy
estimates for this flare made using SDO/AIA data. Using the differential
emission measure estimation method of Cheung et al. (2015), we find that the
thermal energy in the flaring pixels rises from 1.3$\times 10^{27}$ to
2.4$\times 10^{28}$ erg between 18:50 and 18:53 UT, reaching a peak at
2.7$\times 10^{28}$ at 18:54 UT when considering only the emission measure
above 8 MK (corresponding to the plasma measured by NuSTAR). When considering
all temperatures, the thermal energy peaks at 3.6$\times 10^{28}$ erg. This
means that the ratio of non-thermal electron energy (over the three minutes
chosen for analysis) to maximum observed thermal energy (occurring shortly
after this interval) is a factor of 10. This is a larger non-thermal ratio
than any of the C class flares studied using RHESSI and GOES data by Warmuth,
& Mann (2016), but is a very typical ratio for the larger flares in that study
(see, e.g., Figure 7 of that work). Aschwanden et al. (2016) found an average
non-thermal-to-thermal energy ratio of 0.15-6.7, depending on the method used
for the non-thermal calculation. (That study relied on two different ways to
estimate the non-thermal electron energy. Their particular application of the
warm target method resulted in artificially high non-thermal energy estimates
as explained in Kontar et al. (2019), but their other, “cross-over” method
provided a conservative lower limit on the non-thermal energies, so the true
ratio probably lies in between.)

RHESSI flares of this thermal energy are close to RHESSI’s sensitivity limit;
this is evident from Figures 13 and 18 of Hannah et al. (2008). Figure 18
shows the rollover in the RHESSI flare frequency distribution below
$\sim$$10^{29}$ erg, and Figure 13 shows that this flare, with a temperature
of $\sim$10 MK and an EM of $\sim$$5\times 10^{45}$ cm-3, would not have been
included in that previous study. The reason for this is that the selection
criteria used in that study to produce reliable automated fits would have
excluded this microflare. So future work could revisit the RHESSI data to
search for similar microflares, during times the instrument was performing
optimally. However, it might be that with NuSTAR’s sensitivity to faint
emission, flares with higher non-thermal-to-thermal energy ratios can be
observed than in the past. This implies that the ratios of non-thermal to
thermal energies in Warmuth, & Mann (2016) might be much more consistent
across flare size if high HXR sensitivity were available for more small
flares. Due to NuSTAR’s limited solar observing time and limited throughput
(which restricts its spectral dynamic range), a thorough study of this point
must wait for a future, space-based direct HXR imager optimized for the Sun.

Table 2: Energetics values Spectral index | 6.3
---|---
Low energy cutoff | 6.5 keV
Average non-thermal power, |  $2\times 10^{27}$ erg s-1
18:50–18:53 |
Energy deposited over 3 min |  $4\times 10^{29}$ erg
NuSTAR thermal energy |  $3\times 10^{27}$ erg
(isothermal, avgd 18:50–18:53) |
AIA thermal energy |  $4\times 10^{28}$ erg
(DEM, peak energy, 18:54) |
Number accelerated electrons |  $2\times 10^{35}$ e- s-1
Average electron energy | 8 keV

The coronal thick target scenario does not rule out the occurrence of
chromospheric evaporation. The higher-energy tail of the accelerated electron
distribution could persist through enough column length to precipitate to the
flare footpoints and deposit energy to inspire evaporation. Alternatively,
footpoint heating via conduction could serve the same function; Warmuth, &
Mann (2016) found conductive loss rates to be dominant in small flares. Here,
the radiative loss rate assuming an average temperature of 10 MK and an
emission measure of $5\times 10^{45}$ cm-3 is $5\times 10^{23}$ erg s-1 and
the conductive loss rate assuming Spitzer conductivity is approximately
$10^{26}$ erg s-1. However, it is difficult to assess the footpoint area for
this flare, which has a strong impact on the conductive losses. If the scaling
relations in Warmuth, & Mann (2016) extend down (in temperature) to this
flare, then the conductive loss rate would be $\sim$$5\times 10^{26}$ erg s-1,
or about 0.25 of the non-thermal power, a significant fraction. The scenario
here is similar to the flares studied by Veronig & Brown (2004), which first
posited coronal thick-target flares; in those cases both electron
precipitation and conduction contributed to chromospheric evaporation. The two
flares described in Veronig & Brown (2004) also exhibited gradual behavior
(rather than impulsive) and rather steep electron power-law distributions, as
does this microflare.

Microflare loop lengths are not, in general, smaller than regular flare loop
lengths (Hannah et al., 2008; Glesener et al., 2017), and microflares tend to
have lower low-energy cutoffs (although this could be an observational bias)
and steeper spectra than larger flares do; this microflare observation
continues that trend. Given these two qualities, it is expected that flares
smaller in energy should exhibit, on average, a greater degree of coronal
thick-target behavior than larger flares do, and this trend should continue to
even fainter flares. It is likely appropriate to model such flares with a
collisional deposition of energy throughout the flare loop rather than only at
the footpoints.

The NuSTAR spectral observations clearly show that the measured RHESSI photons
for this flare mainly emanate from a non-thermal source. This challenges the
usual assumption made in HXR microflare analysis that the lowest energy
($\lesssim$7 KeV) counts are dominated by thermal emission. Further studies
are needed to investigate how commonly flares with non-thermal emission down
past the iron complex energies occur; the iron complex will be an important
disambiguator in those studies, as it is here. For this purpose, it will also
be useful to exploit the database of AIA observations in tandem with the
RHESSI archive, as was shown to be necessary in Inglis, & Christe (2014); Ryan
et al. (2014); and Aschwanden et al. (2015). As the derived non-thermal energy
content depends strongly on the low energy cutoff, a fresh look at the RHESSI
microflare statistics with different spectral assumptions could potentially
provide a significant update to microflare energetics.

##  4 Summary

In summary, we have reported the first evidence of non-thermal HXR emission
from the Sun using direct focusing instruments. This non-thermal emission was
observed in a small, A5.7 class microflare. The flare-accelerated electrons
have an average energy of $\sim$8 keV and extend down to a low-energy cutoff
energy of $\sim$6.5 keV, dominating the X-ray spectrum down to $<$5 keV. A
clear Neupert effect is observed. The non-thermal electrons deposit most of
their energy in the coronal loop, unlike most larger flares, which deposit
collisional energy primarily at the footpoints. The observation confirms that
flare particle acceleration occurs even in the faintest flares observable with
today’s instrumentation, and that the non-thermal energies can be large in
comparison to previous observations of small flares that were studied with
less sensitive instruments. Based on physical arguments, we suggest that
extremely small microflares and nanoflares may be likely to be coronal thick-
target flares, and simulations of such flares would be best served by
depositing flare-accelerated electron energy throughout the corona and
chromosphere, not only at footpoints.

Future direct-focusing HXR instrumentation with greater sensitivity, solar
optimization, and much more solar observing time will allow the measurement of
many more microflares like these, as well as fainter ones. Such
instrumentation has been developed with the FOXSI sounding rocket program and
has been proposed as the Fundamentals of Impulsive Energy Release in the
Corona Explorer (FIERCE) spacecraft concept (Shih et al., 2019). Instruments
based on this technology will allow full exploration of how particle
acceleration scales to small flares, how frequently small flares occur, and
how capable they are of heating the corona.

Support for this work was provided by an NSF Faculty Development Grant
(AGS-1429512), an NSF CAREER award (NSF-AGS-1752268), and the NASA NuSTAR
Guest Observer program (80NSSC18K1744). IGH is supported by a Royal Society
University Fellowship. The authors are indebted to Natasha Jeffrey and
Alexander Warmuth for their perspectives on this work. Some figures were
produced using IDL color-blind-friendly color tables by Wright (2017).

## References

  * Arnaud (1996) Arnaud, K. A. 1996, Astronomical Data Analysis Software and Systems V, 17
  * Aschwanden et al. (2015) Aschwanden, M. J., Boerner, P., Ryan, D., et al. 2015, ApJ, 802, 53
  * Aschwanden et al. (2016) Aschwanden, M. J., Holman, G., O’Flannagain, A., et al. 2016, ApJ, 832, 27
  * Athiray et al. (2020) Athiray, P.S., et al. 2020, ApJ, accepted
  * Cheung et al. (2015) Cheung, M. C. M., Boerner, P., Schrijver, C. J., et al. 2015, ApJ, 807, 143
  * Christe et al. (2008) Christe, S., Hannah, I. G., Krucker, S., McTiernan, J., & Lin, R. P. 2008, ApJ, 677, 1385-1394
  * Dahlin et al. (2016) Dahlin, J. T., Drake, J. F., & Swisdak, M. 2016, Physics of Plasmas, 23, 120704
  * Dahlin et al. (2017) Dahlin, J. T., Drake, J. F., & Swisdak, M. 2017, Physics of Plasmas, 24, 092110
  * Del Zanna (2013) Del Zanna, G. 2013, A&A, 558, A73
  * Duncan et al. (2019) Duncan, J. M., Glesener, L., Hannah, I., et al. 2019, American Astronomical Society Meeting Abstracts, 234, 204.04
  * Fleishman et al. (2016) Fleishman, G. D., Xu, Y., Nita, G. N., et al. 2016, ApJ, 816, 62
  * Gary et al. (2018) Gary, D. E., Chen, B., Dennis, B. R., et al. 2018, ApJ, 863, 83
  * Glesener et al. (2016) Glesener, L., Krucker, S., Christe, S., et al. 2016, Proc. SPIE, 9905, 99050E
  * Glesener et al. (2017) Glesener, L., Krucker, S., Hannah, I. G., et al. 2017, ApJ, 845, 122
  * Grefenstette et al. (2016) Grefenstette, B. W., Glesener, L., Krucker, S., et al. 2016, ApJ, 826, 20
  * Hannah et al. (2008) Hannah, I. G., Christe, S., Krucker, S., et al. 2008, ApJ, 677, 704
  * Hannah et al. (2011) Hannah, I. G., Hudson, H. S., Battaglia, M., et al. 2011, Space Sci. Rev., 159, 263
  * Hannah et al. (2019) Hannah, I. G., Kleint, L., Krucker, S., et al. 2019, ApJ, 881, 109
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103
  * Holman et al. (2011) Holman, G. D., Aschwanden, M. J., Aurass, H., et al. 2011, Space Sci. Rev., 159, 107
  * Inglis, & Christe (2014) Inglis, A. R., & Christe, S. 2014, ApJ, 789, 116
  * Kontar et al. (2015) Kontar, E. P., Jeffrey, N. L. S., Emslie, A. G., & Bian, N. H. 2015, ApJ, 809, 35
  * Kontar et al. (2019) Kontar, E. P., Jeffrey, N. L. S., & Emslie, A. G. 2019, ApJ, 871, 225
  * Krucker et al. (2014) Krucker, S., Christe, S., Glesener, L., et al. 2014, ApJ, 793, L32
  * Kuhar et al. (2018) Kuhar, M., Krucker, S., Glesener, L., et al. 2018, ApJ, 856, L32
  * Neupert (1968) Neupert, W. M. 1968, ApJ, 153, L59
  * Ryan et al. (2014) Ryan, D. F., O’Flannagain, A. M., Aschwanden, M. J., et al. 2014, Sol. Phys., 289, 2547
  * Shih et al. (2019) Shih, A. Y., Glesener, L., Christe, S., et al. 2019, AGU Fall Meeting Abstracts 2019, SH33A-08
  * Veronig et al. (2002) Veronig, A., Vršnak, B., Dennis, B. R., et al. 2002, A&A, 392, 699
  * Veronig & Brown (2004) Veronig, A. M., & Brown, J. C. 2004, ApJ, 603, L117
  * Veronig et al. (2005) Veronig, A. M., Brown, J. C., & Bone, L. 2005, Advances in Space Research, 35, 1683
  * Vievering et al. (2020) Vievering, J.T., et al., ApJ, in prep.
  * Warmuth, & Mann (2016) Warmuth, A., & Mann, G. 2016, A&A, 588, A116
  * Wright et al. (2017) Wright, P. J., Hannah, I. G., Grefenstette, B. W., et al. 2017, ApJ, 844, 132
  * Wright (2017) Wright, P. J. 2017, doi: 10.5281/zenodo.840393, url: https://github.com/PaulJWright/ColourBlind