# A Surrogate Loss Function for Optimization of $F_{\beta}$ Score in Binary
Classification with Imbalanced Data

Namgil Lee<EMAIL_ADDRESS>Department of Information Statistics,
Kangwon National University, Chuncheon, Gangwon 24341, Republic of Korea
Bionsight, Inc., Gangwondaehak-gil 1, Chuncheon, Gangwon 24341, Republic of
Korea  Heejung Yang<EMAIL_ADDRESS>Bionsight, Inc., Gangwondaehak-gil
1, Chuncheon, Gangwon 24341, Republic of Korea  Department of Pharmacy,
Kangwon National University, Gangwondaehak-gil 1, Chuncheon, Gangwon 24341,
Republic of Korea  Hojin Yoo<EMAIL_ADDRESS>Bionsight, Inc.,
Gangwondaehak-gil 1, Chuncheon, Gangwon 24341, Republic of Korea

###### Abstract

The $F_{\beta}$ score is a commonly used measure of classification
performance, which plays crucial roles in classification tasks with imbalanced
data sets. However, the $F_{\beta}$ score cannot be used as a loss function by
gradient-based learning algorithms for optimizing neural network parameters
due to its non-differentiability. On the other hand, commonly used loss
functions such as the binary cross-entropy (BCE) loss are not directly related
to performance measures such as the $F_{\beta}$ score, so that neural networks
optimized by using the loss functions may not yield optimal performance
measures. In this study, we investigate a relationship between classification
performance measures and loss functions in terms of the gradients with respect
to the model parameters. Then, we propose a differentiable surrogate loss
function for the optimization of the $F_{\beta}$ score. We show that the
gradient paths of the proposed surrogate $F_{\beta}$ loss function approximate
the gradient paths of the large sample limit of the $F_{\beta}$ score. Through
numerical experiments using ResNets and benchmark image data sets, it is
demonstrated that the proposed surrogate $F_{\beta}$ loss function is
effective for optimizing $F_{\beta}$ scores under class imbalances in binary
classification tasks compared with other loss functions.

Key words: binary classification, F1 score, imbalanced data.

##  1 Introduction

An imbalanced data problem refers to the case that the distribution of class
labels in the training data is not uniform but skewed (He and Garcia, 2009).
Such problem arises widely in real life applications of machine learning, such
as computer vision (Xiao et al., 2010; Kubat et al., 1998; Beijbom et al.,
2012), biology (Yu et al., 2012), medicine (Grzymala-Busse et al., 2004; Mac
Namee et al., 2002), and fraud detection (Chan and Stolfo, 1998). The
imbalanced data problem can cause machine learning models, such as the
convolutional neural networks, to perform poorly on test data, especially for
minority classes (Japkowicz and Stephen, 2002; He and Garcia, 2009; G. Van
Horn, 2017; Buda et al., 2018; Johnson and Khoshgoftaar, 2019).

In general, machine learning methods for addressing the imbalanced data
problem can be grouped into two categories: data-level methods and algorithm-
level methods (Krawczyk, 2016). In data-level methods, the number of samples
for each class is directly adjusted by over-sampling the minority classes,
under-sampling the majority classes, or generating synthetic data for the
minority classes. However, data-level methods often suffer from several
problems. That is, over-sampling and synthetic data generation can introduce a
large amount of duplicated samples, which can incur an overfitting problem and
increased computational costs. And under-sampling can discard important
samples, so that the test set performance of the models trained on partial
data can be impaired. Due to such potential problems, we focus on algorithm-
level methods in this study. Algorithm-level methods are studied mostly on the
cost-sensitive re-weighting methods, where loss functions are utilized to
assign higher costs to the misclassified minority classes than to the
misclassified majority classes. Gradient-based learning algorithms, such as
the stochastic gradient descent (SGD) method, search for the network
parameters by minimizing the loss functions.

On the other hand, the performance of a classifier is evaluated by performance
measures such as the accuracy. In the context of imbalanced data, the accuracy
has a severe limitation that it can lead to incorrect conclusions which favors
the majority classes over the minority classes. There are a number of
performance measures that address this issue (Johnson and Khoshgoftaar, 2019).
The $F_{\beta}$ score is one of the most widely used performance measures in
imbalanced binary classification tasks. Note that the precision is defined as
the proportion of positive predictions which are actually correct, and the
recall is defined as the proportion of actual positives which are predicted
correctly. The $F_{\beta}$ score is defined as a harmonic mean of precision
and recall by

$F_{\beta}=\frac{(1+\beta^{2})\cdot\text{Precision}\cdot\text{Recall}}{(\beta^{2}\cdot\text{Precision})+\text{Recall}}=\left(\lambda_{\beta}\cdot\text{Recall}^{-1}+(1-\lambda_{\beta})\cdot\text{Precision}^{-1}\right)^{-1},$
(1)

where the parameter $\beta>0$ adjusts the relative weight,
$\lambda_{\beta}=\beta^{2}/(1+\beta^{2})$, between precision and recall.

However, most of the performance measures cannot be directly used as a loss
function by gradient-based learning algorithms because they are usually not
differentiable. In most cases, loss functions which are used for optimizing
neural network parameters are not identical to the performance measures, which
results in sub-optimal model parameters in terms of performance measures. For
instance, one of the most popular loss functions for binary classification is
the binary cross-entropy (BCE) loss. It is well known that minimizing the BCE
loss is equivalent to maximizing the likelihood for Bernoulli distribution,
but a relationship between the BCE loss and performance measures has not been
studied sufficiently.

Our main contribution can be summarized as follows: (1) We investigate a
relationship between classification performance measures and loss functions in
terms of the gradients with respect to the model parameters. We derive
explicit conditions that the gradients should satisfy approximately at
critical points. (2) We propose a differentiable surrogate loss function for
optimization of $F_{\beta}$ score. (3) It is demonstrated that the gradient
paths of the proposed loss function are approximately same to the gradient
paths of the $F_{\beta}$ score through numerical experiments. (4) We present a
generalization of the proposed loss function for robustness to label noise.

The rest of this paper is organized as follows. In Section 2, we provide the
literature of related works. In Section 3, the gradient conditions are
analyzed for performance measures and standard loss functions. In Section 4,
we present the proposed surrogate $F_{\beta}$ loss function and its condition
for gradients. In Section 5, numerical experiments are performed using ResNet
models and benchmark image data sets to compare the proposed loss function
with other loss functions. Conclusion and future works are provided in Section
6.

##  2 Related Work

For binary classification of imbalanced data sets, weights are often set
proportional to inverse class frequencies for cost-sensitive learning
strategies (Wang et al., 2017). Cui et al. (2019) derived a novel formula for
the effective number of samples and used it to propose a class-balanced loss
function for cost-sensitive learning strategies. In this study, we analyze
class-balanced loss functions to investigate their relationship with
classification performance measures.

In most cases, a machine learning model which yields the best value of a
performance measure is desired, however, a learning algorithm is developed to
optimize a loss function instead of the performance measure. Two major reasons
that loss functions are preferred to performance measures by learning
algorithms are their differentiability and convexity (Bartlett et al., 2006).
A direct extension of the $F_{\beta}$ score has been suggested for
optimization of artificial neural networks (Pastor-Pellicer et al., 2013), but
it is a heuristically formulated loss function without mathematical
justification. There have been developed several algorithms to optimize
differentiable lower bounds of performance measures (Joachims, 2005; Kar et
al., 2014; Narasimhan et al., 2015), but the performance of deep neural
networks trained with those lower bounds should be further improved (Sanyal et
al., 2018). Eban et al. (2017) proposed scalable optimization algorithms based
on lower bounds for the performance measures, but its gradient path and
performance were not analyzed sufficiently.

In this study, we propose a surrogate loss function for $F_{\beta}$ score, and
the same idea can be easily extended to define surrogate loss functions for
other performance measures. We present mathematical conditions for the
gradients of performance measures and the proposed loss function at critical
points. The presented gradient conditions guarantee that the gradient paths of
the proposed loss approximate those of the $F_{\beta}$ score asymptotically.
Moreover, the proposed loss function can be used for scalable optimization by
standard gradient-based optimization methods.

##  3 Gradient Conditions of Performance Measures and Standard Loss Functions

###  3.1 Preliminaries

We consider a $C$-class classification problem. Let ${\mathcal{X}}$ denote the
feature space and ${\mathcal{Y}}=\\{0,1,\ldots,C-1\\}$ denote the label space.
A sample is a set $\\{({\mathbf{x}}_{i},y_{i})\\}_{i=1}^{n}$ with each
$({\mathbf{x}}_{i},y_{i})\in{\mathcal{X}}\times{\mathcal{Y}}$. A classifier is
a vector-valued function
${\mathbf{f}}:{\mathcal{X}}\rightarrow{\mathbb{R}}^{C}$, where each component
function is denoted as $f_{y}$, $y\in{\mathcal{Y}}$. We assume that a
classifier is a neural network with the softmax function as the output layer.
A loss function ${\mathcal{L}}$ is a real valued function,
${\mathcal{L}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)$, where
${\boldsymbol{\theta}}$ is a set of parameters. In the special case of binary
classification problem, $C=2$, a classifier ${\mathbf{f}}=(f_{0},f_{1})$ with
the softmax output layer can be represented by a single real-valued function
$f$ defined by $f\equiv f_{1}$ and $f_{0}=1-f$. In this case, a loss function
can be denoted by ${\mathcal{L}}(f({\mathbf{x}};{\boldsymbol{\theta}}),y)$.

A class-balanced loss function (Cui et al., 2019) is defined by a loss
function multiplied by a class-dependent weight, $w_{y}$, as

${\mathcal{L}}_{\text{CB}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)=w_{y}{\mathcal{L}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y).$
(2)

In this study, we suppose that the class-dependent weight is set to the
inverse class frequency, $w_{y}=n/n_{y}$ for $y\in{\mathcal{Y}}$, where
$n_{y}$ is the number of training samples with label $y$ and
$\sum_{y\in{\mathcal{Y}}}n_{y}=n$.

For classification tasks, most of the performance measures can be represented
in terms of the entries of the confusion matrix. In the case of binary
classification, the predicted class label $\hat{y}_{i}$ is defined by
$\hat{y}_{i}=1$ if $f({\mathbf{x}}_{i};{\boldsymbol{\theta}})\geq 0.5$ and
$\hat{y}_{i}=0$ otherwise. The entries of the confusion matrix can be written
as

$\begin{split}TP({\boldsymbol{\theta}})=\sum_{i=1}^{n}y_{i}\hat{y}_{i},\quad
TN({\boldsymbol{\theta}})=\sum_{i=1}^{n}(1-y_{i})(1-\hat{y}_{i}),\\\
FP({\boldsymbol{\theta}})=\sum_{i=1}^{n}(1-y_{i})\hat{y}_{i},\quad
FN({\boldsymbol{\theta}})=\sum_{i=1}^{n}y_{i}(1-\hat{y}_{i}).\end{split}$ (3)

Let the sample proportions of actual positives, false negatives, and false
positives be denoted by

$p_{n}=\frac{n_{1}}{n},\quad
q_{0n}=\frac{FN({\boldsymbol{\theta}})}{n_{1}},\quad
q_{1n}=\frac{FP({\boldsymbol{\theta}})}{n-n_{1}}.$ (4)

If we assume that the samples are independent and identically distributed,
then the entries of the confusion matrix have the multinomial distribution
with the total number of trials
$n=TP({\boldsymbol{\theta}})+TN({\boldsymbol{\theta}})+FP({\boldsymbol{\theta}})+FN({\boldsymbol{\theta}})$.
Hence, the expected values of the sample proportions can be defined by the
true proportions as

$p={\mathbb{E}}\left[p_{n}\right]=\text{Pr}(y=1),\quad
q_{0}={\mathbb{E}}\left[q_{0n}\right]=\text{Pr}(\hat{y}=0|y=1),\quad
q_{1}={\mathbb{E}}\left[q_{1n}\right]=\text{Pr}(\hat{y}=1|y=0).$ (5)

Based on the weak law of large numbers, we can derive that the sample
proportions converge to the true proportions, i.e., $p_{n}\rightarrow p$,
$q_{0n}\rightarrow q_{0}$, and $q_{1n}\rightarrow q_{1}$, as the sample sizes
$n$, $n_{1}$, and $n-n_{1}$ tend to infinity.

In addition, for a real-valued function $g$ on ${\mathcal{X}}$ and
$y_{0}\in\\{0,1\\}$, we define the conditional sample mean of
$g({\mathbf{x}})$ given $y=y_{0}$ by
$\widehat{{\mathbb{E}}}\left[g({\mathbf{x}})|y=y_{0}\right]=n_{y_{0}}^{-1}\sum_{i:y_{i}=y_{0}}g({\mathbf{x}}_{i})$,
which can be re-written as

$\widehat{{\mathbb{E}}}\left[g({\mathbf{x}})|y=1\right]=\frac{1}{n_{1}}\sum_{i=1}^{n}y_{i}g({\mathbf{x}}_{i}),\quad\widehat{{\mathbb{E}}}\left[g({\mathbf{x}})|y=0\right]=\frac{1}{n_{0}}\sum_{i=1}^{n}(1-y_{i})g({\mathbf{x}}_{i}).$
(6)

Then, the sample proportions can be expressed as
$q_{0n}=\widehat{{\mathbb{E}}}\left[1-\hat{y}|y=1\right]$ and
$q_{1n}=\widehat{{\mathbb{E}}}\left[\hat{y}|y=0\right]$. Recall that the
predicted class label $\hat{y}$ is determined by the classifier output
$f({\mathbf{x}};{\boldsymbol{\theta}})$, and it satisfies that
$|\hat{y}-f({\mathbf{x}};{\boldsymbol{\theta}})|\leq 0.5$. By replacing
$\hat{y}$ with $f({\mathbf{x}};{\boldsymbol{\theta}})$, we define the smoothed
sample proportions by

$\tilde{q}_{0n}=\widehat{{\mathbb{E}}}\left[1-f({\mathbf{x}};{\boldsymbol{\theta}})|y=1\right],\quad\tilde{q}_{1n}=\widehat{{\mathbb{E}}}\left[f({\mathbf{x}};{\boldsymbol{\theta}})|y=0\right].$
(7)

###  3.2 Gradient Condition of the $F_{\beta}$ Score

Most of the performance measures, such as the accuracy, precision, $F_{\beta}$
score, and Jaccard score, can be represented as a combination of the entries
of the confusion matrix (Koyejo et al., 2014). In this study, we focus on the
$F_{\beta}$ score, which is defined in (1). Note that the precision and the
recall are defined by

$\text{Precision}=\frac{TP({\boldsymbol{\theta}})}{TP({\boldsymbol{\theta}})+FP({\boldsymbol{\theta}})},\quad\text{Recall}=\frac{TP({\boldsymbol{\theta}})}{TP({\boldsymbol{\theta}})+FN({\boldsymbol{\theta}})}.$
(8)

By using the sample proportions in (4), the $F_{\beta}$ score can be expressed
as

$F_{\beta}=\frac{(1+\beta^{2})\cdot
p_{n}(1-q_{0n})}{p_{n}(\beta^{2}+1-q_{0n}-q_{1n})+q_{1n}}.$ (9)

Even though the $q_{0n}$ and $q_{1n}$ are functions of the parameters
${\boldsymbol{\theta}}$, they cannot be differentiated with respect to the
parameters because the entries of the contingency matrix are not smooth
functions of ${\boldsymbol{\theta}}$. On the other hand, the limit of the
$F_{\beta}$ score as $n\rightarrow\infty$ can be expressed in terms of the
true proportions as

$\lim_{n\rightarrow\infty}F_{\beta}=\tilde{F}_{\beta}=\frac{(1+\beta^{2})\cdot
p(1-q_{0})}{p(\beta^{2}+1-q_{0}-q_{1})+q_{1}}.$ (10)

Since the true proportions, $q_{0}$ and $q_{1}$, are smooth functions of
${\boldsymbol{\theta}}$, the limit of $F_{\beta}$ is a smooth function.
Considering a gradient-based learning algorithms, taking the gradient of the
logarithm of the limit function yields

$\frac{\partial}{\partial{\boldsymbol{\theta}}}\log\tilde{F}_{\beta}=z^{-1}\cdot\left(-q_{0}^{\prime}(\beta^{2}\cdot
p-pq_{1}+q_{1})+q_{1}^{\prime}(p-pq_{0}-1+q_{0})\right),$ (11)

where $q_{0}^{\prime}=\nabla_{\boldsymbol{\theta}}q_{0}$ and
$q_{1}^{\prime}=\nabla_{\boldsymbol{\theta}}q_{1}$ are the gradients of the
true proportions, and $z=(1-q_{0})\cdot(p(\beta^{2}+1-q_{0}-q_{1})+q_{1})$.
From the condition at the critical point,
$\partial\log\tilde{F}_{\beta}/\partial{\boldsymbol{\theta}}={\bf 0}$, we
obtain the following condition of the gradients:

$-\frac{q_{0}^{\prime}}{1-q_{0}}=\frac{q_{1}^{\prime}}{\beta^{2}\cdot
p/(1-p)+q_{1}}.$ (12)

###  3.3 Standard Loss Functions

We consider the following two standard loss functions for binary
classification:

  * •

The binary cross-entropy (BCE) loss can be expressed as

${\mathcal{L}}_{\text{BCE}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)=-\log
f_{y}({\mathbf{x}};{\boldsymbol{\theta}})=-y\log
f({\mathbf{x}};{\boldsymbol{\theta}})-(1-y)\log(1-f({\mathbf{x}};{\boldsymbol{\theta}})).$
(13)

  * •

The mean absolute error (MAE) loss, also called as the $l_{1}$ loss, can be
written as

${\mathcal{L}}_{\text{MAE}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)=\|{\mathbf{e}}_{y}-{\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}})\|_{1}=2-2f_{y}({\mathbf{x}};{\boldsymbol{\theta}})=2y(1-f({\mathbf{x}};{\boldsymbol{\theta}}))+2(1-y)f({\mathbf{x}};{\boldsymbol{\theta}}),$
(14)

where ${\mathbf{e}}_{y}$ is the vector whose $y$th element is 1 and the other
elements are zeros.

###  3.4 Gradient Condition of the BCE Loss

With the class weight $w_{y}=n/n_{y}$, the sample mean of the class-balanced
BCE losses can be expressed as

$\begin{split}\widehat{{\mathcal{R}}}^{w}_{\text{BCE}}({\boldsymbol{\theta}})&=\frac{1}{n}\sum_{i=1}^{n}w_{y_{i}}{\mathcal{L}}_{\text{BCE}}({\mathbf{f}}({\mathbf{x}}_{i};{\boldsymbol{\theta}}),y_{i})\\\
&=-\frac{1}{n}\sum_{i=1}^{n}\left(\frac{ny_{i}}{n_{1}}\log
f({\mathbf{x}}_{i};{\boldsymbol{\theta}})+\frac{n(1-y_{i})}{n-n_{1}}\log(1-f({\mathbf{x}}_{i};{\boldsymbol{\theta}}))\right)\\\
&=\widehat{{\mathbb{E}}}\left[\left.-\log
f({\mathbf{x}};{\boldsymbol{\theta}})\right|y=1\right]+\widehat{{\mathbb{E}}}\left[\left.-\log(1-f({\mathbf{x}};{\boldsymbol{\theta}}))\right|y=0\right]\\\
&\geq-\log\widehat{{\mathbb{E}}}\left[\left.f({\mathbf{x}};{\boldsymbol{\theta}})\right|y=1\right]-\log\widehat{{\mathbb{E}}}\left[\left.(1-f({\mathbf{x}};{\boldsymbol{\theta}}))\right|y=0\right],\end{split}$
(15)

where the inequality at the last line follows from the Jensen’s inequality.
After all, from (7) and (15), the sample mean of the class-balanced BCE losses
can be approximately represented by

$\widehat{{\mathcal{R}}}^{w}_{\text{BCE}}({\boldsymbol{\theta}})\approx-\log(1-\tilde{q}_{0n})-\log(1-\tilde{q}_{1n}),$
(16)

where the approximation error is determined by the Jensen’s inequality in
(15). Note that, unlike the sample proportions $q_{0n}$ and $q_{1n}$, the
smoothed sample proportions $\tilde{q}_{0n}$ and $\tilde{q}_{1n}$ are
differentiable with respect to the parameters ${\boldsymbol{\theta}}$. From
$\partial\widehat{{\mathcal{R}}}^{w}_{\text{BCE}}({\boldsymbol{\theta}})/\partial{\boldsymbol{\theta}}={\bf
0}$, we can obtain the following condition for the gradients:

$-\frac{\tilde{q}_{0n}^{\prime}}{1-\tilde{q}_{0n}}=\frac{\tilde{q}_{1n}^{\prime}}{1-\tilde{q}_{1n}},$
(17)

where $\tilde{q}_{0n}^{\prime}=\nabla_{\boldsymbol{\theta}}\tilde{q}_{0n}$ and
$\tilde{q}_{1n}^{\prime}=\nabla_{\boldsymbol{\theta}}\tilde{q}_{1n}$.

###  3.5 Gradient Condition of the MAE Loss

The sample mean of the class-balanced MAE losses can be approximately
represented by

$\begin{split}\widehat{{\mathcal{R}}}^{w}_{\text{MAE}}({\boldsymbol{\theta}})&=\frac{1}{n}\sum_{i=1}^{n}w_{y_{i}}{\mathcal{L}}_{\text{MAE}}({\mathbf{f}}({\mathbf{x}}_{i};{\boldsymbol{\theta}}),y_{i})\\\
&=\frac{2}{n}\sum_{i=1}^{n}\left(\frac{ny_{i}}{n_{1}}(1-f({\mathbf{x}}_{i};{\boldsymbol{\theta}}))+\frac{n(1-y_{i})}{n-n_{1}}f({\mathbf{x}}_{i};{\boldsymbol{\theta}})\right)\\\
&=2\widehat{{\mathbb{E}}}\left[\left.1-f({\mathbf{x}};{\boldsymbol{\theta}})\right|y=1\right]+2\widehat{{\mathbb{E}}}\left[\left.f({\mathbf{x}};{\boldsymbol{\theta}})\right|y=0\right]\\\
&=2\tilde{q}_{0n}+2\tilde{q}_{1n}.\end{split}$ (18)

From
$\partial\widehat{{\mathcal{R}}}^{w}_{\text{MAE}}({\boldsymbol{\theta}})/\partial{\boldsymbol{\theta}}={\bf
0}$, we can obtain the following condition for the gradients:

$-\tilde{q}_{0n}^{\prime}=\tilde{q}_{1n}^{\prime}.$ (19)

##  4 Proposed Surrogate Loss for the $F_{\beta}$ Score

Motivated from the gradient conditions of the BCE loss and the MAE loss in
(17) and (19), we propose the surrogate loss function for the $F_{\beta}$
score, which is defined by

${\mathcal{L}}_{F_{\beta}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)=-y\log
f({\mathbf{x}};{\boldsymbol{\theta}})+(1-y)\log\left(\beta^{2}\cdot\frac{p}{1-p}+f({\mathbf{x}};{\boldsymbol{\theta}})\right),$
(20)

where $\beta>0$ controls the balance between the precision and recall as in
the $F_{\beta}$ score, and $0<p<1$ represents the proportion of the positive
samples, $p=n_{1}/n$. Note that the positive part ($y=1$) of the surrogate
loss function is equal to that of the BCE loss, but the negative part ($y=0$)
is different. The reason can be found by analyzing the gradient condition of
the surrogate loss function as follows. The sample mean of the class-balanced
surrogate loss function is

$\begin{split}\widehat{{\mathcal{R}}}^{w}_{F_{\beta}}({\boldsymbol{\theta}})&=\frac{1}{n}\sum_{i=1}^{n}w_{y_{i}}{\mathcal{L}}_{F_{\beta}}({\mathbf{f}}({\mathbf{x}}_{i};{\boldsymbol{\theta}}),y_{i})\\\
&=\frac{1}{n}\sum_{i=1}^{n}\left(-\frac{ny_{i}}{n_{1}}\log
f({\mathbf{x}}_{i};{\boldsymbol{\theta}})+\frac{n(1-y_{i})}{n-n_{1}}\log\left(\beta^{2}\cdot\frac{p}{1-p}+f({\mathbf{x}}_{i};{\boldsymbol{\theta}})\right)\right)\\\
&=\widehat{{\mathbb{E}}}\left[\left.-\log
f({\mathbf{x}};{\boldsymbol{\theta}})\right|y=1\right]+\widehat{{\mathbb{E}}}\left[\left.\log\left(\beta^{2}\cdot\frac{p}{1-p}+f({\mathbf{x}}_{i};{\boldsymbol{\theta}})\right)\right|y=0\right]\\\
&\approx-\log(1-\tilde{q}_{0n})+\log\left(\beta^{2}\cdot\frac{p}{1-p}+\tilde{q}_{1n}\right).\end{split}$
(21)

From
$\partial\widehat{{\mathcal{R}}}^{w}_{F_{\beta}}({\boldsymbol{\theta}})/\partial{\boldsymbol{\theta}}={\bf
0}$, we can derive that

$-\frac{\tilde{q}_{0n}^{\prime}}{1-\tilde{q}_{0n}}=\frac{\tilde{q}_{1n}^{\prime}}{\beta^{2}\cdot
p/(1-p)+\tilde{q}_{1n}}.$ (22)

Note that the gradient property for the $F_{\beta}$ score in (12) is
equivalent to that for the surrogate loss in (22) except that one employs the
true proportions and the other employs the smoothed sample proportions. Hence,
the additional parameters $\beta$ and $p$ for the surrogate loss possess the
same roles as in the $F_{\beta}$ score, and they can be utilized to control
the $F_{\beta}$ score of the trained neural network models.

From the expressions in (16), (18), and (21), we can find that every
expression for the sample mean of the loss functions consists of two terms: a
function of $\tilde{q}_{0n}$ and a function of $\tilde{q}_{1n}$. Figure 1(a)
illustrates the first term, i.e., the function of $\tilde{q}_{0n}$, for each
of the BCE, MAE, and surrogate $F_{\beta}$ loss functions. In specific, the
functions corresponding to the first term can be written as $z=-\log(1-q)$,
$z=q$, and $z=-\log(1-q)$ with $q=\tilde{q}_{0n}$. Likewise, Figure 1(b)
compares the second term, i.e., the function of $\tilde{q}_{10}$, for the BCE,
MAE, and surrogate $F_{\beta}$ loss functions. The three functions
corresponding to the second term can be written as $-\log(1-q)$, $q$, and
$\log(1+q)$ with $q=\tilde{q}_{1n}$, respectively, when $\beta=1$ and $p=0.5$.
Note that the first term of the surrogate $F_{\beta}$ loss is equal to the
first term of the BCE loss, but the second term of the surrogate $F_{\beta}$
loss is quite unique and different from those of the other losses.

|
---|---
(a) | (b)
Figure 1:  Illustration of (a) the functions of $\tilde{q}_{0n}$ and (b) the
functions of $\tilde{q}_{1n}$ comprising of the the expressions in (16), (18),
and (21) for the approximate sample mean of the loss functions.

##  5 Numerical Experiments

In this section, we conduct numerical experiments for evaluating the gradient
paths and the performances of the proposed surrogate $F_{\beta}$ loss
function. We compare several standard loss functions listed as follows:

  1. (a)

The macro soft $F_{\beta}$ loss is the direct extension of the $F_{\beta}$
score (Pastor-Pellicer et al., 2013), which can be written as

${\mathcal{L}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)=1-\tilde{F}_{\beta},$
(23)

where $\tilde{F}_{\beta}$ is the smoothed $F_{\beta}$ score defined by

$\tilde{F}_{\beta}=\frac{(1+\beta^{2})\cdot
p_{n}(1-\tilde{q}_{0n})}{p_{n}(\beta^{2}+1-\tilde{q}_{0n}-\tilde{q}_{1n})+\tilde{q}_{1n}},$
(24)

and the $\tilde{q}_{0n}$ and $\tilde{q}_{1n}$ are the smoothed sample
proportions defined in (7).

  2. (b)

The BCE loss in (13).

  3. (c)

The MAE loss in (14).

For the experiments, ResNet models (He et al., 2016) were trained and tested
on the CIFAR-10111https://www.cs.toronto.edu/ kriz/cifar.html and Fashion-
MNIST222https://github.com/zalandoresearch/fashion-mnist data sets. The
experiments were implemented based on PyTorch (Paszke et al., 2019) running on
Ubuntu 18 with Intel(R) Xeon CPU of 2.30GHz and GeForce RTX 2080 Ti. See the
following section for the detailed experimental settings.

###  5.1 Experimental Settings

ResNet-18 was used for binary classification of T-shirt/top (the first class)
and the other classes in the Fashion-MNIST data set. ResNet-34 was used for
binary classification of airplane images (the first class) and the other
classes in the CIFAR-10 data set. In each of the data sets, the first class
accounts for 10 percent of the whole data set. ResNet models employed ReLU in
hidden layers and softmax layer at output layer. For training ResNet models,
we used a mini-batch size of 100 and stochastic gradient descent (SGD) with
0.9 momentum, a weight decay of $10^{-4}$, and learning rate of 0.01. Each
experiment consisted of 120 epochs in total. We repeated the experiments five
times independently. In each experiment, network weights were randomly
initialized and 10 percent of the training set was randomly separated into
validation set and the rest were remained as training set and used for
training ResNet models. For data augmentation and preprocessing, each image
was transformed by $32\times 32$ random crop after padding 4 pixels on each
side, random horizontal flip, and normalization by per-pixel mean subtraction.
Since ResNet takes three channels as input but the images in the Fashion-MNIST
are gray scale single channel images, each image was repeated into three
channels.

###  5.2 Results on Fashion-MNIST Data Set

For the Fashion-MNIST data set, the ResNet-18 models were trained by the
proposed surrogate $F_{\beta}$ losses and the other standard loss functions.
Figure 2 illustrates the convergence of the ResNet-18 models trained by using
the proposed loss functions and the other loss functions, where all the loss
functions were class-balanced by the inverse class frequency. We present three
performance measures, which are the F1 score on the training set, F1 score on
the validation set, and accuracy on the test set. Each performance values are
the median of the five repeated experiments. We removed the BCE loss function
because the models trained by the BCE loss function failed to achieve any
improvement in performance, maybe because the class-balanced BCE loss requires
a larger number of epochs greater than 120 for training ResNets. We note that
the surrogate $F_{\beta}$ loss with a large $\beta$ value, i.e., $\beta=3$ or
$\beta=2$, requires a relatively larger number of epoch for convergence, which
is partly due to the relatively small increment rate of the loss function,
which is illustrated in Figure 1(b).

|
---|---
(a) F1 score on training set | (b) F1 score on validation set
(c) Accuracy on test set
Figure 2: The convergence of the ResNet-18 models trained by using the
proposed loss functions and the other loss functions on the Fashion-MNIST data
set. Model performances are measured by (a) F1 score on the training set, (b)
F1 score on the validation set, and (c) accuracy on the test set. The
performance values are the median from the five repeated experiments.

Figure 3 shows the scatter plot between the surrogate $F_{\beta}$ loss and the
$F_{\beta}$ score on the training set at the first repetition of the
experiment, which compares the gradient paths of the surrogate $F_{\beta}$
loss and the those of the $F_{\beta}$ score. It is clear that the surrogate
$F_{\beta}$ loss is linearly correlated with the $F_{\beta}$ score over all
choices of the $\beta$ values. It implies that, even if the gradient paths of
the surrogate $F_{\beta}$ loss function may fluctuate over epochs (see, e.g.,
Figure 2(a)), they are perfectly aligned with the gradient paths of the actual
$F_{\beta}$ score. The figure demonstrates that the surrogate $F_{\beta}$ loss
can be an effective method for optimization of the $F_{\beta}$ score.

|
---|---
(a) Surrogate $F_{\beta}$ loss v.s. $F_{0.5}$ score | (b) Surrogate $F_{\beta}$ loss v.s. $F_{1.0}$ score
|
(c) Surrogate $F_{\beta}$ loss v.s. $F_{2.0}$ score | (d) Surrogate $F_{\beta}$ loss v.s. $F_{3.0}$ score
Figure 3: Scatter plot between the surrogate $F_{\beta}$ loss and the
$F_{\beta}$ score on the training set of the Fashion-MNIST data set at the
first repetition of the experiment. The $\beta$ value for the $F_{\beta}$
score were set to (a) 0.5, (b) 1.0, (c) 2.0, and (d) 3.0.

For comparison with the macro soft $F_{\beta}$ loss, Figure 4 shows the
scatter plot between the macro soft $F_{\beta}$ loss and the $F_{\beta}$ score
on the training set at the first repetition of the experiment. In this figure,
we can find that the macro soft $F_{\beta}$ loss is not linearly correlated
with the $F_{\beta}$ score except a few cases of $\beta$ values. It implies
that optimization of ResNet models by using the macro soft $F_{\beta}$ loss
may not yield optimal $F_{\beta}$ scores.

|
---|---
(a) Macro soft $F_{\beta}$ loss v.s. $F_{0.5}$ score | (b) Macro soft $F_{\beta}$ loss v.s. $F_{1.0}$ score
|
(c) Macro soft $F_{\beta}$ loss v.s. $F_{2.0}$ score | (d) Macro soft $F_{\beta}$ loss v.s. $F_{3.0}$ score
Figure 4: Scatter plot between the macro soft $F_{\beta}$ loss and the
$F_{\beta}$ score on the training set of the Fashion-MNIST data set at the
first repetition of the experiment. The $\beta$ value for the $F_{\beta}$
score were set to (a) 0.5, (b) 1.0, (c) 2.0, and (d) 3.0.

###  5.3 Results on CIFAR-10 Data Set

ResNet-34 is considered a deeper model with more parameters than ResNet-18.
Figure 5 illustrates the convergence of the ResNet-34 models trained by using
the proposed loss functions and the other loss functions, where all the loss
functions were class-balanced by the inverse class frequency. In Figure 5, the
ResNet models trained by using the surrogate $F_{\beta}$ loss functions
converge relatively fast than the other loss functions, even if they require
more epochs to find good initial weights. On the other hand, other loss
functions except for the macro soft $F_{1}$ loss converge relatively slowly
with respect to the $F_{1}$ score, which is partly due to the mismatch between
loss function and performance measure. Moreover, Figure 5(c) shows that the
macro soft $F_{\beta}$ loss achieves the best accuracy on test set, but the
accuracy of 0.9 is the base accuracy value for the imbalanced data set in this
binary classification task. Instead, in Figure 5(b), we can see that the
surrogate $F_{\beta}$ loss achieves the best $F_{1}$ score on validation set,
which implies that the accuracy measure can be largely misleading.

|
---|---
(a) F1 score on training set | (b) F1 score on validation set
(c) Accuracy on test set
Figure 5: The convergence of the ResNet-34 models trained by using the
proposed loss functions and the other loss functions. Model performances are
measured by (a) F1 score on the training set, (b) F1 score on the validation
set, and (c) accuracy on the test set. The performance values are the median
from the five repeated experiments.

Figure 6 shows the scatter plot between the surrogate $F_{\beta}$ loss and the
$F_{\beta}$ score on the training set at the first repetition of the
experiment, which compares the gradient paths of the surrogate $F_{\beta}$
loss and those of the $F_{\beta}$ score. It is clear that the surrogate
$F_{\beta}$ loss is linearly correlated with the $F_{\beta}$ score over all
choices of the $\beta$ values. Figure 7, on the other hand, shows the scatter
plot between the macro soft $F_{\beta}$ loss and the $F_{\beta}$ score. After
comparing the figures, it is clear that the surrogate $F_{\beta}$ loss
function can generate gradient paths which effectively approximates the
gradient paths of the actual $F_{\beta}$ score.

|
---|---
(a) Surrogate $F_{\beta}$ loss v.s. $F_{0.5}$ score | (b) Surrogate $F_{\beta}$ loss v.s. $F_{1.0}$ score
|
(c) Surrogate $F_{\beta}$ loss v.s. $F_{2.0}$ score | (d) Surrogate $F_{\beta}$ loss v.s. $F_{3.0}$ score
Figure 6: Scatter plot between the surrogate $F_{\beta}$ loss and the $F_{\beta}$ score on the training set at the first repetition of the experiment. The $\beta$ value for the $F_{\beta}$ score were set to (a) 0.5, (b) 1.0, (c) 2.0, and (d) 3.0.  |
---|---
(a) Macro soft $F_{\beta}$ loss v.s. $F_{0.5}$ score | (b) Macro soft $F_{\beta}$ loss v.s. $F_{1.0}$ score
|
(c) Macro soft $F_{\beta}$ loss v.s. $F_{2.0}$ score | (d) Macro soft $F_{\beta}$ loss v.s. $F_{3.0}$ score
Figure 7: Scatter plot between the macro soft $F_{\beta}$ loss and the
$F_{\beta}$ score on the training set at the first repetition of the
experiment. The $\beta$ value for the $F_{\beta}$ score were set to (a) 0.5,
(b) 1.0, (c) 2.0, and (d) 3.0.

##  6 Conclusion and Extensions

In this study, we investigated the conditions of the gradients for the
$F_{\beta}$ score, which is one of the most popular performance measure in
imbalanced data classification. The derived gradient conditions were compared
with those of the standard loss functions, and we proposed a loss function
which is surrogate to the $F_{\beta}$ score. The proposed loss function, which
is called the surrogate $F_{\beta}$ loss, can effectively approximate the
gradient paths of the $F_{\beta}$ score when the network parameters are
optimized by gradient-based learning algorithms such as the stochastic
gradient descent method. The experimental results demonstrated that the
proposed surrogate $F_{\beta}$ loss can be used to effectively optimize
$F_{\beta}$ score of deep neural networks such as ResNet models.

The approach adopted in this study for analyzing the gradient conditions of
the $F_{\beta}$ score can be extended to other performance measures which are
computed based on contingency matrix to build surrogate loss functions. The
proposed surrogate loss function can be used for a scalable optimization
method where mini-batch or parallel optimization approaches are adopted.

Moreover, the proposed loss function can be extended to a generalized cross
entropy to cope with noisy labels in multiclass classification problems (Zhang
and Sabuncu, 2018). For instance, the generalized surrogate $F_{\beta}$ loss
can be written as

${\mathcal{L}}_{F_{\beta},q}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)=y\cdot\frac{1-f({\mathbf{x}};{\boldsymbol{\theta}})^{q}}{q}+(1-y)\cdot\frac{\left(\beta^{2}\cdot\frac{p}{1-p}+f({\mathbf{x}};{\boldsymbol{\theta}})\right)^{q}-1}{q},$
(25)

where the parameter $0<q\leq 1$ controls the level of noise robustness between
MAE loss and BCE loss. That is, we can show that the limit $\lim_{q\rightarrow
0}{\mathcal{L}}_{F_{\beta},q}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)={\mathcal{L}}_{F_{\beta}}({\mathbf{f}}({\mathbf{x}};{\boldsymbol{\theta}}),y)$
approaches the surrogate $F_{\beta}$ loss, and setting $q=1,\beta=1,p=0.5$
leads to the MAE loss.

## References

  * Bartlett et al. (2006) P. L. Bartlett, M. I. Jordan, and J. D. McAuliffe.  Convexity, classification, and risk bounds.  _Journal of the American Statistical Association_ , 101(473):138–156, 2006.  doi: 10.1198/016214505000000907.
  * Beijbom et al. (2012) O. Beijbom, P. J. Edmunds, D. I. Kline, B. G. Mitchell, and D. Kriegman.  Automated annotation of coral reef survey images.  In _2012 IEEE Conference on Computer Vision and Pattern Recognition_ , pages 1170–1177, 2012.  doi: 10.1109/CVPR.2012.6247798.
  * Buda et al. (2018) M. Buda, A. Maki, and M. A. Mazurowski.  A systematic study of the class imbalance problem in convolutional neural networks.  _Neural Networks_ , 106:249–259, 2018.  doi: 10.1016/j.neunet.2018.07.011.
  * Chan and Stolfo (1998) P. K. Chan and S. J. Stolfo.  Toward scalable learning with non-uniform class and cost distributions: A case study in credit card fraud detection.  In _Proceedings of the Fourth International Conference on Knowledge Discovery and Data Mining_ , KDD’98, pages 164–168. AAAI Press, 1998\.
  * Cui et al. (2019) Y. Cui, M. Jia, T.-Y. Lin, Y. Song, and S. Belongie.  Class-balanced loss based on effective number of samples.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019.
  * Eban et al. (2017) E. Eban, M. Schain, A. Mackey, A. Gordon, R. Rifkin, and G. Elidan.  Scalable learning of non-decomposable objectives.  In A. Singh and J. Zhu, editors, _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ , volume 54 of _Proceedings of Machine Learning Research_ , pages 832–840, Fort Lauderdale, FL, USA, 2017. PMLR.
  * G. Van Horn (2017) P. Perona G. Van Horn.  The devil is in the tails: Fine-grained classification in the wild.  _arXiv Preprint_ , arXiv:1709.01450, 2017.
  * Grzymala-Busse et al. (2004) J.W. Grzymala-Busse, L.K. Goodwin, W.J. Grzymala-Busse, and X. Zheng.  An approach to imbalanced data sets based on changing rule strength.  In S. K. Pal, L. Polkowski, and A. Skowron, editors, _Rough-Neural Computing: Cognitive Technologies_ , chapter 21, pages 543–553. Springer, Berlin, Heidelberg, 2004.  doi: https://doi.org/10.1007/978-3-642-18859-6_21.
  * He and Garcia (2009) H. He and E. A. Garcia.  Learning from imbalanced data.  _IEEE Transactions on Knowledge and Data Engineering_ , 21(9):1263–1284, 2009.  doi: 10.1109/TKDE.2008.239.
  * He et al. (2016) K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, 2016.  doi: 10.1109/CVPR.2016.90.
  * Japkowicz and Stephen (2002) N. Japkowicz and S. Stephen.  The class imbalance problem: A systematic study.  _Intelligent Data Analysis_ , 6(5):429–449, 2002\.  doi: 10.3233/IDA-2002-6504.
  * Joachims (2005) T. Joachims.  A support vector method for multivariate performance measures.  In _Proceedings of the 22nd International Conference on Machine Learning_ , ICML ’05, pages 377–384, New York, NY, USA, 2005. Association for Computing Machinery.  doi: 10.1145/1102351.1102399.
  * Johnson and Khoshgoftaar (2019) J. M. Johnson and T. M. Khoshgoftaar.  Survey on deep learning with class imbalance.  _Journal of Big Data_ , 6:27, 2019.  doi: 10.1186/s40537-019-0192-5.
  * Kar et al. (2014) P. Kar, H. Narasimhan, and P. Jain.  Online and stochastic gradient methods for non-decomposable loss functions.  In _Proceedings of the 27th International Conference on Neural Information Processing Systems_ , NIPS’14, pages 694–702, Cambridge, MA, USA, 2014\. MIT Press.
  * Koyejo et al. (2014) O. O. Koyejo, N. Natarajan, P. K. Ravikumar, and I. S. Dhillon.  Consistent binary classification with generalized performance metrics.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc., 2014.
  * Krawczyk (2016) B. Krawczyk.  Learning from imbalanced data: open challenges and future directions.  _Progress in Artificial Intelligence_ , 5:221–232, 2016\.  doi: 10.1007/s13748-016-0094-0.
  * Kubat et al. (1998) M. Kubat, R. C. Holte, and S. Matwin.  Machine learning for the detection of oil spills in satellite radar images.  _Machine Learning_ , 30:195–215, 1998.  doi: 10.1023/A:1007452223027.
  * Mac Namee et al. (2002) B. Mac Namee, P. Cunningham, S. Byrne, and O.I. Corrigan.  The problem of bias in training data in regression problems in medical decision support.  _Artificial Intelligence in Medicine_ , 24(1):51–70, 2002.  doi: https://doi.org/10.1016/S0933-3657(01)00092-6.
  * Narasimhan et al. (2015) H. Narasimhan, P. Kar, and P. Jain.  Optimizing non-decomposable performance measures: A tale of two classes.  In _Proceedings of the 32nd International Conference on International Conference on Machine Learning_ , ICML’15, pages 199–208. JMLR.org, 2015.
  * Pastor-Pellicer et al. (2013) J. Pastor-Pellicer, F. Zamora-Martínez, S. España-Boquera, and M. J. Castro-Bleda.  F-measure as the error function to train neural networks.  In I. Rojas, G. Joya, and J. Gabestany, editors, _Advances in Computational Intelligence_ , pages 376–384, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.
  * Paszke et al. (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  PyTorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Sanyal et al. (2018) A. Sanyal, P. Kumar, P. Kar, S. Chawla, and F. Sebastiani.  Optimizing non-decomposable measures with deep networks.  _Machine Learning_ , pages 1597–1620, 2018.  doi: 10.1007/s10994-018-5736-y.
  * Wang et al. (2017) Y.-X. Wang, D. Ramanan, and M. Hebert.  Learning to model the tail.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * Xiao et al. (2010) J. Xiao, J. Hays, K. A. Ehinger, A. Oliva, and A. Torralba.  SUN database: Large-scale scene recognition from abbey to zoo.  In _2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , pages 3485–3492, 2010.  doi: 10.1109/CVPR.2010.5539970.
  * Yu et al. (2012) H. Yu, J. Ni, Y. Dan, and S. Xu.  Mining and integrating reliable decision rules for imbalanced cancer gene expression data sets.  _Tsinghua Science and Technology_ , 17(6):666–673, 2012.  doi: 10.1109/TST.2012.6374368.
  * Zhang and Sabuncu (2018) Z. Zhang and M. Sabuncu.  Generalized cross entropy loss for training deep neural networks with noisy labels.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.


# Non-Markovian process with variable memory functions

Athokpam Langlen Chanu, Jyoti Bhadana and R.K. Brojen Singh<EMAIL_ADDRESS>(corresponding author) School of Computational & Integrative Sciences,
Jawaharlal Nehru University, New Delhi-110067, India.

###### Abstract

We present a treatment of non-Markovian character of memory by incorporating
different forms of Mittag-Leffler (ML) function, which generally arises in the
solution of fractional master equation, as different memory functions in the
Generalized Kolmogorov-Feller Equation (GKFE). The cross-over from the short
time (stretched exponential) to long time (inverse power law) approximations
of the ML function incorporated in the GKFE is proven. We have found that the
GKFE solutions are the same for negative exponential and for upto frst order
expansion of stretched exponential function for very small $\tau\rightarrow
0$. A generalized integro-differential equation form of the GKFE along with an
asymptotic case is provided.

Keywords: Kolmogorov-Feller equation; Mittag-Leffler function; Stochastic;
Markov; Non-Markov

##  1 Introduction

Stochastic treatment in complex systems study has been a very powerful
technique to study inherent properties in such systems because modelling with
this technique, generally, closely mimics with actual experimental situations
and can reveal dynamical patterns and behaviors, especially the importance of
noise in the system dynamics. Kolmogorov-Feller (KF) equation in stochastic
theory is an integro-differential equation governing the transition
probability density (field) of stochastic processes with jump changes in the
state of the system under study kolmo ; feller1 ; Feller , and has been
applied to various areas of immense interest. In exact , the authors have
produced the exact stationary solutions to the classical KF equation for some
examples of probability distribution analytically. The probability density
function of bounded jump processes having saturation function and Poisson
white noise is deduced and its analytical stationary solutions are found
denisov . KF equation is also used in some topics of astrophysics petro1 ;
petro2 ; galaxy ; comet . This classical KF equation is extended to
generalized Kolmogorov-Feller equation for one and two dimensional cases in
operator representation and is solved using successive approximations
methodszas . The treatment of fractional generalization of the KF equation is
done and a solution for the asymptotic time case is provided saichev .

The Kolmogorov-Feller equation is Markovian or memoryless in character.
However, any natural dynamical process is not truly Markovian. Non-Markovian
character, may be in the form of time delay, exists everywhere whether big or
small. This necessitates a revisit to the classical Markovian Kolmogorov-
Feller equation. In volatility ; iaeng ; stochastic , the classical KF
equation is analyzed to reach a generalized form of the KF equation with the
treatment of memory function in the equation. The Generalised Kolmogorov
Feller Equation (GKFE) is a very important formalism to study stochastic
volatility, for instance in financial system such as financial time series
analysis of stock market iaeng ; volatility . However, the authors in
volatility ; iaeng ; stochastic , even though they study computationally the
solution of GKFE, have not properly discusssed how they take, for a specific
$m(t)$ memory function, the category of functions of $n(t)$ form. In the
present paper, we attempt to use different forms of memory function and
analyse the solution of GKFE analytically. We will discuss how to find, for an
$m(t)$ memory function, the type $n(t)$ particularly the Mittag-Leffler
function.

Mittag-Leffler (ML) function, discovered by Swedish Mathematician,
G$\ddot{o}$sta Mittag-Leffler in 1903 ml ; erdel , is a very important
function in fractional calculus and rightly, known as fractional calculus’s
Queen function. This functional form arises generally in the solution of
master equation which has a fractional form, like, generalized fractional
kinetic equation, the study of super diffusive transport, studies related to
random walk and L$\acute{e}$vy flights and also many complex systems saxena
study. During recent years, ML function has been studied in several applied
problems of physics, chemistry, biology, engineering, earth sciences, etc.
like electrical networks, fluid dynamics, rheology, diffusion transport and
also in theories of probability and statistical distribution. The function ML
is used in the fractional modeling of some stochastic processes like the model
of Continuous Time Random Walk (CTRW) in finance mainardi ; gorenflo . This
fractional calculus has been used to study properties of viscoelastic
materials since it is able to model phenomena with longer memory (hereditary)
vis1 ; vis2 ; vis3 . Kenneth S. Cole (1933) can be quoted for first applying
ML function in connection with nerve conduction; the solution to the
fractional order differential equations has an exponential functional
generalized to the function ML and to the experimental data regarding cell
membrane, this gives a better fit vis1 ; role . Recently, ML function and its
properties have also been studied in neuroscience neuro1 ; neuro2 ; neuro3
and cosmologycosmo ; cosmo1 . In the present paper, we will discuss different
approximations of Mittag-Leffler function and why we can use this function as
a memory function for our GKFE.

Our work is organized in the following way. Section 1 provides a brief
introduction to Classical Kolmogorov-Feller Equation, Generalized Kolmogorov-
Feller Equation and Mittag-Leffler function. Section 2 deals with the
Classical Kolmogorov-Feller equation which is derived from the Brownian motion
and kinetic theory of gas. Generalized Kolmogorov-Feller equation with memory
function is also presented in section 2. Section 3 discusses Mittag-Leffler
function, followed by two common asymptotic approximations of the ML function
given in section 3.1. Section 4 is concerned with results of the analytical
solutions of GKFE using different memory functions. Discussions based on the
results obtained are given in section 5. Section 6 presents conclusion.

##  2 Theory of Kolmogorov-Feller equation

As early as 1905, Albert Einstein has used stochastic modelling of a natural
dynamical process known as Brownian motion albert  which is a kind of random
walk process, based on the idea of Markov chain developed around the same time
by A.A. Markov markov . In albert , Einstein has studied the disordered
Brownian motion of particles in a liquid arising out of thermal molecular
motion that leads to diffusion. In such motion, one can consider a very small
time interval $\tau$ in comparison to the observable time intervals. An
important assumption in this study is that the motions carried out by a single
particle during any two adjacent intervals of time are taken as independent
events. Consider that $N$ total number of Brownian particles is inside the
liquid having volume $V$ such that density function is denoted by
$\mathcal{U}(t,x)$. Although the measured value per unit volume can be
different, the distribution of the particles will almost remain the same
throughout the volume. Thus, using Markov property, one can determine what is
the particle distribution at some time $(t+\tau)$ from the knowledge of its
previous distribution at time, $t$. In one-dimensional scenario with position
$x$, we can thus write the particle number at any $x$ at time $(t+\tau)$ as
albert ,

$\mathcal{U}(t+\tau,x)=\int\limits_{\delta=-\infty}^{\delta=+\infty}\
\mathcal{U}(t,x-\delta)\ f(\delta)\ d\delta$ (1)

We can easily observe the Markov assumption in equation (1). The function,
$f(\delta)$ is any probability density function (p.d.f) for one-dimensional
distance $\delta$ traversed by the particle during $\tau$ time interval. One
can choose $f(\delta)$ depending on the properties and behavior of the system
dynamics, for instance, a Gaussian p.d.f. The definition of convolution of any
two functions, $G(a)$ and $H(a)$ gives,

$\displaystyle G(a)\otimes_{a}H(a)=\int\limits_{0}^{a}G(z)\ H(a-z)\
dz=\int\limits_{0}^{a}G(a-z)\ H(z)\ dz$

Using the above definition and rewriting equation (1), one can arrive at
volatility ; iaeng ,

$\mathcal{U}(t+\tau,x)=\mathcal{U}(t,x)\otimes_{x}f(x)$ (2)

Here, $\mathcal{U}(t+\tau,x)$ can be regarded as a measured value given by
expectation value of $\mathcal{U}(t,x)$ over the range of $x-$values. If we
consider a Markovian (memoryless) process in the time evolution of the system,
such that $\tau$ is negligibly small ($\tau\ll 1$), one can Taylor expand
R.H.S. of equation (2) in powers of $\tau$ and take only up to first order
term as,

$\displaystyle\mathcal{U}(t+\tau,x)\approx\mathcal{U}(t,x)+\tau\frac{\partial\mathcal{U}(t,x)}{\partial
t}$ (3)

From equations (2) and (3),

$\displaystyle\tau\frac{\partial\mathcal{U}(t,x)}{\partial t}$
$\displaystyle=$
$\displaystyle-\mathcal{U}(t,x)+\mathcal{U}(t,x)\otimes_{x}f(x)$ (4)

This evolution equation of the field is generally known as the Classical
Kolmogorov-Feller (KF) Equation volatility . It represents a Markovian process
of the system with no memory. It assumes that only events with short term
memory affect the field’s evolution (dynamics) whereas long term memory (which
depends on histories) events do not affect the behavior of the field.

Around 1930s, Kolmogorov and Feller generalized the transition probability of
Markov stochastic processes with jump changes in the state of the system under
study kolmo ; feller1 . Thus, for a stochastic jump process continuous in
time, Kolmogorov kolmo  wrote the integro-differential equation of the Markov
function $F(t,x)$ as,

$\frac{\partial F(t,x)}{\partial t}=-\int\limits_{-\infty}^{\infty}\ F(t,x)\
G(t,x,z)\ dz+\ \int\limits_{-\infty}^{\infty}F(t,z)\ G(t,x,z)\ dz$

According to Feller feller1 , the forward and backward integro-differential
equations are respectively,

$\frac{\partial P(\tau,x;t,\Lambda)}{\partial
t}=-\int\limits_{\Lambda}p(t,y)P(\tau,x;t,dE_{y})+\int\limits_{E}p(t,y)\Pi(t,y,\Lambda)P(\tau,x;t,dE_{y})$

and

$\frac{\partial
P(\tau,x;t,\Lambda)}{\partial\tau}=p(\tau,x)[P(\tau,x;t,\Lambda)-\int\limits_{E}P(\tau,y;t,\Lambda)\Pi(\tau,x,dE_{y})]$

Now, in equation (4), the investigation of the behavior of $\mathcal{U}(t,x)$
due to longer term events can be done by incorporating the role of a memory
function, $m(t)$ in the KF equation (4) as given in volatility ; iaeng  as

$\displaystyle m(t)\otimes_{t}\frac{\partial\mathcal{U}(t,x)}{\partial
t}=-\frac{\mathcal{U}(t,x)}{\tau}+\frac{1}{\tau}\
\mathcal{U}(t,x)\otimes_{x}f(x)$ (5)

If one takes $m(t)=\delta(t)$, because of the convolution property of the
delta function $\delta\otimes h=h$, R.H.S. in the equation (5) becomes,
$\delta(t)\otimes_{t}\frac{\partial\mathcal{U}(t,x)}{\partial
t}=\frac{\partial\mathcal{U}(t,x)}{\partial t}$, and hence, equation (5)
reduces to classical KF equation (4). The classical KF equation is classical
in the sense that when the delta function is used as a memory function, one
can trace the exact position and time (memory) of the particle which gives us
a classical insight. Now, let us consider that any arbitrary $m(t)$ memory
function has a category of functions $n(t)$ given by the condition
$n(t)\otimes_{t}m(t)=\delta(t)$ volatility , then the equation (5) can be
written in another format as given in volatility ; iaeng  as

$\displaystyle\tau\frac{\partial\mathcal{U}(t,x)}{\partial t}\
=-n(t)\otimes_{t}\mathcal{U}(t,x)+n(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x)$
(6)

Equation (6) then becomes the Generalized Kolmogorov-Feller equation (GKFE).
Again, if $n(t)=\delta(t)$ in GKFE, then one can easily prove that GKFE yields
classical KF equation given by equation (4). We will be following the form of
GKFE solution as given by Blackledge et al. volatility ; iaeng  . After adding
$\mathcal{U}(t,x)$ on both sides of equation (6)and solving using Green
function technique, then the GKFE solution is volatility ; iaeng

$\displaystyle\mathcal{U}(t,x)=-g(t)\otimes_{t}n(t)\otimes_{t}\mathcal{U}(t,x)+g(t)\otimes_{t}n(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x)+g(t)\otimes_{t}\mathcal{U}(t,x)$
(7)

Here, the function $g(t)$ denotes the Green function which is solution of the
following equation,

$\displaystyle\tau\frac{\partial g(t-t_{0})}{\partial
t}+g(t-t_{0})=\delta(t-t_{0})$ (8)

Its solution is $\ g(t)=\frac{1}{\tau}\mathrm{e}^{({\frac{-t}{\tau}})}$, for
time, $t>0$. The proof is simple, and can be proceeded by taking $t-t_{0}=z$.
Then, $g(z)$ is given by $g(z)=\int\limits_{t_{0}}^{\infty}\mathcal{G}(t,z)\
f(z)\ dz$. The $f(z)$ is a known function and $\mathcal{G}(t,z)$ is its
fundamental solution. We have,

$\displaystyle\frac{\partial\mathcal{G}(t,z)}{\partial
t}+\frac{\mathcal{G}(t,z)}{\tau}=\frac{\delta(t-z)}{\tau}$ (9)

Its homogeneous solution is
$\mathcal{G}_{hom}(t,z)=\mathrm{e}^{\left(-\int\limits_{z}^{t}\
\frac{1}{\tau}\ dt\right)}$. We know that,
$\int\limits_{t_{0}}^{\infty}\delta(t-z)\ dz=\textbf{H}(t-z)$, where Heaviside
function, $\textbf{H}(t-z)$ has the following definition,

$\textbf{H}(t-z)=\begin{cases}0&\text{; $t<z$}\\\ 1&\text{; $t>z$}\\\
\end{cases}$

Therefore, the fundamental solution is
$\mathcal{G}(z,t)=\frac{1}{\tau}\mathrm{e}^{\left(-\int\limits_{z}^{t}\frac{1}{\tau}\
dt\right)}\times\textbf{H}(t-z)$. Thus, the above solution to equation (8) is
proved to be $g(t)=\frac{1}{\tau}e^{(\frac{-t}{\tau})}$ for time, $t>0$ with
intial $t_{0}=0$. Now, taking simultaneous Fourier (k-space) and Laplace
(s-space) transforms of equation (7) as discussed in volatility ; iaeng , one
can arrive at,

$\displaystyle\widetilde{\overline{\mathcal{U}}}(s,k)=-\overline{g}(s)\
\overline{n}(s)\ \widetilde{\overline{\mathcal{U}}}(s,k)+\overline{g}(s)\
\overline{n}(s)\ \widetilde{\overline{\mathcal{U}}}(s,k)\
\widetilde{f}(k)+\overline{g}(s)\ \widetilde{\overline{\mathcal{U}}}(s,k)$
(10)

where, by definition,

$\displaystyle\widetilde{\overline{\mathcal{U}}}(s,k)$ $\displaystyle=$
$\displaystyle\int\limits_{t=0}^{\infty}\int\limits_{x=-\infty}^{+\infty}\
e^{-ikx}\ e^{-st}\ \mathcal{U}(t,x)\ dx\ dt$ $\displaystyle\overline{g}(s)$
$\displaystyle=$ $\displaystyle\int\limits_{t=0}^{\infty}\ e^{-st}\ g(t)\
dt;~{}~{}\overline{n}(s)=\int\limits_{t=0}^{\infty}\ e^{-st}\ n(t)\
dt;~{}~{}\widetilde{f}(k)=\int\limits_{x=-\infty}^{+\infty}\ e^{-ikx}\ f(x)\
dx$

Now, solving equation (10), one can get volatility ,

$\displaystyle\widetilde{\overline{\mathcal{U}}}(s,k)=-\frac{\overline{g}(s)\
\overline{n}(s)\
\widetilde{\overline{\mathcal{U}}}(s,k)}{[1-\overline{g}(s)]}+\frac{\overline{g}(s)\
\overline{n}(s)\ \widetilde{\overline{\mathcal{U}}}(s,k)\
\widetilde{f}(k)}{[1-\overline{g}(s)]}$ (11)

The functional form of $\overline{g}(s)$ can be obtained by taking Laplace
transform of the equation (8), where, $y=t$, with the initial condition
$g(0)=0$ i.e., one can get as
$\tau\left[s\overline{g}(s)-g_{0}\right]+\overline{g}(s)=1$. From this, one
can get the expression, $\overline{g}(s)=\frac{1}{1+\tau s}$. Now, use
$\frac{\overline{g}(s)}{1-\overline{g}(s)}=\frac{1}{\tau s}$ in equation (11),
and solve as follows volatility ; iaeng ,

$\displaystyle\tau s\ \widetilde{\overline{\mathcal{U}}}(s,k)$
$\displaystyle=-\overline{n}(s)\
\widetilde{\overline{\mathcal{U}}}(s,k)+\overline{n}(s)\
\widetilde{\overline{\mathcal{U}}}(s,k)\ \widetilde{f}(k)$

Solving the equation, we get,

$\displaystyle\widetilde{\overline{\mathcal{U}}}(s,k)$
$\displaystyle=\overline{h}(s)\ \widetilde{\overline{\mathcal{U}}}(s,k)\
\widetilde{f}(k)$ (12)

where,

$\overline{h}(s)=\frac{\overline{n}(s)}{\tau s+\overline{n}(s)}$ (13)

Then, one can obtain the functional form of $h(t)$ by the application of
Inverse Laplace transform to exoression of $\overline{h}(s)$ of equation (13),
. This function $h(t)$ relates to the category memory function $n(t)$
indirectly via equation (13) and this category function $n(t)$ again is
related to the main memory function $m(t)$ via condition volatility .

Further, by applying both Inverse Fourier and Inverse Laplace transforms to
$\widetilde{\overline{\mathcal{U}}}(s,k)$ of equation (12), we finally get a
self-consistent form of solution $\mathcal{U}(t,x)$ of GKFE as Blackledge et
al. derived volatility ; iaeng  ,

$\displaystyle\mathcal{U}(t,x)=h(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x)$
(14)

From equation (14), we can see that the concept of memory is incorporated in
the solution of GKFE through $h(t)$. Different functional forms of $h(t)$ will
give different GKFE solutions, $\mathcal{U}(t,x).$ This signifies the role of
memory functions in the GKFE solutions. In the subsequent sections, we will
study these GKFE solutions with respect to different memory functional forms
derived from one particular Mittag-Leffler function.

##  3 Memory functional forms derived from Mittag-Leffler Function,
$E_{\beta}$

The Mittag-Leffler (ML) function of order $\beta$, $E_{\beta}$ represents a
complete transcendental function defined in the $\mathbb{Z}$-complex plane as
mainardi ; gorenflo ; saxena ,

$\displaystyle
E_{\beta}(z):=\sum\limits_{n=0}^{\infty}\frac{z^{n}}{\Gamma(\beta n+1)},\ \
where\ \ \beta>0$

Let us denote, for time, $t>0$,

$\displaystyle e_{\beta}(t):=E_{\beta}(-t^{\beta})=\sum\limits_{n=0}^{\infty}\
(-1)^{n}\ \frac{t^{\beta n}}{\Gamma(1+\beta n)}\ ,\ \ where\ \ \beta\ is\
such\ that\ 0<\beta<1$ (15)

A negative exponential function which is characteristic of a memoryless
Markovian system is obtained in the limit $\beta\rightarrow 1$.

###  3.1 Asymptotic approximations of Mittag-Leffler function

From equation (15) (the representation of power series which is convergent),
the ML function’s short time approximation is obtained as follows.

$\displaystyle e_{\beta}^{(0)}(t)=1-\frac{1}{\Gamma(1+\beta)}t^{\beta}+...\sim
e^{-\frac{t^{\beta}}{\Gamma(\beta+1)}},\ \ \ \ t\rightarrow 0$ (16)

Thus, for small time $t\rightarrow 0$, the power series ML function of
equation (15) approximates to stretched exponential of equation (16). Again,
from the asymptotic power series representation of equation (15), we can get
the ML function’s long term approximation as erdel ; mainardi ; gorenflo ,

$\displaystyle e_{\beta}^{(\infty)}(t)\sim\sum\limits_{n=1}^{\infty}\
(-1)^{n-1}\ \frac{t^{-\beta n}}{\Gamma(1-\beta
n)}\sim\frac{t^{-\beta}}{\Gamma(1-\beta)}\ ,\ \ \ \ t\rightarrow\infty$ (17)

Thus, for large time $t\rightarrow\infty$, the power series ML function of
equation (15) approximates to inverse power law of equation (17). From
equations (16) and (17), we find that the Mittag-Leffler function,
$e_{\beta}(t)$ of equation (15) is the interpolation of intermediate time
between two limiting behaviours, $t\rightarrow 0$ and $t\rightarrow\infty.$
The stretched exponential of equation (16) represents a model of the decay
which is very fast (or memoryless) for small time. On the otherhand, the
inverse power law of equation (17) is because of the decay very slow in nature
(having some memory) for large time. Thus, the Mittag-Leffler function is a
good candidate for being a complete memory function since it can represent
short term and long term memory in its limits. The Mittag-Leffler function,
$e_{\beta}(t)$ still maintains its “complete monotonic” character mainardi
for $0<\beta<1$, even though the function loses its behaviour of very fast
(exponential) decay by exibiting power law tails at larger time. The “complete
monotonic” character of the function, $e_{\beta}(t)$ means erdel

$(-1)^{m}\ \frac{d^{m}}{dt^{m}}e_{\beta}(t)\ \geq 0,\ \ \ \ \ where\ \
time,t\geq 0,\ \ 0\leq\beta\leq 1,\ and\ \ \ n=0,1,...$

##  4 Results: GKFE solutions driven by different forms of Mittag-Leffler
memory function

This section presents analytical results regarding different forms of memory
functions. In Case-I, upto first order expansion of stretched exponential of
equation (16) is taken as memory function $m(t).$ In Case-II, the memory
function $m(t)$ is taken to be the negative power law functional form of
equation (17). In Case-III, the memory function $m(t)$ is the intermediate
functional form of the Mittag-Leffler function of equation (15). In Case-IV,
negative exponential is taken as a memory function $m(t)$. In all the cases,
respective final GKFE solutions are presented. Detailed calculations are given
below.

Case-I: For time, $t\rightarrow 0$ (short memory)
Consider the memory function in this case to be
$m(t)=1-\frac{t^{\beta}}{\Gamma(1+\beta)}$. To find GKFE solution associated
with this memory function, we need to find the expression for $\bar{n}(s)$ to
get the functional form of $h(t)$. In order to find the form of $\bar{n}(s)$
to get the solution of GKFE, we use the condition
$n(t)\otimes_{t}m(t)=\delta(t)$ volatility ; stochastic  as follows,

$\displaystyle\int\limits_{0}^{t}n(\tau)\left[1-\frac{(t-\tau)^{\beta}}{\Gamma(1+\beta)}\right]d\tau=\delta(t)$
(18)

Now, applying the change of variables i.e., $t\rightarrow\tau$ and
$\tau\rightarrow s$ in equation (18) and then multiplying by
$(t-\tau)^{-\beta-1}$ and again integrating w.r.t $\tau$, it gives

$\displaystyle\int\limits_{0}^{t}d\tau\
(t-\tau)^{-\beta-1}\int\limits_{0}^{\tau}n(s)ds\left[1-\frac{(\tau-s)^{\beta}}{\Gamma(1+\beta)}\right]$
$\displaystyle=$ $\displaystyle\int\limits_{0}^{t}d\tau\ \delta(\tau)\
(t-\tau)^{-(\beta+1)}=t^{-(\beta+1)}$ (19)

Then, applying Fubini’s theorem of change of integrating variables to the
L.H.S., we get,

$\displaystyle\int\limits_{0}^{t}n(s)ds\int\limits_{s}^{t}d\tau\
(t-\tau)^{-\beta-1}\left[1-\frac{(\tau-s)^{\beta}}{\Gamma(1+\beta)}\right]=t^{-(\beta+1)}$
(20)

Next, in order to solve the integration, the change of variables,
$Y=\frac{\tau-s}{t-s}$, is introduced where, $\tau=s+Y(t-s)$ and
$d\tau=dY(t-s)$. We then integrate the second integration in equation (20),
and after simplification and arranging the terms, we arrive at,

$\displaystyle\int\limits_{0}^{t}ds\
n(s)\left[\frac{-(t-s)^{-\beta}}{\beta}-\frac{1}{\Gamma(1+\beta)}\
\textbf{\large{B}}(\beta+1,-\beta)\right]=t^{-\beta-1}$ (21)

where, Beta function
$\textbf{\large{B}}(1+\beta,-\beta)=\int\limits_{0}^{1}dY\
Y^{\beta}(1-Y)^{-1-\beta}$. Using the relation
$\textbf{\large{B}}(y,z)=\frac{\Gamma(y)\ \Gamma(z)}{\Gamma(y+z)}$, we get,

$\displaystyle\int\limits_{0}^{t}dsn(s)\left[\frac{-(t-s)^{-\beta}}{\beta}-\Gamma(-\beta)\right]=t^{-(\beta+1)}$
(22)

Now, on both sides of the above equation (22), we take Laplace transform and
after simplification, we get,

$\displaystyle\bar{n}(s)\left[\frac{-L\left\\{t^{-\beta}\\}\right]}{\beta}-\Gamma(-\beta)L\left\\{1\right\\}\right]=L\\{t^{-(\beta+1)}\\}$
(23)

By the definition of Laplace transform, we know,
$L\left\\{t^{-\beta}\right\\}=\frac{\Gamma(1-\beta)}{s^{1-\beta}}$,
$L\left\\{1\right\\}=1/s$, and $L\\{t^{-\beta-1}\\}=s^{\beta}\
\Gamma(-\beta)$. Then using the property, $\Gamma(1-\beta)=-\beta\
\Gamma(-\beta)$, we obtain the expression for $\bar{n}(s)$ as follows,

$\displaystyle\bar{n}(s)=\frac{s}{1-s^{-\beta}}$ (24)

Now, the expression for $\bar{n}(s)$ is put in the equation (13), and for the
condition $\tau\rightarrow 0$, it can be shown that $\bar{h}(s)=1$. Then
taking inverse Laplace transform, we can obtain the expression for $h(t)$ as,
$h(t)\approx\delta(t)$, which is characteristic of classical Markovian
process. Then putting the expression for $h(t)$ in the equation (14), we
obtain the self-consistent GKFE solution for this short term memory function
as given by,

$\displaystyle\mathcal{U}(t,x)$ $\displaystyle=$
$\displaystyle\mathcal{U}(t,x)\otimes_{x}f(x)$ (25)

Iteratively, this solution can be represented as
$\mathcal{U}_{n+1}(t,x)=\mathcal{U}_{n}(t,x)\otimes_{x}f(x)$, with convolution
over $x$. Here, the function $f(x)$ can be any known function according to
one’s choice. f(x) is system dependent and one can analyse the resultant
equation (25) according to it. For example, let us take a Gaussian
distribution given by, $f(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$.
Hence, equation (25) becomes,
$\mathcal{U}(t,x)=\mathcal{U}(t,x)\otimes_{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}.$

Case-II: For large time, $t\rightarrow\infty$ (long memory)
Consider memory function to be $m(t)=\frac{t^{-\beta}}{\Gamma(1-\beta)}$. In
this case also, in order to find out the form of $\bar{n}(s)$, we use the same
condition $n(t)\otimes_{t}m(t)=\delta(t)$,

$\displaystyle\int\limits_{0}^{t}n(\tau)\frac{d\tau}{\Gamma(1-\beta)(t-\tau)^{\beta}}=\delta(t)$
(26)

We now replace $(1-\beta)=\alpha$ in the equation (26). We proceed the change
of the variable from $t$ to $\tau$ and $\tau$ to $s$. Further, in the
equation, the term $(t-\tau)^{-\alpha}$ is multiplied on both sides and with
reference to $\tau$, it is integrated. Then, the other of the integration is
exchanged using Fubini’s theorem and after rearranging the terms, we get,

$\displaystyle\int\limits_{0}^{t}n(s)\ ds\
\int\limits_{s}^{t}\frac{1}{(t-\tau)^{\alpha}}\times\frac{1}{(\tau-s)^{1-\alpha}}d\tau=\int\limits_{0}^{t}\frac{\Gamma(\alpha)\
\delta(\tau)\ d\tau}{(t-\tau)^{\alpha}}$ (27)

Now, we put, $\tau=s+Y(t-s)$ in equation (27), such that $d\tau=dY(t-s)$ in
L.H.S. We then use the Beta function property,
$\int\limits_{0}^{1}Y^{\alpha-1}(1-Y)^{-\alpha}dY=\textbf{B}(\alpha,1-\alpha)=\Gamma(\alpha)\Gamma(1-\alpha)$,
we get,

$\displaystyle\int\limits_{0}^{t}\ \Gamma(1-\alpha)\ n(s)\
ds=\int\limits_{0}^{t}\frac{1}{(t-\tau)^{\alpha}}\delta(\tau)\ d\tau$ (28)

From this, we can get the expression for $n(t)$, given by,

$\displaystyle n(t)$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(1-\alpha)}\times\frac{d}{dt}\int\limits_{0}^{t}\frac{\delta(\tau)\
d\tau}{(t-\tau)^{\alpha}}=\frac{1}{\Gamma(\beta)}\frac{d}{dt}t^{-\alpha}=\frac{1}{\Gamma(\beta-1)\
t^{2-\beta}}$ (29)

Using the standard integral, $\int\limits_{0}^{\infty}\frac{1}{\Gamma(\beta)\
t^{1-\beta}}\exp(-st)dt=\frac{1}{s^{\beta}}$, we get,
$\overline{n}(s)=\frac{1}{s^{\beta-1}}$. Now, by substituting this form of
$\overline{n}(s)$ to the equation (13), we get
$\overline{h}(s)=\frac{1}{1+\tau s^{\beta}}$. If we consider a non-Markovian
process i.e.,$\tau\gg$ 1, then $\overline{h}(s)=\frac{1}{\tau s^{\beta}}$ and
taking inverse Laplace transform of $\overline{h}(s)$, we can write,
$h(t)\sim\frac{1}{\tau\ \Gamma(\beta)\ t^{1-\beta}}$. Thus, the GKFE solution,
which incorporates long time memory is,

$\displaystyle\mathcal{U}(t,x)$ $\displaystyle=$ $\displaystyle\frac{1}{\tau\
\Gamma(\beta)\ t^{1-\beta}}\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x)$ (30)

Iteratively, this solution can also be represented as
$\mathcal{U}_{n+1}(t,x)=\frac{1}{\tau\ \Gamma(\beta)\
t^{1-\beta}}\otimes_{t}\mathcal{U}_{n}(t,x)\otimes_{x}f(x)$, with convolution
over $t$ and $x$. We observe that this final GKFE solution for long term
memory depends on the parameter, $\beta.$ Here, the function $f(x)$ is system
dependent and any known function can be taken according to one’s choice for
analysis. For example, let us consider a Gaussian distribution,
$f(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$. Hence, equation (30)
becomes, $\mathcal{U}(t,x)=\frac{1}{\tau\ \Gamma(\beta)\
t^{1-\beta}}\otimes_{t}\mathcal{U}(t,x)\otimes_{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}.$

Case-III: For intermediate time, $0<t<\infty$

In this case, the memory function is taken to be the Mittag-Leffler function
or, $m(t)=e_{\beta}(t)=E_{\beta}(-t^{\beta})$, given by equation(15). We then
substitute this function directly to the relation
$n(t)\otimes_{t}m(t)=\delta(t)$ volatility  such that
$n(t)\otimes_{t}E_{\beta}(-t^{\beta})=\delta(t)$. Taking on both sides of this
equation the Laplace transform, and employing the relation
$L\\{E_{\beta}(-t^{\beta})\\}=\frac{s^{\beta-1}}{1+s^{\beta}}$ saxena ; feller
, we get the expression for $\bar{n}(s)$ given by,
$\bar{n}(s)=\frac{1+s^{\beta}}{s^{\beta-1}}$. Now, putting this expression of
$\bar{n}(s)$ to the equation (13), we get the expression for $\bar{h}(s)$
given by, $\bar{h}(s)=\frac{1+s^{\beta}}{1+s^{\beta}+\tau s^{\beta}}$ and
$h(t)$ will be its inverse Laplace transform which also depends on $\beta$.
The form of the GKFE solution in this case will be given by
$h(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x).$ Iteratively, this solution
can also be represented as
$\mathcal{U}_{n+1}(t,x)=h(t)\otimes_{t}\mathcal{U}_{n}(t,x)\otimes_{x}f(x)$,
with convolution over $x$ and $t$ and depends on $\beta$. Here, the function
$f(x)$ is system dependent and any known function can be taken according to
one’s choice for analysis. For example, let us consider a Gaussian
distribution, $f(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$. Hence,
$\mathcal{U}(t,x)=h(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}.$

Case-IV: For $\beta\rightarrow 1$, using equation (15), the Mittag-Leffler
function approximates to,

$\displaystyle e(t)_{(\beta\rightarrow 1)}=\sum\limits_{n=0}^{\infty}\ \
\frac{(-1)^{n}\ t^{n}}{\Gamma(n+1)}=e^{-t}$

Now, consider the memory function to be exponentially decaying $m(t)=e^{-t}$.
To find the GKFE solution for this memory function, we use the condition
$n(t)\otimes_{t}m(t)=\delta(t)$ volatility ; stochastic  such that after
taking Laplace transform on both sides of the condition, $L[n(t)]\
L[e^{-t}]=1$, such that, $\bar{n}(s)\ \frac{1}{1+s}=1$, we get,
$\bar{n}(s)=1+s$. Now, equation (13) becomes,

$\displaystyle\bar{h}(s)=\frac{1+s}{\tau s+(1+s)}=\left(1+\frac{\tau
s}{1+s}\right)^{-1}=1-\frac{\tau}{1+\frac{1}{s}}=1-\tau+\frac{\tau}{s}$

For memoryless process $(\tau\rightarrow 0)$, $\bar{h}(s)=1$ and,
$h(t)=\delta(t).$ Hence, from equation (14), the self-consistent GKFE solution
for this memory function becomes,

$\mathcal{U}(t,x)=\delta(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x)=\mathcal{U}(t,x)\otimes_{x}f(x)$
(31)

Iteratively, this solution can be represented as
$\mathcal{U}_{n+1}(t,x)=\mathcal{U}_{n}(t,x)\otimes_{x}f(x)$, with convolution
over $x$. Here, the function $f(x)$ can be any known p.d.f according to
system. For example, let us take a Gaussian distribution given by,
$f(x)=f(x|\mu=0,\sigma^{2}=1)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$.
Hence,
$\mathcal{U}(t,x)=\mathcal{U}(t,x)\otimes_{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}.$

##  5 Discussions

Theorem 1: The Generalised Kolmogorov-Feller equation solutions are the same
for negative exponential function and for upto first order expansion of
stretched exponential function for very small $\tau$.

Proof: In Case-IV, for negative exponential memory function $m(t)=e^{-t}$,
$\bar{h}(s)=1-\tau+\frac{\tau}{s}.$ For very small $\tau$ or,
$(\tau\rightarrow 0)$, $\bar{h}(s)=1,or\ h(t)=\delta(t).$ Now, in case-I, for
first order expansion of stretched exponential function as memory function
i.e., $m(t)=1-\frac{t^{\beta}}{\Gamma(\beta+1)}$, we get,
$\bar{h}(s)=\frac{1}{1+\tau-\frac{\tau}{s^{\beta}}}.$ Again, for
$(\tau\rightarrow 0),\ \bar{h}(s)=1,or\ h(t)=\delta(t).$ From equation, (14),
the form of Generalised Kolmogorov Feller Equation solution will be given by
$\mathcal{U}(t,x)=h(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x).$ Hence, in
both the cases, the final solution is
$\mathcal{U}(t,x)=\delta(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x)=\mathcal{U}(t,x)\otimes_{x}f(x)$,
for any arbitrary p.d.f $f(x)$, where we have used the property of
$\delta$-function that $\delta(t)\otimes_{t}f(t)=f(t)$. Hence proved.

Theorem 2: Mittag-Leffler function for intermediate memory function cross-
overs from $t\rightarrow 0$ to $t\rightarrow\infty$.

Proof: The expression of $\bar{n}(s)$ for intermediate memory function
$(0<t<\infty)$ is $\bar{n}(s)=\frac{1+s^{\beta}}{s^{\beta-1}}$. The expression
for $\bar{h}(s)$ is obtained by substituting the expression of $\bar{n}(s)$ to
the equation (13), and is given by,
$\bar{h}(s)=\frac{1+s^{\beta}}{1+s^{\beta}+\tau s^{\beta}}$. Now, if $\tau\ll
1$, this equation for $\bar{h}(s)$ becomes $\bar{h}(s)\approx 1$. Taking its
Inverse Laplace transform, we get the same result which was obtained in the
Case-I, $h(t)\approx\delta(t)$. On the other hand, rewriting the expression
for $\bar{h}(s)$ in the form, $\bar{h}(s)=\frac{1+s^{\beta}}{1+s^{\beta}+\tau
s^{\beta}}=\frac{1}{1+\frac{s^{\beta}\tau}{1+s^{\beta}}}=\frac{1}{1+\frac{\tau}{1+s^{-\beta}}}$,
and, if $\tau\gg 1$, then the expression of $\bar{h}(s)$ becomes the same
expression which we got in the Case-II,
$\bar{h}(s)\approx\frac{1}{\frac{\tau}{1+s^{-\beta}}}\approx\frac{1+s^{-\beta}}{\tau}\approx\frac{s^{-\beta}}{\tau}+\frac{1}{\tau}\approx\frac{1}{\tau
s^{\beta}}$. Taking its Inverse Laplace transform, we get the same result
which was obtained in the Case-II, $h(t)\sim\frac{1}{\tau\ \Gamma(\beta)\
t^{1-\beta}}$. Thus, we can conclude the Mittag-Leffler function in
intermediate case provides us the cross-over from short to long term memory
results as discussed in the Case-I and Case-II.

Theorem 3: A generalized integro-differential equation of the GKFE is,

$\displaystyle\frac{\partial}{\partial
t}\left[\int\limits_{0}^{x}\int\limits_{0}^{t}\ dx_{1}\ dt_{1}\ h(t_{1})\
\mathcal{U}(x_{1},t-t_{1})\ f(x-x_{1})\right]$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau}\int\limits_{0}^{x}\int\limits_{0}^{t}\ dx_{1}\
dt_{1}\ n(t_{1})\ \mathcal{U}(x_{1},t-t_{1})\ f(x-x_{1})$
$\displaystyle+\frac{1}{\tau}\int\limits_{0}^{x}\int\limits_{0}^{x}\int\limits_{0}^{t}dx_{2}dx_{1}dt\
n(t_{1})\ \mathcal{U}(x_{2},t-t_{1})\ f(x_{1})f(x-x_{2}-x_{1})$

Proof: In order to write a generalized integro-differential equation of the
GKFE of equation (6), we can write the left hand side of the equation (6),

$\displaystyle\tau\frac{\partial\mathcal{U}(t,x)}{\partial t}$
$\displaystyle=$ $\displaystyle\tau\frac{\partial}{\partial
t}\left[h(t)\otimes_{t}\mathcal{U}(t,x)\otimes_{x}f(x)\right]=\tau\frac{\partial}{\partial
t}\left[h(t)\otimes_{t}\int\limits_{0}^{x}\mathcal{U}(t,x_{1})\ f(x-x_{1})\
dx_{1}\right]$ $\displaystyle=$ $\displaystyle\tau\frac{\partial}{\partial
t}\left[\int\limits_{0}^{x}\int\limits_{0}^{t}\ dx_{1}\ dt_{1}\ h(t_{1})\
\mathcal{U}(x_{1},t-t_{1})\ f(x-x_{1})\right]$

Again, the right hand side of the equation (6) becomes,

$\displaystyle-n(t)\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)+n(t)\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)\otimes_{x}f(x)$
$\displaystyle=$
$\displaystyle-n(t)\otimes_{t}\int\limits_{0}^{x}\mathcal{U}(x_{1},t)\
f(x-x_{1})\
dx_{1}+n(t)\otimes_{t}\mathcal{U}(x,t)\otimes_{x}\int\limits_{0}^{x}f(x_{1})\
f(x-x_{1})dx_{1}$ $\displaystyle=$
$\displaystyle-\int\limits_{0}^{x}n(t)\otimes_{t}\mathcal{U}(x_{1},t)\
f(x-x_{1})\
dx_{1}+n(t)\otimes_{t}\int\limits_{0}^{x}\mathcal{U}(x_{2},t)\int\limits_{0}^{x}f(x_{1})\
f(x-x_{2}-x_{1})dx_{2}dx_{1}$ $\displaystyle=$
$\displaystyle-\int\limits_{0}^{x}\int\limits_{0}^{t}n(t_{1})\
\mathcal{U}(x_{1},t-t_{1})\ dt_{1}\ f(x-x_{1})\
dx_{1}+\int\limits_{0}^{t}n(t_{1})\int\limits_{0}^{x}\mathcal{U}(x_{2},t-t_{1})\
dt_{1}\int\limits_{0}^{x}f(x_{1})\ f(x-x_{2}-x_{1})\ dx_{2}\ dx_{1}$
$\displaystyle=$ $\displaystyle-\int\limits_{0}^{x}\int\limits_{0}^{t}\
dx_{1}\ dt_{1}\ n(t_{1})\ \mathcal{U}(x_{1},t-t_{1})\
f(x-x_{1})+\int\limits_{0}^{x}\int\limits_{0}^{x}\int\limits_{0}^{t}n(t_{1})\
\mathcal{U}(x_{2},t-t_{1})\ f(x_{1})\ f(x-x_{2}-x_{1})$

Hence, the generalized integro-differential equation of the GKFE in equation
(5) is proved.

Lemma 1: The integro-differential form of GKFE for the asymptotic case (i.e.,
$t\rightarrow\infty$) is

$\displaystyle\frac{1}{\beta}\frac{\partial}{\partial
t}\left[\int\limits_{0}^{x}\int\limits_{0}^{t}\ dx_{1}\ dt_{1}\
t_{1}^{\beta-1}\ \mathcal{U}(x_{1},t-t_{1})\ f(x-x_{1})\right]$
$\displaystyle=$ $\displaystyle-\frac{\Gamma(\beta)}{\tau\
\Gamma(2\beta-1)}\left[\int\limits_{0}^{x}\int\limits_{0}^{t}dx_{1}\ dt_{1}\
t_{1}^{2\beta-2}\ \mathcal{U}(x_{1},t-t_{1})\ f(x-x_{1})\right.$
$\displaystyle+\left.\int\limits_{0}^{x}\int\limits_{0}^{x}\int\limits_{0}^{t}dx_{2}dx_{1}dt_{1}\
t_{1}^{2\beta-2}\ \mathcal{U}(x_{2},t_{1}-t)f(x_{1})f(x-x_{2}-x_{1})\right]$

Proof: For the case of large time, $t\rightarrow\infty$, we get,

$\displaystyle\mathcal{U}(x,t)=\frac{1}{\tau\Gamma(\beta)}t^{\beta-1}\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)$

Putting this to the equation (6), L.H.S. of equation (6) becomes,

$\displaystyle\tau\frac{\partial}{\partial
t}\mathcal{U}(t,x)=\tau\frac{\partial}{\partial
t}\left[\frac{t^{\beta-1}}{\tau\
\Gamma(\beta)}\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)\right]=\frac{1}{\Gamma(\beta)}\frac{\partial}{\partial
t}\left[\int\limits_{0}^{x}\int\limits_{0}^{t}\ dx_{1}\ dt_{1}\
t_{1}^{\beta-1}\ \mathcal{U}(x_{1},t-t_{1}\ f(x-x_{1})\right]$

For $n(t)=\frac{t^{\beta-2}}{\Gamma(\beta-1)}$, we get,

$\displaystyle R.H.S.~{}of~{}(6)$ $\displaystyle=$
$\displaystyle-\frac{t^{\beta-2}}{\Gamma(\beta-1)}\otimes_{t}\frac{t^{\beta-1}}{\tau\
\Gamma(\beta)}\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)+\frac{t^{\beta-2}}{\Gamma(\beta-1)}\otimes_{t}\frac{t^{\beta-1}}{\tau\
\Gamma(\beta)}\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)\otimes_{x}f(x)$

Now,

$\displaystyle t^{\beta-2}\otimes_{t}t^{\beta-1}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{t}\tau^{\beta-2}\ (t-\tau)^{\beta-1}\
d\tau=\int\limits_{0}^{t}\tau^{\beta-2}\ \ t^{\beta-1}\
\left(1-\frac{\tau}{t}\right)^{\beta-1}\ d\tau$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{1}(tX)^{\beta-2}\ \ t^{\beta-1}\
\left(1-X\right)^{\beta-1}\ t\ dX\ \ \ \ \ \ \ \ \ \ \ \ \ (Putting\ \tau=tX)$
$\displaystyle=$ $\displaystyle t^{2\beta-2}\int\limits_{0}^{1}\ \
\left(1-X\right)^{\beta-1}\ X^{\beta-1-1}\ dX=t^{2\beta-2}\
\textbf{B}(\beta-1,\beta)=\frac{t^{2\beta-2}\ \Gamma(\beta-1)\
\Gamma(\beta)}{\Gamma(2\beta-1)}$ $\displaystyle\therefore R.H.S.~{}of~{}(6)$
$\displaystyle=$ $\displaystyle-\frac{\Gamma(\beta)}{\tau\
\Gamma(2\beta-1)}\left[t^{2\beta-2}\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)+t^{2\beta-2}\otimes_{t}\mathcal{U}(x,t)\otimes_{x}f(x)\otimes_{x}f(x)\right]$
$\displaystyle=$ $\displaystyle-\frac{\Gamma(\beta)}{\tau\
\Gamma(2\beta-1)}\left[\int\limits_{0}^{x}\int\limits_{0}^{t}dx_{1}\ dt_{1}\
t_{1}^{2\beta-2}\ \mathcal{U}(x_{1},t-t_{1})\ f(x-x_{1})\right.$
$\displaystyle\left.+\int\limits_{0}^{x}\int\limits_{0}^{x}\int\limits_{0}^{t}dx_{2}\
dx_{1}\ dt_{1}\ t_{1}^{2\beta-2}\ \mathcal{U}(x_{2},t_{1}-t)\ f(x_{1})\
f(x-x_{2}-x_{1})\right]$

Hence, the generalized integro-differential equation of the GKFE in equation
(5) is proven.

##  6 Conclusion

In conclusion, we have found the self-consistent solutions of the GKFE for
four cases i.e., when the memory functions are Mittag-Leffler function,
derived negative exponential, and ML function’s approximations at short time
and long time. We have found that the GKFE solutions are the same for negative
exponential and for upto frst order expansion of stretched exponential
function for very small $\tau\rightarrow 0$. The cross-over from short time to
long time approximations of the ML function in GKFE is presented in this
paper. We have also found a generalized integro-differential equation of the
GKFE and have expressed the same for an asymptotic case. We conclude by saying
that modelling complex dynamics of natural systems with non-Markovian GKFE can
make the complex systems study more realistic in nature by considering past
event histories. One important area of our study can be in modelling short
term and long term memory in neuroscience. We would like to add that the GKFE
with memory function is a very important formalism which can give valuable
insights about stochastic volatility or extreme events in many complex systems
study such as complex biological systems (like disease mechanisms in
epilepsy); chemical systems such as important information processing
properties of complex molecular events driven at short time and at a longer
time; financial system such as stock market and also social systems such as
outbreak of an epidemic.

Authors Contribution
The present work is conceptualized by RKBS and ALC. ALC and JB produced
analytical results. The manuscript is written by ALC and RKBS. All authors
discussed, finalised and approved the presented manuscript.

Competing financial interests:
The authors declare no competing financial interests.

## References

  * (1) Kolmogoroff, A. (1931) Ueber die analytischen Methoden in der Wahrscheinlichkeitsrechnung. Mathematische Annalen. Vol. 104. No.1. pp.415-458.
  * (2) Feller, W. (1940) On the Integro-Differential Equations of Purely Discontinuous Markoff Processes. Transactions of the American Mathematical Society. Vol. 48. No. 3. pp.488-515.
  * (3) Feller, W. (1957) On Boundaries and Lateral Conditions for the Kolmogorov Differential Equations. Ann. Maths. 65\. pp.527-570 .
  * (4) Rudenko, O. V., Dubkov, A. A., Gurbatov, S. N. (2016) On Exact Solutions to the Kolmogorov-Feller Equation. Doklady Mathematics. Vol. 94. No. 1. pp.476-479.
  * (5) Denisov, S.I., Bystrik, Yu.S. (2019) Exact stationary solutions of the Kolmogorov-Feller equation in a bounded domain. Communications in Nonlinear Science and Numerical Simulation. Volume 74. Pp.248-259
  * (6) Petrovskaya,I.V. (1970) Purely Discontinuous Random Processes in an Irregular Force Field. I. Soviet Astronomy Vol. 13. No.4. pp.647
  * (7) Petrovskaya,I.V. (1970) Purely Discontinuous Random Processes in an Irregular Force Field. II. Velocity Variation of a Star. Soviet Astronomy Vol. 13. No.4. pp.957
  * (8) Petrovskaya, I.V. (1992) The evolution of groups of very light stars in the galaxy. Celestial Mechanics and Dynamical Astronomy. Vol. 54. pp.267-269
  * (9) Petrovskaya, I.V. (1996) Evolution of the cometary cloud under the action of stellar perturbations. Earth, Moon, and Planets Vol. 72. pp.31-34.
  * (10) Tarasov, V.E., Zaslavsky, G.M. (2008) Fokker-Planck equation with fractional coordinate derivatives. Physica A. Vol. 387. pp.6505-6512
  * (11) Saichev, A.I., Zaslavsky G.M. (1998) Fractional kinetic Equations: solutions and applications. Chaos Vol.7. No. 4. pp.753-764.
  * (12) Blackledge,J.,Lamphiere,M.,Murphy,K.,Overton,S.,Panahi,A. (2012): Stochastic Volatility Analysis using the Generalised Kolmogorov-Feller Equation. The 2012 International Conference of Financial Engineering, London, 2012. doi:10.21427/D7R637
  * (13) Blackledge, J., Lamphiere, M., Murphy, K., Overton S. (2013) Financial Forecasting Using the Kolmogorov–Feller Equation. In: Yang GC., Ao S., Gelman L. (eds) IAENG Transactions on Engineering Technologies. Lecture Notes in Electrical Engineering, vol 229. Springer, Dordrecht.
  * (14) Mittag-Leffler, G.M. (1903) Sur la nouvelle fonction $E_{\alpha}$(x) C.R.Acad. Sci. Paris. Vol. 137. pp.554-558
  * (15) Erdelyi, A. (1955). Higher Transcendental Functions. McGraw-Hill, New York. Vol. 3
  * (16) Haubold,H.J., Mathai,A.M., Saxena,R.K. Mittag-Leffler Functions and Their Applications. (2011) Journal of Applied Mathematics. Vol. 2011. Article ID 298628. 51 pages.
  * (17) Mainardi,F.,Raberto,M.,Gorenflo,R.,Scalas,E. (2000) Fractional Calculus and continuous-time Finance II: the waiting time distribution. Physica A. Vol. 287, 00. 468-481
  * (18) Gorenflo,R.,Mainardi,F.,Raberto,M.,Scalas,E. (2000) Fractional Diffusion in Finance: Basic Theory. Modelli Dinamici in Economia e Finanza, Italy, September 28-30,2000
  * (19) Mainardi, F. (2010) Fractional Calculus and Waves in Linear Viscoelasticity. Imperial College Press and World Scientific, Singapore
  * (20) Koeller, R.C. (2010) A theory relating creep and relaxation for linear materials with memory. J. Appl. Mech. Vol. 77. Issue 031008. pp. 1-9.
  * (21) Mainardi, F. (2012) An historical perspective on fractional calculus in linear Viscoelasticity. Fract. Calc. Appl. Anal. Vol. 15. Issue 4. pp. 712–717. E-print: http://arxiv.org/abs/1007.2959.45.
  * (22) Rogosin, S. (2015) The Role of the Mittag-Leffler Function in Fractional Modeling. Mathematics. Vol. 3. pp. 368-381. doi:10.3390/math3020368
  * (23) Pillai, R.N. (1990) On Mittag-Leffler Functions and Related Distributions. Ann. Inst. Statist. Math. Vol. 42. Issue No. I. pp.157-161
  * (24) Bianco S, Ignaccolo M, Rider MS, Ross MJ, Winsor P, Grigolini P. (2007) Brain, music, and non-Poisson renewal processes. Phys. Rev. E. Stat. Nonlin. Soft Matter Phys. Vol 75. Issue 061911.
  * (25) David Papo. (2014) Functional significance of complex fluctuations in brain activity: from resting state to cognitive neuroscience. Front. Syst. Neurosci. Vol.11 https://doi.org/10.3389/fnsys.2014.00112
  * (26) Jiejie Chen, Zhigang Zeng, Ping Jiang. (2014) Global Mittag-Leffler stability and synchronization of memristor-based fractional-order neural networks. Neural Networks Vol. 51. pp. 1-8 https://doi.org/10.1016/j.neunet.2013.11.016
  * (27) El-Nabulsi, R.A. (2016) Int J Theor Phys Vol.55. Issue 2. pp-625-635. https://doi.org/10.1007/s10773-015-2700-5
  * (28) Shchigolev, V.K. (2011) Cosmological Models with Fractional Derivatives and Fractional Action Functional Commun. Theor. Phys. Vol. 56. pp.389-396.
  * (29) A.Einstein.(1905) On the Motion of Small Particles Suspended in Liquids at Rest Required by the Molecular-Kinetic Theory of Heat. Annalen der Physik, Vol. 17, pp. 549-560.
  * (30) Markov, A. A. Extension of the limit theorems of probability theory to a sum of variables connected in a chain. Reprinted in Appendix B of R. Howard. Dynamic Probabilistic Systems. Volume 1: Markov Chains. John Wiley and Sons, New York, 1971.
  * (31) Blackledge,J.,Lamphiere,M.,S.,Panahi,A. (2012): Simulation and Analysis of Stochastic signals using the Kolmogorov-Feller Equation. ISSC 2012, NUI Maynooth, June 28-29.
  * (32) Feller, W. (1966). An Introduction to Probability Theory and Its Applications.Wiley, New York. Vol. II.