11institutetext: University of Edinburgh

# Active and sparse methods in smoothed model checking

Paul Piho  Jane Hillston

###### Abstract

Smoothed model checking based on Gaussian process classification provides a
powerful approach for statistical model checking of parametric continuous time
Markov chain models. The method constructs a model for the functional
dependence of satisfaction probability on the Markov chain parameters. This is
done via Gaussian process inference methods from a limited number of
observations for different parameter combinations. In this work we consider
extensions to smoothed model checking based on sparse variational methods and
active learning. Both are used successfully to improve the scalability of
smoothed model checking. In particular, we see that active learning-based
ideas for iteratively querying the simulation model for observations can be
used to steer the model-checking to more informative areas of the parameter
space and thus improve sample efficiency. Online extensions of sparse
variational Gaussian process inference algorithms are demonstrated to provide
a scalable method for implementing active learning approaches for smoothed
model checking.

##  1 Introduction

Stochastic modelling coupled with verification of logical properties via model
checking has provided useful insights into the behaviour of the stochastic
models from epidemiology, systems biology and networked computer systems. A
large number of interesting models in these fields are too complex for the
application of exact model checking methods [11]. To improve scalability of
model checking there has been significant work on statistical model checking
that aims to estimate the satisfaction probability of logical properties based
on independently sampled trajectories of the stochastic model [2].

This paper considers statistical model checking in the context of parametrised
continuous time Markov chain models. Statistical model checking methods have
generally considered single parametrisations of a model. Based on a large
number of independent sample trajectories one can estimate the probability of
the model satisfying a specified logical property defined over individual
sample trajectories. In order to gain insight across the entire parameter
space associated with a model, it can be necessary to repeat the estimation
procedures with different parametrisations to cover the whole space, which
leads to poor scalability.

As an alternative, a model checking approach based on Gaussian process
classification, named smoothed model checking, was proposed in [4]. The main
result of that paper was to show that under mild conditions the function
mapping parameter values to satisfaction probabilities is smooth. Thus the
problem can be solved as a Gaussian process classification problem where the
aim is to estimate the function describing the satisfaction probability over
the parameter space. Model checking results, returning a label true or false,
of individual simulation trajectories are used as the training data to infer
how the satisfaction probability depends on the model parameters. This method
can be used to greatly reduce the number of simulation trajectories needed to
estimate the satisfaction probability in exchange for some accuracy.

There are two aspects that limit the speed of such model checking procedures.
Firstly, the computational cost of gathering the individual trajectories and
secondly, the cost of (approximate) Gaussian process inference itself. Both
benefit from keeping the number of gathered trajectories as low as possible
while minimising the impact a smaller set of training data has on the accuracy
of the methods. In order to keep the gathered sample size small we propose a
method based on active learning. In particular, we make the observation that
the parameter space of models is usually constrained to physically reasonable
ranges. However, even when constrained to such ranges there can be large parts
of the parameter space where the probability of satisfying a formula exhibits
stiff behaviour. Adaptively identifying stiff and non-stiff parts of the
parameter space in order to decide where to concentrate the computational
effort leads to improved algorithms for smoothed model checking.

Our approach is based on moving the smoothed model checking approach in the
context of continuous time Markov chains to an online setting. This allows us
to implement effective active sampling strategies for the parameter space of
the model that take into account the already gathered information. In order to
make this approach scalable we leverage and demonstrate the usage of state of
the art sparse variational Gaussian process inference methods. In particular,
we consider streaming variational inference with inducing points [6].

##  2 Related work

A wealth of literature exists on statistical model checking of stochastic
systems. The use of statistical methods in the domain of formal verification
is motivated by the fact that in order to perform statistical model checking
it is only necessary to be able to simulate the model. Thus these methods can
be used for systems where exact verification methods are infeasible including
black-box systems [12]. In its classical formulation, this involves hypothesis
testing [18] with respect to the desired (or undesired) property based on
independent trials, or in this case, stochastic simulations.

In addition to the frequentist approaches based on hypothesis testing, there
have been Bayesian approaches [10] to estimate the satisfaction probability of
a given logical formula. Our work follows the approach presented in [4] where
the dependence of the satisfaction probability on model parameters is modelled
as a Gaussian process classification problem.

The problem of deciding where to concentrate the model checking efforts is
closely related to optimal experimental design. Experimental design problems
are commonly treated as optimisation problems where the goal is to allocate
resources in a way that allows the experimental goals to be reached more
rapidly and thus with smaller costs [17]. This idea is also known in the
machine learning literature as active learning [20]. The idea is to design
learning algorithms that interactively query an oracle to label new data
points.

In the context of model checking, active learning was used in [5] to solve a
closely related threshold synthesis problem. That approach used a base grid on
the parameter space for initial estimation. The estimates were then refined
around values where the satisfaction probability was close to a defined
threshold. However, the threshold for synthesis has to be defined a priori
making the introduced active step not applicable when we are interested in the
satisfaction probability. We further address the scalability of the ideas
presented by the authors of [5] by considering sparse approximation results
for Gaussian.

##  3 Background

###  3.1 Continuous time Markov chains

Stochastic models are widely used to model a variety of phenomena in natural
and engineered systems. We focus on a type of stochastic model commonly used
in biological modelling, epidemiology and performance evaluation domains.
Specifically, we consider continuous time Markov chain models (CTMCs). To
define a CTMC we start by noting that it is a continuous-time stochastic
process and thus defined as an indexed collection of random variables
$\\{\mathbf{X}\\}_{t\in\mathbb{R}_{\geq 0}}$. We consider CTMCs defined over a
finite state space $S$ with an $|S|\times|S|$ matrix $Q$ whose entries
$q(i,j)$ satisfy

  1. 1.

$0\leq-q(i,i)<\infty$

  2. 2.

$0\leq q(i,j)$ for $i\neq j$

  3. 3.

$\sum_{j}q(i,j)=0$

A CTMC is then defined by the following: for time indices
$t_{1}<t_{2}\ldots<t_{n+1}$ and states $i_{1},i_{2},\ldots,i_{n+1}$ we have

$\mathbb{P}(\mathbf{X}_{t_{n+1}}=i_{n+1}|\mathbf{X}_{t_{n}}=i_{n},\ldots,\mathbf{X}_{t_{1}}=i_{1})=p(i_{n+1};t_{n+1}|i_{n};t_{n})$

where $p(j;t|i;s)$ is the solution to the following Kolmogorov forward
equation

$\frac{\partial}{\partial t}p(j;t|i;s)=\sum_{k}p(k;t|i;s)q(k,j),\qquad\mbox{on
$(s,\infty)$ with $p(j;s|i;s)=\delta_{ij}$}$

with $\delta_{ij}$ being the Kronecker delta taking the value $1$ if $i$ and
$j$ are equal and the value $0$ otherwise. By convention the sample
trajectories of CTMCs are taken to be right-continuous.

In the rest of the paper we consider parametrised models
$\mathcal{M}_{\mathbf{x}}$ and assume that the model $\mathcal{M}$ for a fixed
parametrisation $\mathbf{x}\in\mathbb{R}^{k}$ defines a CTMC. Thus, the model
$\mathcal{M}$ specifies a function mapping parameters $\mathbf{x}$ to
generator matrix $Q$ of the underlying CTMC. A commonly studied special class
of CTMC models are population CTMCs where each state of the CTMC corresponds
to a vector of counts. These counts are used to model the aggregate counts of
groups of indistinguishable agents in a system. In biological modelling and
epidemiology such models are often defined as chemical reaction networks
(CRN).

###### Example 1

Let us consider the following SIR model defined as a CRN

$\displaystyle S+I$ $\displaystyle\xrightarrow{k_{I}}I+I\qquad\qquad
I\xrightarrow{k_{R}}R$

where $S$ gives the number of susceptible, $I$ the infected and $R$ the
recovered individuals in the system. The first type of transition corresponds
to an infected and susceptible individuals interacting, resulting in the
number of infected individuals increasing and the number of susceptible
decreasing. The second type of transition corresponds to recovery of an
infected individual and results in the number of infected decreasing and the
number of recovered increasing. The states of the underlying CTMC keep track
of the counts of different individuals in the system. For the example let us
set the initial conditions to $(95,5,0)$ — at time $0$ there are $95$
susceptible, $5$ infected and $0$ recovered individuals in the system. The
parameters $k_{I}$ and $k_{R}$ give the infection and recovery rates
respectively. We revisit this example throughout the paper to illustrate the
presented concepts.

###  3.2 Smoothed model checking

Smoothed model checking was introduced in [4] as a scalable method for
statistical model checking where Gaussian process classification methods were
used to infer the functional dependence between a parametrisation of a model
and the satisfaction probability given a logical specification.

As described in Section 3.1, suppose we have a model
$\mathcal{M}_{\mathbf{x}}$ parametrised by vector of values
$\mathbf{x}\in\mathbb{R}^{d}$ such that the model $\mathcal{M}$ for a fixed
parametrisation $\mathbf{x}$ defines a CTMC. Additionally assume we have a
logical property $\varphi$ we want to check against. The logical properties we
consider here are defined as a mapping from the time trajectories over the
states of $\mathcal{M}_{\mathbf{x}}$ to $\\{0,1\\}$ corresponding to whether
the property holds for a given sample trajectory of $\mathcal{M}_{\mathbf{x}}$
or not. One way to define such mappings would be, for example, to specify the
properties in metric interval temporal logic (MiTL) [13] or signal temporal
logic (STL) [8] and map the paths satisfying the properties to $1$ and those
not satisfying the properties to $0$. Through sampling multiple trajectories
for the same parametrisation we gain an estimate of the satisfaction
probability corresponding to the parametrisation.

With that in mind, a logical property $\varphi$ with respect to
$\mathcal{M}_{\mathbf{x}}$ can be seen to give rise to a Bernoulli random
variable. The binary outcomes of the random variable correspond to whether or
not a randomly sampled trajectory of $\mathcal{M}_{\mathbf{x}}$ satisfies the
property $\varphi$. We introduce the notation $f_{\varphi}(\mathbf{x})$ for
the parameter of the said Bernoulli random variable given the model parameters
$\mathbf{x}$. In particular, samples from the distribution
$\mathit{Bernoulli}(f_{\varphi}(\mathbf{x}))$ model whether a randomly sampled
trajectory of $\mathcal{M}_{\mathbf{x}}$ satisfies $\varphi$ — for a parameter
value $\mathbf{x}$ the logical property is said to be satisfied with
probability $f_{\varphi}(\mathbf{x})$.

A naive approach for estimating $f_{\varphi}(\mathbf{x})$ at a given
parametrisation $\mathbf{x}$ is to gather a large number $N$ of sample
trajectories and give simple Monte Carlo estimate for the
$f_{\varphi}(\mathbf{x})$ by dividing the number of trajectories where the
property holds by the total number of sampled trajectories $N$. An accurate
estimate requires are large number of samples. However, having such estimate
at a set of given parametrisations does no provide us with a rigorous way to
estimate the satisfaction function at a nearby point.

In [4] the authors considered population CTMCs. It was shown that the
introduced satisfaction probability $f_{\varphi}(\mathbf{x})$ is a smooth
function of $\mathbf{x}$ under the following conditions: the transition rates
of the CTMC $\mathcal{M}_{\mathbf{x}}$ depend smoothly on the parameters
$\mathbf{x}$; and the transition rates depend polynomially on the state vector
$\mathbf{X}$ of the CTMC.

The result was exploited by treating the estimation of the satisfaction
function $f_{\varphi}(\mathbf{x})$ as a Gaussian process classification
problem. The main benefit of this approach is that, based on sampled model
checking results, we can reconstruct an approximation for the functional
dependence between the parameters and satisfaction probability. This makes it
easy to make predictions about the satisfaction probability at previously
unseen parametrisations.

Simulating $\mathcal{M}_{\mathbf{x}}$ we gather a finite set of observations
$\mathcal{D}=\\{(\mathbf{x}_{i},y_{i})|i=1,\cdots,n\\}$ where $\mathbf{x}_{i}$
are the parametrisations of the model and $y_{i}$ correspond to model checking
output over single trajectories. For classification problems, a Gaussian
process prior with mean $m$ and kernel $k$ is placed over a latent function

$g_{\varphi}(\mathbf{x})\sim\mathit{GP}(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime}))$

Here, let us consider the standard squared exponential kernel defined by

$k(\mathbf{x},\mathbf{x}^{\prime})=\exp\left(-\frac{|\mathbf{x}-\mathbf{x}^{\prime}|^{2}}{2\ell}\right)$

where $\ell$ is the length scale parameter governing how far two distinct
points have to be in order to be considered uncorrelated.

The function $g_{\varphi}$ is then squashed through the standard logistic or
probit transformation $\sigma$ so that the composition
$\sigma(g_{\varphi}(\mathbf{x}))$ takes values between 0 and 1. The quantity
$\sigma(g_{\varphi}(\mathbf{x}))$ is interpreted as the probability that
$\varphi$ holds given model parametrisation $\mathbf{x}$ and thus estimates
the probability $f_{\varphi}(\mathbf{x})$ that a simulation trajectory for
parameters $\mathbf{x}$ satisfies the property $\varphi$.

The general aim of Gaussian process inference is to find the distribution
$p(g_{\varphi}(\mathbf{x}^{*})|\mathcal{D})$ over the values $g_{\varphi}$ at
some test point $\mathbf{x}^{*}$ given the set of training observations
$\mathcal{D}$. This distribution is then used to produce a probabilistic
prediction at parameter $\mathbf{x}^{*}$ of
$\sigma(g_{\varphi}(\mathbf{x}^{*}))\approx f_{\varphi}(\mathbf{x}^{*}$). We
present details of inference in the next section. This section is ended by
returning to the running SIR example.

###### Example 2

The property we consider is the following: there always exists an infected
agent in the population in the time interval $(0.0,100.0)$ and in the time
interval $(100.0,120.0)$ the number of infected becomes $0$. Constraining the
parameters to the ranges $k_{I}\in[0.005,0.3]$ and $k_{R}\in[0.005,0.3]$ gives
satisfaction probabilities as depicted in the Figure 1. There each estimate on
the $20$-by-$20$ grid is calculated based on $2000$ stochastic simulation
sample runs of the model. For comparison, Figure 2 gives the results of the
smoothed model checking where $10$ sample trajectories are drawn for each
parameter on the $12$-by-$12$ grid. The smoothed model checking approximation
for the model checking problem shows good agreement with the baseline surface
and is much faster to perform.

Figure 1: Baseline satisfaction probability surface.

###  3.3 Variational inference with inducing points

In order to infer the latent Gaussian process $g_{\varphi}$ based on training
data $\mathcal{D}$ we have to deal with two problems. Firstly the inference is
analytically intractable due to the non-Gaussian likelihood model provided by
Bernoulli observations. To counter this there exist a wealth of approximate
inference schemes like Laplace approximation, expectation propagation [14, 9],
and variational inference methods [21]. Here we consider variational
inference. The second problem is that the methods for inference in Gaussian
process models have cubic complexity in the number of training cases. To
address that there exist sparse approximations based on inducing variables.
Sparse variational methods [7, 21] are popular methods for reducing the
complexity of Gaussian process inference by constructing an approximation
based on a small set of inducing points that are typically selected from
training data. In this section we detail the inference procedure.

Figure 2: Smoothed model checking satisfaction probability surface. $144$ test
points.

Variational inference methods choose a parametric class of variational
distribution for the posterior and minimising the KL-divergence between the
real posterior and the approximate posterior. To accommodate large training
data sets we work with sparse variational methods. We start by defining the
prior distribution

$p(\mathbf{g}_{\varphi},\mathbf{g}_{\mathbf{u}})=\mathcal{N}\left(\begin{bmatrix}\mathbf{g}_{\varphi}\\\
\mathbf{g}_{\mathbf{u}}\end{bmatrix}\Bigg{|}\mathbf{0},\begin{bmatrix}\mathbf{K}_{nn}&\mathbf{K}_{nm}\\\
\mathbf{K}_{nm}^{T}&\mathbf{K}_{mm}\\\ \end{bmatrix}\right)$

where $\mathbf{g}_{\varphi}$ is a vector of $n$ latent function values at
$g_{\varphi}^{i}=g_{\varphi}(x_{i})$ evaluated at
$\mathbf{x}=[\mathbf{x}_{1},\cdots,\mathbf{x}_{n}]$. Similarly,
$g_{\mathbf{u}}^{i}=g_{\varphi}(\mathbf{u}_{i})$ is a vector of $m$ latent
function values evaluated at chosen _inducing points_
$\mathbf{u}=[\mathbf{u}_{1},\cdots,\mathbf{u}_{m}]$. The matrices
$\mathbf{K}_{nn}$, $\mathbf{K}_{nm}$ and $\mathbf{K}_{mm}$ are defined by the
kernel function. In particular, the $(i,j)$-th element of the matrix
$\mathbf{K}_{nn}$ is given by $k(\mathbf{x}_{i},\mathbf{x}_{j})$. Similarly,
$\mathbf{K}_{nm}$ gives the kernel matrix between the training points
$\mathbf{x}$ and the inducing points $\mathbf{u}$ and $\mathbf{K}_{mm}$ gives
the kernel matrix between the locations of inducing points. We then fit the
variational posterior at those points rather than the whole set of training
data points. The assumption we are making is that
$p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{g}_{\varphi},\mathbf{u})=p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{u})$.
That is, the inducing points $\mathbf{u}$ are a sufficient statistic for a
function value at a test point $\mathbf{x}_{*}$.

Under this assumption we make predictions at a test point $\mathbf{x}^{*}$ as
follows:

$p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{y})=\int
p(g_{\varphi}(\mathbf{x}_{*}),\mathbf{u}|\mathbf{y})d\mathbf{u}=\int
p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{u})p(\mathbf{u}|\mathbf{y})d\mathbf{u}$

Thus, we need posterior distribution $p(\mathbf{u}|\mathbf{y})$ at the
inducing points. From Bayes rule we know that
$p(\mathbf{u}|\mathbf{y})=\frac{p(\mathbf{y}|\mathbf{u})p(\mathbf{u})}{p(\mathbf{y})}$.
The terms $p(\mathbf{y})$ and $p(\mathbf{y}|\mathbf{u})$ are called the
likelihood and the marginal likelihood respectively. The above integral can be
seen as the weighted average of the prior model
$p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{u})$ with weight being determined by
the posterior distribution. Here, as mentioned, we consider variational
approximations where $p(\mathbf{u}|\mathbf{y})$ is approximated by a
multivariate Gaussian $q(\mathbf{u})$ making the expression for
$p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{y})$ tractable.

Finding the parameters of $q(\mathbf{u})$ is done by minimising the KL
divergence between $q(\mathbf{u})$ and $p(\mathbf{u}|\mathbf{y})$.

$\displaystyle\mathcal{D}_{\mathit{KL}}(q(\mathbf{u}),p(\mathbf{u}|\mathbf{y}))$
$\displaystyle=\int
q(\mathbf{u})\log\frac{q(\mathbf{u})}{p(\mathbf{u}|\mathbf{y})}d\mathbf{u}=-\big{\langle}\log\frac{p(\mathbf{y},\mathbf{u})}{q(\mathbf{u})}\big{\rangle}_{q(\mathbf{u})}+\log
p(\mathbf{y})$ (1)

The term $\log p(\mathbf{y})$ is known as the log marginal likelihood. In the
following we use the well-known Jensen’s inequality111For a concave function
$f$ and a random variable $X$ we have the following well-known inequality:
$f\langle X\rangle\geq\langle f(X)\rangle$. to derive a lower lower bound for
the log marginal likelihood. As log function in concave we get the following.

$\displaystyle\log p(\mathbf{y})$ $\displaystyle=\log\int
p(\mathbf{y},\mathbf{u})d\mathbf{u}=\log\Big{\langle}\frac{p(\mathbf{y},\mathbf{u})}{q(\mathbf{u})}\Big{\rangle}_{q(\mathbf{u})}\geq\Big{\langle}\log\frac{p(\mathbf{y},\mathbf{u})}{q(\mathbf{u})}\Big{\rangle}_{q(\mathbf{u})}$
$\displaystyle=\langle\log
p(\mathbf{y}|\mathbf{u})\rangle_{q(\mathbf{u})}-\mathcal{D}_{\mathit{KL}}(q(\mathbf{u}),p(\mathbf{u}))$

The last line of the above equation is known as the evidence-based lower bound
or ELBO. Now note that the first term in the expression for KL-divergence in
Equation 1 is exactly the derived ELBO. As the KL-divergence is always non-
negative then maximising the ELBO minimises the KL-divergence between the
approximate and true posteriors $q(\mathbf{u})$ and
$p(\mathbf{u}|\mathbf{y})$.

Choosing our approximating family of variational distribution to be
multivariate Gaussian makes the KL term in ELBO easy to evaluate. The integral
in the expectation term can be computed via numerical approximation schemes
making in possible to use ELBO as a utility function for optimising the
parameters of the approximate posterior $q(\mathbf{u})$. When $q(\mathbf{u})$
is chosen to be a multivariate Gaussian these parameters are the mean
$\bm{\mu}$ and covariance matrix $\bm{\Sigma}$.

The final part of this section briefly discusses making predictions based on
the approximate posterior $q(\mathbf{u})$. With the approximation
$p(\mathbf{u}|\mathbf{y})\approx q(\mathbf{u})$ the predictions are given by
the following integral

$p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{y})\approx\int
p(g_{\varphi}(\mathbf{x}_{*})|\mathbf{u})q(\mathbf{u})d\mathbf{u}$

This can be shown [16] to be a probability density function of a Normal
distribution with the following mean and covariance matrix

$\displaystyle k(\mathbf{x}_{*},\mathbf{u})K_{mm}^{-1}\bm{\mu}$ $\displaystyle
k(\mathbf{x}_{*},\mathbf{x}_{*})-k(\mathbf{x}_{*},\mathbf{u})K_{mm}^{-1}\left[K_{mm}-\bm{\Sigma}\right]^{-1}\left[k(\mathbf{x}_{*},\mathbf{u})K_{mm}^{-1}\right]^{T}$

###  3.4 Inducing points

Previously we introduced the idea of sparse variational inference where the
posterior distribution is fitted to a selection of inducing points
$\mathbf{u}$ such that $|\mathbf{u}|$ is much smaller than the size of the
whole training data set. In this section we discuss how to choose the set of
inducing points $\mathbf{u}$. A common approach is to use cluster centers from
k-means++ clustering as the inducing points [3]. In this papers we simply take
a regular grid of inducing points.

###  3.5 Active learning

Active learning methods in machine learning are a family of methods which may
query data instances to be labelled for training by an _oracle_ [19]. The
fundamental question asked by active learning research is whether or not these
methods can achieve higher accuracy than passive methods with fewer labelled
examples. This is closely related to the established area of optimal
experimental design, where the goal is to allocate experimental resources in a
way that reduces uncertainty about a quantity or function of interest [17,
20].

In the case of Gaussian process classification problems like smoothed model
checking an active learning procedure can be set up as follows. An _active
learner_ consists of a classifier learning algorithm $\mathcal{A}$ and a
_query function_ $\mathcal{Q}$. The query function is used to select an
unlabelled sample $u$ from the pool of unlabelled samples $\mathcal{U}$. This
sample is then labelled by an oracle. In the case of stochastic model checking
the pool of unlabelled samples $\mathcal{U}$ corresponds to a subset of the
possible parametrisations for the model. An oracle is implemented by running
the stochastic simulation for the selected parametrisation and model checking
the resulting trajectory.

The above describes a pool-based active learner. Common formulations of such
pool-based learners select a single unlabelled sample at each iteration to be
sampled. However, in many applications it is more natural to acquire labels
for multiple training instances at once. In particular, the query function
$\mathcal{Q}$ selects a subset $U\subset\mathcal{U}$. We see in the next
section that the sparse inference methods can be extended to a setting where
batches of training data become available over time making it natural to
decide on a query function that selects batches of queries. The main
difficulty of selecting a batch of queries instead of a single query is that
the instances in the subset $U$ need to be both informative and diverse in
order to make the best use of the available labelling resources.

##  4 Active model checking

The shape and properties of the functional dependence of satisfaction for a
logical specification with respect to parameters are generally not known a
priori and can exhibit a variety of properties. For example, in the running
example the satisfaction probability is non-stiff with respect to parameter
changes in one direction while being stiff with respect changes in another. In
addition, much of the sampling was performed in completely flat regions of the
parameter space. Thus the key challenge addressed in this section is where to
sample to make the posterior estimates as informative as possible about the
underlying mechanics. We aim to decide on the regions where the satisfaction
probability surface is not flat and concentrate most of our model checking
effort there.

In this section we introduce the main contribution of this paper — active
model checking. The general outline of the procedure is given by Algorithm 1.

Algorithm 1 Active smoothed model checking

1:procedure Model checking(model $\mathcal{M}$, property $\varphi$, parameter
space $\mathcal{X}$)

2:
$\mathcal{D}_{\mathit{old}}\leftarrow\mathit{generate\\_initial\\_data}(\mathcal{M},\varphi,\mathcal{X})$

3: $\mathcal{D}_{\mathit{new}}\leftarrow\mathcal{D}_{\mathit{old}}$

4:
$\mathbf{u}\leftarrow\mathit{inducing\\_points}(\mathcal{D}_{\mathit{old}})$

5: $\mathbf{v}\leftarrow\mathbf{u}$

6: $q(\mathbf{v})\leftarrow\mathit{initialise\\_posterior}(\mathbf{v})$

7: while $\mathit{true}$ do

8:
$q(\mathbf{v})\leftarrow\mathit{update\\_variational}(q(\mathbf{u}),\mathbf{v},\mathcal{D}_{\mathit{new}},\mathcal{D}_{\mathit{old}})$

9:
$\mathcal{D}_{\mathit{old}}\leftarrow\mathcal{D}_{\mathit{old}}\cup\mathcal{D}_{\mathit{new}}$

10:
$\mathcal{D}_{\mathit{new}}\leftarrow\mathit{query\\_new}(q(\mathbf{v}),\mathcal{M},\varphi,\mathcal{X})$

11: $\mathbf{u}\leftarrow\mathbf{v}$

12:
$\mathbf{v}\leftarrow\mathit{update\\_inducing}(\mathcal{D}_{\mathit{new}})$

13: end while

14: return $q(\mathbf{v})$

15:end procedure

The first step, given by the procedure $\mathit{generate\\_initial\\_data}$,
is to simulate the initial data set $\mathcal{D}_{\mathit{old}}$ via
stochastic simulation of the CTMC model $\mathcal{M}$ for a sample of the
parameter space $\mathcal{X}$ and checking whether or not the individual
trajectories satisfy the property $\varphi$ or not. The initial set of
parameter samples can for example be a regular grid or sampled uniformly from
the parameter space. We are going to experiment with both initialisations.

Based on the results $\mathcal{D}_{\mathit{old}}$ we choose the inducing
points $\mathbf{v}$, captured by the procedure $\mathit{inducing\\_points}$,
$\mathbf{v}$ based on the approach in Section 3.4 and initialise the
variational posterior. Here we are going to initialise the posterior as a
multivariate Gaussian $q(\mathbf{v})=\mathcal{N}(\mathbf{0},\mathbf{I}))$ with
$0$ mean and identity covariance matrix.

Each iteration of the model checking loop will update the variational
posterior $q(\mathbf{v})$ and use the fitted approximate posterior to query
new points in the parameter space to perform model checking. At the end of
each iteration, based on the updated data set we update the set of inducing
points.

There are two issues to be resolved before the procedure can be implemented.
First is that the direct use of ELBO as introduced in Section 3.3 does not
suffice in the online setting where new data becomes available in batches.
Second is the challenge of choosing an appropriate query function that is
going to suggest more points in the parameter space at which to gather more
model checking data. These will be addressed in the following sections.

###  4.1 Streaming setting

In order to incorporate active learning ideas into the Gaussian process based
model checking approach we need to address the problem that not all of the
training data is available a priori. For our purposes it is important to be
able to conduct inference in a streaming setting where data is gradually added
to the model. A naive approach would refit a Gaussian process from scratch
every time a new batch of data arrives. However, with potentially large data
sets this becomes infeasible. To perform sparse variational inference in a
scalable way the method needs to avoid revisiting previously considered data
points. In particular, we consider the method proposed in [6] that derives a
correction to ELBO that allows us to incorporate streaming data incrementally
into the posterior estimate.

The main question is how to update the variational approximation to the
posterior at time step $n$, denoted $q_{\mathit{old}}(\mathbf{u})$, to form an
approximation at the time step $n+1$, denoted $q_{\mathit{new}}(\mathbf{v})$.
In the following we note the variational posteriors $q_{\mathit{old}}$ and
$q_{\mathit{new}}$ at $\mathbf{g}_{\varphi}$ and inducing points $\mathbf{u}$
and $\mathbf{v}$, respectively, are approximations to the true posteriors
given observations $\mathbf{y}_{\mathit{old}}$ and
$\mathbf{y}_{\mathit{new}}$. It was shown in [6] that the lower bound $\log
p(\mathbf{y}_{\mathit{new}}|\mathbf{y}_{\mathit{old}})$ becomes

$\displaystyle\mathcal{F}(q_{\mathit{new}}(\mathbf{g}_{\varphi},\mathbf{v}))=\int
q_{\mathit{new}}(\mathbf{g}_{\varphi},\mathbf{v})\log
p(\mathbf{y}_{\mathit{new}}|\mathbf{g}_{\varphi},\mathbf{v})d(\mathbf{g}_{\varphi},\mathbf{v})-\mathcal{D}_{KL}(q_{\mathit{new}}(\mathbf{v}),p(\mathbf{v}))$
$\displaystyle-\mathcal{D}_{KL}(q_{\mathit{new}}(\mathbf{u}),q_{\mathit{old}}(\mathbf{u}))+\mathcal{D}_{KL}(q_{\mathit{new}}(\mathbf{u}),p(\mathbf{u}))$

The above can be interpreted as follows: the first two terms give the ELBO
under the assumption that the new data seen at iteration $n+1$ is the whole
data set; the final two terms take into account the old likelihood through the
approximate posteriors at old inducing points and the prior $p(\mathbf{u})$.
This allows us to implement an online version of the smoothed model checking
where observation data arrives in batches.

###  4.2 Query strategies

As discussed in Section 3.5, in order to implement an active learning method
for model checking we need to decide which new parameters are tested based on
the existing information. In the following we consider two query strategies
for active model checking.

####  4.2.1 Predictive variance

The first approach is a commonly used experimental design strategy which aims
to minimise the predictive variance. Recall that in smoothed model checking
for a property $\varphi$ we fit a latent Gaussian processes $g_{\varphi}$. The
posterior satisfaction probability for parameter $\mathbf{x}_{*}$ given the GP
$g_{\varphi}$ is then calculated via

$p(y_{*}=1|\mathcal{D},\mathbf{x}_{*})=\int\sigma(g_{\varphi}(\mathbf{x}_{*}))p(g_{\varphi}(\mathbf{x}_{*})|\mathcal{D})dg_{\varphi}(\mathbf{x}_{*})$

The above can also be seen as the expectation of
$\sigma(g_{\varphi}(\mathbf{x}_{*})$ with respect to the distribution
$g_{\varphi}(\mathbf{x}_{*})$, denoted
$\mathbb{E}\left[\sigma(g_{\varphi}(\mathbf{x}_{*})\right]$. Similarly, we can
consider the variance of this estimate

$\mathbb{E}\left[\sigma(g_{\varphi}(\mathbf{x}_{*})^{2}\right]-\mathbb{E}\left[\sigma(g_{\varphi}(\mathbf{x}_{*})\right]^{2}$

Our aim is then to iteratively train the Gaussian process model so that
predictive variance over the parameter space is minimised.

Before giving the outline of the proposed procedure we address the issue of
redundancy in the query points. As pointed out in Section 3.5 simply taking a
set of top points with respect to the utility function, in this case
predictive variance, leads to querying parameters that are clustered together.
We can overcome this problem by clustering the pool of unlabelled samples
$\mathcal{U}$ from which the query choice is made. In particular, the top
points with respect to predictive variance are chosen from a pool of samples
where the redundancy is already reduced. Informally, this leads to the
following basic outline of the procedure:

  1. 1.

Sample an initial set of training points or parametrisations $\mathbf{x}$ of
the model (via uniform, Latin hypercube sampling or taking points on a regular
grid) and conduct model checking based on sampled trajectories. These points
are used to choose the set of inducing points via the k-means++ clustering
algorithm and fit the first iteration of the Gaussian process model.

  2. 2.

For the next iteration we randomly sample another set of points $U$ and
cluster them via regular kmeans. From the set of cluster centres $U_{k}$ the
query function $\mathcal{Q}$ selects a set of points for model checking. The
query function is simply defined by taking the subset $U^{*}$ of cluster
centres $U_{k}$ where the predictive variance, as defined above, is the
highest. This concentrates the sampling to points where the model is most
uncertain about its prediction.

  3. 3.

The points in $U^{*}$ are labelled by simulating the model for the
parametrisations in $U^{*}$ and checking the resulting trajectories against
the logic specification $\varphi$. The results are incorporated into the
Gaussian process model via the streaming method discussed in Section 4.1.

  4. 4.

Repeat points 2 and 3 until a set computational budget is exhausted.

####  4.2.2 Predictive gradient

The second strategy we consider is based on the predictive mean
$\bar{g}_{\varphi}(\mathbf{x})$ of the Gaussian process. Our aim is to
concentrate the sampling at the locations where the predictive mean undergoes
the most rapid change. This requires gradients of the predictive mean.

We recall from Section 3.3 that for a variational posterior $q(\mathbf{u})$
with mean $\bm{\mu}$ and covariance $\bm{\Sigma}$, the posterior mean at a
point $\mathbf{x}$ is given by

$\bar{g}_{\varphi}(\mathbf{x})=k(\mathbf{x},\mathbf{u})K_{mm}^{-1}\bm{\mu}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}k(\mathbf{x},\mathbf{u})\bm{\alpha}$ (2)

Only the first part, the kernel function, depends on $\mathbf{x}$. Thus, in
order to get the derivative of the predictive mean we need to differentiate
$k(\mathbf{x},\mathbf{u})$. Recall, that in this paper we chose to work with
the squared exponential kernel given by

$k(\mathbf{x},\mathbf{u}_{i})=\exp\left(-\frac{|\mathbf{x}-\mathbf{u}_{i}|^{2}}{2\ell}\right)$
(3)

We have used $\mathbf{u}_{i}$ to denote a single inducing point in the set of
inducing points $\mathbf{u}$. Thus, the derivative of
$k(\mathbf{x},\mathbf{u}_{i})$ with respect to $\mathbf{x}$ is given by

$\frac{dk(\mathbf{x},\mathbf{u}_{i})}{d\mathbf{x}}=-\frac{\mathbf{x}-\mathbf{u}_{i}}{\ell}\exp\left(-\frac{|\mathbf{x}-\mathbf{u}_{i}|^{2}}{2\ell}\right)=-\frac{\mathbf{x}-\mathbf{u}_{i}}{\ell}k(\mathbf{x},\mathbf{u}_{i})$
(4)

As the above is for a single inducing point $\mathbf{u}_{i}$, to compute the
derivative of the posterior mean we need to concatenate this derivative for
all $m$ inducing points. Thus, we get

$\displaystyle\frac{d\bar{g}_{\varphi}(\mathbf{x})}{d\mathbf{x}}=-\ell^{-1}\begin{bmatrix}\mathbf{x}-\mathbf{u}_{1}\\\
\cdots\\\
\mathbf{x}-\mathbf{u}_{m}\end{bmatrix}(k(\mathbf{x},\mathbf{u})\odot\bm{\alpha})$

where $\odot$ denotes element-wise multiplication. Given this we can proceed
as in the case of the predictive variance. The only change is that instead of
considering the predictive variance for each sampled set of parameters we
calculate the norm of $\frac{d\bar{g}_{\varphi}(\mathbf{x})}{d\mathbf{x}}$ and
define the query function to choose a subset $U^{*}$ of cluster centres
$U_{k}$ with the highest norms.

Figure 3: Mean satisfaction probability surface. Sparse smoothed MC.

###  4.3 Results

In this section we evaluate the proposed active learning method for model
checking on the running SIR example. The methods are compared to the baseline
naive stochastic simulation-based model checking and smoothed model checking
without the sparse approximation and the active step. We present several
metrics for comparing the smoothed model checking results with the empirical
mean based on stochastic simulation. The first is the mean and standard
deviation of the difference between the mean probability predicted by the
fitted Gaussian processes and the empirical mean from the stochastic
simulation results at each of the points on the $20\times 20$ grid. Secondly,
we consider the maximum difference between the predicted mean probability and
the naive empirical mean. Finally, we give the root-mean-square error (RMSE)
$\sqrt{\frac{1}{N}\sum{\left(\bar{g}_{\varphi}(\mathbf{x}_{i})-\bar{f}_{\varphi}(\mathbf{x}_{i})\right)^{2}}}$
where $\bar{g}_{\varphi}(\mathbf{x}_{i})$ is the predicted mean satisfaction
probability for parametrisation $\mathbf{x}_{i}$. We denote by
$\bar{f}_{\varphi}(\mathbf{x}_{i})$ the empirical estimate of the satisfaction
probability at $\mathbf{x}_{i}$ given 2000 sample trajectories.

In the active learning experiments we start with a coarser $10\times 10$ grid
followed by an active iteration where an additional 300 points are chosen to
refine the approximation for a total of 400 training points. The inducing
points are initialised by choosing the initial grid of $100$ points as
inducing points and kept constant for the remainder of the fitting procedure.
Similarly we present the results for sparse smoothed model checking for a
$20\times 20$ grid with the $10\times 10$ grid of inducing points as well as
smoothed model checking where inducing points are not chosen. Parameters of
the kernel function are kept constant in each of the smoothed model checking
experiments.

####  4.3.1 SIR

Table 1: Comparison of accuracy for smoothed model checking and the sparse and active learning extensions. SIR model. Naive threshold greater than 0.02. Method | error mean/var | maximum | RMSE
---|---|---|---
Smoothed MC | $(0.042,0.032)$ | $0.13$ | $0.578$
Sparse smoothed MC | $(0.049,0.036)$ | $0.151$ | $0.664$
Active sparse smoothed MC |  |  |
Predictive variance | $(0.044,0.038)$ | $0.16$ | 0.629
Predictive gradient | $\mathbf{(0.044,0.033)}$ | $\mathbf{0.136}$ | $\mathbf{0.6}$
Random sampling | $(0.049,0.04)$ | 0.21 | 0.684
Table 2: Comparison of computation times for smoothed model checking and the sparse and active learning extensions. SIR model. Method | SSA | Inference | Query | Total
---|---|---|---|---
Naive statistical MC | 59.5 | N/A | N/A | 59.5
Smoothed MC | 1.6 | 24.6 | N/A | 26.2
Sparse smoothed MC | 1.6 | 4.6 | N/A | $\mathbf{6.2}$
Active sparse smoothed MC |  |  |  |
Predictive variance | 1.7 | 6.2 | 1.9 | 8.6
Predictive gradient | 1.8 | 5.0 | 0.4 | 8.5
Random sampling | 1.6 | 5.7 | 0.00 | 7.3

We again consider the SIR model. The results for accuracy are summarised in
Table 1. Table 1 gives the comparisons for each point on the $20\times 20$
grid where the naive model checking was conducted. We restrict the view to
those points where the naive model checking-based satisfaction probability
estimates exceed $0.02$ in order to consider locations where the surface is
not completely flat. Note that we would expect the full smoothed MC without
inducing points to offer the best accuracy due to the fact that fewer
approximations are made in the inference algorithm. However, the active
methods with predictive variance and gradient-based query functions provide a
better approximation than the sparse model checking without an active step.
The benefit of the sparse methods comes from significant reductions in
computation costs. From Table 2 we see that there is some overhead associated
with implementing the active learning procedures compared to simply exhausting
the chosen computational budget of 400 points. However, as expected the sparse
methods offer a significant speed-up compared to the smoothed model checking
with no inducing points.

Figure 3 gives the mean satisfaction probability surface based on the fitted
sparse Gaussian process. Figures 4 and 5 present the evolution of the
predictive mean surface through two active learning iterations. All figures
are accompanied by the scatter plots showing where the samples were drawn.

Figure 4: Mean satisfaction probability surface. Active and sparse smoothed
model checking with 2 iterations. Predictive variance-based query function.

Figure 5: Mean satisfaction probability surface. Active and sparse smoothed
model checking with 2 iterations. Predictive mean gradient-based query
function.

###  4.4 Implementation

The prototype implementation in written in the Julia programming language and
makes use of the tools provided as part of Julia Gaussian Processes
repositories [1] to set up the Gaussian process models. The CTMC models are
defined as chemical reaction networks with tools provided as part of the SciML
ecosystem for scientific simulations [15]. The simulations were carried out on
a laptop with Intel i7-10750H CPU.

##  5 Conclusions

In this paper we applied sparse approximation and active learning to smoothed
model checking. By leveraging existing sparse approximations, we improved
scalability of the inference algorithms for Gaussian process classification
corresponding to the smoothed model checking problem. Additionally, we showed
that by concentrating the sampling to high variance or high predictive
gradient areas of the parameter space, we improved the resulting approximation
compared to sparse models with uniform or grid-based sampling of model
parameters. When compared to the standard smoothed model checking approach
with no inducing point approximation and no active step, our method
significantly speeds up the inference procedure while attempting to reduce
errors inherent in sparse approximations. This aligns with the pre-existing
results from active learning literature which aim to construct learning
algorithms that actively query for observations in order to improve accuracy
while keeping the number of observations needed to a minimum.

As further work, we aim to refine our query methods and make a comparison with
other existing methods in the active learning literature. Secondly, we plan to
link the choice of inducing points to the active query methods more directly.
In particular, we will test if the inducing points, and perhaps the underlying
kernel parameters, can be effectively reconfigured through active iterations.
This would further improve the approximation to the satisfaction probability
surface at the non-stiff parts of the parameter space. Finally we will
consider alternative kernel functions. The kernel function chosen in this
paper is a standard first approach in many settings but is best suited for
modelling very smooth functions — not necessarily the case with satisfaction
probability surfaces for parametric CTMCs.

## References

  * [1] Gaussian processes for machine learning in julia. https://github.com/JuliaGaussianProcesses, accessed: 2021-01-14
  * [2] Agha, G., Palmskog, K.: A survey of statistical model checking. ACM Trans. Model. Comput. Simul. 28(1), 6:1–6:39 (2018)
  * [3] Arthur, D., Vassilvitskii, S.: K-means++: The advantages of careful seeding. In: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms. p. 1027–1035. SODA, Society for Industrial and Applied Mathematics, USA (2007)
  * [4] Bortolussi, L., Milios, D., Sanguinetti, G.: Smoothed model checking for uncertain continuous-time markov chains. Inf. Comput. 247, 235–253 (2016)
  * [5] Bortolussi, L., Silvetti, S.: Bayesian statistical parameter synthesis for linear temporal properties of stochastic models. In: TACAS 2018,. Lecture Notes in Computer Science, vol. 10806, pp. 396–413. Springer (2018)
  * [6] Bui, T.D., Nguyen, C.V., Turner, R.E.: Streaming sparse gaussian process approximations. In: Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017. pp. 3299–3307 (2017)
  * [7] Csato, L., Opper, M.: Sparse online gaussian processes. Neural Computation 14, 641–668 (2002)
  * [8] Donzé, A., Maler, O.: Robust satisfaction of temporal logic over real-valued signals. In: FORMATS 2010. Proceedings. Lecture Notes in Computer Science, vol. 6246, pp. 92–106. Springer (2010)
  * [9] Hernandez-Lobato, D., Hernandez-Lobato, J.M.: Scalable gaussian process classification via expectation propagation. In: Gretton, A., Robert, C.C. (eds.) Proceedings of the 19th International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research, vol. 51, pp. 168–176. PMLR (2016)
  * [10] Jha, S.K., Clarke, E.M., Langmead, C.J., Legay, A., Platzer, A., Zuliani, P.: A bayesian approach to model checking biological systems. In: Degano, P., Gorrieri, R. (eds.) Computational Methods in Systems Biology, 7th International Conference, CMSB. Lecture Notes in Computer Science, vol. 5688, pp. 218–234. Springer (2009). https://doi.org/10.1007/978-3-642-03845-7_15
  * [11] Kwiatkowska, M.Z., Norman, G., Parker, D.: Stochastic model checking. In: SFM. Lecture Notes in Computer Science, vol. 4486, pp. 220–270. Springer (2007)
  * [12] Legay, A., Lukina, A., Traonouez, L., Yang, J., Smolka, S.A., Grosu, R.: Statistical model checking. In: Steffen, B., Woeginger, G.J. (eds.) Computing and Software Science - State of the Art and Perspectives, Lecture Notes in Computer Science, vol. 10000, pp. 478–504. Springer (2019). https://doi.org/10.1007/978-3-319-91908-9_23
  * [13] Maler, O., Nickovic, D.: Monitoring temporal properties of continuous signals. In: FORMATS/FTRTFT 2004, Grenoble, France, September 22-24, Proceedings. Lecture Notes in Computer Science, vol. 3253, pp. 152–166. Springer (2004)
  * [14] Minka, T.P.: Expectation propagation for approximate bayesian inference. In: Breese, J.S., Koller, D. (eds.) UAI: Proceedings of the 17th Conference in Uncertainty in Artificial Intelligence. pp. 362–369. Morgan Kaufmann (2001)
  * [15] Rackauckas, C., Nie, Q.: Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia. Journal of Open Research Software 5(1) (2017)
  * [16] Rasmussen, C.E., Williams, C.K.I.: Gaussian processes for machine learning. Adaptive computation and machine learning, MIT Press (2006)
  * [17] Santner, T.J., Williams, B.J., Notz, W.I.: The Design and Analysis of Computer Experiments. Springer series in statistics, Springer (2003)
  * [18] Sen, K., Viswanathan, M., Agha, G.: Statistical model checking of black-box probabilistic systems. In: Alur, R., Peled, D.A. (eds.) Computer Aided Verification, 16th International Conference, CAV. Lecture Notes in Computer Science, vol. 3114, pp. 202–215. Springer (2004). https://doi.org/10.1007/978-3-540-27813-9_16
  * [19] Settles, B.: From theories to queries. In: Active Learning and Experimental Design workshop, In conjunction with AISTATS. JMLR Proceedings, vol. 16, pp. 1–18. JMLR.org (2011), http://proceedings.mlr.press/v16/settles11a/settles11a.pdf
  * [20] Settles, B.: Active Learning. Synthesis Lectures on Artificial Intelligence and Machine Learning, Morgan & Claypool Publishers (2012). https://doi.org/10.2200/S00429ED1V01Y201207AIM018
  * [21] Titsias, M.K.: Variational learning of inducing variables in sparse gaussian processes. In: Dyk, D.A.V., Welling, M. (eds.) Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics, AISTATS. JMLR Proceedings, vol. 5, pp. 567–574 (2009)


# An Experiment of Randomized Hints on an Axiom of Infinite-Valued Łukasiewicz
Logic

Ruo Ando  National Institute of Informatics
2-1-2 Hitotsubashi, Chiyoda-ku
Tokyo, Japan
<EMAIL_ADDRESS>Yoshiyasu Takefuji  Musashino University, , Tokyo135-8181,
Japan
3-3-3 Ariake, Koto-Ku
Tokyo, Japan

###### Abstract

In this paper, we present an experiment of our randomized hints strategy of
automated reasoning for yielding Axiom(5) from Axiom(1)(2)(3)(4) of Infinite-
Valued Łukasiewicz Logic. In the experiment, we randomly generated a set of
hints with size ranging from 30 to 60 for guiding hyper-resolution based
search by the theorem prover OTTER. We have successfully found the most useful
hints list (with 30 clauses) among 150 * 6 hints lists. Also, we discuss a
curious non-linear increase of generated clauses in deducing Axiom(5) by
applying our randomized hints strategy.

###### Index Terms:

Randomized hints, Łukasiewicz Logic, hyper-resolution, condensed detachment,
OTTER

##  I Inifite-valued Łukasiewicz logic

Infinite-valued logic is the classic and promising field in automated
reasoning. Recently, Pykacz discussed whether quantum mechanics should be
described by many-valued logic [1]. Infinite-valued logic can cope with the
world without the excluded middle of which principle Kolmogorov discussed in
detail [2]. Pablo et al. [3] point out that Kleene’s many-valued logic is used
to cope with sentences such as Liar or the Truth-Teller. In the experiment, we
use the automated reasoning program OTTER [4]. We apply a randomized hints
strategy [10] for directing a hyper-resolution based search.

###  I-A The semantics of infinite value logic

Łukasiewicz firstly presented the semantics of the infinite-valued logic
(sentential calculus). His semantics aims to interpret the sentential
variables that appeared in Principia Mathematics by Whitehead and Russell [5].
We call the Łukasiewicz axioms (1)-(5) as $\L_{{\aleph}_{0}}$ corresponding to
Harris et al. [6].

We introduce a predicate $v$ for assigning rational numbers with the closed
unit interval [0,1]. The interval [0,1] is applied to each atomic sentence in
$\L_{{\aleph}_{0}}$. Using $v$, we can cope with the sentences determined
recursively, as follows.

$\displaystyle\rightarrow:v(C_{pq})=p\rightarrow q=\begin{cases}{0\ \text{if
$p=>q$}}\\\ {q-p\ \text{if $p<q$}}\end{cases}$ (1)
$\displaystyle\&:v(K_{pq})=max(p,q)$ (2)
$\displaystyle\vee:v(A_{pq})=min(p,q)$ (3) $\displaystyle\neg:1-p$ (4)

Here in (1), the predicate $v$ represents a conditional of the $C_{pq}$, which
means p implies q. In (2), v assigns a conjunction $K_{pq}$ which yields the
maximum of the values of p and q. In (3), $A_{pq}$ which is $\vee$ means a
disjunction by the minimum of the values of p and q.

Material implecation of Łukasiewicz logic truth table is as follows:

TABLE I: Łukasiewicz logic truth table. A $\rightarrow$ B | B
---|---
| F | U | T
| F | T | T | T
A | U | U | T | T
| T | F | U | T

The Łukasiewicz logic differs from Kleene logic in its definition of the
implication that “unknown implies unknown” is true.

###  I-B The Axiomatics of Hilbert-style system $\L_{{\aleph}_{0}}$

Łukasiewicz has conjectured that the set of five axiom schemes (A1-A5) and the
single rule of condensed detachment (CD) [11]. With A1-A5 and CD, the
semantics of $\L_{{\aleph}_{0}}$ are completed.

$\displaystyle(A1):A\rightarrow(B\rightarrow A)$ (5)
$\displaystyle(A2):((A\rightarrow B)->(B\rightarrow
C))\rightarrow(A\rightarrow C)$ (6) $\displaystyle(A3):((A\rightarrow
B)\rightarrow B)\rightarrow((B\rightarrow A)\rightarrow A)$ (7)
$\displaystyle(A4):(\neg A\rightarrow\neg B)\rightarrow(B\rightarrow A)$ (8)
$\displaystyle(A5):((A\rightarrow B)\rightarrow(B\rightarrow
A))\rightarrow(B\rightarrow A)$ (9)

    (CD) Cpq : {p->q, p} -> q
    (A1) CpCqp
    (A2) CCpqCCqrCpr
    (A3) CCCpqqCCqpp
    (A4) CCNpNqCqp
    (A5) CCCpqCqpCqp

It is said that Wajsberg gave proof of this conjecture in the 1930s. But it
has been lost (not published). Afterward, in the 1950s, Rose and Rosser [7]
provided the proofs. Also, the completeness of the proofs concerning more
general classes of the algebraic structure was presented in [8].

##  II Methodology

###  II-A Hyper-Resolution

Hyper-Resolution takes a non-positive clause called the nucleus and
simultaneously infers each of its negative literals. Those negative literals
are called satellites.

The general scheme is:

$\displaystyle{K_{1,l},.,K{1,n}}$ $\displaystyle...$
$\displaystyle{K{m,l},.,K{m,n}}$ $\displaystyle\underline{\\{\lnot
L_{1},.,\lnot L_{m+1},L_{l}\\}\exists\sigma.\sigma}$
$\displaystyle\underline{=mgu([|K_{1}|,.,|K_{m,1}|],[|L_{1}|,.,|L_{m}|])}$
$\displaystyle\\{K_{1,2},.,K_{1,n},K_{m,2},.K_{m,n},L_{m+1},.,L_{l}\\}\sigma$

$K_{i}$ denotes the clauses in the list above, and $L_{i}$ denotes the
literal. Hyper-Resolution is applied to a set of m unit clauses K1 … Km and a
single nucleus L1, …, Lm+1 consisting of m + 1 literals.

###  II-B Condensed detachment

Condensed detachment [9] enables us to combine substitutions and modus ponens
into a single rule.

    i(s,t) : major premise
    r : minor premise

yields $t_{\sigma}$, where $\sigma$ is a most general unifier for terms r and
s.

For the expression of condensed detachment in the language of OTTER, we apply
hyper-resolution. The single rule of condensed detatchment is a particluar
instance of a nucleus of the form:

    list(usable).
    -P(i(x,y)) | -P(x) | P(y).
    end_of_list.

Here the predicate P (provable) is a unary predicate symbol. If we apply
condensed detachment to equation (1), we obtain clause 61 as follows.

    1 [] -P(i(x,y))| -P(x)|P(y).
    2 [] P(i(x,i(y,x))).
    61 [hyper,1,2,2] P(i(x,i(y,i(z,y)))).

###  II-C Hints: a proof sketch

Hins in OTTER can be described as a proof sketch. The proof sketch for
yielding A5 from A1-A4 in $\L_{{\aleph}_{0}}$ can be thought of as a set of
conditions which is lemmas to prove sufficient for the proof in
$\L_{{\aleph}_{0}}$. Ideally, a proof sketch composed of a hints list is
already proof. Then the conditions are also met.

The hints strategy is closely related to the weighting strategy. Both
strategies provide the desired control over the reasoning .process.

Weighting focus on each clause for assigning a weight by upper-supplied weight
templates. Templates provide a recursive mapping of weight from term and atoms
corresponding the values of weight as follows:

    weight_list(pick_given).
    weight(P(i(x,i(y,x))), 2).
    weight(P(i(i(x,y),i(i(y,z),i(x,z)))), 2).
    end_of_list.

In contrast, the hints strategy copes with the identification of key clauses
in the reasoning process instead of the general calculation of weights.

    list(hints).
    P(i(i(x,y),i(i(z,x),i(z,y))))
    end_of_list.

The hint strategy enhances the weighting strategy in two points.

  1. 1.

Hints strategy focuses on entire facts (clauses) rather than terms and
subterms.

  2. 2.

Hints strategy uses subsumption for determining the value of a generated
clause.

The hints strategy provides both semantic and logical components for the
evaluation of a generated clause by leveraging subsumption.

##  III Experiment

In experiment, we use workstation with Intel(R) Xeon(R) CPU E5-2620 v4
(2.10GHz) and 251G RAM. We set the list of set of support corresponding to the
equations (1)-(4) as follows:

    list(sos).
    P(i(x,i(y,x))).
    P(i(i(x,y),i(i(y,z),i(x,z)))).
    P(i(i(i(x,y),y),i(i(y,x),x))).
    P(i(i(n(x),n(y)),i(y,x))).
    end_of_list.

We set the passive list to search unit conflict to yield the proof from
(1)-(4) to (5).

    list(passive).
    -P(i(i(i(a,b),i(b,a)),i(b,a))).
    end_of_list.

Figure 1 shows the average size of the set of support. In the experiment, we
repeated the random extraction of hints 150 times. Average is calculated for
each size of the hints list. Also, Figure 2 shows the average number of
generated clauses.

$\displaystyle average\,size=\frac{\sum_{n=1}^{150}Xn}{150}$ (10)

The X-axis of Figure 1 and Figure 2 is the size of the hints list, ranging
from 30 to 60. The size of set support shown in Figure 1 is relatively stable
compared with Figure 2. The number of generated clauses is drastically
increased between hints*40 to hints*45.

Figure 1: average size of set of support.

Figure 2: average number of generated clauses.

Figure 3 depicts the size of the set of support for every 150 times of
iterations. We sort the X-axis by the size of the set of support. Figure 4
depicts the number of generated clauses for every 150 times of iterations. We
sort the X-axis by the number of generated clauses. In Figure 3, some change
points are observed in areas A (x=95-100) and B (x=126-129) Curiously, there
is no difference between hints*30 and hints*45 from x=1 to 95. In Figure 4, we
observe two change points, A (x=84-87) and B (x=119-120). During these
periods, the number of generated clauses with hints*30 and hints*45 are
drastically increased.

Figure 3: size of set support.

Figure 4: the number of generated clauses.

## References

We use the automated reasoning program OTTER [4]. Veroff present the
effectiveness of hints strategy in automated reasoning [10]. The axiomatics of
L0 is presented in [6]. The use of condensed detachment is discussed in [11].

##  IV Conclusion

In this paper, we present an experiment of our randomized hints strategy of
automated reasoning for yielding Axiom 5 on Łukasiewicz axiom
$\L_{{\aleph}_{0}}$. We have successfully found the most useful 30 hints among
900 hints lists. Also, we have observed some curious non-linear increases in
the size of set of support and the number of generated clauses in deducing
Axiom (5) by comparing the number of hints from 30 to 60.

## References

  * [1] Jaroslaw Pykacz, “Can many-valued logic help to comprehend quantum phenomena?”, International Journal of Theoretical Physics 54(12), August 2014.
  * [2] A. N. Kolmogorov, On the principle ”tertium non datur”, Mathematicheskil Sbornik, vol. 32 (1924/25), pp. 646-667; English translation, On the principle of excluded middle, in [Hei 67], pp. 416-437
  * [3] Pablo Cobreros, Paul Egre, David Ripley and Robert van Rooij: Foreword: Three-valued logics and their applications. J. Appl. Non Class. Logics 24(1-2): 1-11 (2014)
  * [4] McCune, W.: 1994b, Otter 3.0 Reference Manual and Guide, Technical Report ANL-94/6 Technical report, Argonne National Laboratory, Argonne, Illinois.
  * [5] Bertrand Russell, Alfred North Whitehead, Principia Mathematica, 1910
  * [6] Kenneth Harris, Branden Fitelson: Distributivity in [Lstrok][aleph]0 and Other Sentential Logics. J. Autom. Reason. 27(2): 141-156 (2001)
  * [7] Alan Rose and J. Barkley Rosser. Fragments of many-valued statement calculi. Trans. Am. Math. Soc., 87:1?53, 1958. 9
  * [8] Chang, C. C., A new proof of the completeness of the Łukasiewicz axioms, Transactions of the American Mathematical Society, vol. 93 (1959), pp. 74?80.
  * [9] William McCune, Larry Wos: Experiments in Automated Deduction with Condensed Detachment. CADE 1992: 209-223
  * [10] Robert Veroff: Using Hints to Increase the Effectiveness of an Automated Reasoning Program: Case Studies. J. Autom. Reason. 16(3): 223-239 (1996)
  * [11] J. A. Kalman. Condensed detachment as a rule of inference. Studia Logica, 42(4):443-451, 1983.

TABLE II: The most effective 30 hints in proof in Łukasiewicz axioms $\L_{{\aleph}_{0}}$ . clauses given | 4403
---|---
clauses generated | 256223
hyper res generated | 256223
demod & eval rewrites | 771375
clauses forward subsumed | 1053187
(subsumed by sos) | 94343
clauses kept | 30781
usable size | 4404
sos size | 26382
kbytes malloced | 38085

    list(hints).
    P(i(i(i(x,y),z),i(i(n(y),n(x)),z))).
    P(i(i(x,y),i(n(y),i(x,z)))).
    P(i(i(x,y),i(i(i(y,x),x),y))).
    P(i(i(i(x,n(i(x,y))),n(i(y,x))),y)).
    P(i(i(i(i(x,y),y),n(i(x,y))),n(y))).
    P(i(i(x,i(i(i(y,z),z),n(i(y,z)))),i(x,n(z)))).
    P(i(i(n(x),y),i(n(y),x))).
    P(i(i(x,i(n(y),n(z))),i(x,i(z,y)))).
    P(i(i(x,i(i(i(y,z),z),n(i(z,y)))),i(x,n(y)))).
    P(i(i(x,i(n(n(y)),z)),i(x,i(y,z)))).
    P(i(i(n(x),n(i(i(y,z),i(i(z,w),i(y,w))))),x)).
    P(i(x,i(i(i(y,x),i(x,y)),y))).
    P(i(i(i(n(x),i(y,n(i(y,x)))),z),z)).
    P(i(i(x,i(i(y,z),w)),i(x,i(z,w)))).
    P(i(x,i(i(y,z),i(i(x,y),z)))).
    P(i(i(x,i(i(y,n(i(y,z))),n(i(i(y,n(i(y,z))),n(z))))),i(x,z))).
    P(i(i(x,y),i(n(y),n(x)))).
    P(i(i(n(x),n(i(y,i(z,y)))),x)).
    P(i(i(i(i(x,y),i(i(z,x),y)),w),i(z,w))).
    P(i(i(i(x,n(i(x,y))),n(i(i(x,n(i(x,y))),n(y)))),y)).
    P(i(i(i(x,n(i(x,y))),n(i(i(x,n(i(x,y))),n(y)))),n(n(y)))).
    P(i(i(i(n(x),n(y)),n(y)),i(i(x,y),n(x)))).
    P(i(i(x,i(y,z)),i(x,i(i(z,w),i(y,w))))).
    P(i(i(i(x,y),z),i(n(x),z))).
    P(i(i(i(n(n(x)),y),z),i(i(x,y),z))).
    P(i(i(i(x,y),i(y,x)),i(y,x))).
    P(i(i(x,i(y,z)),i(y,i(x,z)))).
    P(i(i(x,n(i(y,z))),i(x,n(i(i(z,w),n(y)))))).
    P(i(i(i(i(x,y),y),z),i(i(i(y,x),x),z))).
    P(i(n(n(x)),x)).
    end_of_list.


[a,b]Alessio Golzio

# A study on UV emission from clouds with Mini-EUSO

Matteo Battisti  Mario Bertaina  Karl Bolmgren  Giorgio Cambiè  Marco Casolino
Claudio Cassardo  Roberto Cremonini  Silvia Ferrarese  Christer Fuglesang
Massimiliano Manfrin  Laura Marcelli  Lech Piotrowski  Kenji Shinozaki

###### Abstract

Mini-EUSO is the first mission of the JEM-EUSO program located on the
International Space Station. One of the main goals of the mission is to
provide valuable scientific data in view of future large missions devoted to
study Ultra-High Energy Cosmic Rays (UHECRs) from space by exploiting the
fluorescence emission generated by Extensive Air Showers (EAS) developing in
the atmosphere. A space mission like Mini-EUSO experiences continuous changes
in atmospheric conditions, including the cloud presence. The influence of
clouds on space-based observation is, therefore, an important topic to
investigate as it might alter the instantaneous exposure for EAS detection or
deteriorate the quality of the EAS images with consequences on the
reconstructed EAS parameters. For this purpose, JEM-EUSO is planning to have
an IR camera and a lidar as part of its Atmospheric Monitoring System. At the
same time, it would be extremely beneficial if the UV camera itself would be
able to detect the presence of clouds, at least in some specific conditions.
For this reason, we analyze a few case studies by comparing the pixel count
rates from Mini-EUSO during orbits with the cloud cover (as cloud fraction).
This quantity is retrieved from the Global Forecast System (GFS) model at
different height levels over the Mini-EUSO trajectory. The results of this
analysis are reported.

##  1 Introduction

The cloud cover determination is an open issue, and affects several
experimental fields, from the physics of the atmosphere (meteorology and
climatology), to astronomy, and to cosmic-ray detection system. The Mini-EUSO
telescope observes the Earth in the ultraviolet band (UV) [1], and is devoted
to the study of ultra high energy cosmic rays coming from the space through
their extensive air showers in the Earth’s atmosphere. From a first inspection
of Mini-EUSO UV maps it was possible to state the capability of this telescope
to detect tropospheric clouds (see Figure 1). In this particular event a
tropical cyclone was detected in the Indian Ocean north-east of Madagascar
Island (Figure 1a, and a zoom in Figure 1c). The spiral-form of the cloud is
also visible in the RGB Dust composite Image of METEOSAT8, the meteorological
geostationary satellite, reported in Figure 1b.

The Mini-EUSO telescope is on board of International Space Station (ISS) that
completes an orbit around the Earth in $\approx 90$ minutes. The covered area
is a 350-km wide line oscillating between $None$ North-South. The size of the
considered area makes the cloud data retrieval from satellites very complex
due to the fragmentation of the archives and the availability of images. So we
decided, as a first step, to compare the UV maps with the cloud cover
predicted by a Numeric Weather Prediction model (NWP) of the United States’
National Weather Service, the Global Forecast System (GFS), that is available
24h-7days at a resolution of $None$ and with a 1-hour forecast step.

The actual NWPs compute the cloud fraction at each model level from pressure,
cloud condensation mixing ratio, relative humidity and temperature data [2].
NWPs predict the place (grid point), the amount and the thickness of the cloud
on a discrete vertical and horizontal grid. The accuracy of NWPs is in general
high (more than 70%), with a forecast error less than 30% in the first hours
[3]. The cloud mask, resulting in a 2D, latitude-longitude, true/false fields,
may be computed at high horizontal resolution using the NWPs cloud fraction
data, following the same methodology applied to evaluate the cloud top height
in [4]. In the present work, in order to obtain more accurate results,
different thresholds are applied for high, medium and low clouds, and the
cloud mask is compared with the UV maps from Mini-EUSO, to detect the cloud
cover from UV data, a new application of this dataset. The results are
presented and discussed.

Figure 1: a) composite image of the Mini-EUSO pixel counts over the Indian
Ocean at approximately 18.38 UTC of December 5, 2019; b) a RGB Dust Composite
image from the METEOSAT 8 meteorological satellite, showing the tropical
cyclone seen by Mini-EUSO; c) a single image of the Mini-EUSO pixel matrix
where the banded clouds of the cyclone are clearly visible. In the first part
of the image, panel a), the dotted points are middle-low cumulus clouds.

##  2 Data and analysis

Even though 41 Mini-EUSO sessions have been registered up to now, in this
paper only 14 sessions are used, considering the availability of GFS high-
resolution data (available only when the session was taken) and the UV maps.
The available data for each of these sessions correspond to $\approx 44$
minutes, coming from the first and/or the last part of the acquisition.

The analysis is conducted both on all Mini-EUSO available trajectories (over
land and sea regions) and on the sea regions only, to avoid the contamination
coming from antropogenic light. The main light-source in space is due to the
reflectance of the Moon (its phase is reported in Figure 3 as percentage of
full-Moon) and the nightglow/airglow signals, while the sources from the Earth
are the cloud reflection (of space light), lightnings (not visible in the data
because they were triggered by the sensor), boat lights or other atmospheric
events. The session considered span from 15 September 2020 to 6 May 2021. The
area compared in each session over sea span from
$2.61\text{\times}{10}^{6}\text{\,}{\mathrm{km}}^{2}$ to
$6.60\text{\times}{10}^{6}\text{\,}{\mathrm{km}}^{2}$, while over sea and land
span from $5.73\text{\times}{10}^{6}\text{\,}{\mathrm{km}}^{2}$ to
$7.30\text{\times}{10}^{6}\text{\,}{\mathrm{km}}^{2}$.

##  3 Methods

###  3.1 UV maps and regridding

Before the UV observations can be used, the possibility of inter-pixel
differences (e.g. differences in pixel efficiency) must be accounted for. The
different pixels of the Mini-EUSO detector are therefore subjected to a
relative calibration using the lowest count of each pixel during the entire
session as a baseline. This is possible because the physical environment (e.g.
a cloudless ocean) producing this measurement should be the same for each
pixel.

The original pixel resolution of the Mini-EUSO telescope at ground is about
$6.3\text{\,}\mathrm{km}$ on the ground, and the matrix is composed by 48 by
48 pixels with a combined field of view of about $350\times 350$
$\mathrm{km}$. The photo-sensor registers the light in counts per pixel per
time step (the acquisition is every $41\text{\,}\mathrm{ms}$), and to make
possible the comparison between the coarser GFS resolution a regridding
operation is necessary.

This was done following three steps:

  1. 1.

associate the center of each matrix-pixel-time observation with the
corresponding latitude and longitude coordinates. This is done by projecting
the shape of the detector pixel array onto the surface of the Earth by knowing
the position of the ISS and the detector geometry. These calculations are made
in polar coordinates, with the radial distance $\rho_{obs}$ between an
observation on the ground and the ISS subpoint as

$\rho_{obs}=h_{ISS}\tan(FoV_{p})$ (1)

where $rho_{p}$ is the radial distance to the pixel from the center of the
detector, $h_{ISS}$ is the altitude of the ISS and $FoV_{p}$ is the part of
the field of view taken up by one pixel. The polar coordinates are
subsequently transformed geographic coordinates.

  2. 2.

collect overlapping pixels by making spatial averages (in
$41\text{\,}\mathrm{ms}$ the ISS moves of $\approx 300\text{\,}\mathrm{m}$)

  3. 3.

regrid on a $None$ longitude-latitude grid at the same center-points of GFS,
considering the average of the pixel-counts within each bin formed by the
grid.

###  3.2 Cloud mask computation

####  3.2.1 GFS cloud masks

The GFS data111The GFS operative data are real-time retrived from
https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25_1hr.pl used in the
present work consists of the cloud fraction fields with a grid spacing of
$None$ and 39 vertical levels (from
$1\text{\times}{10}^{5}\text{\,}\mathrm{Pa}$ to
$5\text{\times}{10}^{3}\text{\,}\mathrm{Pa}$), resulting in 40491360 points
for each session.

Cloud fraction values span from 0 to 1, where 0 is clear sky and 1 overcast
sky. In order to calculate the cloud mask for low, middle and high layers, the
cloud fraction maximum is identified in the three layers and its value is
compared with the following thresholds: Low clouds ($z<2$ km a.s.l.) are
present if maximum low-cloud fraction overpass $0.2$; Medium clouds ($2<z<6$
km a.s.l.) are present if maximum medium-cloud fraction overpass $0.4$; High
clouds ($6<z<18$ km a.s.l.) are present if maximum high-cloud fraction
overpass $0.8$. This set of thresholds was choosen after studies [e.g. 4] that
shows the different ability of the global model to detected the three levels
of clouds. Three cloud cover masks are computed respectively for low, medium
and high layers, the total cloud mask is the union of these three masks.

####  3.2.2 UV cloud mask

The UV cloud mask is computed on the UV mean counts re-sampled at the same
grid spacing of GFS, using three different thresholds: 1, 5 and 10 counts per
grid point.

###  3.3 Contingency tables

To assess the differences in the cloud masks a contingency table is built
between each cloud level and UV threshold and for every session. The cloud
mask from the global model is considered the truth, while the cloud mask from
Mini-EUSO UV band is considered the test, so the contingency table is built
following Figure 2.

Figure 2: The contingency table used in this study. Letters refer to the
calculated indexes of 3.3

The two considered indexes are the Accuracy (Proportion Correct, PC) and the
Heidke Skill Score [HSS, 5], the first indicates the hits (true positive and
true negative) with respect the overall cases, while the HSS indicates the
quality of the forecast, where the PC measure is scaled with the reference
value from correct forecasts due to chance.

PC is defined as (refering to Figure 2):

$PC=\frac{a+d}{a+b+c+d}=\frac{a+d}{n}$ (2)

The perfect value for PC is the unity (or here a percentage if PC is
multiplied by 100). Its maximum-likelihood estimate is

$E=\bigg{(}\frac{a+c}{n}\bigg{)}\left(\frac{a+b}{n}\right)+\left(\frac{b+d}{n}\right)\left(\frac{c+d}{n}\right)$
(3)

HSS is defined as:

$HSS=\frac{PC-E}{1-E}=\frac{2(ad-bc)}{(a+b)(b+d)+(a+c)(c+d)}$ (4)

The HSS ranges between $-\infty$ and 1, where 1 is the perfect forecast, 0
indicates no-skill and negative values indicates that the chance forecast is
better.

##  4 Results and discussion

The general statistics (Table 1) show the higher performances with the UV
threshold set at one count per pixel, with slightly better results in the
accuracy of low clouds over the sea (61.0%), while over the entire dataset,
middle (total) cloud cover obtain the best value (60.8%). The HSS shows
similar results regarding the UV thresholds but with a more significant gap.
Again the low-clouds prediction obtained the best results: $0.12$ on the sea
and $0.17$ over land and sea together.

Table 1: Averaged results over the 14 sessions considered. Medium-cloud values represent also the Total cloud cover statistics. Points location | UV threshold | Accuracy ($100PC$) | HSS
---|---|---|---
|  | High | Medium | Low | High | Medium | Low
Sea ($6.48\text{\times}{10}^{7}\text{\,}{\mathrm{km}}^{2}$) | 1 | 48.3 | 59.6 | 61.1 | 0.019 | 0.072 | 0.119
5 | 49.0 | 40.1 | 57.0 | 0.009 | 0.047 | 0.088
10 | 50.5 | 36.3 | 56.3 | 0.015 | 0.047 | 0.076
All ($9.48\text{\times}{10}^{7}\text{\,}{\mathrm{km}}^{2}$) | 1 | 47.4 | 60.8 | 56.7 | 0.050 | 0.108 | 0.167
5 | 48.8 | 42.2 | 58.3 | 0.009 | 0.052 | 0.091
10 | 50.5 | 38.1 | 58.0 | 0.015 | 0.048 | 0.080

###  4.1 Sea-data results

We considered only sea points to improve the analysis and remove the pixel
counts related to land lights. The first question was to find a possible
correlation with the Moon phase, as it is the most important light source
during the night. Linear regression Pearson correlation coefficient is
calculated between the accuracy and the HSS indexes and the Moon phase. The
considered cases are the three cloud levels with the three UV thresholds on
the two indexes (Accuracy and HSS). The correlation shows low values mainly
due to the limited statistics, never higher than $0.17$, and the most
correlated series with the moon phase is the low cloud HSS with the higher UV
threshold ($R^{2}=0.166$). The poor correlation values and the flatness of
linear regression lines indicate that the moon does not affect the detection
of clouds. Moreover, considering some single situation, the light curve from
Mini-EUSO suggests that the most appropriate threshold is 1 (as also deducted
from beginning of Section 4 and Table 1).

The accuracy values of UV threshold 1 are between 30% and 82%, but sometimes
very high accuracy values present very low (near zero) HSS, indicating a poor
agreement between observation and prediction. This is the case of session 37
(Figures 3 and 4).

Analyzing the two indexes together, two sessions appear especially
interesting: 32 and 36. In these two cases, the Moon phase was respectively 6%
and 41%, confirming the slight influence of the quantity of moon light, and
accuracy is 70.14% and 64.06%, while HSS is 0.40 and 0.28. In particular HSS
have very high and promising values.

Figure 3: The histogram of HSS with the UV threshold 1 and the sea data only.
Figure 4: The histogram of Accuracy with the UV threshold 1 and the sea data
only.

In Figure 5 the cloud masks calculated for session 36 is shown. It is possible
to see that when the threshold of 1 is overpassed on pixel counts clouds are
present, while in clear sky zones the UV map shows no pixel counts (white is
below the threshold).

Figure 5: Session 36, 6 February 2021. a) The cloud mask computed on the UV
counts from Mini-EUSO, the visible points are the pixels that pass the
threshold of 1 count per pixel; b) the cloud mask from the GFS cloud fraction,
as the union of low, medium and high clouds. In both cloud masks 1 represents
overcast sky and 0 clear sky.

##  5 Conclusions

This study, considering the complexity of cloud detection and forecast, shows
the ability of Mini-EUSO UV sensor to detect cloud cover during its operative
time. The best UV threshold to compute the cloud mask appears to be one count
per pixel, and the overall cloud detection was good (mainly in the two
reported examples of session 32 and 36), also considering the degradation of
the UV maps to reach the same grid spacing of GFS. From the first inspection
showed in Figure 1, and from the complete statistics over the 14 session
considered, it is possible to state that the Mini-EUSO sensor is able to
capture the middle and low cloud better than high clouds, that usually are
thinner and have a smaller optical depth. However they are able able to modify
the detection of the EAS.

Future improvements and investigations will increase the reference NWP
resolution and reach a greater detail in cloud cover to be compared.

## Acknowledgments

This work was supported by State Space Corporation ROSCOSMOS, by the Italian
Space Agency through the ASI INFN agreement n. 2020-26-HH.0 and contract n.
2016-1-U.0, by the French space agency CNES, National Science Centre in Poland
grant 2017/27/B/ST9/02162.

This research has been supported by the Interdisciplinary Scientific and
Educational School of Moscow University "Fundamental and Applied Space
Research". The article has been prepared based on research materials carried
out in the space experiment "UV atmosphere".

## References

  * [1] S. Bacholle, P. Barrillon, M. Battisti, A. Belov, M. Bertaina, F. Bisconti et al., _Mini-EUSO mission to study earth UV emissions on board the ISS_ , _The Astrophysical Journal Supplement Series_ 253 (2021) .
  * [2] K.M. Xu and D.A. Randall, _A semiempirical cloudiness parameterization for use in climate models_ , _Journal of Atmospheric Sciences_ 53 (1996) 3084.
  * [3] Q.-Z. Ye and S.-S. Chen, _The ultimate meteorological question from observational astronomers:how good is the cloud cover forecast?_ , _Monthly Notices of the Royal Astronomical Society_ 428 (2013) 3288.
  * [4] A. Anzalone, M.E. Bertaina, S. Briz, C. Cassardo, R. Cremonini, A.J. de Castro et al., _Methods to retrieve the cloud-top height in the frame of the JEM-EUSO mission_ , _IEE Transctions on geoscience and remote sensing_ 57 (2019) 304.
  * [5] P. Heidke, _Berechnung des erfolges und der güte der windstärkevorhersagen im sturmwarnungsdienst_ , _Geografiska Annaler_ 8 (1926) 301.

Full Authors List: The JEM-EUSO Collaboration

G. Abdellaouiah, S. Abefq, J.H. Adams Jr.pd, D. Allardcb, G. Alonsomd, L.
Anchordoquipe, A. Anzaloneeh,ed, E. Arnoneek,el, K. Asanofe, R. Attallahac, H.
Attouiaa, M. Ave Pernasmc, M. Bagheriph, J. Balázla, M. Bakiriaa, D.
Barghiniel,ek, S. Bartocciei,ej, M. Battistiek,el, J. Bayerdd, B.
Beldjilaliah, T. Belenguermb, N. Belkhalfaaa, R. Bellottiea,eb, A.A. Belovkb,
K. Benmessaiaa, M. Bertainaek,el, P.F. Bertonepf, P.L. Biermanndb, F.
Biscontiel,ek, C. Blaksleyft, N. Blancoa, S. Blin-Bondilca,cb, P. Bobikla, M.
Bogomilovba, K. Bolmgrenna, E. Bozzoob, S. Brizpb, A. Brunoeh,ed, K.S.
Caballerohd, F. Cafagnaea, G. Cambiéei,ej, D. Campanaef, J-N. Capdeviellecb,
F. Capelde, A. Carameteja, L. Carameteja, P. Carlsonna, R. Carusoec,ed, M.
Casolinoft,ei, C. Cassardoek,el, A. Castellinaek,em, O. Catalanoeh,ed, A.
Cellinoek,em, K. Černýbb, M. Chikawafc, G. Chiritoija, M.J. Christlpf, R.
Colalilloef,eg, L. Contien,ei, G. Cottoek,el, H.J. Crawfordpa, R. Cremoniniel,
A. Creusotcb, A. de Castro Gónzalezpb, C. de la Tailleca, L. del Peralmc, A.
Diaz Damiancc, R. Diesingpb, P. Dinaucourtca, A. Djakonowia, T. Djemilac, A.
Ebersoldtdb, T. Ebisuzakift, J. Eserpb, F. Fenuek,el, S. Fernández-Gonzálezma,
S. Ferrareseek,el, G. Filippatospc, W.I. Finchpc C. Fornaroen,ei, M. Foukaab,
A. Franceschiee, S. Franchinimd, C. Fuglesangna, T. Fujiifg, M. Fukushimafe,
P. Galeottiek,el, E. García-Ortegama, D. Gardiolek,em, G.K. Garipovkb, E.
Gascónma, E. Gazdaph, J. Gencilb, A. Golzioek,el, C. González Alvaradomb, P.
Gorodetzkyft, A. Greenpc, F. Guarinoef,eg, C. Guépinpl, A. Guzmándd, Y.
Hachisuft, A. Haungsdb, J. Hernández Carreteromc, L. Hulettpc, D. Ikedafe, N.
Inouefn, S. Inoueft, F. Isgròef,eg, Y. Itowfk, T. Jammerdc, S. Jeonggb, E.
Jovenme, E.G. Juddpa, J. Jochumdc, F. Kajinoff, T. Kajinofi, S. Kalliaf, I.
Kanekoft, Y. Karadzhovba, M. Kasztelania, K. Katahiraft, K. Kawaift, Y.
Kawasakift, A. Kedadraaa, H. Khalesaa, B.A. Khrenovkb, Jeong-Sook Kimga, Soon-
Wook Kimga, M. Kleifgesdb, P.A. Klimovkb, D. Kolevba, I. Kreykenbohmda, J.F.
Krizmanicpf,pk, K. Królikia, V. Kungelpc, Y. Kuriharafs, A. Kusenkofr,pe, E.
Kuznetsovpd, H. Lahmaraa, F. Lakhdariag, J. Licandrome, L. López Campanoma, F.
López Martínezpb, S. Mackovjakla, M. Mahdiaa, D. Mandátbc, M. Manfrinek,el, L.
Marcelliei, J.L. Marcosma, W. Marszałia, Y. Martínme, O. Martinezhc, K.
Masefa, R. Matevba, J.N. Matthewspg, N. Mebarkiad, G. Medina-Tancoha, A.
Menshikovdb, A. Merinoma, M. Meseef,eg, J. Meseguermd, S.S. Meyerpb, J.
Mimouniad, H. Miyamotoek,el, Y. Mizumotofi, A. Monacoea,eb, J.A. Morales de
los Ríosmc, M. Mastafapd, S. Nagatakift, S. Naitamorab, T. Napolitanoee, J. M.
Nachtmanpi A. Neronovob,cb, K. Nomotofr, T. Nonakafe, T. Ogawaft, S. Ogiofl,
H. Ohmorift, A.V. Olintopb, Y. Onelpi G. Osteriaef, A.N. Otteph, A.
Pagliaroeh,ed, W. Painterdb, M.I. Panasyukkb, B. Panicoef, E. Parizotcb, I.H.
Parkgb, B. Pastircakla, T. Paulpe, M. Pechbb, I. Pérez-Grandemd, F.
Perfettoef, T. Peteroc, P. Picozzaei,ej,ft, S. Pindadomd, L.W. Piotrowskiib,
S. Pirainodd, Z. Plebaniakek,el,ia, A. Pollinioa, E.M. Popescuja, R.
Preveteef,eg, G. Prévôtcb, H. Prietomc, M. Przybylakia, G. Puehlhoferdd, M.
Putisla, P. Reardonpd, M.H.. Renopi, M. Reyesme, M. Ricciee, M.D. Rodríguez
Fríasmc, O.F. Romero Matamalaph, F. Rongaee, M.D. Sabaumb, G. Saccáec,ed, G.
Sáez Canomc, H. Sagawafe, Z. Sahnouneab, A. Saitofg, N. Sakakift, H.
Salazarhc, J.C. Sanchez Balanzarha, J.L. Sánchezma, A. Santangelodd, A. Sanz-
Andrésmd, M. Sanz Palominomb, O.A. Saprykinkc, F. Sarazinpc, M. Satofo, A.
Scagliolaea,eb, T. Schanzdd, H. Schielerdb, P. Schovánekbc, V. Scottief,eg, M.
Serrame, S.A. Sharakinkb, H.M. Shimizufj, K. Shinozakiia, J.F. Sorianope, A.
Sotgiuei,ej, I. Stanja, I. Strharskýla, N. Sugiyamafj, D. Supanitskyha, M.
Suzukifm, J. Szabelskiia, N. Tajimaft, T. Tajimaft, Y. Takahashifo, M.
Takedafe, Y. Takizawaft, M.C. Talaiac, Y. Tamedafp, C. Tenzerdd, S.B.
Thomaspg, O. Tibollahe, L.G. Tkachevka, T. Tomidafh, N. Toneft, S. Toscanoob,
M. Traïcheaa, Y. Tsunesadafl, K. Tsunoft, S. Turrizianift, Y. Uchihorifb, O.
Vaduvescume, J.F. Valdés-Galiciaha, P. Vallaniaek,em, L. Valoreef,eg, G.
Vankova-Kirilovaba, T. M. Venterspj, C. Vigoritoek,el, L. Villaseñorhb, B.
Vlcekmc, P. von Ballmooscc, M. Vrabellb, S. Wadaft, J. Watanabefi, J. Watts
Jr.pd, R. Weigand Muñozma, A. Weindldb, L. Wienckepc, M. Willeda, J. Wilmsda,
D. Winnpm T. Yamamotoff, J. Yanggb, H. Yanofm, I.V. Yashinkb, D. Yonetokufd,
S. Yoshidafa, R. Youngpf, I.S Zguraja, M.Yu. Zotovkb, A. Zuccaro Marchift

aa Centre for Development of Advanced Technologies (CDTA), Algiers, Algeria
ab Dep. Astronomy, Centre Res. Astronomy, Astrophysics and Geophysics (CRAAG),
Algiers, Algeria
ac LPR at Dept. of Physics, Faculty of Sciences, University Badji Mokhtar,
Annaba, Algeria
ad Lab. of Math. and Sub-Atomic Phys. (LPMPS), Univ. Constantine I,
Constantine, Algeria
af Department of Physics, Faculty of Sciences, University of M’sila, M’sila,
Algeria
ag Research Unit on Optics and Photonics, UROP-CDTA, Sétif, Algeria
ah Telecom Lab., Faculty of Technology, University Abou Bekr Belkaid, Tlemcen,
Algeria
ba St. Kliment Ohridski University of Sofia, Bulgaria
bb Joint Laboratory of Optics, Faculty of Science, Palacký University,
Olomouc, Czech Republic
bc Institute of Physics of the Czech Academy of Sciences, Prague, Czech
Republic
ca Omega, Ecole Polytechnique, CNRS/IN2P3, Palaiseau, France
cb Université de Paris, CNRS, AstroParticule et Cosmologie, F-75013 Paris,
France
cc IRAP, Université de Toulouse, CNRS, Toulouse, France
da ECAP, University of Erlangen-Nuremberg, Germany
db Karlsruhe Institute of Technology (KIT), Germany
dc Experimental Physics Institute, Kepler Center, University of Tübingen,
Germany
dd Institute for Astronomy and Astrophysics, Kepler Center, University of
Tübingen, Germany
de Technical University of Munich, Munich, Germany
ea Istituto Nazionale di Fisica Nucleare - Sezione di Bari, Italy
eb Universita’ degli Studi di Bari Aldo Moro and INFN - Sezione di Bari, Italy
ec Dipartimento di Fisica e Astronomia "Ettore Majorana", Universita’ di
Catania, Italy
ed Istituto Nazionale di Fisica Nucleare - Sezione di Catania, Italy
ee Istituto Nazionale di Fisica Nucleare - Laboratori Nazionali di Frascati,
Italy
ef Istituto Nazionale di Fisica Nucleare - Sezione di Napoli, Italy
eg Universita’ di Napoli Federico II - Dipartimento di Fisica "Ettore
Pancini", Italy
eh INAF - Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo, Italy
ei Istituto Nazionale di Fisica Nucleare - Sezione di Roma Tor Vergata, Italy
ej Universita’ di Roma Tor Vergata - Dipartimento di Fisica, Roma, Italy
ek Istituto Nazionale di Fisica Nucleare - Sezione di Torino, Italy
el Dipartimento di Fisica, Universita’ di Torino, Italy
em Osservatorio Astrofisico di Torino, Istituto Nazionale di Astrofisica,
Italy
en Uninettuno University, Rome, Italy
fa Chiba University, Chiba, Japan
fb National Institutes for Quantum and Radiological Science and Technology
(QST), Chiba, Japan
fc Kindai University, Higashi-Osaka, Japan
fd Kanazawa University, Kanazawa, Japan
fe Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Japan
ff Konan University, Kobe, Japan
fg Kyoto University, Kyoto, Japan
fh Shinshu University, Nagano, Japan
fi National Astronomical Observatory, Mitaka, Japan
fj Nagoya University, Nagoya, Japan
fk Institute for Space-Earth Environmental Research, Nagoya University,
Nagoya, Japan
fl Graduate School of Science, Osaka City University, Japan
fm Institute of Space and Astronautical Science/JAXA, Sagamihara, Japan
fn Saitama University, Saitama, Japan
fo Hokkaido University, Sapporo, Japan
fp Osaka Electro-Communication University, Neyagawa, Japan
fq Nihon University Chiyoda, Tokyo, Japan
fr University of Tokyo, Tokyo, Japan
fs High Energy Accelerator Research Organization (KEK), Tsukuba, Japan
ft RIKEN, Wako, Japan
ga Korea Astronomy and Space Science Institute (KASI), Daejeon, Republic of
Korea
gb Sungkyunkwan University, Seoul, Republic of Korea
ha Universidad Nacional Autónoma de México (UNAM), Mexico
hb Universidad Michoacana de San Nicolas de Hidalgo (UMSNH), Morelia, Mexico
hc Benemérita Universidad Autónoma de Puebla (BUAP), Mexico
hd Universidad Autónoma de Chiapas (UNACH), Chiapas, Mexico
he Centro Mesoamericano de Física Teórica (MCTP), Mexico
ia National Centre for Nuclear Research, Lodz, Poland
ib Faculty of Physics, University of Warsaw, Poland
ja Institute of Space Science ISS, Magurele, Romania
ka Joint Institute for Nuclear Research, Dubna, Russia
kb Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University,
Russia
kc Space Regatta Consortium, Korolev, Russia
la Institute of Experimental Physics, Kosice, Slovakia
lb Technical University Kosice (TUKE), Kosice, Slovakia
ma Universidad de León (ULE), León, Spain
mb Instituto Nacional de Técnica Aeroespacial (INTA), Madrid, Spain
mc Universidad de Alcalá (UAH), Madrid, Spain
md Universidad Politécnia de madrid (UPM), Madrid, Spain
me Instituto de Astrofísica de Canarias (IAC), Tenerife, Spain
na KTH Royal Institute of Technology, Stockholm, Sweden
oa Swiss Center for Electronics and Microtechnology (CSEM), Neuchâtel,
Switzerland
ob ISDC Data Centre for Astrophysics, Versoix, Switzerland
oc Institute for Atmospheric and Climate Science, ETH Zürich, Switzerland
pa Space Science Laboratory, University of California, Berkeley, CA, USA
pb University of Chicago, IL, USA
pc Colorado School of Mines, Golden, CO, USA
pd University of Alabama in Huntsville, Huntsville, AL; USA
pe Lehman College, City University of New York (CUNY), NY, USA
pf NASA Marshall Space Flight Center, Huntsville, AL, USA
pg University of Utah, Salt Lake City, UT, USA
ph Georgia Institute of Technology, USA
pi University of Iowa, Iowa City, IA, USA
pj NASA Goddard Space Flight Center, Greenbelt, MD, USA
pk Center for Space Science & Technology, University of Maryland, Baltimore
County, Baltimore, MD, USA
pl Department of Astronomy, University of Maryland, College Park, MD, USA
pm Fairfield University, Fairfield, CT, USA