∎

11institutetext:  🖂Jian Li
<EMAIL_ADDRESS>
${}^{\textbf{1}}$School of Artificial Intelligence, Beijing University of
Posts and Telecommunications, Beijing 100876, China.
${}^{\textbf{2}}$College of Information Science and Engineering, ZaoZhuang
University, ZaoZhuang Shandong 277160, China.
${}^{\textbf{3}}$School of Cyberspace Security Security, Beijing University of
Posts Telecommunications, Beijing 100876, China.
${}^{\textbf{4}}$Information Security Center, State Key Laboratory of
Networking and Switching Technology, Beijing University of Post and
Telecommunications, Beijing 100876, China.
${}^{\textbf{5}}$GuiZhou University, Guizhou Provincial Key Laboratory of
Public Big Data, Guizhou Guiyang, 550025, China.

# Quantum adversarial metric learning model based on triplet loss function

Yan-Yan Hou ${}^{\textbf{1,2}}$ Jian Li ${}^{\textbf{3}}$ Xiu-Bo Chen
${}^{\textbf{4,5}}$ Chong-Qiang Ye ${}^{\textbf{1}}$

(Received: date / Accepted: date)

###### Abstract

Metric learning plays an essential role in image analysis and classification,
and it has attracted more and more attention. In this paper, we propose a
quantum adversarial metric learning (QAML) model based on the triplet loss
function, where samples are embedded into the high-dimensional Hilbert space
and the optimal metric is obtained by minimizing the triplet loss function.
The QAML model employs entanglement and interference to build superposition
states for triplet samples so that only one parameterized quantum circuit is
needed to calculate sample distances, which reduces the demand for quantum
resources. Considering the QAML model is fragile to adversarial attacks, an
adversarial sample generation strategy is designed based on the quantum
gradient ascent method, effectively improving the robustness against the
functional adversarial attack. Simulation results show that the QAML model can
effectively distinguish samples of MNIST and Iris datasets and has higher
$\epsilon$-robustness accuracy over the general quantum metric learning. The
QAML model is a fundamental research problem of machine learning. As a
subroutine of classification and clustering tasks, the QAML model opens an
avenue for exploring quantum advantages in machine learning.

###### Keywords:

Metric learning hybrid quantum-classical algorithm quantum machine learning

##  1 Introduction

Machine learning has developed rapidly in recent years and is widely used in
artificial intelligence and big data fields. Quantum computing can efficiently
process data in exponentially sizeable Hilbert space and is expected to
achieve dramatic speedups in solving some classical computational problems.
Quantum machine learning, as the interplay between machine learning and
quantum physics, brings unprecedented promise to both disciplines. On the one
hand, machine learning methods have been extended to quantum world and applied
to the data analysis in quantum physicscong2019quantum . On the other hand,
quantum machine learning exploits quantum properties, such as entanglement and
superposition, to revolutionize classical machine learning algorithms and
achieves computational advantages over classical
algorithmsbenedetti2019parameterized . Metric Learning is the core problem of
some machine learning taskschen2018adversarial , such as $k$-nearest neighbor,
support vector machines, radial basis function networks, and $k$-means
clustering. Its core work is to construct an appropriate distance metric that
maximizes the similarities of samples of the same class and minimizes the
similarities of samples from different classes. Linear and nonlinear methods
can be used to implement metric learning. In general, linear models have a
limited number of parameters and are unsuitable for characterizing high-order
features of samples. Recently, neural networks have been adopted to establish
nonlinear metric learning models, and promising results have been achieved in
face recognition and feature matching.

Classical metric learning models usually extract low-dimensional
representations of samples, which will lose some details of samples. Quantum
states are in high-dimensional Hilbert spaces, and their dimensions grow
exponentially with the number of qubits. This quantum enables quantum models
to learn high-dimensional representations of samples without explicitly
invoking a kernel function. A parameterized quantum circuit is used to map
samples in high-dimensional Hilbert space. The optimal metric model is
obtained by optimizing the loss function based on Hilbert-Schmidt distances.
With the increase of the the dimension, this speed-up advantage will become
more and more pronounced, and it is expected to achieve exponential growth in
computing speeds. In recent years, researchers began to study how to adopt
quantum methods to implement metric learning. Lloydlloyd2020quantum  firstly
proposed a quantum metric learning model based on hybrid quantum-classical
algorithms. A parameterized quantum circuit is used to map samples in high-
dimensional Hilbert space. The optimal metric model is obtained by optimizing
the loss function based on Hilbert-Schmidt distances. This model achieves
better effects in classification tasks. Nhatnghiem2020unified  introduced
quantum explicit and implicit metric learning approaches from the perspective
of whether the target space is known or not. The research establishes the
relationship between quantum metric learning and other quantum supervised
learning models. The above two algorithms mainly focus on classification
tasks. Metric learning is a fundamental problem in machine learning, which can
be applied not only to classification but also to clustering, face
recognition, and other issues. In our research, we are devoted to constructing
a quantum metric learning model that can serve various machine learning tasks.

Angular distance is a vital metric that quantifies the included angle between
normalized samplesmao2019metric . Angular distance focuses on the difference
in the direction of samples and is more robust to the variation of local
featurewang2017deep ,duan2018deep . Considering the similarities between
angular distances of classical data and inner products of quantum states, we
design a quantum adversarial metric learning (QAML) model based on inner
product distances, which is more suitable for image-related tasks. Unlike
other quantum metric learning models, the QAML model maps samples from
different classes into quantum superposition states and utilizes simple
interface circuits to compute metric distances for multiple sample pairs in
parallel. Furthermore, quantum systems in high-dimensional Hilbert space have
counter-intuitive geometrical propertiesliu2020vulnerability . The QAML model
using only natural samples is vulnerable to adversarial attacks, under which
some samples are closer to the false class, so the model is easy to make wrong
judgementsmadry2017towards . To solve this issue, we construct adversarial
samples based on natural samples. The model’s robustness is improved by the
alternative train of natural and adversarial samples. Our work has two main
contributions:(i) We explore a quantum method to compute the triplet loss
function, which utilizes quantum superposition states to calculate sample
distances in parallel and reduce the demand for quantum resources. (ii) We
design an adversarial samples generation strategy based on the quantum
gradient ascent, and the robustness of the QAML model is significantly
improved by alternatively training generated adversarial samples and natural
samples. Simulation results show that the QAML model separates samples by a
larger margin and has better robustness for functional adversarial attacks
than general quantum metric learning models.

The paper is organized as follows. Section 2 gives the basic method of the
QAML model. Section 3 verifies the performances of the QAML model. Finally, we
get a conclusion and discuss the future research directions.

##  2 Quantum adversarial metric learning

###  2.1 Preliminary theory

Triplet loss function is a widely used strategy for metric
learningsalakhutdinov2007learning , commonly used in image retrieval and face
recognition. A triplet set ${(x_{i}^{a},x^{p}_{i},x^{n}_{i})}$ consists of
three samples from two classes, where anchor sample $x_{i}^{a}$ and positive
sample $x^{p}_{i}$ belong to the same class, and negative sample $x^{n}_{i}$
comes from another class. The goal of metric learning based on triplet loss
function is to find the optimal embedded representation space, in which
positive sample pairs ${(x_{i}^{a},x^{p}_{i})}$ are pulled together and
negative sample pairs ${(x_{i}^{a},x^{n}_{i})}$ are pushed away. Fig.1 shows
sample space change in the metric learning process. As we can see, samples
from different classes become linearly separable through metric learning.
Fig.2 shows the schematic of the metric learning model based on triplet loss
function. Firstly, the model prepares multiple triplet sets, and one triplet
set ${(x_{i}^{a},x^{p}_{i},x^{n}_{i})}$ is sent to convolutional neural
networks (CNN), where three CNN with the same structure and parameters are
needed. Each CNN acts on one sample of the triplet set to extract its
features. The triplet loss function is obtained by computing metric distances
for multiple sample pairs of triplet sets. In the learning process, the
optimal parameters of CNN are obtained by minimizing the triplet loss
function.

Figure 1: Sample space change in metric learning process. Before metric
learning, the distances between negative sample pairs are smaller, and samples
from different classes are difficult to separate by linear functions. After
metric learning, the distances between negative sample pairs become larger,
and a large margin separates samples from different classes. Linear functions
can easily separate positive and negative samples. Figure 2: The schematic of
the metric learning model based on triplet loss function. A triplet set
includes an anchor sample, a positive sample, and a negative sample. The input
consists of a batch of triplet sets, and only one triplet set serves as input
in each iteration. Three CNN with the same structure and parameters are used
to map the triplet set into the embedded representation space. CNN, consisting
of multiple convolutions, pooling, and fully connected layers, is responsible
for extracting the features of samples. The triplet loss function is further
constructed based on the extracted features.

Let one batch samples include $N_{1}$ triplet sets. The triplet loss function
is

$\displaystyle\begin{array}[]{ccc}L=\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}[D(g(x^{a}_{i}),g(x^{p}_{i}))-D(g(x^{a}_{i}),g(x^{n}_{i}))+\mu]_{+},\end{array}$
(1)

where $g(\cdot)$ represents the function mapping input samples to the embedded
representation space, $D(\cdot,\cdot)$ denotes the distance between a sample
pair in the embedded representation space, and $[\ \cdot\ ,\ \cdot\
]_{+}=max(0,\ \cdot\ )$ represents the hinge loss function. The goal of metric
learning is to learn a metric that makes the distances between negative sample
pairs greater than the distance between the corresponding positive sample
pairs and satisfies the specified margin $\mu\in\mathbb{R}^{+}$mao2019metric .
In the triplet loss function, $D(g(x^{a}_{i}),g(x^{p}_{i}))$ penalizes the
positive sample pair $(x^{a}_{i},x^{p}_{i})$ that is too far apart, and
$D(g(x^{a}_{i}),g(x^{n}_{i}))$ penalizes the negative sample pair
$(x^{a}_{i},x^{p}_{i})$ whose distance is less than the margin $\mu$.

Metric learning can adopt various distance metric methods. Angular distance
metric is robust to image illumination and contrast variation wang2017deep ,
which is an efficient way for metric learning tasks. In this method, samples
need to be normalized to unit vectors in advance. The distance between a
positive sample pair is

$\displaystyle\begin{array}[]{ccc}D(g(x^{a}_{i}),g(x^{p}_{i}))=1-\frac{|g(x^{a}_{i})\cdot
g(x^{p}_{i})|}{||g(x^{a}_{i})||_{2}||g(x^{p}_{i})||_{2}},\end{array}$ (2)

where $|\ |$ and $||\ ||_{2}$ represent $l_{1}$-norm and $l_{2}$-norm,
respectively, and $\cdot$ denotes the inner product operation for two vectors.
The distance between negative sample pairs can be calculated in the same way.

###  2.2 Framework of quantum metric learning model

For most machine learning tasks, it is often challenging to adopt simple
linear functions to distinguish samples of different classes. According to
kernel theoryblank2020quantum , samples in high-dimensional feature space have
better distinguishability. Classical machine learning algorithms usually adopt
kernel methods to map samples to high-dimensional feature space, where the
mapped samples can be separated by simple linear functions. Quantum states
with $n$-qubits are in $2^{n}$-dimensional Hilbert space, where quantum
systems characterize the nonlinear features of data and efficiently process
data through a series of linear unitary operations.

In the QAML model, samples should be firstly mapped into quantum systems by
qubit encoding. The Hilbert space after encoding usually does not correspond
to the optimal space for separating samples of different classes. To search
for the optimal Hilbert space, the QAML model performs parameterized quantum
circuits $W(\theta)$ on the encoded statesgrant2018hierarchical . As different
variable parameters $\theta$ correspond to different mapping spaces, we can
search the optimal space by modifying parameters
$\theta=(\theta_{1}^{1},...,\theta_{i}^{j})$. As long as $W(\theta)$ has
strong expressivity, we can find the optimal Hilbert space by optimizing the
loss function of metric learningperez2020data ; schuld2021effect . $W(\theta)$
with different structures and layers have different expressivity. The more
layers $W(\theta)$ has, the stronger the expressivity, and the easier it is to
find the optimal metric space.

The classical metric learning model based on triplet loss function requires
three identical CNN to map triplet sets $(x_{i}^{a},x_{i}^{p},x_{i}^{n})$ into
the novel Hilbert space. To reduce the demand for quantum resources, we
construct a quantum superposition state to represent one triplet set so that a
triplet set only needs one $W(\theta)$ to transform it into Hilbert space. The
core work of the building loss function is to compute inner products between
sample pairs, but $W(\theta)$ and subsequent conjugate operation
$W^{\dagger}(\theta)$ counteract each other’s effects. To solve this issue, we
add a repeated encoding operation after $W(\theta)$. It is worth mentioning
that the repeated encoding operation is also conducive to the construction of
high-dimensional features of samples.

The QAML model is mathematically represented as the minimization of the loss
function with respect to the parameters $\theta$. The triplet loss function
consists of metric distances for positive and negative sample pairs, so the
kernel work of the QAML model is constructing the metric distances for sample
pairs in the transformed Hilbert space. The mapping samples
$h(x_{i}^{a})/||h(x_{i}^{a})||_{2}$ and $h(x_{i}^{p})/||h(x_{i}^{p})||_{2}$ of
Equ.2 are replaced by the quantum states of $x_{i}^{a}$ and $x_{i}^{p}$, then
the second term of Equ.2 is converted to the inner product between quantum
states of the positive sample pair $(x_{i}^{a},x_{i}^{p})$, which can be got
by the method of the Hadamard classifierblank2020quantum . The triplet loss
function can be viewed as the weighted sum of the inner product of sample
pairs $(x_{i}^{a},x_{i}^{p})$ and the inner product of sample pairs
$(x_{i}^{a},x_{i}^{n})$. With the help of ancilla registers, the triplet set
can be prepared in superposition states form. According to the entanglement
property of superposition states, the triplet loss function can be implemented
with one parameterized quantum circuit. Then, the triplet loss function value
is transmitted to a classical optimizer, and parameters are optimized until
the optimal metric is obtained. The QAML model constructs adversarial samples
according to the gradient of natural samples and trains alternatively natural
and adversarial samples to improve the model’s robustness against adversarial
attacks. The schematic of the QAML model is shown in Fig.3.

Figure 3: Overview of quantum adversarial metric learning (QAML) model. Panel
(a) shows the framework of quantum adversarial metric learning. $Reg.s$ is the
sample register that stores triplet sets, and $Reg.1$ and $Reg.2$ are ancilla
registers used distinguishing different samples. The model firstly adopts
principal component analysis (PCA) to reduce the input dimension.
Subsequently, anchor, negative and positive samples are encoded into a quantum
superposition state by controlled qubit encoding. The transformation of
Hilbert space is implemented by parameterized quantum circuit $W(\theta)$ and
the subsequent qubit encoding $U_{1}(x_{i})$. Finally, Hadamard and
measurement operations act on ancilla registers to simultaneously compute the
inner products for the positive and negative sample pairs, and the triplet
loss function is further obtained. In each iteration, the parameters $\theta$
are updated by optimizing the triplet loss function with a classical
optimizer. Panel (b) shows the quantum dimension reduction circuit to reduce
the number of output qubits. In each module, only one qubit is measured, and
the controlled unitary based on its measurement result acts on another qubit.
Panel (c) shows another case of the QAML model, where adversarial samples are
built and added to the training process.
$V^{{}^{\prime}}(\lambda\nabla_{i}^{a})$ is the unitary operation based on the
gradient of anchor sample $x_{i}^{a}$ and acts on the encoded quantum states
to produce its adversarial sample. In the QAML model training process, natural
and adversarial samples alternatively serve as input.

## References

  * (1) I. Cong, S. Choi, M.D. Lukin, Nature Physics 15(12), 1273 (2019)
  * (2) M. Benedetti, E. Lloyd, S. Sack, M. Fiorentini, Quantum Science and Technology 4(4), 043001 (2019)
  * (3) S. Chen, C. Gong, J. Yang, X. Li, Y. Wei, J. Li, arXiv preprint arXiv:1802.03170 (2018)
  * (4) S. Lloyd, M. Schuld, A. Ijaz, J. Izaac, N. Killoran, arXiv preprint arXiv:2001.03622 (2020)
  * (5) N.A. Nghiem, S.Y.C. Chen, T.C. Wei, arXiv preprint arXiv:2010.13186 (2020)
  * (6) C. Mao, Z. Zhong, J. Yang, C. Vondrick, B. Ray, Advances in Neural Information Processing Systems 32 (2019)
  * (7) J. Wang, F. Zhou, S. Wen, X. Liu, Y. Lin, in _Proceedings of the IEEE international conference on computer vision_ (2017), pp. 2593–2601
  * (8) Y. Duan, W. Zheng, X. Lin, J. Lu, J. Zhou, in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ (2018), pp. 2780–2789
  * (9) N. Liu, P. Wittek, Physical Review A 101(6), 062331 (2020)
  * (10) A. Madry, A. Makelov, L. Schmidt, D. Tsipras, A. Vladu, arXiv preprint arXiv:1706.06083 (2017)
  * (11) R. Salakhutdinov, G. Hinton, in _Artificial Intelligence and Statistics_ (PMLR, 2007), pp. 412–419
  * (12) C. Blank, D.K. Park, J.K.K. Rhee, F. Petruccione, npj Quantum Information 6(1), 1 (2020)
  * (13) E. Grant, M. Benedetti, S. Cao, A. Hallam, J. Lockhart, V. Stojevic, A.G. Green, S. Severini, npj Quantum Information 4(1), 1 (2018)
  * (14) A. Pérez-Salinas, A. Cervera-Lierta, E. Gil-Fuster, J.I. Latorre, Quantum 4, 226 (2020)
  * (15) M. Schuld, R. Sweke, J.J. Meyer, Physical Review A 103(3), 032430 (2021)
  * (16) C. Zoufal, A. Lucchi, S. Woerner, npj Quantum Information 5(1), 1 (2019)
  * (17) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J.M. Chow, J.M. Gambetta, Nature 549(7671), 242 (2017)
  * (18) T. Miyato, S.i. Maeda, M. Koyama, S. Ishii, IEEE transactions on pattern analysis and machine intelligence 41(8), 1979 (2018)
  * (19) A. Kurakin, I. Goodfellow, S. Bengio, arXiv preprint arXiv:1611.01236 (2016)
  * (20) J.R. McClean, M.E. Kimchi-Schwartz, J. Carter, W.A. De Jong, Physical Review A 95(4), 042308 (2017)
  * (21) G.E. Crooks, arXiv preprint arXiv:1905.13311 (2019)
  * (22) M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, N. Killoran, Physical Review A 99(3), 032331 (2019)
  * (23) K. Mitarai, M. Negoro, M. Kitagawa, K. Fujii, Physical Review A 98(3), 032309 (2018)
  * (24) V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, M.S. Alam, S. Ahmed, J.M. Arrazola, C. Blank, A. Delgado, S. Jahangiri, et al., arXiv preprint arXiv:1811.04968 (2018)
  * (25) M.C. Mukkamala, M. Hein, in _International conference on machine learning_ (PMLR, 2017), pp. 2545–2553
  * (26) J. Guan, W. Fang, M. Ying, CoRR (2020)


11footnotetext: * Corresponding author.
E-mail address<EMAIL_ADDRESS>(Fang Yang)

# Quantifying model uncertainty for the observed non-Gaussian data by the
Hellinger distance

Yayun Zheng1, Fang Yang ∗2, Jinqiao Duan3, Jürgen Kurths4,5

1 School of Mathematical Sciences, Jiangsu University, Zhenjiang 212013,
China.
2 School of Mathematics and Statistics and Center for Mathematical Science,
Huazhong University of Science and Technology, Wuhan, 430074, China.
3 Department of Applied Mathematics, Illinois Institute of Technology,
Chicago, IL 60616, USA.
4 Potsdam Institute for Climate Impact Research, Potsdam, 14473, Germany.
5Department of Physics, Humboldt University, Berlin, 12489, Germany

###### Abstract

Mathematical models for complex systems under random fluctuations often
certain uncertain parameters. However, quantifying model uncertainty for a
stochastic differential equation with an $\alpha$-stable Lévy fluctuation is
still lacking. Here, we propose an approach to infer all the uncertain non-
Gaussian parameters and the other system parameter by minimizing the Hellinger
distance over the parameter space. The Hellinger distance investigate the
similarity between an empirical probability density from non-Gaussian
observations and a solution of the nonlocal Fokker-Planck equation. Numerical
experiments verify that the method is feasible for estimating single and
multiple parameters. Meanwhile, we could find an optimal estimation interval
of the estimated parameters. This method is beneficial for parameter
estimation of data-driven dynamical system models with non-Gaussian
distribution, such as abrupt climate changes in the Dansgaard-Oeschger events.

Keywords: Non-Gaussian observations; Parameters estimation; Hellinger
distance; Probability density

##  1 Introduction

Complex systems under influences of random fluctuations also have uncertain
parameters [1]. An important problem in modeling such random processes by
stochastic differential equations (SDEs) is to estimate uncertain parameters
from observations of the stochastic paths.

A Brownian motion has properties of continuous sample paths, normal diffusion
and light tail (probability density decays exponentially), theoretical results
on parametric estimations for SDEs driven by Brownian motion are relatively
well developed. The Gaussian kernel density estimator [2, 3] and the Bayesian
estimator[4, 5] are well-known approaches for parameters estimation of a drift
function when the observations of complete paths are available. The
nonparametric estimation method based on Kramers-Moyal coefficients [6], the
statistical definitions of conditional first and second moments [7] and
variational formulation of the stationary Fokker-Planck equation [8] could
provide an expression for the drift function and the diffusion one.

However, various complex phenomena involve non-Gaussian fluctuations, with
properties such as intermittent jumps, anomalous diffusion, and heavy tail
(probability density decays with power law) distribution. A heavy-tailed
distribution, like the Lévy distribution, is characterized by a high
likelihood for extreme events, compared to a normal distribution. For
instance, Ditlevsen shows that the paleoclimatic records for Dansgaard-
Oeschger events have a strong non-Gaussian distribution [9]. The protein
production occurs in bursts which are observed during a genetic regulation
[10]. Meanwhile, experimental studies find that Lévy flights are an optimal
pattern when the prey is sparsely and randomly distributed for open-ocean
predatory fish[11]. A Lévy process is also used in other scientific domains,
for example, it has been shown that certain stock price has statistical
properties that are compatible with a Lévy random walk[12]. Additionally, in
the field of cognitive research, a few studies provide evidence of Lévy
processes, e.g., to search and cluster in semantic memory [13] and human
decision making[14].

An $\alpha$-stable Lévy process is thought to be an appropriate model for a
non-Gaussian heavy-tailed process. The problem of parameter estimation for a
stochastic system driven by $\alpha$-stable Lévy process has caught broad
attention. In general, the $p$-th moment of an $\alpha$-stable Lévy random
variable is finite if and only if $p<\alpha$ ($0<\alpha<2$), so it does not
have second moments. Meanwhile, an $\alpha$-stable Lévy process lacks of an
explicit expression of the stationary probability density. Subject to these
disadvantages, unfortunately, the existing research results for parameters
estimation of Brownian motion can not be used in non-Gaussian situations. In
some special cases, it is possible to infer parameters for the Ornstein-
Uhlenbeck processes driven by an $\alpha$-stable Lévy process, i.e., the drift
function is known to be linear. Hu and Long _et al._ addresses a trajectory
fitting estimator and a least squares estimator for the this parameter
assuming other parameters are known[15, 16]. However, in reality, the drift
function is seldom known, Long _et al._ focus on estimating the drift function
by using the nonparametric Nadaraya-Watson estimator [17].

The non-Gaussian index $\alpha$ plays a decisive role in the construction of a
Lévy process. An $\alpha$-stable Lévy process has larger jumps with lower jump
probabilities when $\alpha$ is small ($0<\alpha<1$), while it has smaller
jumps with higher jump frequencies for large $\alpha$ values ($1<\alpha<2$).
The special cases for $\alpha=1$ and $\alpha=2$ correspond to the Cauchy
process and the Brownian motion, respectively. Therefore, the estimation of
the parameter $\alpha$ is extremely important. There are some simple and
straightforward approaches to learn this $\alpha$ from the path observation,
such as the slope of the log-log linear regression [18] or the Hill estimator
[19]. These methods do not assume a parametric form for the entire
distribution function, but focus only on the tail behavior. However, the true
tail behavior of Lévy distribution is visible only for extremely large data
sets, or it is a challenge to choose the right value of the largest order
statistics.

An alternative method for estimating $\alpha$ is relied on the characteristic
function[20]. Based on the ergodic theory and sample characteristic functions,
Cheng _et al._ estimated $\alpha$ and the other parameters by matching the
empirical characteristic function with the corresponding theoretical one[21].
Furthermore, the generalized Lévy characteristic function could also define a
Hurst exponent to take the expression of $\alpha=\frac{1}{H}$ in the
probability space[22]. We note that another method for an $\alpha$-stable
Ornstein-Uhlenbeck process is the maximum likelihood estimation. Chen _et al._
choose a mixture of Cauchy and Gaussian distribution to approximate the
probability density of the $\alpha$-stable Ornstein-Uhlenbeck distribution
[23]. Besides, a method of numerical optimization is devised in [24], where
the mean exit time or the escape probability is observed to estimate the
uncertain parameter and other system parameters. It is based on solving an
inverse problem for a deterministic, nonlocal partial differential equation.

In response to the existing challenge of the parameters estimation for a SDE
driven by $\alpha$-stable Lévy process, we propose an approach to estimate
$\alpha$ and other system parameters simultaneously for a general
$\alpha$-stable Lévy system. The method is based on the probability density
from an observation data set with heavy-tailed distribution. The probability
density function for SDE driven by an $\alpha$-stable Lévy process satisfies a
deterministic, nonlocal differential equation with an initial condition, i.e.,
nonlocal Fokker-Planck equation. Fokker-Planck equations are deterministic
equations describing how the probability density functions propagate and
evolve. Recently, Gao _et al._ developed a fast and accurate numerical
algorithm to simulate nonlocal Fokker-Planck equations under either absorbing
or natural conditions [25]. In theory, we derived the Fokker-Planck equations
for Marcus SDEs driven by Lévy processes in high dimensional[26].

In present paper, we consider the parameter estimation problem of an
$\alpha$-stable Lévy stochastic dynamical system containing uncertain
parameters. In Section 2, we propose a method of estimating the uncertain
parameters based on the Hellinger distance of the probability densities. In
Section 3, we present some simulation results of estimation for single and
multiple parameters by minimizing the Hellinger distance. Finally we give some
concluding and future works in Section 4.

##  2 Methods

We consider a dynamical system with heavy-tailed uncertainty, which could be
modeled by a stochastic process $X(t)$

${\rm d}X(t)=f(X(t),\theta){\rm d}t+\epsilon{\rm d}L^{\alpha}(t),\quad
X(0)=x_{0}\in\mathbb{R}^{1},$ (1)

where the drift function $f(x,\theta)$ has the uncertain parameter $\theta$,
and a scalar symmetric $\alpha$-stable Lévy process $L^{\alpha}(t)$ with the
non-Gaussian index $0<\alpha<2$ is defined in a probability space
$(\Omega,\mathcal{F},\mathbb{P})$. The parameter $\epsilon$ is the non-
negative $\alpha$-stable Lévy noise intensity.

A scalar symmetric $\alpha$-stable Lévy process is characterized by a
generating triplet ($b,Q,\nu_{\alpha}$), a linear coefficient $b$, a diffusion
parameter $Q$, and a nonnegative Borel measure $\nu_{\alpha}$. This jump
measure $\nu_{\alpha}$ is defined on
$(\mathbb{R}^{1},\mathfrak{B}(\mathbb{R}^{1}))$[27] by:

$\nu_{\alpha}=\frac{C_{\alpha}dy}{|y|^{1+\alpha}},$

with $0<\alpha<2$ and
$C_{\alpha}=\frac{\alpha}{2^{1-\alpha}\sqrt{\pi}}\frac{\Gamma(\frac{1+\alpha}{2})}{\Gamma(1-\frac{\alpha}{2})}$.
In this paper, we consider an $\alpha$-stable Lévy process with a triplet
$(0,0,\nu)$, i.e., a pure jump process.

For $0<\alpha<2$, the $\alpha$-stable Lévy process $L^{\alpha}(t)$ has a
heavy-tailed distribution[28]

$\mathbb{P}(|L^{\alpha}(t)|>y)\sim\frac{1}{y^{\alpha}},$

as the tail estimate decays in a power law. Therefore $\alpha$ is also called
the power parameter. The tail behavior is different from the Brownian motion
with light tail, as the tail decays exponentially.

We assume that the drift term $f$ is local Lipschitz continuous. Then the SDE
(1) has a unique solution [27]. The conditional probability density
$p(x,t|x_{0},0)\triangleq p(X(t)=x|X(0)=x_{0})$ represents the density of the
$X(t)$ given a value $x_{0}$ at initial time. For convenience, we drop the
initial condition and simply denote it by $p(x,t)$. There exists sufficient
condition for the existence and regularity of the probability density $p(x,t)$
for some SDEs driven by Lévy processes. The existence is based on Malliavin
calculus with jumps under Hörmander’s condition, see Refs. [29, 30, 31] and
the references therein for more details.

We see that the stochastic process $X_{t}$ in Eq.(1) under Lévy noise depends
on the following parameters. The first one is an uncertain system parameter
$\theta$. In general, the estimated parameter $\theta$ plays a key role in the
system model, which could be a bifurcation parameter inducing a transition
between states. The control parameter $\theta$ could be a greenhouse factor in
the case of the energy balance model [32], or a freshwater forcing strength in
the thermohaline circulation one[33]. Besides, there are uncertain Lévy
parameters: the non-Gaussian index $\alpha$ and the Lévy noise intensity
$\epsilon$.

Let us assume that we have access to a set of observations
$y=(y_{1},y_{2},\cdots,y_{n})$, which are the version of the process $X_{t}$
with a non-Gaussian distribution sampled at discrete times $t_{k}\in[0,T]$ for
$k=1,2,\cdots,n$, i.e., $y_{k}=X_{t_{k}}$ for $k=1,2,\cdots,n$. In this paper,
we will discuss the problem of estimating the parameters $\alpha$, $\epsilon$
and $\theta$ simultaneously using the observations
$y=(y_{1},y_{2},\cdots,y_{n})$.

(a)

(b)

Figure 1: (a) The stochastic sample path of the process $X_{t}$ (Eq.(4)) with
$\alpha=1.7$, $\epsilon=0.3$ and $\theta=1$. (b) The empirical probability
density $p_{d}$ of the stochastic trajectory (a).

To reach this purpose, we would like to introduce the Hellinger distance. It
is used to quantify the similarity between two probability distributions [34].
The Hellinger distance between two probability density functions $p(x)$ and
$p_{d}(x)$ is

$H^{2}(p,p_{d})=\frac{1}{2}\int_{\mathbb{R}^{1}}\left(\sqrt{p(x)}-\sqrt{p_{d}(x)}\right)^{2}dx.$

The Hellinger distance $H$ satisfies the property: $0\leqslant
H(p,p_{d})\leqslant 1$. Here, $p_{d}$ is the empirical probability density
from an observation data set $y=(y_{1},y_{2},\cdots,y_{n})$. The probability
density function $p(x)$ is a solution of the nonlocal Fokker-Planck equation
$p(x,t)$ at time $t$.

$\frac{\partial}{\partial t}p(x,t)=-\frac{\partial}{\partial
x}(f(x,\theta)p(x,t))+\epsilon^{\alpha}\int_{\mathbb{R}^{1}\backslash\\{0\\}}\left[p(x+y,t)-p(x,t)-I_{|y|<1}(y)\;y\frac{\partial}{\partial
x}p(x,t)\right]\nu_{\alpha}(dy).$ (2)

The integral part in the right hand side is actually the nonlocal Laplacian
operator, reflecting the non-Gaussian Lévy fluctuations [35]. The equation
fulfills an initial condition

$\lim_{t\rightarrow 0}p(x,t|x_{0},0)=\delta(x-x_{0}).$

We consider that the observation set comes from an $\alpha$-stable Lévy
distribution $p(x,t)$. Associated with each probability density is the
parameters set $\lambda=\\{\theta,\alpha,\epsilon\\}\in\Theta$, where $\Theta$
is called the parameter space, a finite-dimensional subset of the Euclidean
space. Evaluating the Hellinger distance at the observed data set $y$ gives an
objective function

$G(\lambda)=H^{2}(p(x,\lambda),p_{d}(x)).$

The Hellinger distance estimation aims to find the value of the model
parameters that minimize the objective function over the parameter space
$\Theta$, that is

$\hat{\lambda}=\arg\min_{\lambda\in\Theta}{G}(\lambda).$

To address the probability density $p(x,\lambda)$, we use the numerical
algorithm of Gao _et al._[25] to solve the nonlocal differential equation in
Eq.(2) under the absorbing condition. This absorbing condition means that the
probability of finding “partical” $X_{t}$ outside the finite interval
$D=(a,b)$ is zero. We decompose the integral part of Eq.(2) into three parts
$\int_{\mathbb{R}^{1}}=\int^{a-x}_{-\infty}+\int^{a-x}_{b-x}+\int^{\infty}_{b-x}$
in $\mathbb{R}^{1}$ and analytically evaluate the first and third integrals,
then Eq.(2) changes to

$\displaystyle\frac{\partial}{\partial t}p(x,t)$
$\displaystyle=-\frac{\partial}{\partial
x}(f(x)p(x,t))-\frac{\epsilon^{\alpha}C_{\alpha}}{\alpha}\left[\frac{1}{(x-a)^{\alpha}}+\frac{1}{(b-x)^{\alpha}}\right]p(x,t)$
$\displaystyle+\epsilon^{\alpha}C_{\alpha}\int_{a-x}^{b-x}\frac{p(x+y,t)-p(x,t)-I_{|y|<1}(y)y\frac{\partial}{\partial
x}p(x,t)}{|y|^{1+\alpha}}dy,$ (3)

for $x\in(a,b)$. The non-Gaussian index $\alpha\in(0,2)$ and the Lévy
intensity $\epsilon\in(0,1]$.

(a)

(b)

Figure 2: (a) The Hellinger distance estimation of $\alpha$. The estimated
value $\hat{\alpha}=1.738$ corresponds to the minimum objective function
$G_{\alpha}=0.001399$ (red $\star$). (b) The estimated density function $p(x)$
with $\hat{\alpha}=1.738$, $\theta=1$ and $\epsilon=0.3$ compared with the
empirical density $p_{d}$ in Fig.1(b).

##  3 Numerical Experiments

We now explore how the Hellinger distance can be used to estimate the
parameters of SDE driven by a symmetric $\alpha$-stable Lévy process. We
consider the following example

${\rm d}X(t)=(-\theta X(t)^{3}+X(t)){\rm d}t+\epsilon{\rm
d}L^{\alpha}(t),\quad X(0)=0\in\mathbb{R}^{1},$ (4)

In this example, a nonlinear drift term is $f(x,\theta)=-\theta x^{3}+x$ with
uncertain parameter $\theta$. We start with training data from numerical
simulations of Eq.(4). The stochastic trajectory can be regarded as a heavy-
tailed time series with $6\times 10^{6}$ data points with the parameters
$\alpha=1.7$, $\theta=1$ and $\epsilon=0.3$ shown in the Fig. 1(a). The
empirical probability density function $p_{d}$ for observations (Fig. 1(a))
could be determined by the normal kernel method. In the simulation, we use the
MATLAB function _ksensity_ to evaluate the $p_{d}(x)$ for $x\in(-3,3)$ as
shown in Fig. 1(b). A selected bandwidth is $h=1.8s/n^{1/5}$, where $n=6\times
10^{6}$ is the number of observed data points and $s$ is the standard
deviation of the data set.

Now we provide the details of the computation to infer parameters from the
observations. The value of the parameters are estimated by minimizing the
Hellinger distance over the parameter space $\Theta$. We shall first estimate
single parameter assuming that the values of the other parameters are known,
and then estimate multiple parameters

###  3.1 Estimation for single parameter

We want to find out an estimation of $\alpha$ by achieving a numerical
optimization of the objective function $G$ of the Hellinger distance. We
consider the parameter set $\lambda=\\{\alpha\\}$ is in the parameter space
$\Theta=(0,2)\subset\mathbb{R}^{1}$, and assume that the other parameters are
known, i.e. $\epsilon=0.3$ and $\theta=1$. Then the objective function $G$ of
the Hellinger distance is

$G(\alpha)=\frac{1}{2}\int_{\mathbb{R}^{1}}\left(\sqrt{p(x,\alpha)}-\sqrt{p_{d}(x)}\right)^{2}dx.$

The probability density $p(x,\alpha)$ is solved by the nonlocal differential
equation (2) given $\alpha\in(0,2)$ for $x\in D=(-3,3)$ at $t=50$. In the
numerical simulations, the probability profile of its initial position is
Gaussian $p(x,0)=\sqrt{\frac{40}{\pi}}e^{-40x^{2}}$. We have chosen the
spatial resolution $h=0.003$ and the time step size $\Delta t=0.5h^{2}$.

In Fig. 2(a), we employ a discretization step of $\Delta\alpha=0.035$ and use
55 grid points for $\alpha\in(0,2)$. Then the estimation of
$\hat{\alpha}=1.738$ is obtained with the minimum value of the Hellinger
distance $G(\hat{\alpha})=0.001399$ over the parameter space $\Theta=(0,2)$.
Furthermore, we restrict $\alpha$ on a small region $[1.5,1.9]$ for accurately
estimation. The result illustrates that the estimated values of
$\hat{\alpha}\in[1.684,1.764]$ contains the true value of $\alpha=1.7$ with
the Hellinger distance $G(\hat{\alpha})<0.00145$ (inset figure in Fig. 2(a)) .
It means that we could find an optimal interval for the estimated parameter
$\alpha$. As an illustration, we show the results of the probability density
$p(x,\alpha)$ of SDE (4) with estimated value $\hat{\alpha}=1.738$ (dashed)
and the empirical density $p_{d}$ from the observed data (dotted) in Fig.
1(b). We can see that the estimated probability density presents the goodness-
of-fit to the empirical one.

(a)

(b)

Figure 3: (a) The Hellinger distance estimation of $\epsilon$. The estimated
value $\hat{\epsilon}=0.3268$ corresponds to the minimum objective function
$G(\hat{\epsilon})=0.000951$ (red $\star$). (b) The estimated density function
$p(x)$ with $\hat{\epsilon}=0.3268$, $\theta=1$ and $\alpha=1.7$ compared with
the empirical density $p_{d}$ in Fig.1(b).

Similarity, we explore the dependence of the objective function $G$ on the
value of $\epsilon$ keeping the other parameters fixed. Fig. 3(a) shows the
minimized the $G(\hat{\epsilon})=0.000951$ of the Hellinger distance that
corresponds to the estimated value of the Lévy noise intensity
$\hat{\epsilon}=0.3268$. Meanwhile, we could get the optimal estimation
interval $\hat{\epsilon}\in[0.292,0.348]$ for the Hellinger distance
$G(\hat{\epsilon})<0.0019$. This domain includes the true value of
$\epsilon=0.3$. In Figure 3(b), the empirical density $p_{d}$ is well fitted
by the probability density $p(x,\epsilon)$ with $\hat{\epsilon}=0.3268$.

We have inferred the parameters $\alpha$ and $\epsilon$ by considering the
Hellinger distance, respectively. Next, we would like to compare the Hellinger
distance with other commonly used distances, such as the $L^{2}$ norm and the
absolute approximation error, to quantify the similarity between two
probability distributions. The objective function of the $L^{2}$ norm is
defined as

$G(\lambda)=\frac{\|p(\lambda,x)-p_{d}(x)\|^{2}_{2}}{\|p_{d}(x)\|^{2}_{2}},$

The estimation of the uncertain parameters set $\lambda$ could be achieved by
minimizing $G$, i.e., $\hat{\lambda}=\arg\min G(\lambda)$ for
$\lambda\in\Theta$. We could also give an expression of the objective function
for the absolute error

$G(\lambda)=|p(\lambda,x)-p_{d}(x)|,$

while the estimation of the parameters are taken to be maximum of the absolute
error distance, i.e., $\hat{\lambda}=\arg\max G(\lambda)$ for
$\lambda\in\Theta$. The reason is that the minimum of the absolute error at a
point does not guarantee the similarity of the two probability densities.

Next, let us examine the effect of three kinds of distances on the estimation
of $\alpha$ and $\epsilon$, respectively. We keep the other parameters and the
divided subintervals the same as those in Figs. 2 and 3. In Table 1, the
results on these distance show that the Hellinger distance gives a better
estimation for $\alpha$ than the other one. In contrast, all three show a good
fit to the true $\epsilon$. In this example, the Hellinger distance is the
most effective method to estimate parameters.

(a)

(b)

Figure 4: (a) Estimation of $\alpha$ and $\theta$ by optimization of the Hellinger distance. The estimated values are $\hat{\alpha}=1.707$ and $\hat{\theta}=0.973684$ with the minimum $G(\hat{\alpha},\hat{\theta})=0.0014$ (red $\star$). (b) The optimal estimation domain of $\hat{\alpha}$ and $\hat{\theta}$ is $[1.6867,1.78]\times[0.98,1]$. Distance | Estimated $\hat{\alpha}(True1.7)$ | Estimated $\hat{\epsilon}(True0.3)$ | $G_{\alpha}$ | $G_{\epsilon}$
---|---|---|---|---
Hellinger distance | 1.738 | 0.3268 | 0.0014 | 0.0009
$L^{2}$ norm | 1.8460 | 0.3070 | 0.0028 | 0.0032
Absolute error | 1.9180 | 0.3070 | 0.0427 | 0.062
Table 1: Compared with different distance: Hellinger distance, $L^{2}$ norm
and absolute error.

###  3.2 Estimation for multiple parameters

The above example has verified that our method is feasible for estimating a
single parameter by minimizing the Hellinger distance. Next, we will apply
this approach to estimate the multiple unknown parameters. First, we simplify
our model by assuming that one parameter $\epsilon=0.3$ is known and then
estimate $\alpha$ and $\theta$, while keeping the other factors the same as in
the section 3.1. The objective function $G$ is given by Hellinger distance

$G(\alpha,\theta)=\frac{1}{2}\int_{\mathbb{R}^{1}}\left(\sqrt{p(x,\alpha,\theta)}-\sqrt{p_{d}(x)}\right)^{2}dx.$

The probability density $p(x,\alpha,\theta)$ is a solution of the nonlocal
differential equation (2) given values of $\alpha\in(0,2)$ and
$\theta\in[0.5,1.5]$ at $t=50$. Fig. 4(a) shows that the objective function
$G$ changes with the values of $\alpha$ and $\theta$ in the parameter space
$\Theta=(0,2)\times[0.5,1.5]\subset\mathbb{R}^{2}$. The minimum value of
$G(\hat{\alpha},\hat{\theta})=0.0014$ is identified with $\hat{\alpha}=1.707$
and $\hat{\theta}=0.9828$. In the same manner, the optimal estimation domain
of $\hat{\alpha}$ and $\hat{\theta}$ is $[1.6867,1.78]\times[0.98,1]$ (orange
rectangular frame) by further restricting the range of parameters as shown in
Fig. 4(b).

(a)

(b)

Figure 5: (a) Estimation of $\epsilon$ and $\theta$. The estimated values are
$\hat{\epsilon}=0.3083$ and $\hat{\theta}=0.973684$ with the minimum
$G(\hat{\epsilon},\hat{\theta})=0.0012$ (red $\star$). (b) The optimal
estimation domain is $[0.3,0.3267]\times[0.9,1.02]$.

(a)

(b)

Figure 6: (a) Estimation of $\alpha$ and $\epsilon$. The estimated values are $\hat{\alpha}=1.707$ and $\hat{\epsilon}=0.304$ with the minimum $G(\hat{\alpha},\hat{\epsilon})=0.0011$ (red $\star$). (b) The optimal estimation domain is $[1.6733,1.78]\times[0.3133,0.3267]$. Parameter | True value | Estimated $\hat{\lambda}_{1}=\\{\hat{\epsilon},\hat{\theta}\\}$ | Estimated $\hat{\lambda}_{2}=\\{\hat{\alpha},\hat{\epsilon}\\}$
---|---|---|---
$\theta$ | 1.0 | 1.0185 | –
$\alpha$ | 1.7 | – | 1.707
$\epsilon$ | 0.3 | 0.3303 | 0.307
$G$ | 0 | 0.0012 | 0.0011
Table 2: Estimation of $\lambda_{1}=\\{\epsilon,\theta\\}$ and
$\lambda_{2}=\\{\alpha,\epsilon\\}$.

Second, we take into account the estimation of the other two combinations of
all three parameters, $\lambda_{1}=\\{\epsilon,\theta\\}$ and
$\lambda_{2}=\\{\alpha,\epsilon\\}$ corresponding to the parameter spaces
$\Theta_{1}=(0,1]\times[0.5,1.5]$, $\Theta_{2}=(0,2)\times(0,1]$,
respectively. The estimated results are found by the minimized the Hellinger
distance as shown in Table 2. Meanwhile, we can also determine the optimal
domains of the estimated parameters sets $\hat{\lambda}_{1}$ and
$\hat{\lambda}_{2}$ as shown in Figs. 5(b) and 6(b). The result shows that the
Hellinger distance $G(\hat{\alpha},\hat{\epsilon})\leqslant 0.001$ if the
estimated parameters set $\hat{\lambda}_{1}$ belongs to the domain
$[1.6733,1.78]\times[0.3,0.3267]$. Meanwhile, the estimation domain of
$\hat{\lambda}_{2}$ is $[0.3,0.34]\times[0.98,1]$ if the Hellinger distance
$G(\hat{\epsilon},\hat{\theta})\leqslant 0.0015$.

Finally, we seek all these parameters $\lambda=\\{\alpha,\epsilon,\theta\\}$
such that the Hellinger distance reaches the minimum value in the parameter
spaces $\Theta=(0,1]\times(0,2]\times(0,2)\subset\mathbb{R}^{3}$.

$G(\theta,\epsilon,\alpha)=\frac{1}{2}\int_{\mathbb{R}^{1}}\left(\sqrt{p(x,\theta,\epsilon,\alpha)}-\sqrt{p_{d}(x)}\right)^{2}dx.$

The values of the model parameters $\hat{\epsilon}=0.307$, $\hat{\theta}=1.05$
and $\hat{\alpha}=1.7380$ are achieved by minimizing the Hellinger distance
$G(\hat{\alpha},\hat{\theta},\hat{\epsilon})=0.0027$ over the parameter space
$\Theta\subset\mathbb{R}^{3}$. Estimated results $\hat{\epsilon}=0.307$,
$\hat{\theta}=1.05$ and $\hat{\alpha}=1.7380$ defines a slice plane in the
$\epsilon$-axis, $\theta$-axis, or $\alpha$-axis direction as shown in Fig. 7.

Figure 7:  The Hellinger distance estimation of $\alpha$, $\epsilon$ and
$\theta$. The estimated $\hat{\epsilon}=0.307$, $\hat{\theta}=1.05$ and
$\hat{\alpha}=1.7380$ for the minimum
$G(\hat{\alpha},\hat{\theta},\hat{\epsilon})=0.0027$.

##  4 Conclusion

In summary, we consider a non-Gaussian dynamical system containing uncertain
parameters. An approach of parameter estimation is proposed by numerical
optimization of the Hellinger distance between two probability distributions.
The one probability density $p(x)$ is a solution of the nonlocal Fokker-Planck
equation at time $t$ for a stochastic dynamical system $X(t)$ driven by an
$\alpha$-stable Lévy process. The other one is the empirical probability
density $p_{d}$ from observations data of discrete version of the process. The
approach is used to find all out the estimation of single parameter and
multiple parameters, by a numerical optimization of the Hellinger distance
over the parameters space. The results of an example verified that this method
is feasible for estimating non-Gaussian parameters $\alpha$, $\epsilon$ and
other system parameter by the Hellinger distance. Compared with $L^{2}$ norm
and absolute error distances, the Hellinger distance is the most effective
method to estimate parameters in this example. Meanwhile, we could find an
optimal interval for the estimated parameters.

This approach can be used to establish parameter estimations for a data-driven
dynamical system, the observations data with jumps and heavy-tailed
distribution. A very important future work will be a model study of the abrupt
climate changes in the Dansgaard-Oeschger events with non-Gaussian
distribution. The approach would be applied to estimate the system parameters
and non-Gaussian parameters in this model.

## References

  * [1] J. Bishwal.  Parameter estimation in stochastic differential equations.  Springer-Verlag Berlin Heidelberg, 2008.
  * [2] S. Iacus.  Simulation and inference for stochastic differential equations: With R examples.  Springer-Verlag New York, 2008.
  * [3] Z. Botev, J. Grotowski, and D. Kroese.  Kernel density estimation via diffusion.  Annals of Statistics, 38(5):2916–2957, 2010.
  * [4] O. Papaspiliopoulos, Y. Pokern, G. Roberts, and A. Stuart.  Nonparametric estimation of diffusions: a differential equations approach.  Biometrika, 99(3):511–531, 2012.
  * [5] M. Beaumont, W. Zhang, and D. Balding.  Approximate bayesian computation in population genetics.  Genetics, 162(4):2025–2035, 2002.
  * [6] C. Honisch and R. Friedrich.  Estimation of kramers-moyal coefficients at low sampling rates.  Physical Review E, 84(2):066701, 2011.
  * [7] S. Siegert, R. Friedrich, and J. Peinke.  Analysis of data sets of stochastic systems.  Physics Letters A, 243(5-6):275–280, 1998.
  * [8] P. Batz, A. Ruttor, and M. Opper.  Variational estimation of the drift for stochastic differential equations from the empirical density.  Journal of Statistical Mechanics-Theory and Experiment, 2016(8):083404, 2016.
  * [9] P. Ditlevsen.  Observation of $\alpha$-stable noise induced millennial climate changes from an ice-core record.  Geophysical Research Letters, 26(10):1441–1444, 1999.
  * [10] L. Cai, N. Friedman, and X. Xie.  Stochastic protein expression in individual cells at the single molecule level.  Nature, 440(7082):358–362, 2006.
  * [11] N. Humphries, N. Queiroz, J. Dyer, N. Pade, M. Musyl, K. Schaefer, D. Fuller, J. Brunnschweiler, T. Doyle, J. Houghton, G. Hays, C. Jones, L. Noble, V. Wearmouth, E. Southall, and D. Sims.  Environmental context explains lévy and brownian movement patterns of marine predators.  Nature, 465(7301):1066–1069, 2010.
  * [12] R. Mantegna.  Lévy walks and enhanced diffusion in milan stock exchange.  Physica A: Statistical Mechanics and its Applications, 179(2):232–242, 1991.
  * [13] P. Montez, G. Thompson, and C. Kello.  The role of semantic clustering in optimal memory foraging.  Cognitive science, 39(8):1925–1939, 2015.
  * [14] A. Voss, V. Lerche, U. Mertens, and J. Voss.  Sequential sampling models with variable boundaries and non-normal noise: A comparison of six models.  Psychonomic bulletin and review, 26(3):813–832, 2019.
  * [15] Q. Yu, G. Shen, and M. Cao.  Parameter estimation for Ornstein-Uhlenbeck processes of the second kind driven by $\alpha$-stable lévy motions.  Communications in Statistics-Theory Methods, 46(21):10864–10878, 2017.
  * [16] Y. Hu and H. Long.  Least squares estimator for Ornstein-Uhlenbeck processes driven by $\alpha$-stable motions.  Stochastic Processes and Their Applications, 119(8):2465–2480, 2009\.
  * [17] H. Long and L. Qian.  Nadaraya-watson estimator for stochastic processes driven by stable Lévy motions.  Electronic Journal of Statistics, 7:1387–1418, 2013.
  * [18] P. Gopikrishnan, V. Plerou, L. Amaral, M. Meyer, and H. Stanley.  Scaling of the distribution of fluctuations of financial market indices.  Physical Review E, 60(5):5305–5316, 1999.
  * [19] B. Hill.  A simple general approach to inference about the tail of a distribution.  Annals of Statistics, 3(5):1163–1174, 1975.
  * [20] J. Press.  Estimation in univariate and multivariate stable distributions.  Journal of the American Statistical Association, 67(340):842–846, 1972.
  * [21] Y. Cheng, Y. Hu, and H. Long.  Generalized moment estimators for $\alpha$-stable Ornstein-Uhlenbeck motions from discrete observations.  Statistical Inference for Stochastic Processes, 23(1):53–81, 2020\.
  * [22] L. Mello.  Fat-tailed distributions and Lévy processes.  arXiv preprint arXiv:0808.2393, 2008.
  * [23] Z. Chen and X. Chen.  Maximum likelihood estimation for symmetric $\alpha$-stable Ornstein-Uhlenbeck processes.  Stochastics and Dynamics, 2150018, 2020.
  * [24] T. Gao and J. Duan.  Quantifying model uncertainty in dynamical systems driven by non-Gaussian Lévy stable noise with observations on mean exit time or escape probability.  Communications in Nonlinear Science and Numerical Simulation, 39:1–6, 2016.
  * [25] T. Gao, J. Duan, and X. Li.  Fokker-planck equations for stochastic dynamical systems with symmetric Lévy motions.  Applied Mathematics and Computation, 278:1–20, 2016.
  * [26] X. Sun, X. Li, and Y. Zheng.  Governing equations for probability densities of Marcus stochastic differential equations with Lévy noise.  Stochastics and Dynamics, 17(05):1750033, 2017.
  * [27] D. Applebaum.  Lévy processes and stochastic calculus.  Cambridge University Press, 2009.
  * [28] K. Sato.  Lévy processes and infinitely divisible distributions.  Cambridge University Press, 1999.
  * [29] X. Zhang.  Densities for sdes driven by degenerate $\alpha$-stable processes.  Annals of Probability, 42(5):1885–1910, 2014.
  * [30] Y. Song and X. Zhang.  Regularity of density for SDEs driven by degenerate Lévy noises.  Electronic Journal of Probability, 20:1–27, 2015.
  * [31] Z. Chen, R. Song, and X. Zhang.  Stochastic flows for Lévy processes with Hölder drifts.  Revista Matematica Iberoamericana, 34(4):1755–1788, 2018.
  * [32] H. Kaper and H. Engler.  Mathematics and climate.  SIAM, 2013.
  * [33] P. Cessi.  A simple box model of stochastically forced thermohaline flow.  Journal of physical oceanography, 24(9):1911–1920, 1994.
  * [34] R. Beran.  Minimum Hellinger distance estimates for parametric models.  Annals of Statistics, 5(3):445–463, 1977.
  * [35] J. Duan.  An introduction to stochastic dynamics.  Cambridge University Press, 2015.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Author Contributions

Y. Zheng designed the research. Y. Zheng and F. Yang performed computations
and wrote the first draft of the manuscript. J. Kurths and J. Duan analysed
the results and concepts development. All authors conducted research
discussions and reviewed the manuscript.

## Acknowledgements

We would like to thank Xiaoli Chen, Xiujun Cheng and Yang Liu for discussions
about computation. This work is supported by the National Natural Science
Foundation of China (grants No.11801192 and No. 11771449).