# Bounded-memory adjusted scores estimation in generalized linear models with
large data sets

Patrick Zietkiewicz<EMAIL_ADDRESS>Department of
Statistics, University of Warwick
Gibbet Hill Road, Coventry, CV4 7AL, UK  Ioannis Kosmidis
<EMAIL_ADDRESS>Department of Statistics, University of Warwick
Gibbet Hill Road, Coventry, CV4 7AL, UK

###### Abstract

The widespread use of maximum Jeffreys’-prior penalized likelihood in
binomial-response generalized linear models, and in logistic regression, in
particular, are supported by the results of Kosmidis and Firth (2021,
Biometrika), who show that the resulting estimates are also always finite-
valued, even in cases where the maximum likelihood estimates are not, which is
a practical issue regardless of the size of the data set. In logistic
regression, the implied adjusted score equations are formally bias-reducing in
asymptotic frameworks with a fixed number of parameters and appear to deliver
a substantial reduction in the persistent bias of the maximum likelihood
estimator in high-dimensional settings where the number of parameters grows
asymptotically linearly and slower than the number of observations. In this
work, we develop and present two new variants of iteratively reweighted least
squares for estimating generalized linear models with adjusted score equations
for mean bias reduction and maximization of the likelihood penalized by a
positive power of the Jeffreys-prior penalty, which eliminate the requirement
of storing $O(n)$ quantities in memory, and can operate with data sets that
exceed computer memory or even hard drive capacity. We achieve that through
incremental QR decompositions, which enable IWLS iterations to have access
only to data chunks of predetermined size. We assess the procedures through a
real-data application with millions of observations, and in high-dimensional
logistic regression, where a large-scale simulation experiment produces
concrete evidence for the existence of a simple adjustment to the maximum
Jeffreys’-penalized likelihood estimates that delivers high accuracy in terms
of signal recovery even in cases where estimates from ML and other recently-
proposed corrective methods do not exist.

Keywords: iteratively reweighted least squares, incremental QR decomposition,
mean bias reduction, Jeffreys’-prior penalty, data separation

##  1 Introduction

Albert and Anderson (1984) shows that, in logistic regression, data separation
is a necessary and sufficient condition for the maximum likelihood (ML)
estimate to have at least one infinite component. Data separation occurs when
there is a linear combination of covariates that perfectly predicts the
response values, and results in the log-likelihood having an asymptote in some
direction of the parameter space. When separation occurs, ML fitting
procedures typically result in apparently finite estimates and estimated
standard errors, which are mere artefacts of the numerical optimization
procedure stopping prematurely by meeting the optimizer’s convergence
criteria. If infinite estimates go undetected, then inference about the model
parameters can lead to misleading conclusions. Mansournia et al. (2018) is a
recent overview of the practical issues associated with infinite estimates in
logistic regression. Infinite ML estimates can also occur for other binomial-
response regression models.

The detection of separation and, more generally, identification of which
components of the maximum likelihood estimate are infinite can be done prior
to fitting the model by solving appropriate linear programs. For example, the
detectseparation R package (Kosmidis et al., 2022) provides methods that
implement the linear programs in Konis (2007) for many binomial-response
generalized linear models (GLMs), and the linear program in Schwendinger et
al. (2021) for log-binomial models for relative risk regression. Such
procedures, however, while helpful to identify issues with ML, they provide no
constructive information about what to do when infinite ML estimates are
encountered.

It is important to note that infinite estimates occur for both small and large
data sets. For an illustration, see Section 4, where infinite estimates are
observed for a probit regression model fit on a data set with about 5.5
million observations and 37 covariates. Furthermore, for logistic regression
models with $n$ observations and $p$ covariates, where $p/n\to\kappa\in(0,1)$,
Candès and Sur (2020) prove that for model matrices generated from particular
distributions, there are combinations of $\kappa$ and the variance of the
linear predictor for which the ML estimate has infinite components with
probability one as $p/n\to\kappa$, regardless of the specific values of $p$ or
$n$.

The linear programs for detecting infinite estimates can be slow for large
model matrices, and infeasible in their vanilla implementations if the model
matrix and response vector do not fit in computer memory. For example, solving
the Konis (2007) linear program for the model of Section 4 with the default
settings of detect_infinite_estimates() method of the detectseparation R
package did not return a result even after 12 hours of computation on a 2021
MacBook Pro with an Apple M1 Max chip and 64 GB RAM.

To summarize, on one hand, ML estimation for binomial-response GLMs may lead
to arbitrarily large estimates and estimated standard errors, which are in
reality infinite, and if they go undetected can cause havoc to standard
inferential procedures. On the other hand, linear programs for detecting the
occurrence of infinite ML estimates prior to fitting are not constructive
about the consequences of having infinite estimates in the model fit or the
impact this may have on other parameters in the model, and have either long
run times for large data sets or their standard implementations are infeasible
for data sets that do not fit computer memory.

For those reasons, researchers typically resort to alternative estimators,
which in many settings have the same optimal properties as the ML estimator
has, and are guaranteed to be finite even in cases where the ML estimator is
not, regardless of the data set or its dimension. For example, Kosmidis and
Firth (2021) show that penalization of the likelihood by Jeffreys’ invariant
prior, or by any positive power thereof, always produces finite-valued maximum
penalized likelihood estimates in a broad class of binomial-response GLMs,
only under the assumption that the model matrix is of full rank. Another
approach that has been found to deliver finite estimates for many binomial-
response GLMs are the adjusted score equations for mean bias reduction (mBR)
in Firth (1993), which, for fixed $p$, result in estimators with mean bias of
smaller asymptotic order than what the ML estimator has in general. To date,
there is only empirical evidence that the reduced-bias estimator of Firth
(1993) has always finite components for binomial-response GLMs with arbitrary
link functions. For fixed-$p$ logistic regression, the bias-reducing adjusted
score functions end up being the gradient of the logarithm of the likelihood
penalized by Jeffreys’ invariant prior. Hence, in that case, the results in
Kosmidis and Firth (2021) show that the reduced-bias estimators have also
always finite components. Both Sur and Candès (2019) and Kosmidis and Firth
(2021) also illustrate that maximum Jeffreys’-prior penalized likelihood
(mJPL) results in a substantial reduction in the persistent bias of the ML
estimator in high-dimensional logistic regression problems with
$p/n\rightarrow\kappa\in(0,1)$, when the ML estimator exists. Those results
underpin the increasingly widespread use of mBR and similarly penalized
likelihood estimation in binomial regression models in many applied fields.

The ML estimates for GLMs can be computed through iterative reweighted least
squares (IWLS; Green 1984). It is well-known that due to the fact that the
each component of the working variates vector for the IWLS iteration depends
on the current value of the parameters and the corresponding observation, ML
estimation can be performed even with larger-than-memory data sets using
incremental QR decompositions as in Miller (1992). The incremental IWLS is
implemented in the biglm R package (Lumley, 2020) for all GLMs.

In this work, we begin by presenting a unifying IWLS procedure for mBR and
mJPL in GLMs through a modification of the ML working variates. However,
unlike the working variates for ML, the modified working variates involve
quantities, like the leverage, that depend on the whole dataset (Kosmidis et
al., 2020), and as a result IWLS through incremental QR decompositions is not
directly possible. To overcome this difficulty, we present and analyze the
properties two variants of IWLS for solving the mBR and mJPL adjusted score
equations, which only require access to data blocks of fixed size that the
user can specify in light of any local memory constraints. As a result, they
eliminate the requirement of keeping $O(n)$ quantities in memory, opening the
door for fitting GLMs with adjusted score equations and penalized likelihood
methods on data sets that are larger than local memory or hard drive capacity,
and are stored in a remote database. The procedures operate with either one-
or two-passes through the data set per iteration, and return the exact, not
approximate, mBR and mJPL estimates.

Using the IWLS variants, we conduct a comprehensive simulation experiment in
the high-dimensional logistic regression setting of Sur and Candès (2019),
where $p/n\to\kappa\in[0,1)$. The experiment results in novel evidence for
using mJPL. Specifically, we examine in depth and expand the observations of
Sur and Candès (2019) and Kosmidis and Firth (2021), and illustrate that mJPL
performs well and markedly better than the ML estimator, in the region below
Candès and Sur (2020)’s phase transition curve where the ML estimates
asymptotically exist. In addition, an analysis of the evidence we gather
allows us to derive and propose a remarkably simple scaling factor for the
mJPL estimator that depends only on $\kappa$, the variance of the linear
predictor and the size of the intercept parameter, and enables mJPL to recover
the true signal to high accuracy and for a wide range of settings in the
region above Candès and Sur (2020)’s phase transition curve, where the ML
estimates do no exist with probability approaching one. To our knowledge,
there has been no other proposal for signal recovery in the high-dimensional
logistic regression setting of Sur and Candès (2019) that applies regardless
of whether the ML estimator exists or not. The experiment also allows us to
assess the performance of the two IWLS variants. It is found that, despite the
two-pass implementation is slower than the one-pass, it requires less
iterations to converge. Also, in notable contrast to the observations in Sur
and Candès (2019, Supporting Information document, Section D), where mJPL is
reported to be computationally infeasible, requiring about $2.5$ hours for
$n=2000$ and $p=400$, the maximum average runtime we observed for mJPL with
the two-pass IWLS variant ranged from milliseconds to just above two minutes
in all settings we tried with $n=2000$ and $p$ ranging from $20$ to $1100$,
with no particular care in the choice of starting values; see, Table 2 for one
of those settings.

Section 2 sets notation by introducing GLMs and their ML estimation using
IWLS, and details how that estimation can be done using only data chunks of
fixed size thorough incremental QR decomposition’s. Section 3 presents the
adjusted score equations for mBR and mJPL in GLMs, and introduces and analyzes
the memory and computation complexity of two variants of IWLS that operate in
a chunk-wise manner, avoiding the memory issues that vanilla IWLS
implementations can encounter. Section 4 applies the algorithms for the
modelling of the probability of a diverted flight using probit regression,
based data on all $5,683,047$ commercial flights within the USA in 2000.
Section 5 presents a comprehensive computer experiment to empirically
investigate the limits of mJPL in terms of frequentist performance in the
framework of Candès and Sur (2020) for logistic regression with high-
dimensional data sets. The analysis of the evidence suggests that a rescaled
version of the mJPL estimator can recover the true signal to high accuracy and
for a wide range of settings, where other estimators do not exist. Finally,
Section 6 provides discussion and concluding remarks.

##  2 Generalized linear models

###  2.1 Model

Suppose that $y_{1},\ldots,y_{n}$ are observations on random variables
$Y_{1},\ldots,Y_{n}$ that are conditionally independent given
$x_{1},\ldots,x_{n}$, where $x_{i}$ is a $p$-vector of covariates. A GLM
(McCullagh and Nelder, 1989) assumes that, conditionally on $x_{i}$, $Y_{i}$
has an exponential family distribution with probability density or mass
function of the form

$f_{Y_{i}}(y;\theta_{i},\phi)=\exp\left\\{\frac{y\theta_{i}-b(\theta_{i})-c_{1}(y)}{\phi/m_{i}}-\frac{1}{2}a\left(-\frac{m_{i}}{\phi}\right)+c_{2}(y)\right\\}$

for some sufficiently smooth functions $b(.)$, $c_{1}(.)$, $a(.)$ and
$c_{2}(.)$, and fixed observation weights $m_{1},\ldots,m_{n}$. The expected
value and the variance of $Y_{i}$ are then

$\displaystyle E(Y_{i})$ $\displaystyle=\mu_{i}=b^{\prime}(\theta_{i})$
$\displaystyle Var(Y_{i})$
$\displaystyle=\frac{\phi}{m_{i}}b^{\prime\prime}(\theta_{i})=\frac{\phi}{m_{i}}V(\mu_{i})\,.$

The mean $\mu_{i}$ is linked to a linear predictor $\eta_{i}$ through a a
monotone, sufficiently smooth link function $g(\mu_{i})=\eta_{i}$ with

$\eta_{i}=\sum_{t=1}^{p}\beta_{t}x_{it}$ (1)

where $x_{it}$ can be thought of as the $(i,t)$th component of a model matrix
$X$, assumed to be of full rank, and
$\beta=(\beta_{1},\ldots,\beta_{p})^{\top}$. An intercept parameter is
customarily included in the linear predictor, in which case $x_{i1}=1$ for all
$i\in\\{1,\ldots,n\\}$.

###  2.2 Likelihood, score functions and information

The log-likelihood function for a GLM is $\ell(\beta)=\sum_{i=1}^{n}\log
f_{Y_{i}}(y_{i};g^{-1}(\eta_{i}),\phi)$, where $g^{-1}(.)$ is the inverse of
the link function, and $\eta_{i}$ is as in (1). Temporarily suppressing the
dependence of the various quantities on the model parameters and the data, the
derivatives of the log-likelihood function with respect to the components of
$\beta$ and $\phi$ are

$s_{\beta}=\frac{1}{\phi}X^{T}W(z-\eta)\quad\text{and}\quad
s_{\phi}=\frac{1}{2\phi^{2}}\sum_{i=1}^{n}(q_{i}-\rho_{i})\,,$ (2)

respectively, where $z=\eta+D^{-1}(y-\mu)$ is the vector of working variates
for ML, $\eta=X\beta$, $y=(y_{1},\ldots,y_{n})^{\top}$,
$\mu=(\mu_{1},\ldots,\mu_{n})^{\top}$, $W={\rm
diag}\left\\{w_{1},\ldots,w_{n}\right\\}$ and $D={\rm
diag}\left\\{d_{1},\ldots,d_{n}\right\\}$, with $w_{i}=m_{i}d_{i}^{2}/v_{i}$
being the $i$th working weight, and $d_{i}=d\mu_{i}/d\eta_{i}$,
$v_{i}=V(\mu_{i})$. Furthermore,
$q_{i}=-2m_{i}\\{y_{i}\theta_{i}-b(\theta_{i})-c_{1}(y_{i})\\}$ and
$\rho_{i}=m_{i}a^{\prime}_{i}$ are the $i$th deviance residual and its
expectation, respectively, with $a^{\prime}_{i}=a^{\prime}(-m_{i}/\phi)$,
where $a^{\prime}(u)=da(u)/du$.

The ML estimators $\hat{\beta}$ and $\hat{\phi}$ can be found by solution of
the score equations $s_{\beta}=0_{p}$ and $s_{\phi}=0$, where $0_{p}$ is a
$p$-vector of zeros. Wedderburn (1976) derives necessary and sufficient
conditions for the existence and uniqueness of the ML estimator. Given that
the dispersion parameter $\phi$ appears in the expression for $s_{\beta}$ in
(2) only multiplicatively, the ML estimate of $\beta$ can be computed without
knowledge of the value of $\phi$. This fact is exploited in popular software
like the glm.fit() function in R (R Core Team, 2023). The $j$th iteration of
IWLS updates the current iterate $\beta^{(j)}$ by solving the weighted least
squares problem

$\beta^{(j+1)}:=\left(X^{\top}W^{(j)}X\right)^{-1}X^{\top}W^{(j)}z^{(j)}\,,$
(3)

where the superscript $(j)$ indicates evaluation at $\beta^{(j)}$ (Green,
1984). The updated $\beta$ from (3) is equal to that from the Fisher scoring
step $\beta^{(j)}+\\{i_{\beta\beta}^{(j)}\\}^{-1}s_{\beta}^{(j)}$ where
$i_{\beta\beta}$ is the $(\beta,\beta)$ block of the expected information
matrix about $\beta$ and $\phi$

$\left[\begin{array}[]{cc}i_{\beta\beta}&0_{p}\\\
0_{p}^{\top}&i_{\phi\phi}\end{array}\right]=\left[\begin{array}[]{cc}\frac{1}{\phi}X^{\top}WX&0_{p}\\\
0_{p}^{\top}&\frac{1}{2\phi^{4}}\sum_{i=1}^{n}m_{i}^{2}a^{\prime\prime}_{i}\end{array}\right]\,,$
(4)

with $a^{\prime\prime}_{i}=a^{\prime\prime}(-m_{i}/\phi)$, where
$a^{\prime\prime}(u)=d^{2}a(u)/du^{2}$.

The weighted least squares problem in (3) is typically solved either through
the method of normal equations, which involves a Cholesky decomposition of
$X^{\top}W^{(j)}X$, or through the QR decomposition of $(W^{(j)})^{1/2}X$
(see, e.g. Golub and Van Loan, 2013, Section 5.3). Despite that the QR
approach requires more computations than the method of normal equations, the
former is more appealing in applications and is the default choice in popular
least squares software because it can solve, in a numerically stable manner, a
wider class of least squares problems. In particular, the method of normal
equations can numerically break down with $(W^{(j)})^{1/2}X$ that are not
particularly close to being numerically rank deficient, while the QR approach
solves a “nearby” least squares problem (see, Golub and Van Loan, 2013,
Section 5.3 for an extensive analysis).

ML estimation of $\phi$ can then take place by solving $s_{\phi}=0$ after
evaluating $q_{i}$ at the ML estimates $\hat{\beta}$. This can be done through
the Fisher scoring iteration

$\phi^{(j+1)}:=\phi^{(j)}\left\\{1+\phi^{(j)}\frac{\sum_{i=1}^{n}(\hat{q}_{i}-\rho_{i}^{(j)})}{\sum_{i=1}^{n}m_{i}^{2}a_{i}^{\prime\prime(j)}}\right\\}\,,$
(5)

where $\hat{q}_{i}$ is $q_{i}$ evaluated at $\hat{\beta}$. McCullagh and
Nelder (1989) recommend against estimating $\phi$ using ML, and instead
propose the moment estimator
$(n-p)^{-1}\sum_{i=1}^{n}\hat{w}_{i}(\hat{z}_{i}-\hat{\eta}_{i})^{2}$ where
$\hat{w}_{i}$, $\hat{z}_{i}$, and $\hat{\eta_{i}}$ are $w_{i}$, $z_{i}$, and
$\eta_{i}$, respectively, evaluated at $\hat{\beta}$. The moment estimator of
$\phi$ is considered to have less bias that the ML estimator.

###  2.3 Bounded-memory procedures for least squares

Miller (1992) proposed a procedure to solve least squares problems using the
QR approach, which does not require keeping the full model matrix $X$ and
response vector $Y$ in memory, and operates by sequentially bringing only a
fixed-size chunk of data in memory. This is particularly useful in getting a
numerically stable least squares solution in cases where the model matrix has
too many rows to fit in memory. We briefly review the method in Miller (1992).

Consider the least squares problem $A\psi=b$, where $A$ is an $n\times p$
matrix, with $n>p$, and $b$ is a $n$-dimensional vector. The least squares
solution $\hat{\psi}$ for $\psi$ can be found by first computing the QR
decomposition

$A=QR=\left[\begin{array}[]{cc}Q_{1}&Q_{2}\end{array}\right]\left[\begin{array}[]{c}\bar{R}\\\
0_{(n-p)\times p}\end{array}\right]\,,$

where $Q_{1}$ and $Q_{2}$ are $n\times p$ and $n\times(n-p)$ matrices,
respectively, $Q$ is an orthogonal matrix (i.e. $Q^{\top}=Q^{-1}$), $\bar{R}$
is a $p\times p$ upper triangular matrix, and $0_{u\times v}$ denotes a
$u\times v$ matrix of zeros. Then, $\hat{\psi}$ is found by using back-
substitution to solve $\bar{R}\hat{\psi}=\bar{b}$, where
$\bar{b}=Q_{1}^{\top}b$. To describe the proposal in Miller (1992), let
$A_{k}$ be the $k$th chunk of $c\leq n$ rows of $A$, and $b_{k}$ be the $k$th
chunk of $c$ elements of $b$, and denote $A_{1:k}$ the first $k$ chunks of
$A$. The last chunks may have $c$ or less observations. Suppose that the QR
decomposition of $A_{1:k}$ is $Q_{k}R_{k}$. Then, in light of a new chunk
$A_{k+1}$, the QR decomposition can be updated as

$\left[\begin{array}[]{c}A_{1:k}\\\
A_{k+1}\end{array}\right]=\underbrace{\left[\begin{array}[]{ccc}Q_{1,k}&0_{ck\times
k}&Q_{2,k}\\\ 0_{c\times
p}&I_{c}&0_{c\times(ck-p)}\end{array}\right]G_{1}\cdots
G_{cp}}_{Q_{k+1}}\overbrace{G_{cp}^{\top}\cdots
G_{1}^{\top}\left[\begin{array}[]{c}\bar{R}_{k}\\\ A_{k+1}\\\ 0_{(n-p)\times
p}\end{array}\right]}^{R_{k+1}}\,,$

where $G_{1},\ldots,G_{cp}$ is the set of the Givens rotation matrices (see,
Golub and Van Loan, 2013, Section 5.1.8) required to eliminate $A_{k+1}$ in
the right hand side. By the orthogonality of Givens rotation matrices,
$Q_{k+1}R_{k+1}$ is a valid QR decomposition. The only values that are needed
for computing the least squares estimates $\hat{\psi}$ are $\bar{R}$ and
$\bar{b}$. Hence, their current values after a new chunk arrives are typically
the only objects that are updated and kept in computer memory. This is useful
for large $n$, as the need of having an $n\times p$ matrix in memory ($X$) is
swapped by only keeping $p(p+1)/2+p$ real numbers in memory at any given time.

###  2.4 Iteratively re-weighted least squares in chunks

The approach of Miller (1992) can readily be used for (3) by replacing $A$ by
$(W^{(j)})^{1/2}X$ and $b$ by $z^{(j)}$. This is possible because the diagonal
entries of $W$ and the components of $z$ only depend on the corresponding
components of $\eta=X\beta$, and hence, can be computed in a chunk-wise
manner. Specifically, the $k$th chunk of the diagonal of $W^{(j)}$, and the
$k$th chunk of $z^{(j)}$ can be computed by just bringing in memory $X_{k}$,
$Y_{k}$, and $m_{k}$, computing the current value of the $k$th chunk of the
linear predictor as $X_{k}\beta^{(j)}$, and using that to compute the $k$th
chunks of $\mu^{(j)}$, $d^{(j)}$ and $v^{(j)}$. This is the procedure that the
biglm R package (Lumley, 2020) implements.

##  3 Bias reduction and maximum penalized likelihood

###  3.1 Adjusted score equations

Consider the adjusted score equations

$\displaystyle 0_{p}$ $\displaystyle=\frac{1}{\phi}X^{\top}W\left\\{z+\phi
H(b_{1}\xi+b_{2}\lambda)\right\\}\,,$ (6) $\displaystyle 0$
$\displaystyle=\frac{1}{2\phi^{2}}\sum_{i=1}^{n}(q_{i}-\rho_{i})+c_{1}\left(\frac{p-2}{2\phi}+\frac{\sum_{i=1}^{n}m_{i}^{3}a_{i}^{\prime\prime\prime}}{2\phi^{2}\sum_{i=1}^{n}m_{i}^{2}a_{i}^{\prime\prime}}\right)\,,$
(7)

where $0_{p}$ is a $p$-vector of zeros, $\xi=(\xi_{1},\ldots,\xi_{n})^{\top}$,
and $\lambda=(\lambda_{1},\ldots,\lambda_{n})^{\top}$ with

$\displaystyle\xi_{i}=\frac{d_{i}^{\prime}}{2d_{i}w_{i}}\,,\quad\text{and}\quad\lambda_{i}=\frac{1}{2}\left\\{\frac{d_{i}^{\prime}}{d_{i}w_{i}}-\frac{v_{i}^{\prime}}{m_{i}d_{i}}\right\\}\,,\quad$

and $a^{\prime\prime\prime}_{i}=a^{\prime\prime\prime}(-m_{i}/\phi)$, where
$a^{\prime\prime\prime}(u)=d^{3}a(u)/du^{3}$,
$d_{i}^{\prime}=d^{2}\mu_{i}/d\eta_{i}^{2}$, and
$v^{\prime}_{i}=d^{2}V(\mu_{i})/d\mu_{i}^{2}$. In the above expressions, $H$
is the diagonal matrix, whose diagonal is the diagonal of the “hat” matrix
$X(X^{\top}WX)^{-1}X^{\top}W$. The ML estimators of $\beta$ and $\phi$ are
obtained by solving (6) and (7) with $b_{1}=b_{2}=c_{1}=0$. Kosmidis et al.
(2020) show that mBR estimators of $\beta$ and $\phi$ can be obtained by
solving (6) and (7) with $b_{1}=1$, $b_{2}=0$, and $c_{1}=1$. On the other
hand, direct differentiation of the Jeffreys’-prior penalized likelihood shows
that the mJPL estimators for
$\beta=\arg\max\\{\ell(\beta)+\log|X^{\top}WX|/2\\}$ in Kosmidis and Firth
(2021) can be obtained by solving (6) for $b_{1}=b_{2}=1$.

For binomial-response GLMs, Kosmidis and Firth (2021) show that the mJPL
estimates of $\beta$ has always finite components for a wide range of well-
used link functions, including logit, probit and complementary log-log, even
in cases where the ML estimate have infinite components. The finiteness
property of the mJPL estimator is attractive for applied work, because it
ensures the stability of all numerical and inferential procedures, even when
the data set is separated and infinite ML estimates occur (see, Mansournia et
al., 2018, for a recent review of the practical issues associated with
infinite estimates in logistic regression). The same holds for any power $t>0$
of the Jeffreys’ prior, in which case
$\lambda_{i}=t\\{d_{i}^{\prime}/(d_{i}w_{i})-v_{i}^{\prime}/(m_{i}d_{i})\\}/2$
in (6).

The mJPL estimators, though, do not necessarily have better asymptotic bias
properties than the ML estimator for all GLMs . For GLMs that are full
exponential families, such as logistic regression for binomial data and log
linear models for Poisson counts, the mJPL estimators with $t=1$ and the mBR
estimators coincide. This has been shown in Firth (1993), and is also
immediately apparent by (6). For canonical links $d_{i}=v_{i}$,
$w_{i}=m_{i}d_{i}$, and, by the chain rule
$v_{i}^{\prime}=d_{i}^{\prime}/d_{i}$, hence $\lambda_{i}=0$. Hence the mBR
estimator for logistic regression not only has bias of smaller asymptotic
order than what the ML estimator generally has, but their components also
always take finite values.

We should note that the first-order asymptotic properties expected by the ML,
mJPL and mBR estimators of $\beta$ and $\phi$ are preserved for any
combination of $b_{1}$, $b_{2}$ and $c_{1}$ in (6) and (7). For example, mBR
estimators of $\beta$ can be obtained for $b_{1}=1$, $b_{2}=0$ and $c_{1}=0$,
effectively mixing mBR for $\beta$ with ML for $\phi$. This is the result of
the orthogonality of $\beta$ and $\phi$ (Cox and Reid, 1987); the detailed
argument is a direct extension of the argument in Kosmidis et al. (2020,
Section 4) for mixing adjusted score equations to get mBR for $\beta$ and
estimators of $\phi$ that have smaller median bias. Another option is to mix
the adjusted score equations (6) for $\beta$ with the estimating function

$\phi=\frac{1}{n-p}\sum_{i=1}^{n}w_{i}(z_{i}-\eta_{i})^{2}\,,$ (8)

which gives rise to the moment estimator of $\phi$, once $\mu_{i}$ is
evaluated at an estimator for $\beta$.

###  3.2 Iteratively re-weighted least squares for solving adjusted score
equations

Using similar arguments as for ML estimation, we can define the following IWLS
update to compute the ML, the mBR, and the mJPL estimators of $\beta$
depending on the value of $b_{1}$, $b_{2}$ in (6). That update has the form

$\beta^{(j+1)}:=\left(X^{\top}W^{(j)}X\right)^{-1}X^{\top}W^{(j)}\left(z^{(j)}+\phi^{(j)}H^{(j)}\kappa^{(j)}\right)\,,$
(9)

where $\kappa=b_{1}\xi+b_{2}\lambda$. The value of $\phi$, which is required
for implementing mBR and mJPL, can be found via the quasi-Fisher scoring step
for solving (7), which using (4), takes the form

$\phi^{(j+1)}:=\phi^{(j)}\left[1+\phi^{(j)}\frac{\sum_{i=1}^{n}(q_{i}^{(j)}-\rho_{i}^{(j)})}{\sum_{i=1}^{n}m_{i}^{2}a_{i}^{\prime\prime(j)}}+c_{1}\phi^{(j)}\left\\{\frac{\sum_{i=1}^{n}m_{i}^{3}a_{i}^{\prime\prime\prime(j)}}{(\sum_{i=1}^{n}m_{i}^{2}a_{i}^{\prime\prime(j)})^{2}}+\phi^{(j)}\frac{p-2}{\sum_{i=1}^{n}m_{i}^{2}a_{i}^{\prime\prime(j)}}\right\\}\right]$
(10)

The candidate value for the dispersion parameter when one of $b_{1}$ or
$b_{2}$ is non-zero, can be computed using (8) at $\beta^{(j)}$ in every
iteration. Kosmidis et al. (2020) have already derived the special case of (9)
and (10) for mBR estimation, that is for $b_{1}=1$, $b_{2}=0$ and $c_{1}=1$.
Convergence can be declared if
$\lVert\beta^{(j+1)}-\beta^{(j)}\rVert_{\infty}<\epsilon$ and
$\lVert\phi^{(j^{\prime}+1)}-\phi^{(j)}\rVert_{\infty}<\epsilon$ for some
small $\epsilon>0$.

###  3.3 Adjusted score estimation in chunks

When either (8) or (5) is used, the updates for $\phi$ can be performed in a
chunkwise manner once $\beta^{(j)}$ has been computed. In fact, use of (8) has
the advantage of requiring only $w^{(j)}$ and $z^{(j)}$, which have already
been computed after performing the IWLS update (6).

Nevertheless, the IWLS update (9) for $\beta$ cannot be readily performed in a
chunkwise manner because the diagonal of $H$ cannot be computed at
$\beta^{(j)}$ by just bringing in memory chunks of $X$, $Y$ and $m$, and the
current estimates. In particular, the $i$th diagonal element of $H$ is
$h_{i}=x_{i}^{\top}(X^{\top}WX)^{-1}x_{i}w_{i}$, and, hence, its computation
requires the inverse of the expected information matrix $X^{\top}WX$. In what
follows, we present two alternatives for computing adjusted score estimates
for $\beta$ in a chunkwise manner.

###  3.4 Two-pass implementation

The direct way to make the update (6) for $\beta$ possible in a chunk-wise
manner comes from the considering the projection that is being made. The left
plot of Figure 1 shows how update (6) projects the current value of the vector
$W^{1/2}(z+\phi H\kappa)$ onto the column space of the current value of
$W^{1/2}X$.

Figure 1: Demonstration of the IWLS update (9). All quantities in the figure
should be understood as being pre-multiplied by $W^{1/2}$ and evaluated at
$\beta^{(j)}$ and $\phi^{(j)}$. The left figure shows the addition of $\phi
H\kappa$ to the maximum likelihood working variates vector $z$, and the
subsequent projection onto $\mathcal{C}$ (the column space of $W^{1/2}X$) that
gives the updated value for the mBR or mJPL estimates $\tilde{\beta}$. The
right figure shows the projection of $z$, and the subsequent projection of
$\phi H\kappa$ onto $\mathcal{C}$, and the addition of the projected vectors
that gives the updated value for the mBR or mJPL estimates
$\tilde{\beta}$.$z$$\phi H\kappa$$z+\phi H\kappa$$X\tilde{\beta}$$z$$\phi
H\kappa$$X\hat{\beta}$$Xv$$X\tilde{\beta}$$\mathcal{C}$$\mathcal{C}$

The right of Figure 1, then shows how we can achieve exactly the same
projection in two passes through the data: i) project the current value of
$W^{1/2}z$ onto the column space of the current value of $W^{-1/2}X$ using QR
decomposition, and ii) project the current value of $\phi W^{1/2}H\kappa$ onto
the column space of the current value of $W^{-1/2}X$. Adding up the
coefficient vectors from the two projections returns the required updated
value of $\beta$. Section 2.4 describes how the projection in step i) can be
done in a chunkwise manner. Then, given that the current value of $\bar{R}$ is
available after the completion of the incremental QR decomposition from the
first pass through the data, the current value of $h_{1},\ldots,h_{n}$ in step
ii) can be computed in a chunkwise manner through another pass through the
data.

###  3.5 One-pass implementation

An alternative way to make the update (6) for $\beta$ possible in a chunk-wise
manner is to change the IWLS update into

$\beta^{(j+1)}:=\left(X^{\top}W^{(j)}X\right)^{-1}X^{\top}W^{(j)}(z^{(j)}+\phi^{(j)}H^{(j-1)}\kappa^{(j)})\,.$
(11)

Iteration (11) has the correct stationary point, that is the solution of the
adjusted score equations (6), and can be performed in a chunkwise manner
following the descriptions in Section 2.4. This is because the $i$th diagonal
element of $H^{(j-1)}$,
$h_{i}=x_{i}^{\top}(X^{\top}W^{(j-1)}X)^{-1}x_{i}w_{i}^{(j-1)}$, depends on
the weight $w_{i}^{(j-1)}$ at the previous value of $\beta$, which can be
recomputed using knowledge of $\beta^{(j-1)}$ and $x_{i}$ only, and the
inverse of the expected information $(X^{\top}W^{(j-1)}X)^{-1}$ from the
previous iteration, which can be computed using $\bar{R}^{(j-1)}$ that results
from the incremental QR decomposition at the $(j-1)$th iteration.

In addition to the starting value for $\beta^{(0)}$ that iteration (3)
requires, iteration (11) also requires starting values for $h_{i}$; a good
starting value is $h_{i}^{(0)}=p/n$, which corresponds to a balanced model
matrix $X$.

###  3.6 Memory requirements and computational complexity

The direct implementation of IWLS for mBR and mJPL, as done in popular R
packages such as brglm2 and logistf, requires $O(np+p^{2})$ memory, as ML
does. On the other hand, the chunkwise implementations require $O(cp+p^{2})$
memory, where $c$ is the user-specified chunk size as in Section 2.3. The
computational cost of all implementations using the QR decomposition remains
at $O(np^{2}+p^{3})$.

The two-pass implementation has almost twice the iteration cost of the one-
pass implementation, because two passes through all observations is required
per IWLS iteration. However, the two-pass implementation reproduces exactly
iteration (9) where the adjusted score, rather that just part of it as in
(11), is evaluated at the current parameter values. From our experience with
both implementations the one-pass one tends to require more iterations to
converge to the same solution, with the two implementations having the same
computational complexity per iteration. In addition, the two-pass
implementation requires starting values only for $\beta$, while the one-pass
implementation requires starting values for both $\beta$ and
$h_{1},\ldots,h_{n}$.

##  4 Demonstration: Diverted US flights in 2000

We demonstrate here the one-pass and two-pass implementations using data on
all $5,683,047$ commercial flights within the USA in 2000. The data set is
part of the data that was shared during The Data Exposition Poster Session of
the Graphics Section of the Joint Statistical meetings in 2009, and is
available at the Harvard Dataverse data repository (air, 2008).

Table 1: Estimates and estimated standard errors for selected parameters of model (12) through ML, mBR and mJPL with one- and two-pass IWLS implementations. The table also reports the elapsed time, number of iterations, and average time per iteration for each method and implementation. The parameters corresponding to reference categories are set to $0$. |  | ML | mBR | mJPL
---|---|---|---|---
|  | one-pass | one-pass | one-pass | two-pass | one-pass | two-pass
$\alpha$ |  | $-5.52$ | $-6.38$ | $-4.10$ | $-4.10$ | $-4.12$ | $-4.12$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
⋮ |  |  |  |  |  |  |
$\delta_{1}$ | AA | $2.49$ | $3.36$ | $1.08$ | $1.08$ | $1.09$ | $1.09$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{2}$ | AQ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
$\delta_{3}$ | AS | $2.59$ | $3.45$ | $1.17$ | $1.17$ | $1.18$ | $1.18$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{4}$ | CO | $2.33$ | $3.19$ | $0.91$ | $0.91$ | $0.92$ | $0.92$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{5}$ | DL | $2.36$ | $3.22$ | $0.94$ | $0.94$ | $0.95$ | $0.95$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{6}$ | HP | $2.18$ | $3.05$ | $0.76$ | $0.76$ | $0.78$ | $0.78$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{7}$ | NW | $2.44$ | $3.30$ | $1.02$ | $1.02$ | $1.04$ | $1.04$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{8}$ | TW | $2.36$ | $3.23$ | $0.94$ | $0.94$ | $0.96$ | $0.96$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{9}$ | UA | $2.35$ | $3.22$ | $0.94$ | $0.94$ | $0.95$ | $0.95$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{10}$ | US | $2.54$ | $3.41$ | $1.12$ | $1.12$ | $1.14$ | $1.14$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
$\delta_{11}$ | WN | $2.41$ | $3.28$ | $0.99$ | $0.99$ | $1.01$ | $1.01$
|  | $(4.89)$ | $(53.18)$ | $(0.35)$ | $(0.35)$ | $(0.36)$ | $(0.36)$
⋮ |  |  |  |  |  |  |
Time (sec) | $216.78$ | $276.82$ | $220.52$ | $380.75$ | $217.59$ | $378.03$
Iterations | $15$ | $20$ | $12$ | $12$ | $12$ | $12$
Time/Iteration (sec) | $14.45$ | $13.84$ | $18.38$ | $31.73$ | $18.13$ | $31.50$

Suppose that interest is in modelling the probability of diverted US flights
in 2000 in terms of the departure date attributes, the scheduled departure and
arrival times, the coordinates and distance between departure and planned
arrival airports, and the carrier. Towards this goal, we assume that the
diversion status of the $i$th flight is a Bernoulli random variable with
probability modelled as

$\Phi^{-1}(\pi_{i})=\alpha+\beta^{\top}M_{i}+\gamma^{\top}W_{i}+\delta^{\top}C_{i}+\zeta_{(d)}T_{(d),i}+\zeta_{(a)}T_{(a),i}+\rho
D_{i}+\psi_{(d)}^{\top}L_{(d),i}+\psi_{(a)}^{\top}L_{(a),i}\,,$ (12)

and that diversions are conditionally independent given the covariate
information that appears in the above expression. The covariate information
consists of $M_{i}$, which is a vector of $12$ dummy variables characterizing
the calendar month that the $i$th flight took place, $W_{i}$, which is a
vector of $7$ dummy variables characterizing the week day that the $i$th
flight took place, $C_{i}$, which is a dummy variable with $11$ levels,
characterizing the carrier (out of $11$ carriers) of the $i$th flight,
$T_{(d),i}$ and $T_{(a),i}$, which are the planned departure and arrival times
in a 24 hour format, $D_{i}$, which is the distance between the origin and
planned destination airport, respectively and $L_{(d),i}$ and $L_{(a),i}$,
which are the $(x,y,z)$ coordinates of the departure and arrival airport,
respectively, computed from longitude (lon) and latitude (lat) information as
$x=\cos(\rm lat)\cos(\rm lon)$, $y=\cos(\rm lat)\sin(\rm lon)$, and
$z=\sin(\rm lat)$.

For identifiability reasons, we fix $\beta_{1}=0$ (January as a reference
category), $\gamma_{1}=0$ (Monday as a reference category) and $\delta_{2}=0$
(carrier AQ as a reference category). Ignoring the columns corresponding to
those parameters, the model matrix for the model with linear predictor as in
(12) has dimension $5683047\times 37$, which requires about $1.6$GB of memory.

Despite that the memory requirements for this data set are manageable by the
available memory of modern laptops, testing for data separation by solving,
for example, the Konis (2007) linear programs is computationally demanding.
For example, the detect_infinite_estimates() method of the detectseparation R
package did not complete even after 12 hours of computation on a 2021 MacBook
Pro with an Apple M1 Max chip and 64 GB RAM.

Table 1 shows the ML estimates of the parameters $\alpha$ and $\delta$ of
model (12) after 15 and 20 IWLS iterations using the bounded-memory procedures
in Section 2.3 as implemented in the biglm R package (Lumley, 2020), using
chunks of $c=10,000$ observations and $\epsilon=10^{-3}$. The estimates and
estimated standard errors for $\alpha$ and the components of $\delta$ grow in
absolute value with the number of IWLS iterations, which is typically the case
when the maximum likelihood estimates are infinite (see, for example, Lesaffre
and Albert, 1989, for such behaviours in the case of multinomial logistic
regression). The ML estimates for the other parameters do not change in the
reported accuracy as we move from 15 to 20 IWLS iterations; see Table
LABEL:tab:diverted_model_all of the Supplementary Materials document for
estimates for all parameters. In contrast, mBR and mJPL return finite
estimates for all parameters by declaring convergence before reaching the
limit of allowable iterations, for both their one- and two-pass
implementations. As expected by the discussion in Section 3.6, the one-pass
implementations require about $58\%$ of the time that the two-pass
implementations require. No memory issues have been encountered obtaining the
ML, mJPL and mBR fits, even when the fits were re-computed on a Ubuntu virtual
machine with 2 cores and 2GB and 4GB of RAM, where also no swapping took
place. On the other hand, the mBR and mJPL fits could not be obtained using
the brglm2 R package, because available physical memory was exhausted.

The observations made in this example highlight that, even for large data
sets, use of estimation methods that return estimates in the interior of the
parameter space is desirable, especially so when the computational cost for an
IWLS iteration for ML scales linearly with the number of observations.

##  5 Maximum Jeffreys’ penalized likelihood in high-dimensional logistic
regression

###  5.1 Preamble

Firth (1993) shows that for logistic regression with fixed-$p$ asymptotics
mJPL and mBR coincide. More recently, Sur and Candès (2019) and Kosmidis and
Firth (2021) have empirically found that mBR still achieves a substantial
reduction in the persistent bias of the ML estimator in specific high-
dimensional logistic regression problems under a more extreme asymptotic
framework than Firth (1993) used to develop mBR.

In particular, Candès and Sur (2020, Section 1.5) prove that there exists a
sharp phase transition about when the ML estimate has infinite components in
high-dimensional logistic regression with linear predictors
$\eta_{i}=\alpha+x_{i}^{\top}\beta$, where $x_{i}$ is a $p$-vector generated
from a Normal distribution with mean $0_{p}$ and unknown, non-singular
variance-covariance matrix $\Sigma$, $p/n\to\kappa\in(0,1)$,
$x_{1},\ldots,x_{n}$ are generated independently from each other, and assuming
that ${\mathop{\rm var}}(x_{i}^{\top}\beta)\to\gamma_{0}^{2}$. That phase
transition is useful because it characterizes when the corrective method
developed in Sur and Candès (2019) can be applied. That corrective method re-
scales the ML estimate by a factor from the solution of the system of
nonlinear equations in Sur and Candès (2019, expression (5)). As Sur and
Candès (2019) clearly note, that system of nonlinear equations does not have a
unique solution if $(\kappa,\gamma)$ are in the region where the ML estimate
has finite components, and their approach cannot be used when at least one of
the ML estimates diverges. In contrast, Kosmidis and Firth (2021) prove that
mJPL for logistic regression always delivers finite estimates under the sole
assumption that the model matrix has full rank.

Figure 2: The phase transition curves derived using Candès and Sur (2020,
Theorem 2.1) for $\rho^{2}=0$ (top) and $\rho^{2}=0.75$ (bottom) for the
computer experiment in Section 5, overlaid with scatterplots of ML (red) and
mJPL estimates of $\beta$ (orange) from 5 independent realizations of $y$ and
$X$ vs the true $\beta$, at $30$ $(\kappa,\gamma)$ points, for $n=2000$, and
with $\beta^{*}$ set to an equi-spaced grid of length $p=\lceil n\kappa\rceil$
between $-10$ and $10$. The ML estimates asymptotically exist in the white
regions, and do not asymptotically exist in the grey regions. The black lines
are intercept zero and slope one reference lines, and the green lines are
estimated from a simple linear regression of the estimates on the truth.

Motivated by this fact, in this section, we empirically investigate the limits
of mJPL in terms of frequentist performance in the setting of Candès and Sur
(2020) at either side of the phase transition curve derived therein. We
provide evidence that mJPL performs exceptionally well in terms of bias in the
region where the ML estimate is expected to exist, and tends to over-correct
in the region where the ML estimate has infinite components. The empirical
evidence also suggests that the amount of over-correction can be accurately
characterized in terms of $\kappa$, $\gamma$ and the size of the intercept
parameter, for a wide range of values of $\kappa$ and $\gamma$. This leads to
a fruitful conjecture about adjusting the mJPL estimates post-fit to almost
recover their performance even in the region where the ML estimates do not
exist, where the methods of Sur and Candès (2019) do not apply.

Figure 3: Figure 2, with the mJPL estimator $\tilde{\beta}$ adjusted to
$\kappa\gamma(1-\rho^{2})^{-1/2}\tilde{\beta}$ in the grey region, where the
ML estimates do not exist asymptotically.

###  5.2 Simulation experiment

We use the asymptotic framework used for producing Candès and Sur (2020,
Figure 1(b)). For generating a data set under that asymptotic framework, we
need to specify $\kappa\in(0,1)$, $n$, $\rho^{2}\in[0,1)$, $\gamma>0$, and an
initial parameter vector $\beta^{*}$. We form an $n\times p$ matrix $X$ by
generating its entries from a collection of $np$ independent and identically
distributed standard normal random variables, where $p=\lceil n\kappa\rceil$.
We define $\alpha=\rho\gamma$, and $\gamma_{0}=\gamma\sqrt{1-\rho^{2}}$, and
rescale $\beta^{*}$ to $\beta=\gamma_{0}\beta^{*}/\lVert\beta^{*}\rVert_{2}$,
so that $\lVert\beta\rVert_{2}=\gamma_{0}$. In this way, ${\mathop{\rm
var}}(x_{i}^{\top}\beta)=\gamma_{0}^{2}$, which is a sufficient condition for
the developments in Candès and Sur (2020). We then generate
$y_{1},\ldots,y_{n}$ independently with
$Y_{i}\sim\text{Bernoulli}(1/(1+e^{-\eta_{i}}))$, where
$\eta_{i}=\alpha+x_{i}^{\top}\beta$ $(i=1,\ldots,n)$.

Figure 2 shows the phase transition curves derived using Candès and Sur (2020,
Theorem 2.1) for $\rho^{2}=0$ (top) and $\rho^{2}=0.75$ (bottom). At each from
a set of $30$ points on the $(\kappa,\gamma)$ plane ($\kappa\in(0,0.6)$ and
$\gamma\in(0,20)$), and for $n=2000$, $\beta^{*}$ is set to an equi-spaced
grid of length $p$ between $-10$ and $10$. Then, $5$ samples of $X$ and $y$
are drawn independently as detailed above at each point. The ML estimates of
$\alpha$ and $\beta$ are computed for the samples at the region where the ML
estimates asymptotically exist, and the mJPL estimates are computed everywhere
using the two-pass implementation with $c=1000$, $\epsilon=10^{-3}$ (see
Section 3.2), and allowing for up to $250$ IWLS iterations. The mJPL estimates
are computed starting at zero for all parameters, and are used as starting
values for the ML estimates whenever those are computed.

The scatterplots of the estimates versus the true $\beta$ vectors are overlaid
in Figure 2 at each of the $30$ points on the $(\kappa,\gamma)$ plane. As is
evident, in the region where the ML estimates exist, the mJPL estimates
$\tilde{\beta}$ illustrate excellent performance in recovering the true
signal, even when the ML estimates over-estimate it. Their performance is
similar to that of the re-scaled ML estimator in Sur and Candès (2019). On the
other hand, in the region where the ML estimates do not exist, mJPL appears to
overshrink towards zero moving further away from the true signal as $\kappa$
and $\gamma$ increase.

A simple exploratory regression analysis of mJPL estimates versus the true
values in the region where the ML estimates do not exist results in strong
evidence that a post-hoc multiplicative adjustment of the mJPL estimates to
$\kappa\gamma(1-\rho^{2})^{-1/2}\tilde{\beta}$ can recover the signal. Figure
3 shows the resulting estimates, confirming that the adjusted estimator
recovers the signal to high accuracy, in the region where the ML estimates do
not exist. To our knowledge, there has been no other viable proposal for
signal recovery in that region.

Table 2: Runtime summaries for the two-pass IWLS implementation of mJPL (average computational time in seconds, and minimum, average and maximum required number of iterations) across the $5$ samples at each of the $(\kappa,\gamma)$ points in Figure 3, for $n=2000$ and $\rho=0$. $\kappa$ | $\gamma$ | $p$ | average | minimum | average | maximum
---|---|---|---|---|---|---
|  |  | time (sec) | iterations | iterations | iterations
$0.01$ | $1.0$ | $20$ | $0.04$ | $4$ | $4.0$ | $4$
$0.01$ | $8.0$ | $20$ | $0.10$ | $8$ | $8.6$ | $9$
$0.01$ | $15.0$ | $20$ | $0.14$ | $11$ | $11.2$ | $12$
$0.05$ | $4.5$ | $100$ | $0.87$ | $7$ | $7.2$ | $8$
$0.05$ | $11.5$ | $100$ | $1.29$ | $12$ | $13.0$ | $14$
$0.05$ | $18.5$ | $100$ | $2.58$ | $17$ | $26.2$ | $36$
$0.15$ | $1.0$ | $300$ | $3.17$ | $4$ | $4.0$ | $4$
$0.15$ | $8.0$ | $300$ | $27.49$ | $27$ | $38.4$ | $55$
$0.15$ | $15.0$ | $300$ | $24.22$ | $27$ | $31.8$ | $40$
$0.22$ | $8.0$ | $440$ | $45.82$ | $24$ | $29.6$ | $42$
$0.22$ | $15.0$ | $440$ | $32.11$ | $18$ | $21.4$ | $25$
$0.25$ | $4.5$ | $500$ | $74.22$ | $23$ | $40.0$ | $91$
$0.25$ | $11.5$ | $500$ | $41.39$ | $18$ | $21.2$ | $29$
$0.25$ | $18.5$ | $500$ | $37.15$ | $17$ | $19.0$ | $22$
$0.30$ | $8.0$ | $600$ | $59.52$ | $17$ | $21.4$ | $25$
$0.30$ | $15.0$ | $600$ | $51.22$ | $16$ | $18.0$ | $21$
$0.35$ | $1.0$ | $700$ | $21.58$ | $6$ | $6.6$ | $7$
$0.35$ | $4.5$ | $700$ | $101.74$ | $22$ | $29.0$ | $40$
$0.35$ | $11.5$ | $700$ | $63.17$ | $16$ | $17.0$ | $18$
$0.35$ | $18.5$ | $700$ | $56.39$ | $14$ | $15.4$ | $17$
$0.40$ | $8.0$ | $800$ | $73.69$ | $14$ | $15.6$ | $20$
$0.40$ | $15.0$ | $800$ | $65.75$ | $13$ | $13.8$ | $15$
$0.45$ | $4.5$ | $900$ | $102.83$ | $17$ | $18.0$ | $19$
$0.45$ | $11.5$ | $900$ | $80.31$ | $13$ | $13.0$ | $13$
$0.45$ | $18.5$ | $900$ | $69.78$ | $10$ | $11.8$ | $13$
$0.50$ | $8.0$ | $1000$ | $86.58$ | $11$ | $12.0$ | $13$
$0.50$ | $15.0$ | $1000$ | $83.12$ | $10$ | $11.8$ | $15$
$0.55$ | $4.5$ | $1100$ | $135.47$ | $13$ | $15.6$ | $18$
$0.55$ | $11.5$ | $1100$ | $96.42$ | $10$ | $10.8$ | $11$
$0.55$ | $18.5$ | $1100$ | $92.16$ | $10$ | $10.8$ | $12$

The Supplementary Materials document provides results for
$\rho\in\\{0,0.25,0.5,0.75,0.9\\}$, for both when $\beta^{*}$ is set to an
equi-spaced grid of length $p$ between $-10$ and $10$ (see Section
LABEL:sec:high-dim-set-a), and when $\beta^{*}$ is set to have $20\%$ of its
values set to $-10$, $20\%$ of its values set to $10$, and the remaining set
to $0$ (see Section LABEL:sec:high-dim-set-b), for $n\in\\{1000,2000,3000\\}$.
From the observations made through those computer experiments, and other
settings we have experimented with, not reported here, and for both smaller
and larger $n$, we are confident to state that the estimator
$\beta^{\dagger}=q(\kappa,\gamma,\rho)\tilde{\beta}$ with

$q(\kappa,\gamma,\rho)=\left\\{\begin{array}[]{ll}1\,,&\text{if the ML
estimate exists asymptotically}\\\
\frac{\kappa\gamma}{\sqrt{1-\rho^{2}}}\,,&\text{if the ML estimate does not
exist asymptotically}\end{array}\right.$ (13)

is effective in recovering the true signal, with its effectiveness
deteriorating as $\rho^{2}$ approaches $1$ (and, hence, as $\alpha$ increases
relative of $\lVert\beta\rVert_{2}$). This can be seen in the figures in
Supplementary Materials document for $\rho=0.9$. Figure 2 is a composite of
the top of Figure LABEL:fig:set-a-n-2000-rhosq-0 and Figure LABEL:fig:set-
a-n-2000-rhosq-0.75, and Figure 3 is a composite of the bottom of Figure
LABEL:fig:set-a-n-2000-rhosq-0 and Figure LABEL:fig:set-a-n-2000-rhosq-0.75.

###  5.3 Remarks on computational performance and extensions

The two-pass implementation of mJPL estimation converged in all cases within
$250$ IWLS iterations. We also tried the one-pass implementation, which either
required more iterations to converge than the two-pass, or did not converge
within $250$ iterations.

As a side note, Sur and Candès (2019, Supporting Information document, Section
D) reports that mJPL has been computationally infeasible for high dimensions,
with a runtime of approximately 10 minutes for $n=1000$ and $p=200$, and of
over 2.5 hours for $n=2000$ and $p=400$. Table 2 summarizes the performance of
the two-pass implementation of mJPL both in terms of runtime and number of
iterations for $\rho^{2}=0$. In notable contrast to the observations in Sur
and Candès (2019, Supporting Information document, Section D), the maximum
average runtime for mJPL we observed has been about 2 minutes and 15 seconds
for $n=2000$ and $p=1100$ ($\kappa=0.55$, $\gamma=4.5$), with no particular
care in the choice of starting values. The Supplementary Materials document
provides the summaries in Table 2, for $\rho^{2}=0$,
$n\in\\{1000,2000,3000\\}$, for both when $\beta^{*}$ is set to an equi-spaced
grid of length $p$ between $-10$ and $10$ (see Section LABEL:sec:high-dim-set-
a-perf), and when $\beta^{*}$ is set to have $20\%$ of its values set to
$-10$, $20\%$ of its values set to $10$, and the remaining set to $0$ (see
Section LABEL:sec:high-dim-set-b-perf).

For small $\kappa$, the incremental IWLS implementations of mJPL can be also
used to carry out the computer experiments of the current Section with no
memory problems even for larger $n$, and at least as long as the $O(cp+p^{2})$
memory requirements and $O(np^{2}+p^{3})$ computational complexity do not
become prohibitive for the hardware at hand. The $n=2000$ we used here ensures
that numerical results that support our conjecture can be easily reproduced in
a modern laptop with the code we provide in the Supplementary Material. For
larger $p$, it is worthwhile choosing reasonable starting values for mJPL.
These can be obtained through a few IWLS iterations for ML, potentially after
adjusting the responses to $(1-2\delta)y_{i}+\delta$, for some small
$\delta>0$.

##  6 Concluding remarks

We have developed two variants of IWLS that can estimate the parameters of
GLMs using adjusted score equations for mean bias reduction and maximum
Jeffreys’-prior penalized likelihood, for data sets that exceed computer
memory or even hard-drive capacity and are stored in remote databases. The two
procedures return the exact, not approximate, mBR and mJPL estimates, and they
have been used in Section 4 to obtain finite estimates of the parameters of a
probit regression model with $37$ parameters from about $5.5$ million
observations, where the ML estimates have been found to have infinite
components. The two-pass implementation has been also used in a large-scale
experiment to develop a simple adjustment to the maximum Jeffreys’-penalized
likelihood estimates that delivers excellent performance in terms of signal
recovery even in cases where estimates from ML and other recently-proposed
corrective methods do not exist.

A proof of the conjecture about the performance of
$q(\kappa,\gamma,\rho)\tilde{\beta}$ estimator in high-dimensional logistic
regression may be possible by adapting the framework and results in Salehi et
al. (2019) to the case where the likelihood is penalized by Jeffreys prior and
examining when the solution of the system of nonlinear equations in Salehi et
al. (2019, expression (6)) is close to the scaling factor
$q(\kappa,\gamma,\rho)$ in expression (13) that the numerical evidence points
to.

Current work focuses on reducing the cubic complexity on $p$, without
impacting the finiteness and bias reducing properties of mJPL, or its highly
predictable behaviour in high-dimensional logistic regression settings.
Drineas et al. (2012) on the approximation of the leverage seems relevant.

As Kosmidis et al. (2020) show, median bias reduction for $\beta$ can also be
achieved using an IWLS procedure after modifying the ML working variates. The
IWLS update for median bias reduction has the form

$\beta^{(j+1)}:=\left(X^{\top}W^{(j)}X\right)^{-1}X^{\top}W^{(j)}\left(z^{(j)}+\phi^{(j)}\left\\{H^{(j)}\xi^{(j)}+Xu^{(j)}\right\\}\right)\,.$

The particular form of the $p$-vector $u$ is given in Kosmidis et al. (2020,
expression (10)), and depends on the inverse $i_{\beta\beta}$ and hence on
$\bar{R}$. Since $Xu$ can be computed in a chunkwise manner for any given $u$,
it is possible to develop one- and two-pass implementations of the IWLS
procedure for median bias reduction by the same arguments as those used in
Section 3.5 and Section 3.4. These procedures are computationally more
expensive than the procedures for mBR and mJPL because each component of $u$
requires $O(np^{3})$ operations.

##  7 Supplementary materials

The Supplementary Material provides i) the Supplementary Material document
that is cross-referenced above and contains all numerical results and figures
from the case study of Section 4 and the computer experiment of Section 5, and
ii) R code to reproduce all numerical results and figures in the main text and
in the Supplementary Materials document. The code is organized in the two
directories diverted-flights and high-dim-logistic, for the case study of
Section 4 and the computer experiment of Section 5, respectively. The README
file in each directory provides specific instructions about have the results
can be reproduced, along with the specific versions of the contributed R
packages that have been used to produce the results. The biglm directory has a
port of the biglm R package (Lumley, 2020), which implements the one- and two-
pass IWLS variants for solving the bias-reducing adjusted score equations
(Firth, 1993) and for maximum Jeffreys’-penalized likelihood estimation
(Kosmidis and Firth, 2021). The Supplementary Material is available at
https://github.com/ikosmidis/bigbr-supplementary-material.

## References

  * air (2008) (2008).  2000.csv.bz2.  In Data Expo 2009: Airline on time data. Harvard Dataverse.  https://doi.org/10.7910/DVN/HG7NV7/YGU3TD.
  * Albert and Anderson (1984) Albert, A. and J. Anderson (1984).  On the existence of maximum likelihood estimates in logistic regression models.  Biometrika 71(1), 1–10.
  * Candès and Sur (2020) Candès, E. J. and P. Sur (2020).  The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression.  Annals of Statistics 48(1), 27–42.
  * Cox and Reid (1987) Cox, D. R. and N. Reid (1987).  Parameter orthogonality and approximate conditional inference.  Journal of the Royal Statistical Society, Series B: Methodological 49, 1–18.
  * Drineas et al. (2012) Drineas, P., M. Magdon-Ismail, M. W. Mahoney, and D. P. Woodruff (2012).  Fast approximation of matrix coherence and statistical leverage.  Journal of Machine Learning Research 13(111), 3475–3506.
  * Firth (1993) Firth, D. (1993).  Bias reduction of maximum likelihood estimates.  Biometrika 80(1), 27–38.
  * Golub and Van Loan (2013) Golub, G. H. and C. F. Van Loan (2013).  Matrix computations (Fourth edition ed.).  Johns Hopkins studies in the mathematical sciences. Baltimore: The Johns Hopkins University Press.
  * Green (1984) Green, P. J. (1984).  Iteratively reweighted least squares for maximum likelihood estimation, and some robust and resistant alternatives.  Journal of the Royal Statistical Society. Series B (Methodological) 46(2), 149–192.
  * Konis (2007) Konis, K. (2007).  Linear programming algorithms for detecting separated data in binary logistic regression models.  Ph. D. thesis, University of Oxford.
  * Kosmidis and Firth (2021) Kosmidis, I. and D. Firth (2021).  Jeffreys-prior penalty, finiteness and shrinkage in binomial-response generalized linear models.  Biometrika 108, 71–82.
  * Kosmidis et al. (2020) Kosmidis, I., E. C. Kenne Pagui, and N. Sartori (2020).  Mean and median bias reduction in generalized linear models.  Statistics and Computing (to appear) 30, 43–59.
  * Kosmidis et al. (2022) Kosmidis, I., D. Schumacher, and F. Schwendinger (2022).  detectseparation: Detect and Check for Separation and Infinite Maximum Likelihood Estimates.  R package version 0.3.
  * Lesaffre and Albert (1989) Lesaffre, E. and A. Albert (1989).  Partial separation in logistic discrimination.  Journal of the Royal Statistical Society. Series B (Methodological) 51(1), 109–116.
  * Lumley (2020) Lumley, T. (2020).  biglm: Bounded Memory Linear and Generalized Linear Models.  R package version 0.9-2.1.
  * Mansournia et al. (2018) Mansournia, M. A., A. Geroldinger, S. Greenland, and G. Heinze (2018).  Separation in logistic regression: Causes, consequences, and control.  American Journal of Epidemiology 187(4), 864–870.
  * McCullagh and Nelder (1989) McCullagh, P. and J. A. Nelder (1989).  Generalized Linear Models (2nd ed.).  London: Chapman and Hall.
  * Miller (1992) Miller, A. J. (1992).  Algorithm AS 274: Least Squares Routines to Supplement Those of Gentleman.  Applied Statistics 41(2), 458.
  * R Core Team (2023) R Core Team (2023).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Salehi et al. (2019) Salehi, F., E. Abbasi, and B. Hassibi (2019).  The Impact of Regularization on High-Dimensional Logistic Regression.  Number 1075. Red Hook, NY, USA: Curran Associates Inc.
  * Schwendinger et al. (2021) Schwendinger, F., B. Grün, and K. Hornik (2021).  A comparison of optimization solvers for log binomial regression including conic programming.  Computational Statistics.
  * Sur and Candès (2019) Sur, P. and E. J. Candès (2019).  A modern maximum-likelihood theory for high-dimensional logistic regression.  Proceedings of the National Academy of Sciences 116(29), 14516–14525.
  * Wedderburn (1976) Wedderburn, R. W. M. (1976).  On the existence and uniqueness of the maximum likelihood estimates for certain generalized linear models.  Biometrika 63(1), 27–32.


11institutetext: Helmholtz Institute Mainz, GSI Helmholtzzentrum für
Schwerionenforschung, Darmstadt, Johannes Gutenberg University, 55099 Mainz,
Germany 22institutetext: Institute for Nuclear Physics, Johannes Gutenberg
University, 55099 Mainz, Germany

# Exploring the neutron skin by hyperon–antihyperon production in
antiproton–nucleus interactions

Martin Christiansen 11 Falk Schupp 11 Patrick Achenbach 1122 Michael Bölting
11 Josef Pochodzalla 1122<EMAIL_ADDRESS>Marcell Steinen 11

###### Abstract

In this work we propose a new method to measure the evolution of the neutron
skin thickness between different isotopes. We consider antiproton–nucleus
interactions close to the production threshold of $\Lambda\overline{\Lambda}$
and $\Sigma^{-}\overline{\Lambda}$ pairs. At low energies,
$\Lambda\overline{\Lambda}$ pairs are produced in
$\overline{\text{p}}+\text{p}$ collisions, while
$\Sigma^{-}\overline{\Lambda}$ pairs can only be produced in
$\overline{\text{p}}+\text{n}$ interactions. Measuring these cross sections
provides information on the neutron skin thickness.

##  1 Introduction

Protons and Neutrons are the building blocks of conventional atomic nuclei.
Their distribution in nuclei are related to bulk properties of nuclear matter,
which are often expressed in terms of a nuclear equation of state (EoS).
Indeed, the isospin dependence of the EoS correlates strongly with the
distribution of neutrons in nuclei doi:10.1063/PT.3.4247 .

The neutron skin thickness, ${\Delta}r_{\text{np}}$, is defined as the
difference between the root-mean-squared (rms) point radii of the neutron and
proton density distributions. The focus of their studies often lies on neutron
rich doubly magic nuclei such as 40Ca, 48Ca, and 208Pb. However, it is the
evolution of the neutron skin and proton distributions along isotope chains
that provide important information for our understanding of the nuclear
structure over the whole nuclear chart.

In this work we propose a new method to measure the difference of the neutron
skin thickness between two different isotopes of a given element. We consider
antiproton–nucleus interactions close to the thresholds of
$\Lambda\overline{\Lambda}$ and $\Sigma^{-}\overline{\Lambda}$ pair
production. At low energies, $\Lambda\overline{\Lambda}$ pairs are produced in
$\overline{\text{p}}$+p collisions, while $\Sigma^{-}\overline{\Lambda}$ pairs
can only be produced in $\overline{\text{p}}$+n interactions. Measuring the
probabilities $p_{\Lambda\overline{\Lambda}}$ and
$p_{\Sigma^{-}\overline{\Lambda}}$ for the two processes for a reference
isotope (I) and a second isotope (II), allows to determine the double ratio

$DR=\frac{p^{II}_{\Sigma^{-}\overline{\Lambda}}\Big{/}p^{II}_{\Lambda\overline{\Lambda}}}{p^{I}_{\Sigma^{-}\overline{\Lambda}}\Big{/}p^{I}_{\Lambda\overline{\Lambda}}}.$
(1)

Within a simple geometrical picture we show that this ratio is strongly
related to the difference of the neutron skin thicknesses of the two
considered isotopes (I) and (II). In a second step, we explore this ratio with
the Gießen Boltzmann–Uehling–Uhlenbeck transport model for a xenon isotope
chain.

Figure 1: Illustration of $\Lambda\overline{\Lambda}$ and
$\Sigma^{-}\overline{\Lambda}$ pair production in antiproton–nucleus
interactions close to their thresholds. At low energies,
$\Lambda\overline{\Lambda}$ pairs are produced in $\overline{\text{p}}$+p
collisions, while $\Sigma^{-}\overline{\Lambda}$ pairs can only be produced in
$\overline{\text{p}}$+n interactions. We consider two isotopes (I) and (II)
which differ by an additional outer neutron layer with thickness
$\Delta_{\text{n}}$, indicated in pink.

##  2 Motivation of the method

Because of the strong absorption of antiprotons in nuclei the production of
hyperon–antihyperon pairs happens in the nuclear periphery. For simplicity we
consider a nearly central antiproton–nucleus collision (see left part of Fig.
1). Close to threshold, the probability $p^{I}_{\Lambda\overline{\Lambda}}$ to
produce a $\Lambda\overline{\Lambda}$ pair within the reference nucleus (I)
can be written as

$p^{I}_{\Lambda\overline{\Lambda}}=\kappa_{\Lambda\overline{\Lambda}}\cdot\frac{\rho(p)}{\rho(p)+\rho(n)}\cdot\frac{\sigma_{\Lambda\overline{\Lambda}}}{\sigma_{\text{tot}}}.$
(2)

Here $\sigma_{\text{tot}}$ denotes the total ${\overline{\text{p}}+A}$ cross
section and $\sigma_{\Lambda\overline{\Lambda}}$ is the elementary
$\overline{\text{p}}+{\text{p}}\rightarrow\Lambda\overline{\Lambda}$ cross
section. $\rho(p)$ and $\rho(n)$ denote the effective densities in the
periphery of the target nucleus of protons and neutrons, respectively. The
factor $\kappa_{\Lambda\overline{\Lambda}}$ describes the loss of
$\Lambda\overline{\Lambda}$ pairs by the absorption of the produced $\Lambda$
or $\overline{\Lambda}$ after their production. Similarly, the production
probability of $\Sigma^{-}\overline{\Lambda}$ pairs can be approximated by

$p^{I}_{\Sigma^{-}\overline{\Lambda}}=\kappa_{\Sigma^{-}\overline{\Lambda}}\cdot\frac{\rho(n)}{\rho(p)+\rho(n)}\cdot\frac{\sigma_{\Sigma^{-}\overline{\Lambda}}}{\sigma_{\text{tot}}}.$
(3)

Because of the large annihilation cross section of antibaryons in nuclei,
both, $\kappa_{\Lambda\overline{\Lambda}}$ and
$\kappa_{\Sigma^{-}\overline{\Lambda}}$ are dominated by the
$\overline{\Lambda}$ absorption. Therefore, we assume
$\kappa_{\Lambda\overline{\Lambda}}\approx\kappa_{\Sigma^{-}\overline{\Lambda}}\equiv\kappa_{I}$.

We now turn to a isotope (II) with a larger neutron number and hence a more
extended neutron distribution. For simplicity, we assume that the proton
distribution remains identical to the one of nucleus (I) and that the neutron
distribution is only extended by an additional halo $\Delta_{\text{n}}$ at the
surface, see pink area in the right part of Fig. 1. The production of
$\Lambda\overline{\Lambda}$ pairs is reduced by the absorption probability
$p_{\text{abs}}$ of the incident antiprotons within this additional neutron
skin $\Delta_{\text{n}}$:

$p^{II}_{\Lambda\overline{\Lambda}}=(1-p_{\text{abs}})\cdot\kappa_{II}\cdot\frac{\rho(p)}{\rho(p)+\rho(n)}\cdot\frac{\sigma_{\Lambda\overline{\Lambda}}}{\sigma_{\text{tot}}}.$
(4)

The absorption probability $p_{\text{abs}}$ can be expressed in terms of the
total $\overline{\text{p}}+\text{n}$ reaction cross section
$\sigma_{\overline{\text{p}}\text{n}}$ PhysRevD.98.030001  and the integrated
skin density $\int_{\Delta_{\text{n}}}\rho_{\text{n}}\,\text{d}r_{\text{n}}$
of the additional neutron skin of nucleus (II) with respect to the reference
nucleus (I):

$1-p_{\text{abs}}\approx\exp\Big{\\{}-\sigma_{\overline{\text{p}}\text{n}}\cdot\int_{\Delta_{\text{n}}}\rho_{\text{n}}\text{d}r_{\text{n}}\Big{\\}}$
(5)

Like the incoming antiprotons, the produced $\overline{\Lambda}$ are also
absorbed in the additional neutron layer. We use the factorisation ansatz
$\kappa_{II}=\kappa_{I}\cdot\kappa_{n}$. With this simplification we can
express Eq. (4) as

$p^{II}_{\Lambda\overline{\Lambda}}=\kappa_{n}\cdot(1-p_{abs})\cdot
p^{I}_{\Lambda\overline{\Lambda}}.$ (6)

The production of $\Sigma^{-}\overline{\Lambda}$ pairs gains an additional
component from the additional neutron skin. On the other hand, the
contribution from the inner part of the nucleus (II) is reduced by the loss of
antiprotons in the additional neutron layer:

$p^{II}_{\Sigma^{-}\overline{\Lambda}}=\kappa_{II}\cdot{p_{\text{abs}}}\cdot{\frac{\sigma_{\Sigma^{-}\overline{\Lambda}}}{\sigma_{\text{tot}}}}\\\
+(1-p_{\text{abs}})\cdot\kappa_{II}\frac{\rho(n)}{\rho(p)+\rho(n)}\cdot\frac{\sigma_{\Sigma^{-}\overline{\Lambda}}}{\sigma_{\text{tot}}}.$
(7)

Here, we assume that the loss of outgoing $\Sigma^{-}$ and$/$or
$\overline{\Lambda}$ for pairs produced in the additional neutron skin is the
same as for pairs produced within the core nucleus. We thus obtain for the
double ratio

$DR=\frac{p^{II}_{\Sigma^{-}\overline{\Lambda}}\Big{/}p^{II}_{\Lambda\overline{\Lambda}}}{p^{I}_{\Sigma^{-}\overline{\Lambda}}\Big{/}p^{I}_{\Lambda\overline{\Lambda}}}=\frac{p_{\text{abs}}\cdot\frac{\rho(p)+\rho(n)}{\rho(n)}+\Big{(}1-p_{\text{abs}}\Big{)}}{1-p_{\text{abs}}}$
(8)

With the assumption $\rho(p)$ = $\rho(n)$, we finally find for the double
ratio the simple expression

$DR=\frac{1+p_{\text{abs}}}{1-p_{\text{abs}}}$ (9)

To evaluate Eq. (9), we use neutron skin thicknesses which are predicted
during the initialization stage of the GiBUU simulations described in the next
section (cf. Fig. 3). In the following we will compare the predictions of Eq.
(9) with the result of microscopic GiBUU simulations of
$\overline{\text{p}}+A$ interactions at 2.4$\mathrm{\ GeV/{\mit c}}$.

##  3 Transport calculations for a Xe isotope chain

Figure 2: Proton (red) and neutron (blue) rms radii of the initial
distributions for stable or long-lived xenon isotopes used by the GiBUU
simulations. The red squares show experimental values for the proton radii
ANGELI201369 ; Atanasov_2017 .

Figure 3: Green points show the double ratio predicted by the full GiBUU
simulations. For comparison, the horizontal bars are the results of Eq. (9),
adopting the same proton and neutron radii as shown in Fig. 3.

Of course, this extremely simplified geometrical picture has several
deficiencies. A more realistic description of the hyperon pair production can
be achieved by microscopic transport calculations. As such a model, GiBUU
describes many features of $\mathrm{\overline{p}}$–nucleus interactions in the
FAIR energy range Buss20121  .

In the following we employ the GiBUU model to study
$\overline{\text{p}}+\text{Xe}$ for a series of isotopes in the mass range
from $A=129$ to $A=136$. These isotopes have a sizable abundance and can — in
principle — be used for experimental studies. The simulations were performed
with an incident antiproton momentum of 2.4 GeV/c. For each isotope about 28
million events were generated, requiring a computational time of about one day
on the MOGON2 high performance computing cluster at the University of Mainz.
During the initialization of these simulations, the proton and neutron
distributions of the target nuclei were generated by a self consistent
relativistic mean-field model.

The red and blue points in Fig. 3 show the rms radii of the proton and neutron
distributions as a function of the Xe mass number. In these calculations not
only the neutron radius but — unlike in the simple geometrical picture
presented above — also the proton radius is rising with increasing neutron
number. This increasing charge radius is qualitatively consistent with
experimental data, see open squares in Fig. 3 ANGELI201369 ; Atanasov_2017 ,
though the experimental slope with the mass number is only about half as
large.

The green points in Fig. 3 show the double ratio DR predicted by the GiBUU
simulations. For comparison, the horizontal bars in this figure show the
results of Eq. (9), using the same proton and neutron radii which are used as
starting point for the GiBUU simulations. While the absolute agreement is less
relevant, we find a remarkable correlation between the double ratio evaluated
with our simple analytical expression and the double ratio predicted by the
complete GiBUU simulations. This correlation shows that the double ratio may
be a good measure for the neutron skin difference. In future we will apply
this method to the case of 40Ca and 48Ca.

## Acknowledgements

We thank Horst Lenske and Theo Gaitanos for valuable comments. This project
has received funding from the European Unionâs Horizon 2020 research and
innovation programme under grant agreement No 824093. The presented data were
in part collected within the framework of the bachelor theses of Martin
Christiansen and the PhD thesis of Falk Schupp at the Johannes Gutenberg
University Mainz.

## References

  * (1) J. Piekarewicz, F.J. Fattoyev, Physics Today 72, 30 (2019), https://doi.org/10.1063/PT.3.4247
  * (2) M. Tanabashi, K. Hagiwara, K. Hikasa, K. Nakamura, Y. Sumino, F. Takahashi, J. Tanaka, K. Agashe, G. Aielli, C. Amsler et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018)
  * (3) I. Angeli, K. Marinova, Atomic Data and Nuclear Data Tables 99, 69 (2013)
  * (4) D. Atanasov, D. Beck, K. Blaum, C. Borgmann, R.B. Cakirli, T. Eronen, S. George, F. Herfurth, A. Herlert, M. Kowalska et al., Journal of Physics G: Nuclear and Particle Physics 44, 044004 (2017)
  * (5) O. Buss, T. Gaitanos, K. Gallmeister, H. van Hees, M. Kaskulov, O. Lalakulich, A. Larionov, T. Leitner, J. Weil, U. Mosel, Physics Reports 512, 1 (2012)


# Performance study of a Time of Flight Method used for cosmic ray detection

ALİ YILMAZ
Department of Electrical and Electronics Engineering, Giresun University,
28200, Giresun, Turkey.

###### Abstract

Time of Flight methods have been rapidly developed and used in many
experiments recently for determination of particle direction, identification
of particles and energy resolutions. This paper describes a method of time-
mark determination on the reconstruction algorithm, based on the sampled
signal, used for time-of-flight measurements. This method was developed for
distinguishing the signals by fitting to pulse shape which were received from
scintillator detector with a silicon photomultiplier readout have been
developed for a cosmic ray counter telescope. The method was verified using
experimental data taken in the location $40\degree
54^{\prime}52^{\prime\prime}N$ and $38\degree 19^{\prime}26^{\prime\prime}E$
with the elevation of 30 m above the sea level. The data samples were acquired
by the counters which have a scintillator with dimensions of
20$\times$20$\times 1.4~{}cm^{3}$, optically coupled from one side to silicon
photomultiplier, then the signals readout by fast sampling digitizer board
Domino Ring Sampler Board version 4. The method can reconstruct each pulse
even for multiple events without losing the count within the small time
window. Using this method 4.969 ns time-of-flight value were established and
the rise times for scintillation counters, named Tile 1 and Tile 2, were
measured about $6.27\pm 0.16~{}ns$ and $4.979\pm 0.165~{}ns$, respectively.

##  1 Introduction

Precisely measuring the timing information of a particle in an experiment
would allow one to successfully reconstruct physical events. Time of Flight
(TOF) methods have been developed and used in many fields such as high energy
physics experiments [1, 2], astroparticle detectors [3, 4], TOF cameras [5,
6], TOF-PET detectors [7, 8] for identifying the particle species, determining
the direction of the particles and energy measurements. The counter telescope
is useful in researches which is a relatively a few number of desired
particles has to be counted in a major unwanted background. The analyzing the
coincidence time, passing through the counters telescopes, is practical to
remove the spurious counts registered because of the accidental coincidence of
counts in the each detectors. While comparing the Silicon PhotoMultipliers
(SiPMs) with traditional PhotoMultiplier Tubes (PMTs), which have been widely
used in the past experiments for decades, SiPMs have some advantages such as
advanced photon detection efficiency, compactness, excellent single photon
time resolution, insensitive to magnetic fields and inexpensive price. Due to
these advantages, they have been started to use in a wide range of physics
experiments nowadays. A SiPM signal, is constituted of a few different
avalanche events formed at arbitrary times. The difficulties in achieving high
accuracy on the time-marking of the signal and TOF resolutions of a counter
telescope with existing algorithms consist the subsequent subjects:
Variation of the baseline: If the dark count rate increases considerably then
peaks of the signal are mostly overlying on each other within a time latency.
Therefore, the method requires to detect and evaluate signal starting time for
different heights while working at higher over-voltages.
Electrical noise: The signal to noise ratio needs to be at lowest point 5 in
the worst conditions (minimum over-voltages and higher electronic noise). If
the primary signals will be evaluated mostly which yields an ambiguity of 20%
on the signal size while computing directly thus a fitting algorithm is
required.
Very different signal heights and widths: Due to the proposed method can be
able to adjust itself while measuring various properties without a user
concerning to change and optimize the parameters for any evaluation, very
different signal heights and widths vary strongly with over-voltage and SiPM
type. [9].

The rest of the paper is organized as follows. In the next section is
described used experimental setup for signal reconstruction. In Section 3, a
method developed to extract time information of the signal coming from
scintillation counters without varying difficulties listed above. Finally, an
application of this method to the experimental data is constructed and the
results of its is given.

##  2 Experimental Details

###  2.1 Description of the setup

The hardware of the scintillation counter telescope constructed with two
identical scintillator plates, called $Tile~{}1$ and $Tile~{}2$, separated by
$160~{}cm$ apart as seen in Figure 1. Each counter box consist of a KURARAY
organic scintillator plate ($20\times 20\times 1.4~{}cm^{3}$). The
scintillator has excellent properties, in the view of getting the accurate
time information, like yielding the light in blue region of the spectrum, the
emission peak is around 430 nm [10]. Each scintillator panel wrapped in Tyvek
paper for diffusing the reflection, and one SensL SiPM ($3\times 3$ mm2)
contacted to read the produced signal. This SiPM is able to produce a sharp
output pulse of $<2~{}ns$ at Full Width at Half Maximum (FWHM). The bias
voltage of this device is about $27.5~{}V$ and the dynamic range over the
breakdown voltage is about $2~{}V$. Hence the filling factor of the device is
64%, the gain is $2.3\times 10^{6}$ [11]. The produced signal is digitized by
Domino Ring Sampler Board (DRS4), developed by Stefan Ritt [12]. The Data
Acquisition (DAQ) program is based on the setup schematically shown in Figure
1. It is managed by a shell script that controls the two main C++ programs.
One of the program is controlling the Arduino so that reads temperatures from
the SiPM readout circuit and adjust the operating voltage in order to keep
constant the gain of the SiPM on each tile. The other program manages the
Domino Ring Sampler board $v4$ (DRS4) which digitizes the signal detected by
the SiPM and stores in ROOT binary format for further analysis. The DAQ is
based on waveform sampling at 2 GS/s, covering a $2.5~{}\mu s$ window. This
detector can be used to select horizontal tracks for detecting tau shower
produced by the neutrino interacting in Earth crust [13, 14, 15]. So the TOF
resolution of the selected tracks can be achieved to $0.5~{}ns$.

Figure 1: Schematic view of data acquisition.

###  2.2 Description of the reconstruction method

The reconstruction program, is based on ROOT analysis package [16] using C++
programming language, represents a necessary tool for the analysis of the
data. The method optimized for high sensitivity while determining the peak
level and rise time of the signal are required to effectively use the peak
height readout option. The method explained in this section is developed to
work on a continuous SiPM signal.

####  2.2.1 Least Square Method

A function was written in the analysis program in order to determine the the
starting point of the signal. This point gives the timing information of the
particle so that a mathematical procedure can be applied for estimating the
best-fitting curve to a dedicated set of points by minimizing the sum of the
squares of the offsets (”the residuals”) from the curve. By substituting a set
of $N$ data points ($x_{i},y_{i}$) to a straight-line function

$\displaystyle y(x)=mx+c$ (1)

where $y(x)$ and $x$ are holding the voltage and time informations of the SiPM
signal, respectively. The equation is generally called least square method or
linear regression [17]. Assuming that the uncertainty $\sigma_{i}$ associated
with each measurement $y_{i}$ is known, and that the $x_{i}$’s (values of the
dependent variable) are known exactly. The chi-square, $\chi^{2}$, function is
used to compute how well the model matches up with the data.

$\displaystyle\chi^{2}(m,c)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{N-1}\Big{[}\frac{y_{i}-mx_{i}-c}{\sigma_{i}}\Big{]}^{2}$
(2)

When the errors of a measurement are spread out normally, then this function
will yield maximum likelihood coefficient evaluations of $m$ and $c$; if they
are not spread out normally, then the evaluations are not maximum likelihood
but might be advantageous to use in a practical sense. In order to determine
$m$ and $c$ the Eq.2 is minimized. At its minimum, derivatives of
$\chi^{2}(m,c)$ with respect to $m,c$ become zero:

$\displaystyle\frac{\partial\chi^{2}}{\partial c}$ $\displaystyle=$
$\displaystyle-2\cdot\sum_{i=0}^{N-1}\frac{y_{i}-mx_{i}-c}{\sigma_{i}^{2}}=0$
$\displaystyle\frac{\partial\chi^{2}}{\partial m}$ $\displaystyle=$
$\displaystyle-2\cdot\sum_{i=0}^{N-1}\frac{x_{i}\cdot(y_{i}-mx_{i}-c)}{\sigma_{i}^{2}}=0$
(3)

These conditions can be revised in a suitable form given in the following
sums:

$\displaystyle\xi_{\sigma}\cong\sum_{i=0}^{N-1}\frac{1}{\sigma_{i}^{2}}\qquad\xi_{x}\cong\sum_{i=0}^{N-1}\frac{x_{i}}{\sigma_{i}^{2}}\qquad\xi_{y}\cong\sum_{i=0}^{N-1}\frac{y_{i}}{\sigma_{i}^{2}}$
$\displaystyle\xi_{xx}\cong\sum_{i=0}^{N-1}\frac{x_{i}^{2}}{\sigma_{i}^{2}}\qquad\xi_{xy}\cong\sum_{i=0}^{N-1}\frac{x_{i}y_{i}}{\sigma_{i}^{2}}$
(4)

the Eq.3 can be rewritten by using the definitions in the Eqs. 4;

$\displaystyle m\xi_{x}+c\xi_{\sigma}$ $\displaystyle=$ $\displaystyle\xi_{y}$
$\displaystyle m\xi_{xx}+c\xi_{x}$ $\displaystyle=$ $\displaystyle\xi_{xy}$
(5)

The result of these two equations in two unknowns is computed as

$\displaystyle\triangle$ $\displaystyle\cong$
$\displaystyle\xi_{\sigma}\xi_{xx}-(\xi_{x})^{2}$ $\displaystyle c$
$\displaystyle=$
$\displaystyle\frac{\xi_{xx}\xi_{y}-\xi_{x}\xi_{xy}}{\triangle},\quad
m=\frac{\xi_{\sigma}\xi_{xy}-\xi_{x}\xi_{y}}{\triangle}$ (6)

The solution for the best-fit model parameters $m$ and $c$ is given in the
Eq.6. The expected uncertainties in the calculation of $m$ and $c$ should also
be evaluated for the reason that the measurement errors in the data naturally
bring some of uncertainty in the setting of those parameters. Assuming the
data set are independent, then each value puts its own bit of uncertainty to
the parameters. Taking into account the spread of errors indicates that the
variance $\sigma_{f}^{2}$ in the value of any function is as follow:

$\displaystyle\sigma_{f}^{2}$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{N-1}\sigma_{i}^{2}\cdot\Big{(}\frac{\partial
f}{\partial y_{i}}\Big{)}^{2}$ (7)

The derivatives of $m$ and $c$ in respect of $y_{i}$ might be directly
calculated from the solution for the straight line equation:

$\displaystyle\frac{\partial c}{\partial
y_{i}}=\frac{\xi_{xx}-\xi_{x}\xi_{x_{i}}}{\sigma_{i}^{2}\cdot\triangle},\quad\frac{\partial
b}{\partial
y_{i}}=\frac{\xi_{\sigma}\xi_{x}-\xi_{x_{i}}}{\sigma_{i}^{2}\cdot\triangle}$
(8)

after the summation over the points as in the Eq.7;

$\displaystyle\sigma_{c}^{2}=\frac{\xi_{xx}}{\triangle},\quad\sigma_{m}^{2}=\frac{\xi_{xx}}{\triangle}$
(9)

which are called the variances in the estimates of $m$ and $c$, respectively
[18, 19].

##  3 Results

###  3.1 Application of the Time of Flight Method to the experimental data

The counter telescope is useful in researches where a relatively a few number
of desired particles has to be counted in a major unwanted background. The
analyzing the coincidence time is efficient way of removing the spurious
counts occurred due to the accidental coincidence in the each detectors. TOF
is a method measuring the time that it takes for a particle to travel a
distance of a medium. While a particle hitting one of the scintillator tile as
shown in Figure 2 that starts the time counter and it will be stopped when the
particle hits with the second scintillator tile.

Figure 2: Schematic view of a counter telescope.

TOF method is also useful to discriminate the particle direction (_upward or
downward_). Firstly, if the particle hit the Tile 1, then hits to Tile 2 means
the particle is going to _downward_ direction, the expected result of the time
difference (_$t_{0}ch2-t_{0}ch1$_) should be _positive_ , if not, the particle
is going to _upward_ direction.
Figure 3 shows the algorithm used for calculating the time of flight method
which have mentioned in section section 2.2.1. The signals are registered by
using DRS4 board and and then plotted by _ROOT_ program to observe the time
difference between two tiles. If the signal is greater than the threshold
voltage ($\sim 40~{}mV$), then the program estimates the baseline of each
signal in order to be insensitive to the baseline instability. After that, if
the signal is inside the gate ($180~{}ns-250~{}ns$), the analysis program
firstly finds the peak level (absolute minimum) of the signal is the 100% of
the amplitude, then goes backward until safePoint (which is at least $3\sigma$
below the baseline). Next, the program stores those points for being used in
fit method. The program marks the 10% and 90% percents of the amplitude, which
are not fix number of points, is just below the threshold value and estimated
safePoint. The program fits a line equation using these registered points
using the method explained in section 2.2.1, and crosses its own baseline axis
(time axis). This point registered as _$t_{0}$_. The difference of _$t_{0}$_
gives us the time-of-flight between two counters.

Figure 3: Description of proposed TOF calculation method used in this study.
Tile separation is 160 cm. Magenta and blue lines indicate the fitted lines
according to the method works forward direction from 10% until 90% of the
amplitude.

Figure 3 shows the scheme of this algorithm. The computed TOF is about $5\pm
1.7~{}ns$ for a $downward$ going particle is shown in Figure 4.

Figure 4: Time of Flight distribution of the counter telescope for tile
separation is $160~{}cm$. The peak around $+5~{}ns$ is due to the particle
coming from $Tile~{}1$ to $Tile~{}2$ (downward direction).

###  3.2 Rise Time Test

Time taken by the signal to rise from minimum level to maximum level is named
the rise time, and the time taken by the signal goes from maximum level to
minimum level is named fall time. The nonlinearity of the signal typically
takes place at the bottom and at the top of the signal so that the rise time
is generally determined between the $10\%$ and $90\%$ percent of the amplitude
of the signal. A typical SiPM signal registered by the detector is seen in
Figure 5. It has a final rise and fall times. The results are given in the
Figure 6.

Figure 5: Signal shape of the rising time between the boundaries are $10\%$
percent of the amplitude and $90\%$ percent of the amplitude.

In Figure 6, the mean of the gaussian fit on the rising time distribution are
$6.27\pm 0.16~{}ns$ and $4.979\pm 0.165~{}ns$ for $Tile~{}1$ and $Tile~{}2$,
respectively.

(a) Tile 1

(b) Tile 2

Figure 6: Rising time distribution of the signal is between $10\%$ and $90\%$
percent of the amplitude. (a) is for $Tile~{}1$ and in (b) is for $Tile~{}2$,
respectively.

####  3.2.1 Testing the robustness of the Fit

Goodness of fit, $R^{2}$, is called correlation coefficient, is a quantity
which indicates the quality of a least square fitting to the data. If
$R^{2}=1$ means the fit is perfect but generally expected that this is close
to $1$. Figure 7 depicts that the $R^{2}$ distribution for all registered
events for each tile. The mean value is $0.9752$ and $0.9654$ for $Tile~{}1$
and $Tile~{}2$, respectively.

(a) Tile 1

(b) Tile 2

Figure 7: Goodness of Fit distribution of the signal is between $10\%$ and
$90\%$ percent of the amplitude. (a) is for $Tile~{}1$ and in (b) is for
$Tile~{}2$, respectively.

###  3.3 Multiplicity of the Signal Test

The signal (seen in Figure 3) is a single cosmic ray (generally muon) points
towards the region of the SiPMs. Cosmic rays are large collections of
particles. Because these cosmic rays contain many particles, the particles
will often separate, which could cause multiple hits in the region of the
SiPM. One can see that there are two main hits in the given particular event
(seen in Figure 8).
Multiplicity is defined as a number of hits inside the full sampling window of
the SiPMs had per event which is sketched in Figure 8.

Figure 8: Multiplicity of registered SiPM signal in a full sampling window.

Since the detector is not limited to one characteristic of event, those both
particles would be present, and point towards a nearby region, causing the
total charge to be excessed.

In order to find the multiplicity of the signals, the total charge of the
expected signal in the window gate (180 ns - 250 ns) and the total charge of
the all signal in the full sampling window (512 ns) were integrated and the
baseline of the each signal was removed from the total charge. Figure 9 shows
that the integrated charge for the signals in the full sampling window ((a)
$Tile~{}1$ and (b) $Tile~{}2$) and the signals within the gated window ((c)
$Tile~{}1$ and (d) $Tile~{}2$).

Figure 9: Integrated charge for all signal where the baseline was removed from
each signal. Histograms depicts that the excessive number of particles (a) is
for $Tile~{}1$ and in (b) is for $Tile~{}2$, respectively. The integrated
charge for gated signals (c) is for $Tile~{}1$ and in (d) is for $Tile~{}2$,
respectively. DAQ trigger channel is $Tile~{}2$

.

Figure 10 ((a) Tile 1, (b) Tile 2), depicts that the multiplicity as a
function of the integrated charge which was estimated by subtracting the total
charge of the signal in the defined region from the total charge of the
signals in the full sampling window. The multiplicity number as a function of
the expected signal was also estimated and plotted in the Figure 10 ((c) Tile
1, (d) Tile 2). Those values above 1 are due to more than one particle
interacting with one cell. Since the $Tile~{}2$, was selected as a trigger in
the DAQ, the number of registered signals were statistically more than
$Tile~{}1$ as shown in Figure 9 and Figure 10.

Figure 10: Multiplicity events in terms of integrated charge for the full
sampling window (a) for $Tile~{}1$ and (b) for $Tile~{}2$, and (c) and (d)
give the multiplicity in terms of the number of events. DAQ trigger channel is
$Tile~{}2$

.

##  4 Discussion and Conclusion

In this paper, the performance study of the developed method for time-of-
flight measurements using the scintillation counters are presented. A new DAQ
system, described in section 2.1, working in automatically without needing a
user for interaction. Moreover a new analysis program and a time-of-flight
technique were developed and tested for further analysis. The rise time of the
SiPM signals was accumulated and found to be $6.27\pm 0.16$ and $4.979\pm
0.165~{}ns$ for each tile. The coincidence timing resolution between the tiles
were found about 1.7 ns signal from Figure 4. Multiplicity as a function of
the integrated charge of the signals were estimated around -119.6 pVs (-420
pVs) for $Tile~{}1$ ($Tile~{}2$). The limitation on the results is due to the
present electronics used in the DAQ system. These results proves that the
using this developed method satisfy very well the selection of direction of
cosmic rays while removing the background.

## Acknowledgement

I would like to thank Prof. Dr. Haluk Denizli and Prof. Dr. Maurizio Iori for
their support and encouragement during the course of this work.

## References

  * [1] S. Barwick, J. Beatty, C. Bower, C. Chaput, S. Coutu, G. de Nolfo, D. Ellithorpe, D. Ficenec, J. Knapp, D. Lowder, S. McKee, D. Müller, J. Musser, S. Nutter, E. Schneider, S. Swordy, K. Tang, G. Tarle, A. Tomasch, and E. Torbet, “The high-energy antimatter telescope (heat): An instrument for the study of cosmic-ray positrons,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 400, no. 1, pp. 34 – 52, 1997.
  * [2] A. Alici, “Particle identification with the alice time-of-flight detector at the lhc,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 766, pp. 288 – 291, 2014.  RICH2013 Proceedings of the Eighth International Workshop on Ring Imaging Cherenkov Detectors Shonan, Kanagawa, Japan, December 2-6, 2013.
  * [3] B. N. Ratcliff, “Advantages and limitations of the rich technique for particle identification,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 595, no. 1, pp. 1 – 7, 2008.  RICH 2007.
  * [4] M. Iori, I. Atakisi, G. Chiodi, H. Denizli, F. Ferrarotto, M. Kaya, A. Yilmaz, L. Recchia, and J. Russ, “Sipm application for a detector for uhe neutrinos tested at sphinx station,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 742, pp. 265 – 268, 2014.  4th Roma International Conference on Astroparticle Physics.
  * [5] F. Mufti and R. Mahony, “Statistical analysis of signal measurement in time-of-flight cameras,” ISPRS Journal of Photogrammetry and Remote Sensing, vol. 66, no. 5, pp. 720 – 731, 2011.
  * [6] M. Gilles, H. Fayad, P. Miglierini, N. Boussion, O. Pradier, S. Scheib, L. Cozzi, and D. Visvikis, “Ep-1511: Patient positioning using surface images from time-of-flight (tof) cameras,” Radiotherapy and Oncology, vol. 115, p. S823, 2015.  3rd ESTRO Forum 24-28 April 2015, Barcelona, Spain.
  * [7] A. Ronzhin, M. Albrow, S. Los, M. Martens, P. Murat, E. Ramberg, H. Kim, C.-T. Chen, C.-M. Kao, K. Niessen, A. Zatserklyaniy, M. Mazzillo, B. Carbone, G. Condorelli, G. Fallica, A. Piana, D. Sanfilippo, G. Valvo, and S. Ritt, “A sipm-based tof-pet detector with high speed digital drs4 readout,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 703, pp. 109 – 113, 2013\.
  * [8] C. Joram, “Imaging results and tof studies with axial pet detectors,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 732, pp. 586 – 590, 2013\.  Vienna Conference on Instrumentation 2013.
  * [9] M. Putignano, A. Intermite, and C. Welsch, “A non-linear algorithm for current signal filtering and peak detection in sipm,” Journal of Instrumentation, vol. 7, pp. P08014–P08014, Aug 2012.
  * [10] “http://kuraraypsf.jp/psf/.”
  * [11] “https://www.sensl.com/.”
  * [12] S. Ritt, “The drs chip: cheap waveform digitizing in the ghz range,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 518, no. 1, pp. 470 – 471, 2004.  Frontier Detectors for Frontier Physics: Proceedin.
  * [13] M. Iori, A. Sergi, and D. Fargion, “Test results of a prototype designed to detect horizontal cosmic ray flux,” arXiv preprint astro-ph/0409159, 2004\.
  * [14] M. Iori and A. Sergi, “An orientable time-of-flight detector for cosmic rays,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 588, no. 1, pp. 151 – 154, 2008.  Proceedings of the First International Conference on Astroparticle Physics.
  * [15] A. Yilmaz, M. Iori, H. Denizli, K. Y. Oyulmaz, S. A. Yilmaz, U. Keskin, and J. S. Russ, “Simulation of horizontal tau-neutrino induced shower to optimize the site parameters,” PoS, vol. ICRC2017, p. 471, 2018.
  * [16] R. Brun and F. Rademakers, “Root–an object oriented data analysis framework,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 389, no. 1-2, pp. 81–86, 1997.
  * [17] S. Margulies, “Fitting experimental data using the method of least squares,” Review of Scientific Instruments, vol. 39, pp. 478–480, Apr 1968.
  * [18] W. H. Press, S. A. Teukolsky, B. P. Flannery, and W. T. Vetterling, Numerical Recipices, the Art of Scientific Computing.  Cambridge University Press, 3 ed., 2007.
  * [19] E. W. Weisstein, “Least squares fitting. from mathworld-a wolfram web resource,” 2002.