# Lights out: training RL agents robust to temporary blindness

Nathan Ordonez  Marije Tromp  Pau Marquez Julbe  Wendelin Böhmer
<EMAIL_ADDRESS>
###### Abstract

Agents trained with DQN rely on an observation at each timestep to decide what
action to take next. However, in real world applications observations can
change or be missing entirely. Examples of this could be a light bulb breaking
down, or the wallpaper in a certain room changing. While these situations
change the actual observation, the underlying optimal policy does not change.
Because of this we want our agent to continue taking actions until it receives
a (recognized) observation again. To achieve this we introduce a combination
of a neural network architecture that uses hidden representations of the
observations and a novel $n$-step loss function. Our implementation is able to
withstand location based blindness stretches longer than the ones it was
trained on, and therefore shows robustness to temporary blindness. For access
to our implementation, please email Nathan, Marije, or Pau.

## 1\. Introduction

While modern Reinforcement Learning (RL) agents have successfully learned to
complete many different tasks in multiple environments, their robustness to
superficial changes in state (changes that shouldn’t impact the task at hand)
is not a given when using state-of-the-art methods, even though this might be
important for safety-critical systems. The name of this problem, as described
by a recent survey on RL robustness, is Observation robustness[1]. Here,
disturbances in the input state are handled by modeling the uncertainty of the
environment.

But what should an agent do, when it enters an uncertain state? Specifically,
if a state observation has high uncertainty but the agent knows that the
underlying environment (specifically, the transition function) has not
changed, what should the agent do? A straight-forward to this solution could
be that, given the last certain state, an agent could simply remember from the
other times it was in that state, what the best sequence of actions should be,
and then execute them in an open-loop fashion (meaning, performing actions
without paying attention to the states). One could implement this by simply
looking at past trajectory recordings, however one would then miss the
advantages of using neural networks end-to-end, such as the possibility for
better generalization, and the applicability to more complex state spaces.
Implementing such an end-to-end neural network and a loss function specialized
for this reinforcement learning task is the ultimate goal of this project.

Previous works have explored the open-loop challenge in an offline RL setting.
One recent work proposed encoding the input state into a latent
representation, and using a recurrent neural network to make future
predictions of this latent representation[2]. However, the problem they solved
was about performing high-frequency actions with access to low-frequency
observations, rather than tackling observation robustness. Our approach, and
therefore our contributions, are as follows:

  1. 1.

We contribute a novel loss function based on DQN that successfully trains an
open-loop policy on recorded trajectories using a recurrent neural network, so
that it learns to be robust to open-loop situations during evaluations.

  2. 2.

To the best of our knowledge, we contribute the first open-loop loss that is
applied to the online deep reinforcement learning setting. In doing so, an
agent is able to learn an open-loop policy while exploring and learning to
perform a novel task (i.e. a new maze).

  3. 3.

Finally, experiments are shown to qualitatively understand the strengths,
weaknesses and relevant parameters for agents trained using our loss function
in a gridworld maze setting.

In this project we are interested in how well an agent performs when
navigating a maze while being blinded. This means that the environment
includes a mask, which defines the positions at which the agent does not
receive a state vector, and is therefore required to perform open loop
control. We then train an open-loop policy on a set number of blind steps,
which is a parameter to experiment with. Our research questions are:
Can a DQN agent learn a policy that is robust to temporarily not having access
to the state during evaluation by switching between closed- and open-loop
control?
This main question can be divided into multiple sub-questions:

  1. 1.

Can our loss function train both an open-loop and a closed-loop controller to
be robust to blind state masks without training on them?

  2. 2.

How does length/frequency of blindness during evaluation affect agent’s
ability to complete its task?

  3. 3.

How does the presence or absence of open-loop policy trajectories in the
replay buffer affect robustness to temporary blindness?

  4. 4.

How many successive steps in the future can we use to train an open-loop
policy robust to temporary blindness?

The paper is organized as follows: In Section 2 and 3 we will outline the
basic theory on which our architecture is based. Then, the proposed algorithm
and training procedure is detailed in Section 4, followed by an explanation of
the experiments in Section 5. Finally, our conclusions are presented in
Section 6.

## 2\. Background

First of all, we define a fully-observable MDP by the tuple
$\langle\mathcal{S},\mathcal{A},r,\mathcal{P},\mu,\gamma\rangle$, where
$\mathcal{S}$ denotes the state space, $\mathcal{A}$ the action space,
$r:\mathcal{A}\times\mathcal{S}\mapsto\mathbb{R}$ is the reward function,
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\mapsto[0,1]$ is the
state transition function, $\mu$ is the initial state distribution and
$\gamma\in[0,1]$ is the discount factor. Moreover, $s_{t}\in\mathcal{S}$ is
the state, $a_{t}\in\mathcal{A}$ is the action, and $r_{t}\in\mathbb{R}$ is
the reward at time step $t$.

### Deep Q-Networks

Deep Q-Networks (DQN), as introduced by [3], are a family of deep Q-Learning
algorithms where a value function is estimated for every state-action pair in
the system. The goal of reinforcement learning is to find a stationary policy
$\pi:\mathcal{S}\times\mathcal{A}\to[0,1]$ that maximizes the expected
discounted sum of future rewards:

$\begin{aligned}
J^{\pi}=\mathbb{E}\left[\sum_{t=0}^{\infty}\gamma^{t}r_{t}\Bigg{|}\\!\\!\begin{array}[]{c}\scriptstyle
s_{0}\sim\mu,\;a_{t}\sim\pi(\cdot|s_{t})\\\\[-2.84526pt] \scriptstyle
r_{t}=r(s_{t},a_{t})\\\\[-2.84526pt] \scriptstyle
s_{t+1}\sim\mathcal{P}(\cdot|s_{t},a_{t})\end{array}\\!\\!\right]\end{aligned}\,.$
(1)

Q-learning achieves this by estimating the Q-value of a greedy policy that
selects the action that maximizes that Q-value, yielding the Bellman
optimality equation [4]:

$\begin{aligned}
Q(s,a)=\mathbb{E}\big{[}r(s,a)+\gamma\,\operatorname*{max}_{a^{\prime}\in\mathcal{A}}Q(s^{\prime},a^{\prime})\big{|}\\!\\!\begin{array}[]{c}\scriptstyle
s^{\prime}\sim P(s,a)\end{array}\\!\\!\big{]}\end{aligned}\,,\qquad\forall
s\in\mathcal{S},\forall a\in\mathcal{A}\,.$ (2)

To deal with large state and action spaces, the Q-value function is
implemented as a neural network, and its loss function over transitions
$(s_{t},a_{t},r_{t},s_{t+1})$ experienced in the environment is defined as
follows:

$L(\theta)=\mathbb{E}\Big{[}\big{(}r_{t}+\gamma\max_{a^{\prime}\in\mathcal{A}}Q_{(s_{t+1},a^{\prime};\theta^{\prime})}-Q_{(s_{t},a_{t};\theta)}\big{)}^{2}\,\Big{]}$
(3)

Because the target Q-value function at $s_{t+1}$ evolves at the same speed as
the regular Q-value function, training curves can be quite unstable and
prevent the agent from converging towards an optimal policy. To avoid this, a
soft target update is used where the target network, with parameters
$\theta^{\prime}$, is a copy of the original Q-value function, which is
updated after every gradient upsate step using the rule
$\theta^{\prime}\leftarrow(1-\tau)\theta^{\prime}+\tau\theta$ where $\tau$ is
the soft update parameter. Further instability both early on during training
and after having converged to the optimal policy can happen if an agent only
learns from the last step it has taken. To mitigate this, a replay buffer can
be helpful. Replay buffers are filled with a set number of past recordings of
actions, states and rewards. During training, a batch is sampled from the
replay buffer and used to perform gradient descent steps based on the Q-value
update function.
In this project, both soft target updates and a replay buffer were used as
they improved dramatically the performance of our agent.

### Closed and Open-Loop

The version of DQN described in the previous subsection uses what we call a
closed-loop controller. Closed-loop controllers receive an observation in
every state they are in and take actions based on these observations.

A different type of controller is an open-loop controller. These controllers
receive an observation and then take multiple actions based on that
observation.

A situation in which it is useful to train an open-loop controller is when the
state-observation distribution can change, but the optimal action sequence to
solve a task does not. An example of this could be the wallpaper in a room
changing, or the lights in a hallway turning off. Neither of these will
influence the optimal policy, if its goal is to go from one side of the room
to the other, for example. A closed-loop controller would not generally be
expected to perform well in this scenario, as it is essentially an out of
distribution situation for its training, whereas an open-loop controller would
simply take actions based off of the last in-distribution observation it has
recorded until it reaches the next in-distribution observation. Once the agent
is in an in-distribution state, it can switch back to a closed-loop controller
and continue its path.

### Distribution Shift

Distribution shift is a phenomenon that occurs in offline and off-policy
learning when the policy applies actions that are out-of-distribution, i.e.
not seen during training. To solve this issue, different solutions exist [5],
[6]. This can be a potential issue when we only train the open-loop controller
on decisions made by the closed-loop controller. In our case, the amount of
open-loop transitions that are sampled during training can be controlled by a
parameter to check for this, and the difference in trajectory distributions is
indicated by the difference in performance when the agent is evaluated with
blindness masks.

## 3\. Related work

The task of dealing with a lack of state observations can be seen as that of
building an agent that can deal with closed and open-loop control. [7]
describes a Q-learning algorithm that works in an environment where sensing
comes at a cost by switching between open-loop and closed-loop, and is the
earliest mention of open-loop RL we could find. In our problem sensing does
not come at a cost, instead not being available at times. This paper also does
not use neural networks.

Open-loop control in deep RL was explored by [2], where the agent has access
to a state input only every $n$ steps, thus requiring the agent to make action
decisions on the basis of the potentially sparse state inputs. Instead of
receiving input every $n$ steps our agent will encounter location based
blindness. We therefore test how robust it is to not receiving state inputs.

The concept of encoding the state vector and computing predictions using this
embedding has contributed to new state of the art in multiple RL tasks with
MuZero [8]. Here, the authors use embeddings to build a search tree of
possible strategies to play games, evaluating their value and play
accordingly. In contrast to this, we plan to exploit previous experience to
repeat a known desirable previous behavior and increase agent robustness.
Additionally, they do not employ the MuZero algorithm on open-loop control
tasks.

In [2], state embeddings are also used to make predictions. The goal there is
that the hidden embeddings lead to actions with similar value to actions that
would have been taken if the agent had access to the ground truth state input.
The task is also different in terms of the frequency at which the agent
performs open-loop control. In their work, the agent is mostly open-loop and
regularly gets access to the state, however in our case blindness is not
assumed to be regular and the agent is expected to be merely robust to
occasional loss of access to state observations.

## 4\. Methodology

Instead of implementing our DQN agent with an open-loop controller from
scratch, we used the DQN implementation from the existing codebase cleanrl [9]
as a starting point. To transition this implementation into a controller that
can handle both closed and open-loop we changed the neural network model
architecture, the replay buffer, and loss function. Since our replay buffer is
a standard trajectory buffer, we will only explain the architecture and loss
functions in the following subsections.

### Architecture

Figure 1: A single step in used architecture.

In a closed-loop setting, the controller outputs an action $a_{t}$ for each
state $s_{t}$. Concretely, with the DQN architecture, it predicts the Q-value
of each state-action pair and selects the optimal one. However, when open-loop
control starts, the agent only has access to the first state, and must make an
arbitrary number of actions on the basis of this one state. The architecture
we used to learn the open-loop controller is inspired by MuZero [8]. For each
open-loop timestep, instead of taking the actual state as the input, our
architecture takes a hidden representation of said state. This hidden
representation is not directly regularized, meaning that the neural network
learns by itself what information to store in this embedding. This also means
that this is a model-based algorithm, in the sense that the neural network can
approximate a transition function of the environment.

Figure 2: Architecture in closed-loop mode.

Figure 3: Architecture in open-loop mode.

To get the initial hidden state $h_{0}$, we use an encoder function
$g_{\theta}$, which takes as input the state $s_{t}$ and gives as output the
hidden state $h_{0}$. We then use this hidden state together with a sequence
of actions as input for a function $f_{\theta}$, which is an Long Short-Term
Memory (LSTM) [10]. Using an LSTM allows us to give an action sequence input
of variable length, therefore allowing us to use the same neural network for
both the closed- and open-loop controller, based on Q-values:

$Q(s,a_{0:{n-1}})\;=\;\sum_{k=0}^{n-1}\gamma^{k}\underbrace{r_{\theta}(h_{k+1})}_{r_{k}}+\,\gamma^{n}\underbrace{v_{\theta}(h_{n})}_{v_{n}}\,,\qquad
h_{0}\;=\;g_{\theta}(s)\,,\qquad h_{k+1}\;=\;f_{\theta}(h_{k},a_{k})\,.$ (4)

The LSTM $f_{\theta}$ then returns a new hidden state, which represents the
network’s interpretation of what the (hidden) state $h_{k+1}$ looks like after
taking the action $a_{k}$ from the given (hidden) state $h_{k}$. From this new
hidden state we use a reward function $r_{\theta}$ and a value function
$v_{\theta}$ to predict the value $v_{k+1}$ of the new hidden state and the
reward $r_{k}$ of taking the action in the previous hidden state. A single
step of this can be seen in Figure 1. Then, based on whether or not there is a
new state available, we can use $g$ to transform that state and use it as
input for the next time step, or we can use the output from $f$ at the
previous timestep as the hidden state input.

In Figure 3, we show our architecture applied to the closed-loop setting. If
we recursively apply function $f$ to the hidden state, we get an open-loop
controller as in Figure 3.

### Loss function

Because of the potential distribution shift when training the open-loop
controller on closed-loop trajectories, we can train both policies on
trajectories drawn from both policies. This is done for a certain n-step
target loss $n$ during training by becoming blind with probability $p/n$,
where $p$ is a parameter to control how often the agent is blind during
training. When the agent becomes blind, it then stays blind for $n$ steps.
Then, when performing gradient descent we sample $N$-step trajectories
$(s_{t},a_{t},r_{t})_{t=\tau}^{\tau+N}$ from the replay buffer of past
episodes and apply the following loss to each one:

$\displaystyle\mathcal{L}_{val}(\theta)$
$\displaystyle\;=\;\sum_{n=1}^{N}\mathbb{E}\biggr{[}\Big{(}v_{\theta}(h_{n};s_{t},a_{0:n-1})-y(s_{t+n})\Big{)}^{2}\biggr{]}\,,$
(5) $\displaystyle y(s)$
$\displaystyle\;=\;\max_{a^{\prime}\in\mathcal{A}}\Bigl{(}r_{\theta^{\prime}}\big{(}f_{\theta^{\prime}}(g_{\theta^{\prime}}(s),a^{\prime})\big{)}+\gamma
v_{\theta^{\prime}}\big{(}f_{\theta^{\prime}}(g_{\theta^{\prime}}(s),a^{\prime})\big{)}\Bigl{)}\,.$

The goal of this loss function is that value estimations after $n$ open-loop
steps need to be as close as possible to closed-loop value estimations. For
each time step $t+n$, the optimal action in time step $t+n$ is used as target
$y(s_{t+n})$. This target is based on single-timestep embeddings from the
$f_{\theta^{\prime}}$ function, using slowly changing target parameters
$\theta^{\prime}$, whereas the value estimation function $v_{\theta}$ on the
left uses $n$-timestep embeddings with differentiated parameters $\theta$,
i.e. function $f_{\theta}$ is applied $n$ times to the state embedding
$g_{\theta}(s_{t})$ using the actions from the trajectory buffer that were
performed at those timesteps.
To ground the reward head of the network, we compute the following loss, which
enforces the predictions of the reward head $r_{\theta}$ to represent the
environment rewards for each time step in the trajectory:

$\mathcal{L}_{rew}(\theta)\;=\;\sum_{n=1}^{N}\mathbb{E}\biggr{[}(r_{\theta}(h_{n};s_{t})-r_{t+n-1})^{2}\biggr{]}\,.$
(6)

The complete loss is:
$\mathcal{L}(\theta)=\mathcal{L}_{val}(\theta)+\mathcal{L}_{rew}(\theta)\,.$

Figure 4: Mazes used in the experiments. In light green the agent’s starting
point, in white the goal. The left maze is intended to test different kinds of
challenges, while the right maze is intended to measure performance for an
increasingly long blind area.

## 5\. Experiments and Results

### Experimental Setup

To test our controller in different ways, we designed 2 maze environments, as
shown in Figure 4, and multiple masks for each. In both mazes, the agent
spawns on the light green square, and its goal is to arrive at the white
square. To better understand the strengths and weaknesses of our approach, the
left maze in Figure 4 was designed to contain three kinds of challenges; the
first part is called the open room, where there are many different paths to
the exit. The second is the zigzag path, chosen because early experimentation
indicated it presented a tough challenge for agents. The third is the forks,
where making a wrong turn could bring you far from the end goal, to punish bad
individual decisions. The second maze was designed to empirically check how
long of a blindness mask our agents can still solve. In both mazes, when the
agent finds the goal it receives a reward of $+1$, when it walks into a wall,
and thus stays in the same spot, it receives a reward of $-0.02$, and when it
takes a valid step it receives a reward of $-0.01$.

Figure 5: Mazes used in the experiments. The three masks are overlapped in
purple, cyan and brown, which mask the open room, zigzag and forks parts,
respectively. The light green square on the top left is the starting point of
the agent, and the white square is the goal.

Most importantly, the masks used to determine whether the agent should use
open-loop control or not in the first maze, is indicated in Figure 5. The left
maze shows the masks used together to serve as a benchmark for the agent, so
that it faces all three challenges while being blinded. Results for this are
in Figure 6, which shows the performance of the agent for different values of
N. The masks on the right were used individually in order to understand which
challenges are most difficult for our agents, and the results are shown in the
appendix, Figure 9. For the maze on the left of Figure 5, all of the masks
were combined to evaluate whether the agents are able to switch between closed
and open-loop control. These masks are slightly smaller to facilitate this
switch.

Since we are training a policy to perform on a single maze with a single goal
the trained policy would fail to generalize to different mazes or goals. For
this reason, the observation given to the policy is the one-hot encoded
position of the agent, which is the only dynamic component of the environment.
During both training and testing the actions of the agent are restricted to
the four vertical and horizontal movements.

To test this project’s proposed algorithm, we tried different configurations
of hyperparameters on different experimental settings.

### Hyperparameters

During our experiments, after finding hyperparameters that seemed to give good
performance, we chose two as being primarily relevant to test in this project.

The first is parameter $p$, which influences the probability that an open-loop
trajectory of length $N$ is recorded into the trajectory buffer during
training, which is set to $p/N$ to adjust for different values of N. The
second is parameter $N$, which determines the length of the trajectories
sampled from the replay buffer during gradient descent, and the maximum length
of the $n$-step target losses applied for one neural network update.

### Experiments

To answer our research questions, we tested the performance of our algorithm
in 3 different settings, which use different shape blindness masks. In the
areas covered by the masks, the agent will use the open-loop policy to move
towards the goal in the optimal amount of steps.

The first experimental setting serves as an overall benchmark and answers our
first sub-question: Can our loss function train both an open-loop and a
closed-loop controller to be robust to blind state masks without training on
them? For this experiment, we use the left maze in Figure 4, combined with all
the masks shown in the left maze in Figure 5. Since none of these masks
overlap, the controller is forced to switch between both open and closed loop
control.

The second experimental setting answers our second and third sub-questions:
How does length/frequency of blindness during evaluation affect loss and
reward during testing? and How does the presence or absence of open-loop
policy trajectories in the replay buffer affect robustness to temporary
blindness?. It also attempts to answer our fourth sub-question: How many
successive steps in the future can we use to train an open-loop policy robust
to temporary blindness?. In this experiment we use the right maze in Figure 4.
For multiple values of $N$ of our proposed loss, we find out the largest
number of steps for which it can use the open-loop controller reliably. We do
this by starting with the first block of the optimal path being masked, then
the first 2, the first 3, etc. until all blocks on the optimal path between
the starting position and the goal are covered, or until the agent can no
longer reliably reach the goal.

The third experimental setting takes a closer look at the sub-question Can our
loss function train both an open-loop and a closed-loop controller to be
robust to blind state masks without training on them? by using the left maze
in Figure 4, combined with the masks shown in the right maze in Figure 5. We
report the relationship between the relevant hyperparameters and each
individual mask, to better understand the kinds of challenges that our agents
struggle most with. This can be found in the appendix.

### Results and Discussion

#### Switching between closed and open-loop

Figure 6: Average episode length on left masks of Figure 5 over $3$ seeds with
different $N$ parameter during training, showing that higher N values allow
the agent to solve the maze despite the blindness masks.

In this experiment we trained agents with values of $N$ from 1 to 12 for
$60.000$ global steps in the left maze in Figure 4 with $p$ set to $0.5$. The
optimal path for this maze is $34$ steps. We evaluate whether the controller
is able to switch between closed and open-loop using the left maze in Figure 5
with epsilon set to $0.05$. The results of this experiment are shown in Figure
6.

As we can see in the graph, agents that have been trained for $N<4$ are not
able to find the goal within the given timesteps, which might be due to the
shortest masked path being 6 steps long.

What is perhaps more interesting, is that while all of the others were able to
reach the goal at some point, only the agent with $N=7$ was able to do so in
(close to) the optimal number of steps. The longest optimal path through any
of the masks takes 8 steps. This means that the agent can still generalize at
least somewhat to states that it has not been trained on.

Overall, the Figure shows that the controller is able to successfully switch
between open and closed-loop to solve the masked maze, in some cases even
performing optimally. However, since we only used 3 different draws to compute
the statistics, the results are not statistically reliable enough.

#### Maximum blind length

In this experiment we trained agents with values 1-25 for $N$ for 45000 global
steps in the zigzag maze in Figure 4 with $p$ $0.0$ and $0.5$. The optimal
path for this maze is $40$ steps. As mentioned before, we iteratively test how
much blindness the agent can handle by adding a mask to each block on the
optimal path one by one. If the agent is able to finish the maze within 150
maximum steps it moves on to the next mask size. To evaluate the actual policy
we set epsilon to 0, otherwise the agent would be able to pass through the
next mask size on accident.

The results are shown in Figure 7, which shows the maximum length the agent
was able to pass through for each $N$, calculated from 3 different seeds. The
shaded areas are between the minimum blindness length and the maximum
blindness length that the agent was able to traverse in these 3 seeds.

As we can see in the Figure the controller is able to solve masks that are
many times the size of $N$. Between $N=4$ and $N=10$ the controller is even
able to solve all masked blocks for some of the seeds. This means that the
controller is able to predict the reward and value relatively accurately, even
when it has been blind for much longer than it was trained on.

Figure 7: Maximum length of the blindness mask that an agent was able to
solve, when trained on different values of N, averaged over 3 seeds. It shows
that our method trains an agent robust to blindness lengths multiple times
bigger than the $N$ parameter it is trained with, even without sampling any
blind trajectories (p=0). Figure 8: Mean episode length across $3$ different
seeds on the zigzag maze with no masks. This figure shows how non-blind
performance is affected when varying the $N$ parameter, and explains that the
drop-off in performance for higher $N$ values on Figure 7 is not only because
of the blindness.

However, for $N$ larger than 10 the performance quickly degrades, and after
$N=16$ it is almost entirely unable to solve any blind steps. At first glance,
it seems like the agent simply loses its ability to be robust to blindness.

However, looking at Figure 8, we notice something else. This graph shows the
mean episode length across three seeds for each $N$ for the same zigzag maze
without any masks, where epsilon is also set to $0$. We can see that for N
larger than 10, the closed-loop controller was unable to solve the problem,
which explains why the open-loop controller was not able to solve it either.
For this reason, we cannot reach any conclusion about the open-loop controller
for $N>15$. A reason why this could be happening is that it needs more
training steps to converge, but this has not been tested for.

This, combined with the fact that we only use 3 different seeds to create
these results means that we cannot make any certain conclusions about exactly
how the solvable blind period scales with $N$. All we can say is that the
performance of the controller is highly dependent on $N$, and that it seems to
be able to generalize to blindness lengths far longer than what it was trained
with.

The last point of note for this experiment is the 2 different values of $p$.
As we can see in the graph, the performance of the open-loop policy does not
seem to statistically depend on $p$. One can also see how high the variance
is, and in general the results are difficult to interpret with only an average
over three seeds. Since $p$ is influencing the amount of open-loop
trajectories that we record during training, we can conclude that there is no
distribution shift in this dataset, i.e. the trajectories accessed by both
policies during training are similar. This fact could be due to the size of
the dataset, where random actions can lead to explore the whole state-space,
which would hide the distribution shift.

#### Mask evaluation

By evaluating the three different masks on the left maze one by one, as
depicted on the left of Figure 5, we are able to investigate which qualities
of the maze give more difficulties to the agent’s open-loop controller. Since
the experiments showed that hyperparameter $p$ did not have an impact on
performance on any of the masks, we will focus on the hyperparameter $N$. The
relationship shows that higher values of $N$ from 1 to 12 generally increase
the ability to deal with the three different masks. Then, comparing the
minimum $N$ that achieved the best performance for each mask, we can see that
the forks, room, and zigzag masks required $N$ values of 3, 4, and 9 to reach
optimal performance. This indicates that the zigzag pattern is especially
challenging for the open-loop controller to deal with, and that elongating the
$n$-step loss has a major effect on the kinds of challenges the agent is
robust to. This could be an indication of something we noticed when looking at
episode recordings, that the open-loop controller can get stuck repeating the
last action while being blind. More details can be seen in the appendix.

## 6\. Conclusions

In this paper we presented a method that uses a recurrent neural network and
an $n$-step loss function to learn an open-loop policy using Q-learning. We
have showed its high performance and robustness to different kind of mazes in
a gridworld by masking the state in different experiment settings. The open-
loop policy was successfully able to find the optimal path, even for a large
amount of blinded steps.

Overall it is important to keep in mind that the statistics found in this work
have high variance due to a small sample of trained neural networks with
different seeds.

We have shown that in the small environments we tested, while the trajectory
distribution shift between open and closed loop controllers is present as they
show differences in performance, training the open loop controller on open-
loop sampled trajectories did not have much of an effect on performance.
However, since our mazes have a relatively small state space, further
experiments in larger mazes should be done in future work to further explore
this.

## References

  * [1] Janosch Moos, Kay Hansel, Hany Abdulsamad, Svenja Stark, Debora Clever, and Jan Peters.  Robust reinforcement learning: A review of foundations and recent advances.  Machine Learning and Knowledge Extraction, 4:276–315, 03 2022.
  * [2] Alex Lewandowski, Vincent Zhang, and Dale Schuurmans.  Recurrent open-loop control in offline reinforcement learning.  In NeurIPS 2020 Offline Reinforcement Learning Workshop, 2020.
  * [3] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al.  Human-level control through deep reinforcement learning.  nature, 518(7540):529–533, 2015.
  * [4] Richard S. Sutton and Andrew G. Barto.  Reinforcement Learning: An Introduction.  The MIT Press, second edition, 2018.
  * [5] Seunghyun Lee, Younggyo Seo, Kimin Lee, Pieter Abbeel, and Jinwoo Shin.  Offline-to-Online Reinforcement Learning via Balanced Replay and Pessimistic Q-Ensemble.  arXiv e-prints, page arXiv:2107.00591, July 2021.
  * [6] Vivian Lin, Kuk Jin Jang, Souradeep Dutta, Michele Caprio, Oleg Sokolsky, and Insup Lee.  Take Me Home: Reversing Distribution Shifts using Reinforcement Learning.  arXiv e-prints, page arXiv:2302.10341, February 2023.
  * [7] Eric Hansen, Andrew Barto, and Shlomo Zilberstein.  Reinforcement learning for mixed open-loop and closed-loop control.  In M.C. Mozer, M. Jordan, and T. Petsche, editors, Advances in Neural Information Processing Systems, volume 9. MIT Press, 1996.
  * [8] Julian Schrittwieser, Ioannis Antonoglou, Thomas Hubert, Karen Simonyan, Laurent Sifre, Simon Schmitt, Arthur Guez, Edward Lockhart, Demis Hassabis, Thore Graepel, et al.  Mastering atari, go, chess and shogi by planning with a learned model.  Nature, 588(7839):604–609, 2020.
  * [9] Shengyi Huang, Rousslan Fernand Julien Dossa, Chang Ye, Jeff Braga, Dipam Chakraborty, Kinal Mehta, and João G.M. Araújo.  Cleanrl: High-quality single-file implementations of deep reinforcement learning algorithms.  Journal of Machine Learning Research, 23(274):1–18, 2022.
  * [10] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.

## 7\. Appendix

### A: performance on individual masks

Figure 9: For each value of n, we plot the average and standard deviation of
the lowest episode lengths of 9 agents.

Performance on the three masks experimentally significantly depended on
parameter n, but not on parameter p, so we plotted the per-mask lowest episode
length averaged across the nine values of p that were evaluated between 0 and
0.9. This shows performance for the three qualitatively different challenges
of the room, zigzag and fork masks. We can clearly see that the room mask is
the easiest to solve with the closed loop (n=1) policy, but that the forks
mask can be solved with a lower value of n. The zigzag mask is the one
requiring the largest value of n to solve. 150 is the maximum number of steps
and indicates that the agent could not reach the goal on time.


# Fractional non-homogeneous counting process

Nick Laskin E-mail address<EMAIL_ADDRESS>
(TopQuark Inc.
Toronto, ON, M6P 2P2
Canada)

###### Abstract

A new fractional non-homogeneous counting process has been introduced and
developed using the Kilbas and Saigo three-parameter generalization of the
Mittag-Leffler function. The probability distribution function of this process
reproduces for certain set of the fractality parameters the famous Poisson and
fractional Poisson probability distributions as well as the probability
distribution function of a counting process displaying stretched exponential
interarrival times distribution.

Applications of the developed fractional non-homogeneous counting process
cover fractional compound process, the generalization of combinatorial
polynomials and numbers, statistics of the fractional non-homogeneous counting
process and a new representation of the Kilbas and Saigo function.

PACS numbers: 05.10.Gg; 05.45.Df; 42.50.-p.

Keywords: Fractional non-homogeneous counting process, fractional Poisson
process, generating functions, combinatorial polynomials and numbers, Mittag-
Leffler function, Kilbas and Saigo function

##  1 Introduction

Recently, a discrete probability distribution function

$P_{\sigma}(n,x)=\frac{x^{\sigma n}}{n!}\exp(-x^{\sigma}),\quad x\geq 0,\quad
0<\sigma\leq 1,\quad n=0,1,2,...,$ (1)

has been implemented to introduce and study stretched coherent states [1].
Substitution $x\rightarrow\lambda_{\sigma}^{1/\sigma}t$ leads to the
probability distribution function of the following counting process

$P_{\sigma}(n,t)=\frac{(\lambda_{\sigma}t^{\sigma})^{n}}{n!}\exp(-\lambda_{\sigma}t^{\sigma}),$
(2)

where $\lambda_{\sigma}$ has units of $[\lambda_{\sigma}]=time^{-\sigma}.$
This expression is recognized as a non-homogeneous Poisson process whose
“rate” $r(t)$ changes with time,

$r(t)=\sigma\lambda_{\sigma}t^{\sigma-1}.$ (3)

Then $P_{\sigma}(n,t)$ can be represented as follows,

$P_{\sigma}(n,t)=\frac{(\Lambda_{\sigma}(t))^{n}}{n!}\exp(-\Lambda_{\sigma}(t)),$
(4)

where

$\Lambda_{\sigma}(t)=\int\limits_{0}^{t}d\tau
r(\tau)=\sigma\lambda_{\sigma}\int\limits_{0}^{t}d\tau\tau^{\sigma-1}=\lambda_{\sigma}t^{\sigma}.$
(5)

Given the probability distribution $P_{\sigma}(n,t)$, it is easy to obtain the
probability distribution function of interarrival times, i.e., the time
intervals between two consecutive arrivals. It turns out, that for the
counting process described by the probability distribution function specified
by Eq.(2), interarrival times probability distribution function
$\psi_{\sigma}(\tau)$ has the form

$\psi_{\sigma}(\tau)=\sigma\lambda_{\sigma}\tau^{\sigma-1}\exp(-\lambda_{\sigma}\tau^{\sigma}),\qquad\tau>0,$
(6)

which is the well-known Weibull distribution with the shape parameter $\sigma$
and the scale parameter $\lambda_{\sigma}$. The probability distribution
function $\psi_{\sigma}(\tau)$ is widely used to fit “fat-tail” economic,
social and natural data and is called the stretched exponential distribution
(see, for example, [2]).

The function $P_{\sigma}(n,t)$ satisfies the following system of the equations

$\frac{\partial P_{\sigma}(n,t)}{\partial
t}=\sigma\lambda_{\sigma}t^{\sigma-1}(P_{\sigma}(n-1,t)-P_{\sigma}(n,t)),\qquad
n\geq 1,$ (7)

$\frac{\partial P_{\sigma}(0,t)}{\partial
t}=-\sigma\lambda_{\sigma}t^{\sigma-1}P_{\sigma}(0,t),$ (8)

with some initial condition at $t=0$,

$P_{\sigma}(n,t=0)=p_{0}(n).$ (9)

Usually, it is supposed that at $t=0$ nothing had arrived, that is the initial
condition is specified as $P_{\sigma}(n,t)|_{t=0}=p_{0}(n)$, where
$p_{0}(n)=0,\ \mathrm{if}\ n\neq 0$ and $p_{0}(n)=1,\ \mathrm{if}\ n=0$.

Thus, Eq.(7) is linear differential equation with respect to $t$ and
difference equation with respect to $n$, and is called a differential-
difference equation. When $\sigma=1$, we get the famous Poisson probability
distribution function $P_{1}(n,t)\equiv P_{\sigma}(n,t)|_{\sigma=1}$, which
satisfies the system of equations [3]

$\frac{\partial P_{1}(n,t)}{\partial
t}=\lambda_{1}(P_{1}(n-1,t)-P_{1}(n,t)),\qquad n\geq 1,$ (10)

$\frac{\partial P_{\sigma}(0,t)}{\partial t}=-\lambda_{1}P_{\sigma}(0,t),$
(11)

and has the form

$P_{1}(n,t)=\frac{(\lambda_{1}t)^{n}}{n!}\exp(-\lambda_{1}t),\qquad
n=0,1,2,...,.$ (12)

where $\lambda_{1}$ has units of $[\lambda_{1}]=time^{-1}.$ As can be seen
from Eq.(6), when $\sigma=1$, $\psi_{\sigma}(\tau)|_{\sigma=1}=$
$\psi_{1}(\tau)$ becomes an exponential distribution,

$\psi_{1}(\tau)=\lambda_{1}\exp(-\lambda_{1}\tau),\qquad\tau>0.$ (13)

Therefore, $P_{\sigma}(n,t)$ defined by Eq.(2) is generalization of the
Poisson probability distribution function $P_{1}(n,t)$ defined by Eq.(12) and
$\psi_{\sigma}(\tau)$ is generalization of the exponential distribution
defined by Eq.(13).

Another important and widely used generalization of the Poisson process is
fractional Poisson process developed to model and study a long-memory impact
on the counting processes [4] empirically observed in complex quantum and
classical systems. The quantum system example is the fluorescence
intermittency of single CdSe quantum dots, that is, the fluorescence emission
of single nanocrystals exhibits intermittent behavior, namely, a sequence of
”light on” and ”light off” states departing from Poisson statistics. In fact,
the waiting time distribution in both states displays non-exponential behavior
[5]. As examples of classical systems let’s mention the distribution of
waiting times between two consecutive transactions in financial markets [6]
and another, which comes from network communication systems, where the
duration of network sessions or connections exhibits non-exponential power-law
behavior [7]. The probability distribution function of the fractional Poisson
process satisfies the following system of fractional differential-difference
equations

${}_{0}^{C}D_{t}^{\mu}P_{\mu}(n,t)=\lambda_{\mu}\left(P_{\mu}(n-1,t)-P_{\mu}(n,t)\right),\qquad
n\geq 1,$ (14)

and

${}_{0}^{C}D_{t}^{\mu}P_{\mu}(0,t)=-\lambda_{\mu}P_{\mu}(0,t),\qquad 0<\mu\leq
1,$ (15)

subject to normalization condition

$\sum\limits_{n=0}^{\infty}P_{\mu}(n,t)=1,$ (16)

and the initial condition at $t=0$, $P_{\mu}(n,t=0)=\delta_{n,0},$ where $\mu$
is fractality parameter, the operator ${}_{0}^{C}D_{t}^{\mu}$ is the Caputo
fractional derivative111The basic formulas on fractional calculus can be found
in Refs. [8] \- [10]. whose action on a function $\mathrm{f}(t)$ is defined as
[11]

${}_{0}^{C}D_{t}^{\mu}\mathrm{f}(t)=\frac{1}{\Gamma(1-\mu)}\int\limits_{0}^{t}\frac{d\tau(d\mathrm{f}(\tau)/d\tau)}{(t-\tau)^{\mu}},\qquad
0<\mu\leq 1,$ (17)

with $\Gamma(\mu)$ being the Gamma function, which has the familiar
representation $\Gamma(\mu)=\int\limits_{0}^{\infty}dte^{-t}t^{\mu-1}$,
$\mathrm{Re}\mu>0$, and the parameter $\lambda_{\mu}$ has units
$[\lambda_{\mu}]=time^{-\mu}$.

The solution of the system of the equations (14) and (15) has the form

$P_{\mu}(n,t)=\frac{(\lambda_{\mu}t^{\mu})^{n}}{n!}\sum\limits_{k=0}^{\infty}\frac{(k+n)!}{k!}\frac{(-\lambda_{\mu}t^{\mu})^{k}}{\Gamma(\mu(k+n)+1))},\quad
n=0,1,2,...,\quad 0<\mu\leq 1,$ (18)

which is probability distribution function of the fractional Poisson process
[4]. The $P_{\mu}(n,t)$ gives us the probability that in the time interval
$[0,t]$ we observe $n$ arrivals.

Equations (14) and (15) are fractional generalization of Eqs.(7) and (8) for
probability distribution function of the Poisson process. When $\mu=1$,
Eq.(18) turns into Eq.(12).

The goal of this work is to develop a system of fractional non-homogeneous
differential-difference equations the exact solution of which is new
fractional probability distribution of non-homogeneous fractional counting
process. It turns out that this system reproduces for certain set of the
fractality parameters the systems of equations (7) and (8), (10) and (11),
(14) and (15). This goal was achieved by implementing the Kilbas and Saigo
function [12], which is a three-parameter generalization of the Mittag-Leffler
function [13].

The paper is organized as follows.

In Section 1 a non-homogeneous fractional differential-difference equation for
the counting probability distribution functions has been developed and
discussed. To solve this equation the probability generating function was
calculated in Section 2. The fractional counting probability distribution
function has been obtained in Section 3 using the probability generating
function. Section 4 is devoted to the calculation of the moment generating
function. The mean and variance of the new fractional counting process were
found using the moment generating function. Interarrival time distribution was
obtained in Section 5. Section 6 presents applications of new fractional
counting probability distribution function covering fractional compound
process and the generalization of combinatorial polynomials and numbers.

The results obtained are summarized in the Conclusion.

The Appendix A lists special cases of Kilbas and Saigo three-parameter
function.

The Appendix B presents analytical expressions for the first four central
moments of the introduced fractional non-homogeneous probability distribution.

##  2 Fundamentals

We define a new fractional non-homogeneous counting process by introducing the
following system of fractional differential-difference equations for
probability distribution function $P_{\mu,\beta}(n,t)$ of this process

${}_{0}D_{t}^{\mu}P_{\mu,\beta}(n,t)=\lambda_{\mu+\beta}t^{\beta}\left(P_{\mu,\beta}(n-1,t)-P_{\mu,\beta}(n,t)\right),\quad
n=1,2,3,...,\quad t\geq 0,$ (19)

${}_{0}D_{t}^{\mu}P_{\mu,\beta}(0,t)=-\lambda_{\mu+\beta}t^{\beta}P_{\mu,\beta}(0,t),$
(20)

subject to normalization condition

$\sum\limits_{n=0}^{\infty}P_{\mu,\beta}(n,t)=1,$ (21)

and the initial condition

$P_{\mu,\beta}(n,t=0)=\delta_{n,0}.$ (22)

Here, $\mu$ and $\beta$ are dimensionless fractality parameters subject to the
conditions,

$0<\mu\leq 1,\qquad-\mu<\beta\leq 1-\mu,$ (23)

the rate of arrivals $\lambda_{\mu+\beta}$ has units $[\lambda_{\mu+\beta}]=$
$time^{-(\mu+\beta)}$, ${}_{0}D_{t}^{\mu}$ is the operator of the Caputo
fractional time derivative of order $\mu$ defined by Eq.(17). We will see
later where the conditions for $\mu$ and $\beta$ come from.

The system of Eqs.(19) and (22) introduces the counting process with
probability distribution function $P_{\mu,\beta}(n,t)$ of arriving $n$ items
($n=0,1,2,...$) by time $t,$ when the rate of arrivals depends on time $t$. In
other words, the above defined counting process is generalization to non-
homogeneous case (due to the dependency of the arrival rate on time) of the
fractional Poisson process originally introduced and developed in [4].

When $\mu=1$ and $\beta=0$, the system of Eqs.(19) and (22) turns into the
system of equations (10) and (11), which define the Poisson process with rate
of arrival $\lambda_{1}$.

When $\beta=0$, the system of Eqs.(19) and (22) turns into the system of
equations (14) and (15), which define the fractional Poisson process with rate
of arrival $\lambda_{\mu}$ where $0<\mu\leq 1$.

When $\mu=1$, the system of Eqs.(19) and (22) turns into the the system of
equations (7) and (8), which define the probability distribution function
$P_{\sigma}(n,t)$, with $\sigma=1+\beta$, $0<\sigma\leq 1$.

##  3 Probability generating function

To solve the system of equations (19) - (22) it is convenient to use the
method of the generating function. Introducing the probability generating
function $G_{\mu,\beta}(s,t)$

$G_{\mu,\beta}(s,t)=\sum\limits_{n=0}^{\infty}s^{n}P_{\mu,\beta}(n,t),$ (24)

allows to calculate $P_{\mu,\beta}(n,t)$ as

$P_{\mu,\beta}(n,t)=\frac{1}{n!}\frac{\partial^{n}G_{\mu,\beta}(s,t)}{\partial
s^{n}}|_{s=0}.$ (25)

Then, by multiplying Eqs.(19) and (20) by $s^{n}$, summing over $n$, we obtain
the following fractional differential-difference equation for the probability
generating function $G_{\mu}(s,t)$

${}_{0}D_{t}^{\mu}G_{\mu}(s,t)=\lambda_{\mu+\beta}t^{\beta}\left(\sum\limits_{n=1}^{\infty}s^{n}P_{\mu}(n-1,t)-\sum\limits_{n=0}^{\infty}s^{n}P_{\mu}(n,t)\right)=$
(26)

$\lambda_{\mu+\beta}t^{\beta}(s-1)G_{\mu}(s,t).$

Thus, we come to the following equation

${}_{0}D_{t}^{\mu}G_{\mu,\beta}(s,t)=\lambda_{\mu+\beta}t^{\beta}(s-1)G_{\mu,\beta}(s,t),$
(27)

subject to the initial condition

$G_{\mu,\beta}(s,t=0)=1,$ (28)

due to the definition Eq.(24) and the initial condition for
$P_{\mu,\beta}(n,t)$ given by Eq.(22). The solution to Eq.(27) can be obtained
using the ansatz [12]

$G_{\mu,\beta}(s,t)=\sum\limits_{n=0}^{\infty}\frac{c_{n}(\mu,\beta)\\{\lambda_{\mu+\beta}t^{\mu+\beta}(s-1)\\}^{n}}{\Gamma(n(\mu+\beta)+1)},$
(29)

where $c_{n}(\mu,\beta)$ are still unknown coefficients and $\Gamma(\mu)$ is
the Gamma function. Substituting this ansatz into Eq.(27) and taking into
account that due to the definition of Caputo fractional derivative Eq.(17)

${}_{0}D_{t}^{\mu}t^{n(\mu+\beta)}=\frac{\Gamma(n(\mu+\beta)+1)}{\Gamma(n(\mu+\beta)+1-\mu)}t^{n(\mu+\beta)},$
(30)

we find that Eq.(27) is satisfied if the following recurrence relation holds,

$c_{n+1}(\mu,\beta)=\frac{\Gamma(n(\mu+\beta)+\beta+1)}{\Gamma(n(\mu+\beta)+1)}c_{n}(\mu,\beta).$
(31)

This recurrence relation allows us to represent Eq.(29) as

$G_{\mu,\beta}(s,t)=c_{0}(\mu,\beta)+$ (32)

$\sum\limits_{n=1}^{\infty}\\{\lambda_{\mu+\beta}t^{\mu+\beta}(s-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma(k(\mu+\beta)+\beta+1)}{\Gamma(k(\mu+\beta)+\mu+\beta+1)}c_{0}(\mu,\beta).$

Taking into account the initial condition given by Eq.(28), we come to the
conclusion that $c_{0}(\mu,\beta)=1$. Therefore, the solution of the initial
problem defined by Eqs.(27) and (28) is

$G_{\mu,\beta}(s,t)=1+\sum\limits_{n=1}^{\infty}\\{\lambda_{\mu+\beta}t^{\mu+\beta}(s-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma(k(\mu+\beta)+\beta+1)}{\Gamma(k(\mu+\beta)+\mu+\beta+1)}.$
(33)

It turns out that the probability generating function $G_{\mu,\beta}(s,t)$ can
be expressed in terms of Kilbas and Saigo three-parameter generalization of
the Mittag-Leffler function [12], defined as

$E_{\alpha,m,l}(z)=\sum\limits_{n=0}^{\infty}c_{n}z^{n},\qquad
z\in\mathbb{C},$ (34)

where

$c_{n}=\prod\limits_{j=0}^{n-1}\frac{\Gamma(\alpha(jm+l)+1)}{\Gamma(\alpha(jm+l)+1)+1)},\qquad
n=1,2,...,\qquad c_{0}=1,$ (35)

and

$\alpha,m,l\in\mathbb{R},\quad\alpha>0,\quad
m>0\quad\mathrm{and}\quad\alpha(jm+l)\neq-1,-2,-3,....$ (36)

To show how the Kilbas and Saigo function appears in Eq.(33) we transform the
arguments of the Gamma functions under the product sign in Eq.(33) as follows

$\prod\limits_{k=0}^{n-1}\frac{\Gamma(k(\mu+\beta)+\beta+1)}{\Gamma(k(\mu+\beta)+\mu+\beta+1)}=\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$
(37)

Hence, it can be seen from Eqs.(33) and (37) that the power series holds for
$G_{\mu,\beta}(s,t)$

$G_{\mu,\beta}(s,t)=1+\sum\limits_{n=1}^{\infty}\\{\lambda_{\mu+\beta}t^{\mu+\beta}(s-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$
(38)

On the other hand, using the definitions given by Eqs.(34) and (35), we see
that the right side of Eq.(38) is the following Kilbas and Saigo
function222The special cases of the function
$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)$ are presented in Appendix
A.,

$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)=1+$ (39)

$\sum\limits_{n=1}^{\infty}z^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$

with $z=\lambda_{\mu+\beta}t^{\mu+\beta}(s-1)$. Therefore, the generating
function $G_{\mu,\beta}(s,t)$ has the form

$G_{\mu,\beta}(s,t)=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(\lambda_{\mu+\beta}t^{\mu+\beta}(s-1)),$
(40)

where $\mu$ and $\beta$ are subject to the conditions given by Eq.(23). Thus,
we expressed in terms of Kilbas and Saigo three-parameter function the
probability generating function $G_{\mu,\beta}(s,t)$ of the counting process
introduced by the system of fractional differential-difference equations (19)
and (20).

When $\mu=1$ and $\beta=0$, then according to Eq.(134) the generating function
$G_{1,0}(s,t)$ takes the form

$G_{1,0}(s,t)=\exp(\lambda_{1}t(s-1)),$

which is the generating function of the Poisson probability distribution [3].

When $0<\mu\leq 1$ and $\beta=0$, then due to Eq.(136) we obtain from Eq.(38)

$G_{\mu,0}(s,t)=1+\sum\limits_{n=1}^{\infty}\frac{\\{\lambda_{\mu}t^{\mu}(s-1)\\}^{n}}{\Gamma(\mu
n+1)}=E_{\mu}(\lambda_{\mu}t^{\mu}(s-1)),$ (41)

where $E_{\mu}(z)$ is the Mittag-Leffler function defined as [13]

$E_{\mu}(z)=\sum\limits_{n=0}^{\infty}\frac{z^{n}}{\Gamma(\mu n+1)}.$ (42)

Thus, Eq.(41) gives us the probability generating function of the fractional
Poisson process [4].

When $\mu=1$ and $-1<\beta\leq 0$, then due to Eq.(138) and having introduced
the notations,

$\sigma=1+\beta\qquad\mathrm{and}\qquad\overline{\lambda}_{\sigma}=\frac{\lambda_{\sigma}}{\sigma},$
(43)

we present Eq.(38) in the form

$G_{\sigma}(s,t)=\exp(\overline{\lambda}_{\sigma}t^{\sigma}(s-1)),$ (44)

where $0<\sigma\leq 1$, if $-1<\beta\leq 0$. Thus, Eq.(44) is probability
generating function of the probability distribution function given by Eq.(2).

When $\mu+\beta=1$, then taking into account Eq.(139) we can express Eq.(38)
in the form

$G_{\mu,1-\mu}(s,t)=1+\sum\limits_{n=1}^{\infty}\\{\lambda
t(s-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma(k+2-\mu)}{\Gamma(k+2)},$

or

$G_{\mu,1-\mu}(s,t)=E_{\mu,\frac{1}{\mu},\frac{1-\mu}{\mu}}(\lambda t(s-1)),$

where the following notation was introduced

$\lambda=\lambda_{\mu+\beta}\left|{}_{\mu+\beta=1}\right..$ (45)

###  3.1 Probability distribution function

Using the definition given by Eq.(25) we obtain the probability distribution
function $P_{\mu,\beta}(n,t)$

$P_{\mu,\beta}(n,t)=\frac{1}{n!}\frac{\partial^{n}G_{\mu,\beta}(s,t)}{\partial
s^{n}}|_{s=0}=\frac{(\lambda_{\mu+\beta}t^{\mu+\beta})^{n}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(-\lambda_{\mu+\beta}t^{\mu+\beta}),$
(46)

where the following notation is introduced

$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(z)=\frac{d^{n}}{dz^{n}}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z).$
(47)

The equation (46) can be rewritten as

$P_{\mu,\beta}(n,t)=\frac{(-z)^{n}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(z)\left|{}_{z=-\lambda_{\mu+\beta}t^{\mu+\beta}}\right.,$
(48)

and

$P_{\mu,\beta}(n=0,t)=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(-\lambda_{\mu+\beta}t^{\mu+\beta}).$
(49)

The function $P_{\mu,\beta}(n,t)$ is normalized

$\sum\limits_{n=0}^{\infty}P_{\mu,\beta}(n,t)=\frac{(-z)^{n}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(z)|_{z=-\lambda_{\mu+\beta}t^{\mu+\beta}}=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(0)=1,$
(50)

due to the definition given by Eqs.(34) and (35), and positive due to the
complete monotonicity of the three-parameter Kilbas and Saigo function
$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(-x)$ when $\mu$ and $\beta$ are
subject to the conditions333The conditions  $0<\mu\leq 1,\qquad-\mu<\beta,$
(51) follow directly from the condition of complete monotonicity of the three-
parameter Kilbas and Saigo function [14]. We have added an additional
condition  $\beta\leq 1-\mu,$ to ensure consistency with the parameter
$\sigma$ defined by Eq.(43) and included in Eq.(1). defined by Eq.(23), see,
[14], [15],

$(-1)^{n}\frac{d^{n}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(-x)}{dx^{n}}\geq
0,\qquad x\geq 0,$ (52)

Therefore, $P_{\mu,\beta}(n,t)$ is indeed a probability distribution function.
The power series representation for $P_{\mu,\beta}(n,t)$ is

$P_{\mu,\beta}(n,t)=\frac{(\lambda_{\mu+\beta}t^{\mu+\beta})^{n}}{n!}\times$
(53)

$\sum\limits_{m=0}^{\infty}\frac{(m+n)!}{m!}(-\lambda_{\mu+\beta}t^{\mu+\beta})^{m}\prod\limits_{k=0}^{m+n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$

Replacing $\lambda_{\mu+\beta}t^{\mu+\beta}$ with $x^{\mu+\beta}$ in Eqs.(48)
and (53), we obtain a new family of fractional probability distributions
$P_{\mu,\beta}(x)$,

$P_{\mu,\beta}(n,x)=\frac{(-z)^{n}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(z)|_{z=-\lambda_{\mu+\beta}t^{\mu+\beta}},$
(54)

or in a series form

$P_{\mu,\beta}(n,x)=\frac{x^{(\mu+\beta)n}}{n!}\sum\limits_{m=0}^{\infty}\frac{(m+n)!}{m!}\times$
(55)

$(-x^{(\mu+\beta)})^{m}\prod\limits_{k=0}^{m+n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$

where $\mu$ and $\beta$ are subject to the conditions given by Eq.(23).

When $\mu=1$ and $\beta=0$, using Eq.(134) we get from Eq.(46)

$P_{1,0}(n,t)=\frac{(\lambda_{1}t)^{n}}{n!}\exp(-\lambda_{1}t),$

which is the celebrated Poisson probability distribution.

When $0<\mu\leq 1$ and $\beta=0$, using Eq.(136) we get from Eq.(46)

$P_{\mu,0}(n,t)=\frac{(\lambda_{\mu}t^{\mu})^{n}}{n!}\sum\limits_{m=0}^{\infty}\frac{(m+n)!}{m!}\frac{(-\lambda_{\mu}t^{\mu})^{m}}{\Gamma(\mu
m+1)},$

which is the fractional Poisson probability distribution [4].

When $\mu=1$ and $-1<\beta\leq 0$, then according to Eq.(138) we obtain from
Eq.(46)

$P_{1,\beta}(n,t)\equiv
P_{1+\beta}(n,t)=\frac{(\frac{\lambda_{1+\beta}}{1+\beta}t^{1+\beta})^{n}}{n!}\exp(-\frac{\lambda_{1+\beta}t^{1+\beta}}{1+\beta}).$
(56)

Using the notations given by Eq.(43), we represent Eq.(56) in the form

$P_{\sigma}(n,t)=\frac{(\overline{\lambda}_{\sigma}t^{\sigma})^{n}}{n!}\exp(-\overline{\lambda}_{\sigma}t^{\sigma}),\quad
0<\sigma\leq 1,$

which is recognized as Eq.(1).

With $\mu+\beta=1$ and using Eq.(139) we can present Eq.(46) as

$P_{\mu,1-\mu}(n,t)=\frac{(\lambda
t)^{n}}{n!}E_{\mu,\frac{1}{\mu},\frac{1-\mu}{\mu}}^{(n)}(-\lambda t),$

where $\lambda$ is defined by Eq.(45). The series representation for
$P_{\mu,1-\mu}(n,t)$ has the form

$P_{\mu,1-\mu}(n,t)=\frac{(\lambda
t)^{n}}{n!}\sum\limits_{m=0}^{\infty}\frac{(m+n)!}{m!}(-\lambda
t)^{m}\prod\limits_{k=0}^{m+n-1}\frac{\Gamma(k+2-\mu)}{\Gamma(k+2)}.$ (57)

##  4 Moment generating function

An equation for the moment of any integer order $m$ of the probability
distribution function $P_{\mu,\beta}(n,t)$ introduced by Eq. (48) can be
easily found using the moment generating function $M_{\mu,\beta}(s,t)$ defined
as

$M_{\mu,\beta}(s,t)=\sum\limits_{n=0}^{\infty}e^{-sn}P_{\mu,\beta}(n,t).$ (58)

Then for the moment of $m^{\mathrm{th}}$ order we have

$<n_{\mu,\beta}^{m}>=\sum\limits_{n=0}^{\infty}n^{m}P_{\mu,\beta}(n,t)=(-1)^{m}\frac{\partial^{m}}{\partial
s^{m}}M_{\mu,\beta}(s,t)\left|{}_{s=0}\right..$ (59)

Multiplying Eq. (48) by $e^{-sn}$ and summing over the $n$ we get

$M_{\mu,\beta}(s,t)=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(\lambda_{\mu+\beta}t^{\mu+\beta}(e^{-s}-1)),$
(60)

or in a series form

$M_{\mu,\beta}(s,t)=1+\sum\limits_{n=1}^{\infty}\\{\lambda_{\mu+\beta}t^{\mu+\beta}(e^{-s}-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$
(61)

Replacing $\lambda_{\mu+\beta}t^{\mu+\beta}$ with $x^{\mu+\beta}$ in Eqs.(60)
and (61) we get the moment generating function for the fractional generalized
probability distributions $P_{\mu,\beta}(x)$,

$M_{\mu,\beta}(s,x)=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(x^{\mu+\beta}(e^{-s}-1)),$
(62)

and in series form

$M_{\mu,\beta}(s,x)=1+\sum\limits_{n=1}^{\infty}\\{x^{\mu+\beta}(e^{-s}-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$
(63)

When $\mu=1$ and $\beta=0$, using Eq.(134) we get from Eq.(63)

$M_{1,0}(s,x)=1+\sum\limits_{n=1}^{\infty}\frac{(\lambda_{1}t(e^{-s}-1))^{n}}{n!}=\exp(\lambda_{1}t(e^{-s}-1)),$

which is the moment generating function of the Poisson process.

With $0<\mu\leq 1$ and $\beta=0$, using Eq.(136) we get from Eq.(63)

$M_{\mu,0}(s,x)=E_{\mu}(\lambda_{\mu}t^{\mu}(e^{-s}-1)),$

which is the moment generating function of the fractional Poisson process [4].

For $\mu=1$ and $-1<\beta\leq 0$, Eq.(63) gives us the moment generating
function of the counting process, the probability distribution function of
which is determined by Eq.(2)

$M_{\sigma}(s,x)\equiv
M_{1+\beta}(s,t)=\exp(\overline{\lambda}_{\sigma}t^{\sigma}(e^{-s}-1)),$ (64)

where $\sigma$ and $\overline{\lambda}_{\sigma}$ are defined by Eq.(43).

When $\mu+\beta=1$, Eq.(62) goes to

$M_{\mu,1-\mu}(s,t)=E_{\mu,\frac{1}{\mu},\frac{1-\mu}{\mu}}(\lambda
t(e^{-s}-1),$

and Eq.(63) goes to

$M_{\mu,1-\mu}(s,t)=1+\sum\limits_{n=1}^{\infty}\\{\lambda
t(e^{-s}-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma(k+2-\mu)}{\Gamma(k+2)},$

where $\lambda$ is defined by Eq.(45).

###  4.1 Mean and variance

The mean $<n_{\mu,\beta}>$ can be calculated using the definition (59),

$<n_{\mu,\beta}>=-\frac{\partial}{\partial
s}M_{\mu,\beta}(s,t)\left|{}_{s=0}\right.,$ (65)

with $M_{\mu,\beta}(s,t)$ given either by Eq.(60) or Eq.(61). Thus, we have

$<n_{\mu,\beta}>=\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\lambda_{\mu+\beta}t^{\mu+\beta}.$
(66)

It is easy to see that the moment of the second order $<n_{\mu,\beta}^{2}>$ is

$<n_{\mu,\beta}^{2}>=\frac{\partial^{2}}{\partial
s^{2}}M_{\mu,\beta}(s,t)\left|{}_{s=0}\right.=$ (67)

$2\frac{\Gamma(\beta+1)\Gamma(\mu+2\beta+1)}{\Gamma(\mu+\beta+1)\Gamma(2\mu+2\beta+1)}(\lambda_{\mu+\beta}t^{\mu+\beta})^{2}+\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\lambda_{\mu+\beta}t^{\mu+\beta}.$

Then variance $Var_{\mu,\beta}$ of the probability distribution
$P_{\mu,\beta}(n,t)$ is

$Var_{\mu,\beta}=<n_{\mu,\beta}^{2}>-<n_{\mu,\beta}>^{2}=$ (68)

$\left\\{\frac{2\Gamma(\beta+1)\Gamma(\mu+2\beta+1)}{\Gamma(\mu+\beta+1)\Gamma(2\mu+2\beta+1)}-\left(\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\right)^{2}\right\\}(\lambda_{\mu+\beta}t^{\mu+\beta})^{2}+$

$\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\lambda_{\mu+\beta}t^{\mu+\beta},$

here $\mu$ and $\beta$ are subject to the conditions given by Eq.(23).

##  5 Interarrival times distribution

An important characteristic of any counting process is the time between
successive counts/arrivals - the interarrival time. Interarrival time is a
random variable having the probability distribution function
$\psi_{\mu,\beta}(\tau),$

$\psi_{\mu,\beta}(\tau)=-\frac{d}{d\tau}\mathcal{P}_{\mu,\beta}(\tau),$ (69)

where $\mathcal{P}_{\mu,\beta}(\tau)$ is the probability that a given
interarrival time is greater or equal to $\tau$ [4]

$\mathcal{P}_{\mu,\beta}(\tau)=1-\sum\limits_{n=1}^{\infty}P_{\mu,\beta}(n,\tau)=P_{\mu,\beta}(0,\tau),$
(70)

here $P_{\mu,\beta}(n,\tau)$ is defined by Eq.(53).

It follows from Eqs.(69) and (70) that the probability distribution function
of the interarrival times has the form

$\psi_{\mu,\beta}(\tau)=\frac{\mu+\beta}{\tau}\sum\limits_{l=0}^{\infty}(-1)^{l}(l+1)(\lambda_{\mu+\beta}\tau^{\mu+\beta})^{l+1}\prod\limits_{k=0}^{l}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$
(71)

The Laplace transform $\varphi_{\mu,\beta}(u)$ of the interarrival probability
distribution function is defined as

$\varphi_{\mu,\beta}(u)=\int\limits_{0}^{\infty}d\tau
e^{-u\tau}\psi_{\mu,\beta}(\tau),$ (72)

and has the form

$\varphi_{\mu,\beta}(u)=(\mu+\beta)\sum\limits_{l=0}^{\infty}(-1)^{l}(l+1)(\lambda_{\mu+\beta})^{l+1}\frac{\Gamma[(\mu+\beta)(l+1)]}{u^{(\mu+\beta)(l+1)}}\times$
(73)

$\prod\limits_{k=0}^{l}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$

It easy to see that when $\mu=1$ and $\beta=0$, Eq.(71) transforms into the
exponential distribution defined by Eq.(13). From Eq.(73) we obtain the
Laplace transform of the exponential distribution

$\varphi_{1,0}(u)=\frac{\lambda_{1}}{\lambda_{1}+u}.$

When $0<\mu\leq 1$ and $\beta=0$, Eq.(71) becomes the probability distribution
of interarrival times of the fractional Poisson process

$\psi_{\mu,0}(\tau)=\lambda_{\mu}\tau^{\mu-1}E_{\mu,\mu}(-\lambda_{\mu}\tau^{\mu}),$
(74)

where $E_{\mu,\mu}(z)$ is the two parameter Mittag-Leffler function defined by
the series [16]

$E_{\mu,\mu}(z)=\sum\limits_{n=0}^{\infty}\frac{z^{n}}{\Gamma(\mu n+\mu)}.$
(75)

From and Eq.(73) it follows that the Laplace transform has the form [4]

$\varphi_{\mu,0}(u)=\frac{\lambda_{\mu}}{\lambda_{\mu}+u^{\mu}}.$

When $\mu=1$ and $-1<\beta\leq 0$, Eq.(71) goes to stretched exponential
distribution defined by Eq.(6) if we take into account the notation introduced
by Eq.(43). From Eq.(73) we obtain the Laplace transform of the stretched
exponential distribution

$\varphi_{\sigma}(u)=\sigma\overline{\lambda}_{\sigma}\sum\limits_{l=0}^{\infty}\frac{(-\overline{\lambda}_{\sigma})^{l}}{l!}\frac{\Gamma(\sigma(l+1))}{u^{\sigma(l+1)}},$
(76)

here $\sigma$ and $\overline{\lambda}_{\sigma}$ are defined by Eq.(43).

With $\mu+\beta=1$ Eq.(71) goes to

$\psi_{\mu,1-\mu}(\tau)=\frac{1}{\tau}\sum\limits_{l=0}^{\infty}(-1)^{l}(l+1)(\lambda\tau)^{l+1}\prod\limits_{k=0}^{l}\frac{\Gamma(k+2-\mu)}{\Gamma(k+2)},$
(77)

and Eq.(73) goes to

$\varphi_{\mu,1-\mu}(u)=\frac{\lambda}{u}\sum\limits_{l=0}^{\infty}(-\frac{\lambda}{u})^{l}(l+1)!\prod\limits_{k=0}^{l}\frac{\Gamma(k+2-\mu)}{\Gamma(k+2)},$
(78)

where $\lambda$ is defined by Eq.(45).

##  6 Applications

###  6.1 Fractional compound process

We call stochastic process $\\{X_{\mu,\beta}(t)$, $t\geq 0\\}$ a fractional
compound process if it is represented by

$X_{\mu,\beta}(t)=\sum\limits_{i=1}^{N_{\mu,\beta}(t)}Y_{i},$ (79)

where $\\{N_{\mu,\beta}(t)$, $t\geq 0\\}$ is a fractional counting process
whose probability distribution function is introduced by Eq.(48) or Eq.(53),
and $\\{Y_{i}$, $i=1,2,...\\}$ are independent and identically distributed
random variables with probability distribution function $p(Y_{i})$ for each
$Y_{i}$. It is assumed that the process $\\{N_{\mu,\beta}(t)$, $t\geq 0\\}$
and the sequence of random variables $\\{Y_{i}$, $i=1,2,...\\}$ are
independent.

The moment generating function $J_{\mu,\beta}(s,t)$ of fractional generalized
compound process is defined as follows

$J_{\mu,\beta}(s,t)=<\exp\\{sX_{\mu,\beta}(t)\\}>_{Y,N_{\mu,\beta}(t)},$ (80)

where $<...>_{Y_{i},N_{\mu,\beta}(t)}$ stands for average which includes two
statistically independent averaging procedures:

1\. Averaging $<...>_{Y}$ over independent random variables $Y_{i}$, defined
by

$<...>_{Y}=\int
dY_{1}...dY_{n}p(Y_{1})...p(Y_{n})...=\prod\limits_{i=1}^{n}\int
dY_{i}p(Y_{i})....$ (81)

2\. Averaging over random number of counts defined by

$<...>_{N_{\mu,\beta}(t)}=\sum\limits_{n=0}^{\infty}P_{\mu,\beta}(n,t)...,$
(82)

where $P_{\mu,\beta}(n,t)$ is given by Eq.(48).

One can see from Eq.(80) that $k^{\mathrm{th}}$ order moment of fractional
compound process $X_{\mu,\beta}(t)$ is obtained by differentiating
$J_{\mu,\beta}(s,t)$ $k$ times with respect to $s$ and then setting $s=0$,
that is

$<(X_{\mu,\beta}(t))^{k}>_{Y,N_{\mu,\beta}(t)}=\frac{\partial^{k}}{\partial
s^{k}}J_{\mu,\beta}(s,t)_{|s=0}.$ (83)

To obtain equation for the moment generating function $J_{\mu,\beta}(s,t)$ we
apply Eqs.(81) and (82) to Eq.(80)

$J_{\mu,\beta}(s,t)=\sum\limits_{n=0}^{\infty}<\exp\\{sX_{\mu,\beta}(t)|N_{\mu,\beta}(t)=n\\}>_{Y}P_{\mu,\beta}(n,t)=$
(84)

$\sum\limits_{n=0}^{\infty}<\exp\\{s(Y_{1}+...+Y_{n})\\}>_{Y}\times\frac{(-z)^{n}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(z)\left|{}_{z=-\lambda_{\mu+\beta}t^{\mu+\beta}}\right.,$

where we used the independence of $\\{Y_{1},Y_{2},...\\}$ and
$\\{N_{\mu,\beta}(t)$, $t\geq 0\\}$ and the independence of the $Y_{i}$’s
between themselves. Hence, setting

$g(s)=\int dY_{i}p(Y_{i})\exp\\{sY_{i}\\},$ (85)

for the moment generating function of random variables $Y_{i}$, we obtain from
Eq.(84) the moment generating function $J_{\mu,\beta}(s,t)$ of the fractional
compound process

$J_{\mu,\beta}(s,t)=\sum\limits_{n=0}^{\infty}g^{n}(s)\times\frac{(-z)^{n}}{n!}\frac{d^{n}}{dz^{n}}\left(E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)\right)|_{z=-\lambda_{\mu+\beta}t^{\mu+\beta}}=$
(86)
$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(\lambda_{\mu+\beta}t^{\mu+\beta}(g(s)-1)).$

By differentiating the above with respect to $s$ and setting $s=0$, it is easy
to obtain the mean of the fractional compound process,

$<X_{\mu,\beta}(t)>_{Y_{i},N_{\mu,\beta}(t)}=\frac{\partial}{\partial
s}J_{\mu,\beta}(s,t)_{|s=0}=<Y>_{Y}\lambda_{\mu+\beta}t^{\mu+\beta}\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)},$
(87)

which is a manifestation of independency between fractional counting process
$\\{N_{\mu,\beta}(t)$, $t\geq 0\\}$ and random variables $Y_{i}.$

###  6.2 New polynomials and associated numbers

Based on the probability distribution function defined by Eq.(54) we introduce
polynomials $B_{\mu,\beta}(x,m)$,

$B_{\mu,\beta}(x,m)=\sum\limits_{n=0}^{\infty}n^{m}P_{\mu,\beta}(n,x)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(-x),$
(88)

$B_{\mu,\beta}(x,0)=1,$

where the parameters $\mu$, and $\beta$ satisfy the conditions given by
Eq.(23). From Eq.(88) at $x=1$ we come to fractional numbers
$B_{\mu,\beta}(m)$,

$B_{\mu,\beta}(m)=B_{\mu,\beta}(x,m)|_{x=1}=\sum\limits_{n=0}^{\infty}\frac{n^{m}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(-1).$
(89)

As an example, here are the first two newly introduced polynomials

$B_{\mu,\beta}(x,1)=\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}x,$ (90)

and

$B_{\mu,\beta}(x,2)=\frac{\Gamma(\beta+1)\Gamma(\mu+2\beta+1)}{\Gamma(\mu+\beta+1)\Gamma(2\mu+2\beta+1)}x^{2}+\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}x.$
(91)

The first two fractional numbers associated with the polynomials
$B_{\mu,\beta}(x,m)$ are

$B_{\mu,\beta}(1)=\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)},$ (92)

and

$B_{\mu,\beta}(2)=\frac{\Gamma(\beta+1)\Gamma(\mu+2\beta+1)}{\Gamma(\mu+\beta+1)\Gamma(2\mu+2\beta+1)}+\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}.$
(93)

Series representations of the polynomials $B_{\mu,\beta}(x,m)$ and numbers
$B_{\mu,\beta}(m)$ are

$B_{\mu,\beta}(x,m)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}\sum\limits_{r=0}^{\infty}\frac{(r+n)!}{r!}(-x)^{r}\prod\limits_{k=0}^{r+n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$
(94)

$B_{\mu,\beta}(m)=B_{\mu,\beta}(x,m)|_{x=1}=\sum\limits_{n=0}^{\infty}\frac{n^{m}}{n!}\times$
(95)

$\sum\limits_{r=0}^{\infty}\frac{(r+n)!}{r!}\prod\limits_{k=0}^{r+n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$

where we used Eq.(55).

####  6.2.1 Generating function of polynomials $B_{\mu,\beta}(x,m)$

The generating function of the polynomials $B_{\mu,\beta}(x,m)$ is defined as

$\mathcal{B}_{\mu,\beta}(s,x)=\sum\limits_{m=0}^{\infty}\frac{s^{m}}{m!}B_{\mu,\beta}(x,m).$
(96)

To get the polynomial $B_{\mu,\beta}(x,m)$ we differentiate
$\mathcal{B}_{\mu,\beta}(s,x)$ $m$ times with respect to $s$, and then let
$s=0$. That is,

$B_{\mu,\beta}(x,m)=\frac{\partial^{m}}{\partial
s^{m}}\mathcal{B}_{\mu,\beta}(s,x)|_{s=0}.$ (97)

To find an explicit equation for $\mathcal{B}_{\mu,\beta}(s,x)$, let’s
substitute Eq.(88) into Eq.(96)

$\mathcal{B}_{\mu,\beta}(s,x)=\sum\limits_{m=0}^{\infty}\frac{s^{m}}{m!}\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(-x).$

Summation over $m$ gives us

$\sum\limits_{m=0}^{\infty}\frac{s^{m}}{m!}n^{m}=e^{sn},$

and we obtain

$\mathcal{B}_{\mu,\beta}(s,x)=\sum\limits_{n=0}^{\infty}\frac{x^{n}}{n!}e^{sn}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}^{(n)}(-x)=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(x(e^{s}-1)),$
(98)

where the function $E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)$ is
defined by Eq.(39). It is easy to see that the representation of the
generating function $\mathcal{B}_{\mu,\beta}(s,x)$ in a power series has the
form

$\mathcal{B}_{\mu,\beta}(s,x)=1+\sum\limits_{n=1}^{\infty}\\{x(e^{s}-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$
(99)

where $\mu$ and $\beta$ are subject to the conditions given by Eq.(23).

If we put $x=1$ in Eq.(98), then we immediately come to the generating
function $\mathcal{B}_{\mu,\beta}(s)$ of the numbers $B_{\mu,\beta}(m)$
introduced by Eq.(89)

$\mathcal{B}_{\mu,\beta}(s)=\mathcal{B}_{\mu,\beta}(s,x)|_{x=1}=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(e^{s}-1).$
(100)

To get the numbers $B_{\mu,\beta}(m)$ we differentiate
$\mathcal{B}_{\mu,\beta}(s)$ $m$ times with respect to $s$, and then set
$s=0$. That is,

$B_{\mu,\beta}(m)=\frac{\partial^{m}}{\partial
s^{m}}\mathcal{B}_{\mu,\beta}(s)|_{s=0}.$ (101)

When $\mu=1$ and $\beta=0$, using Eq.(134) we get from Eq.(88)

$B_{1,0}(x,m)=e^{\text{-}x}\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!},$

which are the Bell polynomials [17]. The new formula Eq.(89) is a
generalization of the so-called Dobiński relation known since 1877 [18] for
the Bell numbers444The number of ways a set of $m$ elements can be partitioned
into non-empty subsets.

$B_{{}_{1,0}}(m)=e^{\text{-}1}\sum\limits_{n=0}^{\infty}\frac{n^{m}}{n!}.$

When $0<\mu\leq 1$ and $\beta=0$, using Eq.(136) we get from Eq.(88)

$B_{\mu,0}(x,m)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}E_{\mu,1,0}^{(n)}(-x)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}E_{\mu}^{(n)}(-x),$

or in the form of a power series

$B_{\mu,0}(x,m)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}\sum\limits_{k=0}^{\infty}\frac{(n+k)!}{k!}\frac{(-x)^{k}}{\Gamma(\mu(n+k)+1)},\quad
B_{\mu,0}(x,0)=1,$

which are the fractional Bell polynomials [19]. Fractional generalization of
the Dobiński relation is [19]

$B_{\mu,0}(m)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}E_{\mu,1,0}^{(n)}(-1)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}E_{\mu}^{(n)}(-1).$

With $\mu+\beta=1$ and using Eq.(139) we can present Eq.(88) as

$B_{\mu,1-\mu}(x,m)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}E_{\mu,\frac{1}{\mu},\frac{1-\mu}{\mu}}^{(n)}(-x),$

where $\lambda$ is defined by Eq.(45). The series representation for
polynomials $B_{\mu,1-\mu}(x,m)$ has the form

$B_{\mu,1-\mu}(x,m)=\sum\limits_{n=0}^{\infty}n^{m}\frac{x^{n}}{n!}\sum\limits_{m=0}^{\infty}\frac{(m+n)!}{m!}(-\lambda
t)^{m}\prod\limits_{k=0}^{m+n-1}\frac{\Gamma(k+2-\mu)}{\Gamma(k+2)}.$ (102)

###  6.3 Fractional combinatorial numbers associated with polynomials
$B_{\mu,\beta}(x,m)$

Let us introduce fractional combinatorial numbers $S_{\mu,\beta}(m,l)$
associated with the polynomials $B_{\mu,\beta}(x,m)$ using the following
equations

$B_{\mu,\beta}(x,m)=\sum\limits_{l=0}^{m}S_{\mu,\beta}(m,l)x^{l},$ (103)

$S_{\mu,\beta}(m,0)=\delta_{m,0},\qquad S_{\mu,\beta}(m,l)=0,\ l=m+1,\ m+2,\
...,$

where $B_{\mu,\beta}(x,m)$ is the polynomials given by Eq.(88) and the
parameters $\mu$, and $\beta$ satisfy the conditions (23).

At $x=1$, when the polynomials $B_{\mu,\beta}(x,m)$ become the fractional
numbers, $B_{\mu,\beta}(m)=B_{\mu,\beta}(x,m)\left|{}_{x=1}\right.$, Eq.(103)
goes into the expression for the fractional numbers $B_{\mu,\beta}(m)$ in
terms of fractional combinatorial numbers $S_{\mu,\beta}(m,l)$

$B_{\mu,\beta}(m)=\sum\limits_{l=0}^{m}S_{\mu,\beta}(m,l).$ (104)

To find a generating function of the fractional combinatorial numbers
$S_{\mu,\beta}(m,l)$, let’s expand the generating function
$\mathcal{B}_{\mu,\beta}(s,x)$ given by Eq.(96). Upon substituting
$B_{\mu}(x,m)$ from Eq.(103) we have the following chain of transformations

$\mathcal{B}_{\mu,\beta}(s,x)=\sum\limits_{m=0}^{\infty}\frac{s^{m}}{m!}\left(\sum\limits_{l=0}^{m}S_{\mu,\beta}(m,l)x^{l}\right)=$

$\sum\limits_{m=0}^{\infty}\frac{s^{m}}{m!}\left(\sum\limits_{l=0}^{\infty}\theta(m-l)S_{\mu,\beta}(m,l)x^{l}\right)=\sum\limits_{l=0}^{\infty}\left(\sum\limits_{m=l}^{\infty}S_{\mu,\beta}(m,l)\frac{s^{m}}{m!}\right)x^{l},$
(105)

here $\theta(l)$ is the Heaviside step function, $\theta(l)=1,\ \mathrm{if\
}l\geq 0$ and $\theta(l)=0,\ \mathrm{if\ }l<0$.

On the other hand, from Eq.(98), we have for $\mathcal{B}_{\mu,\beta}(s,x)$

$\mathcal{B}_{\mu,\beta}(s,x)=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(x(e^{s}-1))=1+$
(106)

$\sum\limits_{n=1}^{\infty}\\{x(e^{s}-1)\\}^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$

Upon comparing Eq.(105) and Eq.(106), we conclude that

$\sum\limits_{m=l}^{\infty}S_{\mu,\beta}(m,l)\frac{s^{m}}{m!}=(e^{s}-1)^{l}\prod\limits_{k=0}^{l-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},\quad
l=1,2,....$ (107)

Let us introduce two generating functions $\mathcal{G}_{\mu,\beta}(s,l)$ and
$\mathcal{F}_{\mu,\beta}(s,t)$ of the fractional combinatorial numbers
$S_{\mu,\beta}(m,l)$,

$\mathcal{G}_{\mu,\beta}(s,l)=\sum\limits_{m=l}^{\infty}S_{\mu,\beta}(m,l)\frac{s^{m}}{m!}=(e^{s}-1)^{l}\prod\limits_{k=0}^{l-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$
(108)

and

$\mathcal{F}_{\mu,\beta}(s,t)=\sum\limits_{l=0}^{m}\sum\limits_{m=l}^{\infty}S_{\mu,\beta}(m,l)\frac{s^{m}t^{l}}{m!}=E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(t(e^{s}-1).$
(109)

Then the fractional combinatorial numbers are obtained as follows,

$S_{\mu,\beta}(m,l)=\frac{\partial^{m}}{\partial
s^{m}}\mathcal{G}_{\mu,\beta}(s,l)\left|{}_{s=0}\right.=$ (110)

$\prod\limits_{k=0}^{l-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}\frac{\partial^{m}}{\partial
s^{m}}(e^{s}-1)^{l}\left|{}_{s=0}\right.,$

and

$S_{\mu,\beta}(m,l)=\frac{1}{l!}\frac{\partial^{m+l}}{\partial s^{m}\partial
t^{l}}\mathcal{F}_{\mu,\beta}(s,t)\left|{}_{s=0},_{t=0}\right.=$ (111)

$\frac{1}{l!}\frac{\partial^{m+l}}{\partial s^{m}\partial
t^{l}}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(t(e^{s}-1))|_{s=0},_{t=0},$

where $l\leq m$ and $S_{\mu,\beta}(m,0)=\delta_{m,0}$.

When $\mu=1$ and $\beta=0$ it follows from Eq.(110) that

$S_{1,0}(m,l)=\frac{1}{l!}\frac{\partial^{m}}{\partial
s^{m}}(e^{s}-1)^{l}\left|{}_{s=0}\right..$ (112)

The right side of this equation is generating function of the standard
Stirling numbers of the second kind [20]. Hence, $S_{1,0}(m,l)$ is the
standard Stirling number of the second kind555The number of ways to partition
a set of $m$ elements into $l$ non-empty sets. $S(m,l)\equiv S_{1,0}(m,l)$,
and we conclude that

$S_{\mu,\beta}(m,l)=l!\prod\limits_{k=0}^{l-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}S(m,l).$
(113)

Thus, we have represented fractional combinatorial numbers
$S_{\mu,\beta}(m,l)$ through the the standard Stirling number of the second
kind $S(m,l)$. It is easy to see, for example, that

$S_{\mu,\beta}(m,1)=\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)},$

$S_{\mu,\beta}(m,2)=2(2^{m-1}-1)\prod\limits_{k=0}^{1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$

$S_{\mu,\beta}(m,m-1)=\frac{m!(m-1)}{2}\prod\limits_{k=0}^{m-2}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]},$

$S_{\mu,\beta}(m,m)=m!\prod\limits_{k=0}^{m-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$

In a particular case, for $\mu=1$ and $\beta=0$, equations (108) and (109)
turn into the well-know generating functions for standard Stirling numbers of
the second kind $S(m,l)$,

$\mathcal{G}_{1,0}(s,l)=\sum\limits_{m=l}^{\infty}S(m,l)\frac{s^{m}}{m!}=\frac{(e^{s}-1)^{l}}{l!},\qquad
l=0,1,2,...,$ (114)

and

$\mathcal{F}_{1,0}(s,t)=\sum\limits_{m=0}^{\infty}\sum\limits_{l=0}^{m}S(m,l)\frac{s^{m}t^{l}}{m!}=\exp(t(e^{s}-1)).$
(115)

(see, Eqs.(2.17) and (2.18) in [21]).

####  6.3.1 Statistics of fractional probability distribution

We now use the fractional combinatorial numbers introduced by Eq.(103) to
obtain the moments of the fractional non-homogeneous counting process with the
probability distribution function $P_{\mu,\beta}(n,t)$ given by Eqs.(48) and
(49). By the definition of the $m^{\mathrm{th}}$ order moment, we have

$<n_{\mu,\beta}^{m}>=\sum\limits_{n=0}^{\infty}n^{m}P_{\mu,\beta}(n,t)=\sum\limits_{n=0}^{\infty}n^{m}\frac{(-z)^{n}}{n!}\frac{d^{n}}{dz^{n}}\left(E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)\right)|_{z=-\lambda_{\mu+\beta}t^{\mu+\beta}}.$
(116)

It is easy to see that $<n_{\mu,\beta}^{m}>$ is the fractional polynomial
$B_{\mu,\beta}(\lambda_{\mu+\beta}t^{\mu+\beta},m)$ defined by Eq.(88),

$<n_{\mu,\beta}^{m}>=B_{\mu,\beta}(\lambda_{\mu+\beta}t^{\mu+\beta},m).$

Then using the definition Eq.(103) we obtain

$<n_{\mu,\beta}^{m}>=\sum\limits_{l=0}^{m}S_{\mu,\beta}(m,l)(\lambda_{\mu+\beta}t^{\mu+\beta})^{l},$
(117)

where $S_{\mu,\beta}(m,l)$ is fractional combinatorial numbers defined by
Eq.(113). Therefore, the fractional combinatorial numbers $S_{\mu,\beta}(m,l)$
naturally appear in the power series over $\lambda_{\mu+\beta}t^{\mu+\beta}$
for the $m^{\mathrm{th}}$ order moment of the fractional non-homogeneous
counting process introduced by Eqs.(19) \- (23).

Using the analytical expression given by Eq.(113), let’s list first four
moments of the fractional non-homogeneous probability distribution
$P_{\mu,\beta}(n,t)$. The first and second moments are given by Eqs.(66) and
(67) respectively. Third order moment $<n_{\mu,\beta}^{3}>$ is

$<n_{\mu,\beta}^{3}>=\sum\limits_{n=0}^{\infty}n^{3}P_{\mu,\beta}(n,t)=$ (118)

$6\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}\frac{\Gamma(2\mu+3\beta+1)}{\Gamma(3\mu+3\beta+1)}(\lambda_{\mu+\beta}t^{\mu+\beta})^{3}+$

$6\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}(\lambda_{\mu+\beta}t^{\mu+\beta})^{2}+\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\lambda_{\mu+\beta}t^{\mu+\beta},$

and fourth order moment $<n_{\mu,\beta}^{4}>$ is

$<n_{\mu,\beta}^{4}>=\sum\limits_{n=0}^{\infty}n^{4}P_{\mu,\beta}(n,t)=$

$24\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}\frac{\Gamma(2\mu+3\beta+1)}{\Gamma(3\mu+3\beta+1)}\frac{\Gamma(3\mu+4\beta+1)}{\Gamma(4\mu+4\beta+1)}(\lambda_{\mu+\beta}t^{\mu+\beta})^{4}+$

$36\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}\frac{\Gamma(2\mu+3\beta+1)}{\Gamma(3\mu+3\beta+1)}(\lambda_{\mu+\beta}t^{\mu+\beta})^{3}+$
(119)

$14\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}(\lambda_{\mu+\beta}t^{\mu+\beta})^{2}+\frac{\Gamma(\beta+1)}{\Gamma(\mu+\beta+1)}\lambda_{\mu+\beta}t^{\mu+\beta}.$

When $0<\mu\leq 1$ and $\beta=0$, Eqs.(66), (67), (118) and (119) turn into
equations for the first four moments of the fractional Poisson process [19].

In the case when $\mu=1$ and $\beta=0$, Eqs.(66), (67), (118) and (119) become
the well-know equations for moments of the standard Poisson probability
distribution with the parameter
$\lambda_{1}=\lambda_{\mu+\beta}|_{\mu=1,\beta=0}$ (see, Eqs.(22) - (24) in
Ref. [22]).

When $\mu=1$ and $-1<\beta\leq 0$ Eqs.(66), (67), (118) and (119) become the
equations for the moments of the process, defined by Eq.(1). Taking into
account Eq.(43) we come to the known equations for the moments of the standard
Poisson probability distribution with the parameter
$\overline{\lambda}_{\sigma}$.

###  6.4 Variance, skewness and kurtosis excess of the fractional non-
homogeneous counting process

To find analytical expressions for the variance, skewness and kurtosis excess
of a fractional non-homogeneous counting process determined by the probability
distribution function $P_{\mu,\beta}(n,t)$ given by Eq.(55), let’s introduce
the central moment of the $m^{\mathrm{th}}$ order $\mathrm{M}_{\mu,\beta}(m)$,

$\mathrm{M}_{\mu,\beta}(m)=<(n-<n_{\mu,\beta}>)^{m}>=\sum\limits_{n=0}^{\infty}(n-<n_{\mu,\beta}>)^{m}P_{\mu,\beta}(n,t)=$

$\sum\limits_{n=0}^{\infty}\sum\limits_{r=0}^{m}(-1)^{m-r}\binom{m}{r}n^{r}(<n_{\mu,\beta}>)^{m-r}P_{\mu,\beta}(n,t)=$
(120)

$\sum\limits_{r=0}^{m}(-1)^{m-r}\binom{m}{r}(<n_{\mu,\beta}>)^{m-r}\sum\limits_{l=0}^{r}S_{\mu,\beta}(r,l)(\lambda_{\mu+\beta}t^{\mu+\beta})^{l},$

where $<n_{\mu,\beta}>$ is the first order moment defined by Eq.(66) and
$S_{\mu,\beta}(r,l)$ is given by Eq.(113). The first four central moments
$\mathrm{M}_{\mu,\beta}(m)$, $m=1,2,3,4$, are presented in the Appendix B. In
terms of the central moments $\mathrm{M}_{\mu,\beta}(m)$, $m=1,2,3,4$, the
variance $Var_{\mu,\beta}$, skewness $s_{\mu,\beta}$, and kurtosis excess
$k_{\mu,\beta}$ of the fractional non-homogeneous counting process are

$Var_{\mu,\beta}=\mathrm{M}_{\mu,\beta}(2),$ (121)

${\large
s}_{\mu,\beta}=\frac{\mathrm{M}_{\mu,\beta}(3)}{\mathrm{M}_{\mu,\beta}^{3/2}(2)},$
(122)

${\large
k}_{\mu,\beta}=\frac{\mathrm{M}_{\mu,\beta}(4)}{\mathrm{M}_{\mu,\beta}^{2}(2)}-3.$
(123)

When $0<\mu<1$ and $\beta=0$, Eqs.(121) - (123) turn into equations for the
variance, skewness and kurtosis excess of the fractional Poisson process [19].
In the case where $\mu=1$ and $\beta=0$, Eqs.(121) - (123) become equations
for the variance, skewness and kurtosis excess of the standard Poisson
probability distribution, see, for example, Eqs.(25) - (27) in Ref. [22]. When
$\mu=1$ and $-1<\beta\leq 0$ the variance, skewness and kurtosis excess of the
process defined by Eq.(1) are exactly the same as for the Poisson process if
we take into account that for the process defined by Eq.(1) we have

$<n_{1,\beta}>=<n_{1+\beta}>|_{\mu=1}=\frac{\lambda_{1+\beta}t^{1+\beta}}{1+\beta},$

which can be expressed in Poisson-like form using Eq.(43)

$<n_{1+\beta}>\equiv<n_{\sigma}>=\overline{\lambda}_{\sigma}t^{\sigma}.$

####  6.4.1 New representation of the Kilbas and Saigo function

To obtain a new representation for the Kilbas and Saigo function
$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)$ defined by Eq.(39) we use
the generating function of the Stirling numbers of the first kind666The
Stirling numbers of the first kind are defined as the coefficients $s(m,l)$ in
the expansion
$\frac{x!}{(x-m)!}=x(x-1)...(x-m+1)=\sum\limits_{l=0}^{m}s(m,l)x^{l}.$ (124)

$(1+z)^{n}=\sum\limits_{m=0}^{\infty}\frac{z^{m}}{m!}\sum\limits_{l=0}^{m}s(m,l)n^{l},$
(125)

where $s(m,l)$ stands for the Stirling numbers of the first kind, see Sec.
24.1.3 in [20]. Then, we evaluate expectations of both sides of Eq.(125) using
the probability distribution function of fractional non-homogeneous counting
process $P_{\mu,\beta}(n,\rho)$ given by Eq.(54)

$\sum\limits_{n=0}^{\infty}P_{\mu,\beta}(n,1)(1+z)^{n}=\sum\limits_{n=0}^{\infty}P_{\mu,\beta}(n,1)\sum\limits_{m=0}^{\infty}\frac{z^{m}}{m!}\sum\limits_{l=0}^{m}s(m,l)n^{l}.$
(126)

It is easy to see that the calculations on the left side of Eq.(126) lead to
$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)$, and we

$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)=\sum\limits_{n=0}^{\infty}P_{\mu}(n,1)\sum\limits_{m=0}^{\infty}\frac{z^{m}}{m!}\sum\limits_{l=0}^{m}s(m,l)n^{l}.$

Then, using the definition of the fractional numbers $B_{\mu,\beta}(m)$ given
by Eq.(89) we obtain new representations of the Kilbas and Saigo function

$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)=\sum\limits_{m=0}^{\infty}\frac{z^{m}}{m!}\sum\limits_{l=0}^{m}s(m,l)B_{\mu,\beta}(l)$
(127)

through Stirling numbers of the first kind and the fractional numbers
$B_{\mu,\beta}(l)$. Using Eq.(104) the Kilbas and Saigo function can be
expressed through Stirling numbers of the first kind and fractional
combinatorial numbers $S_{\mu,\beta}(m,l)$

$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)=\sum\limits_{m=0}^{\infty}\frac{z^{m}}{m!}\sum\limits_{l=0}^{m}\sum\limits_{r=0}^{l}s(m,l)S_{\mu,\beta}(l,r).$
(128)

The two equations (127) and (128) are new representations of the Kilbas and
Saigo function.

From Eq.(127) we see that the Kilbas and Saigo function can be considered as
generating function of the sum $\sum\limits_{l=0}^{m}s(m,l)B_{\mu,\beta}(l)$,
that is

$\sum\limits_{l=0}^{m}s(m,l)B_{\mu,\beta}(l)=\frac{\partial^{m}}{\partial
z^{m}}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)|_{z=0}.$ (129)

Using Eq.(39), we obtain the identity

$\sum\limits_{l=0}^{m}s(m,l)B_{\mu,\beta}(l)=m!\prod\limits_{k=0}^{m-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$
(130)

From Eq.(128) we see that the Kilbas and Saigo function can be considered as
generating function of the sum
$\sum\limits_{l=0}^{m}\sum\limits_{r=0}^{l}s(m,l)S_{\mu,\beta}(l,r)$, that is

$\sum\limits_{l=0}^{m}\sum\limits_{r=0}^{l}s(m,l)S_{\mu,\beta}(l,r)=\frac{\partial^{m}}{\partial
z^{m}}E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)|_{z=0}.$

Using Eq.(39), we obtain the identity

$\sum\limits_{l=0}^{m}\sum\limits_{r=0}^{l}s(m,l)S_{\mu,\beta}(l,r)=m!\prod\limits_{k=0}^{m-1}\frac{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu})+1]}{\Gamma[\mu(k(1+\frac{\beta}{\mu})+\frac{\beta}{\mu}+1)+1]}.$
(131)

##  7 Conclusions

A new fractional non-homogeneous counting process has been introduced and
developed using the Kilbas and Saigo three-parameter generalization of the
Mittag-Leffler function. The probability distribution function of this process
reproduces for certain set of the fractality parameters the famous Poisson and
fractional Poisson probability distributions as well as the probability
distribution function of a counting process displaying stretched exponential
interarrival times distribution.

Applications of the developed fractional non-homogeneous counting process
cover fractional compound process, the generalization of combinatorial
polynomials and numbers and their generating functions, variance, skewness and
kurtosis excess of the fractional non-homogeneous counting process expressed
through four central moments and a new representation of the Kilbas and Saigo
function.

The presented results will stimulate further development of the fundamentals
and various applications in a new field of probability theory and stochastic
processes, which can be called fractional stochastic processes and their
applications.

##  8 Appendix A

###  8.1 Special cases of the Kilbas and Saigo function
$E_{\mu,1+\frac{\beta}{\mu},\frac{\beta}{\mu}}(z)$

Defined by Eq.(39) the Kilbas and Saigo function has the following special
cases.

$\bullet$ When $\mu=1$ and $\beta=0$, Eq.(39) goes to

$E_{1,1,0}(z)=1+\sum\limits_{n=1}^{\infty}z^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[k+1]}{\Gamma[(k+1)+1]},$

Then due to

$\prod\limits_{k=0}^{n-1}\frac{\Gamma[k+1]}{\Gamma[(k+1)+1]}=\prod\limits_{k=0}^{n-1}\frac{1}{k+1}=\frac{1}{n!},$
(132)

where the following property of the Gamma function was used

$\Gamma(z+1)=z\Gamma(z),$ (133)

we obtain

$E_{1,1,0}(z)=1+\sum\limits_{n=1}^{\infty}\frac{z^{n}}{n!}=\exp(z).$ (134)

$\bullet$ When $0<\mu\leq 1$ and $\beta=0$, Eq.(39) goes to

$E_{\mu,1,0}(z)=1+\sum\limits_{n=1}^{\infty}z^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu
k+1]}{\Gamma[\mu(k+1)+1]}.$

Then due to

$\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu
k+1]}{\Gamma[\mu(k+1)+1]}=\frac{1}{\Gamma(\mu n+1)},$ (135)

we obtain

$E_{\mu,1,0}(z)=1+\sum\limits_{n=1}^{\infty}\frac{z^{n}}{\Gamma(\mu
n+1)}=E_{\mu}(z),$ (136)

here $E_{\mu}(z)$ is the Mittag-Leffler function defined by Eq.(42).

$\bullet$ When $\mu=1$ and $-1<\beta\leq 0$, Eq.(39) goes to

$E_{1,1+\beta,\beta}(z)=1+\sum\limits_{n=1}^{\infty}z^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma[\mu(k(1+\beta)+\beta)+1]}{\Gamma[\mu(k(1+\beta)+\beta+1)+1]}.$

Then due to

$\prod\limits_{k=0}^{n-1}\frac{\Gamma[(k(1+\beta)+\beta)+1]}{\Gamma[(k(1+\beta)+\beta+1)+1]}=\prod\limits_{k=0}^{n-1}\frac{1}{(k+1)(\beta+1)}=\frac{1}{(1+\beta)^{n}}\frac{1}{n!},$
(137)

we obtain

$E_{1,1+\beta,\beta}(z)=1+\sum\limits_{n=1}^{\infty}\frac{1}{n!}\frac{z^{n}}{(1+\beta)^{n}}=\exp(\frac{z}{1+\beta}).$
(138)

$\bullet$ When $\mu+\beta=1$, Eq.(39) goes to

$E_{\mu,\frac{1}{\mu},\frac{1-\mu}{\mu}}(z)=1+\sum\limits_{n=1}^{\infty}z^{n}\prod\limits_{k=0}^{n-1}\frac{\Gamma(k+2-\mu)}{\Gamma(k+2)}.$
(139)

##  9 Appendix B

###  9.1 Four first central moments

Using Eq.(120), the first four central moments can be calculated. We expressed
the results of the calculation in the form of a power series over the first-
order moment $<n_{\mu,\beta}>$:

$\mathrm{M}_{\mu,\beta}(1)=0,$ (140)

$\mathrm{M}_{\mu,\beta}(2)=\left(\frac{2\Gamma(\mu+\beta+1)\Gamma(\mu+2\beta+1)}{\Gamma(\beta+1)\Gamma(2\mu+2\beta+1)}-1\right)<n_{\mu,\beta}>^{2}+<n_{\mu,\beta}>,$
(141)

$\mathrm{M}_{\mu,\beta}(3)=2{\LARGE(}\frac{3\Gamma(\mu+\beta+1)}{\Gamma(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}\frac{\Gamma(\mu+\beta+1)}{\Gamma(\beta+1)}\frac{\Gamma(2\mu+3\beta+1)}{\Gamma(3\mu+3\beta+1)}-$

$\frac{3\Gamma(\mu+\beta+1)}{\Gamma(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}+1{\LARGE)}<n_{\mu,\beta}>^{3}+$
(142)

$3\left(\frac{2\Gamma(\mu+\beta+1)}{\Gamma(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}-1\right)<n_{\mu,\beta}>^{2}+<n_{\mu,\beta}>,$

and

$\mathrm{M}_{\mu,\beta}(4)=3{\LARGE(}8\frac{\Gamma^{3}(\mu+\beta+1)}{\Gamma^{3}(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}\frac{\Gamma(2\mu+3\beta+1)}{\Gamma(3\mu+3\beta+1)}\frac{\Gamma(3\mu+4\beta+1)}{\Gamma(4\mu+4\beta+1)}-$

$8\frac{\Gamma^{2}(\mu+\beta+1)}{\Gamma^{2}(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}\frac{\Gamma(2\mu+3\beta+1)}{\Gamma(3\mu+3\beta+1)}+$

$4\frac{\Gamma(\mu+\beta+1)}{\Gamma(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}-1{\LARGE)}<n_{\mu,\beta}>^{4}+$
(143)

$6(\frac{6\Gamma^{2}(\mu+\beta+1)}{\Gamma^{2}(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}\frac{\Gamma(2\mu+3\beta+1)}{\Gamma(3\mu+3\beta+1)}-$

$4\frac{\Gamma(\mu+\beta+1)}{\Gamma(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}+1)<n_{\mu,\beta}>^{3}+$

$2\left(7\frac{\Gamma(\mu+\beta+1)}{\Gamma(\beta+1)}\frac{\Gamma(\mu+2\beta+1)}{\Gamma(2\mu+2\beta+1)}-2\right)<n_{\mu,\beta}>^{2}+<n_{\mu,\beta}>.$

## References

  * [1] N. Laskin, Stretched Coherent States, https://arxiv.org/abs/2302.07563
  * [2] J. Laherrere and D. Sornette, Stretched exponential distributions in nature and economy: ”Fat tails” with characteristic scales, https://arxiv.org/pdf/cond-mat/9801293.pdf
  * [3] W. Feller, An introduction to probability theory and its applications, Vol. 1, John Wiley & Sons, Inc., New York, 1968.
  * [4] N. Laskin, Fractional Poisson process, Communications in Nonlinear Science and Numerical Simulation, 8 (2003) 201-213.
  * [5] M. Kuno, D.P. Fromm, H. F. Hamann, A. Gallagher, and D.J. Nesbitt, Nonexponential blinking kinetics of single CdSe quantum dots: A universal power law behavior, J. Chem. Phys. 112 (2000) 3117-3120.
  * [6] L. Sabatelli et al, Waiting time distributions in financial markets, The European Physical Journal B - Condensed Matter and Complex Systems, 27 (2002) 273-275.
  * [7] W. Willinger, and V. Paxson, Where Mathematics Meets the Internet, Notices of the American Mathematical Society, 45 (1998) 961-970.
  * [8] K.B. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York, 1974.
  * [9] S.G. Samko, A.A. Kilbas, and O.I. Marichev, Fractional Integrals and Derivatives and Their Applications, Gordon and Breach Science Publishers, Langhorne, PA, 1993.
  * [10] K.S. Miller and B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, J. Wiley & Sons, New York, 1993\.
  * [11] M. Caputo, Linear model of dissipation whose Q is almost frequency independent-II, Geophys., J. R. Astron. Soc., 13 (1967) 529–539.
  * [12] A. Kilbasi and M. Saigo, On Mittag-Leffler type function, fractional calculus operators and solutions of integral equations, Integral Transforms and Special Functions, 4, (1996) 355-370.
  * [13] G. Mittag-Leffler, Sur la repriésentation analytique d’une branche uniforme d’une fonction monogène, Acta Mathematics, 29 (1905) 101-182.
  * [14] L. Boudabsa, T. Simon,, and P. Valloist, Fractional extreme distributions, Electron. J. Probab., 25 (2020) 1-20 (2020), DOI: 10.1214/20-EJP520
  * [15] L. Boudabsa, T. Simon, Some Properties of the Kilbas-Saigo Function, Mathematics, 9 (2021) 217, 24 pages, https://doi.org/10.3390/math9030217
  * [16] A. Erdélyi, Ed., Higher Transcendental Functions, Vol.3, pp. 206-227 (Chapter 18 Miscellaneous Functions) McGraw-Hill, New York, 1955.
  * [17] E. T. Bell, Exponential Polynomials, The Annals of Mathematics, Second Series, 35 (1934), 258-277.
  * [18] G. Dobiński, Summirung der Reihe $\sum n^{m}/n!$, für $m=1,2,3,4,5,...$, Grunert Archiv (Arch. für M. und Physik) 61 (1877) 333-336.
  * [19] N. Laskin, Some applications of the fractional Poisson probability distribution, J Math Phys, 50 (2009) 113513, 12 pages.
  * [20] M. Abramowitz and I.A. Stegun, (Eds.). Stirling Numbers of the Second Kind. §24.1.4 in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing, pp. 824-825, New York, Dover, 1972.
  * [21] C. A. Charalambides and J. Singh: A review of the Stirling numbers, their generalizations and statistical applications, Communications in Statistics - Theory and Methods, 17 (1988) 2533-2595.
  * [22] E.W. Weisstein, Poisson Distribution, From MathWorld \- A Wolfram Web Resource, http://mathworld.wolfram.com/PoissonDistribution.html