**2.2 Approximate transition density for multidimensional and nonlinear SDME model**

Usually, a formal general solution of the stochastic differential equations as in Eq. (1) cannot be given, which makes the calculation of the transition density for this process more complicated. Moreover, this process has a lot of fluctuations that its exact position cannot be determined but can be known given a region by its probability density; with the FP equation, such a probability density can be determined. The FP equation is a differential equation for the distribution function describing a Brownian motion by which the probability density of the stochastic process can be calculated in a much simpler way by solving this equation. This motion equation is usually used for variables describing a macroscopic but small

system, where the fluctuations are important as for some cases in physics, e.g., the position and the speed of the Brownian motion of a small particle. However, it can be also used for the larger system where, in spite of their small fluctuations, the stochastic description remains necessary when the deterministic equations may not be stable for this type of system.

So, under some assumptions, the transition density *PY y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�1, *bi* , *θ* � � of the process *Y* in the SDME model Eq. (1) is the solution of the following functional partial differential Equation [23, 34]:

$$\frac{\partial P\_Y\left(y\_j^i, \Delta\_j^i | y\_{j-1}^i, b^i, \theta\right)}{\partial t} = L\_{FP} P\_Y\left(y\_j^i, \Delta\_j^i | y\_{j-1}^i, b^i, \theta\right),\tag{10}$$

gradients of the likelihood analytically. Here, we propose to use the genetic algorithm as an optimization tool to maximize the approximate likelihood function

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential…*

^*θ*, <sup>Ψ</sup>^ � � <sup>¼</sup> *argminGA* � *log L*ð Þ *<sup>a</sup>* ð Þ *<sup>θ</sup>*, <sup>Ψ</sup>

The genetic algorithm is a random search technique to look for an exact or approximated optimum points for optimization problems [36–38]. It is based on the concepts of natural genetic evolution which contains the following stages: the reproduction, the crossing, and the mutation of a constantly evolving population. It sets up the evolution of a random population of potential solutions of the N cardinal; then, the N simultaneous iterative trajectories interact with each other by following or imitating the biological evolution, for a convergence of certain elements of the population towards an optimal point of the fitness

The GA can search in multiple directions to explore all the search space by the possibility of jumping across them, so that the seeds spread uniformly over the whole search space. In this algorithm, we have a diversity of initial populations which gives the global optimum faster than other algorithms, where the initial value is very important and should be enough close to the global optimum. All of these features allows the GA to be regarded as a driving tool of evolution giving good results for optimization processes [37, 39]. In the literature, there were many works about the application of GA in optimizing problems specially for likelihood function [40, 41]. For the use of GA, we must first define some parameters of the algorithm: Population size N, EN, SR, CP, MP, fitness function, and convergence criteria. In

> ð Þ 0 <sup>1</sup> , *β* ð Þ 0 <sup>2</sup> , … , *<sup>β</sup>*ð Þ <sup>0</sup> *N*

**4**. Replacement step (by using SR and EN): At the SR rate, individuals with the

**5**. Selection operator by using roulette wheel method, based on the fact that the more the individual has a good result of fitness function, the more likely he will be

**6**. Crossover operator by using CP and mutation operator by using MP: It is a mechanism of perturbation on the candidate individuals (parents) according to CP and MP to generate new groups of individuals, and we obtain a new ð Þ *m* þ 1 *nd*

worst results in step 2 of fitness function are replaced by new ones randomly generated, and a number EN of individuals is selected and accepted for the next

� � � � *:* (13)

n o, *<sup>m</sup>* <sup>¼</sup> 0 via an initialization

Eq. (12) using MATLAB software:

*DOI: http://dx.doi.org/10.5772/intechopen.90751*

the following we present the GA steps:

**1**. Generate initial population *β*

<sup>1</sup> , *<sup>β</sup>*ð Þ *<sup>m</sup>*þ<sup>1</sup>

**7**. Evolution stops; get GA output.

strategy (random generation), in our case *β* ¼ ð Þ *θ*, Ψ .

**3**. While (convergence criteria are not satisfied):

<sup>2</sup> , … , *<sup>β</sup>*ð Þ *<sup>m</sup>*þ<sup>1</sup> *N*

n o.

**<sup>2</sup>**. Evaluate the Fitness function *log* �*L*ð Þ *<sup>a</sup> <sup>θ</sup>*ð Þ *<sup>m</sup>* , <sup>Ψ</sup>ð Þ *<sup>m</sup>* � � � � .

**Steps of GA:**

For *m* ¼ 0:

**Do**:

step.

selected.

population *β*ð Þ *<sup>m</sup>*þ<sup>1</sup>

**8**. *m* ¼ *m* þ 1. End For.

**Else**:

**97**

**2.3 Genetic algorithm**

function.

where *LFP* is the FP operator and *LFP* ¼ �P*<sup>N</sup> k*¼1 *∂ <sup>∂</sup>Yk <sup>μ</sup><sup>k</sup> <sup>Y</sup><sup>i</sup>* , *<sup>t</sup>*, *<sup>θ</sup>*, *<sup>b</sup><sup>i</sup>* � � � � <sup>þ</sup> 1 2 P*<sup>N</sup> k*¼1 P*<sup>N</sup> l*¼1 *∂*2 *<sup>∂</sup>Yk∂Yl* Σ*kl* ½ � and *μ<sup>k</sup>* is the k-th element of the drift vector and Σ*kl* is the klth element of the diffusion matrix. We assume that the diffusion matrix is positive definite so that its inverse exists *Det*ð Þ Σ 6¼ 0.

In [34], H. Risken deals in this book with the derivation of the FP equation and its solution methods with some of its applications especially for problems of Brownian motion. So, here, we propose to characterize the transition density of the N-dimensional process *Y<sup>i</sup> <sup>t</sup>* in Eq. (1) using the Risken approximation based on the formal solution of Eq. (10) proposed in [34], p: 4.109, denoted by *P*ð Þ *<sup>a</sup> <sup>Y</sup> yi j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�<sup>1</sup>, *bi* , *θ* � �:

$$\begin{split} P\_{Y}^{(a)}\left(\mathbf{Y}\_{j}^{i},\Delta\_{j}^{i}|Y\_{j-1}^{i},b^{i},\boldsymbol{\theta}\right) &= \left(2\sqrt{\boldsymbol{\Pi}\Delta\_{j}}\right)^{-N}[\mathbf{Det}\boldsymbol{\Sigma}]^{-\frac{1}{2}} \ast \exp\left(-\frac{1}{4\Delta\_{j}}\sum\_{l=1}^{N}\sum\_{k=1}^{N}\left[\boldsymbol{\Sigma}^{-1}\right]\_{lk} \\ \mathbb{E}\left[\left(Y\_{j}^{i}\right)\_{l}-\left(Y\_{j-1}^{i}\right)\_{l}-\mu\_{l}\left(Y\_{j-1}^{i},t,\theta,b^{i}\right)\Delta\_{j}^{i}\right]\left[\left(Y\_{j}^{i}\right)\_{k}-\left(Y\_{j-1}^{i}\right)\_{k}-\mu\_{k}\left(Y\_{j-1}^{i},t\_{j-1},\theta,b^{i}\right)\Delta\_{j}^{i}\right]\right) \end{split} \tag{11}$$

To test the effectiveness of this approach, we will guide our statistical methodology by simulation studies in order to examine the flexibility of its application to deduce its advantages and disadvantages. So, we substitute *P*ð Þ *<sup>a</sup> <sup>Y</sup> Y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*Yi <sup>j</sup>*�1, *<sup>b</sup><sup>i</sup>* , *θ* � � in Eq. (6), and by solving the integral with respect to the random effects density and maximizing the obtained likelihood function, we get the approximated estimators ^*θ* and Ψ^.

#### *2.2.1 Approximated estimators*

For a nonlinear SDME model with Gaussian random effects using Eqs. (6), (7), and (11), we obtain the following approximated likelihood function:

$$\log L^{(a)}(\theta,\Psi) \simeq \sum\_{i=1}^{M} \left[ -\frac{q}{2} \log \left( (2\pi) \right) - \frac{1}{2} \log \left( \det(\phi) \right) - \frac{n\_i}{2} \log \left( \det(\Sigma) \right) + \left( \sum\_{j=1}^{n\_i} \log \left( \left( 2\sqrt{\Pi \Delta\_{\hat{\mathbb{}}}} \right)^{-N} \right) - \frac{1}{4\Delta\_{\hat{\mathbb{}}}} \right) \right]$$

$$\left[ \Sigma^{-1} \right]\_{\mathbb{H}} \left[ \left( Y^i\_j \right)\_{\mathbb{I}} - \left( Y^i\_{j-1} \right)\_{\mathbb{I}} - \mu\_l \left( Y^i\_{j-1}, t, \theta, \mathring{b}^i \right) \Delta^i\_l \right] \left[ \left( Y^i\_j \right)\_{\mathbb{k}} - \left( Y^i\_{j-1} \right)\_{\mathbb{k}} - \mu\_k \left( Y^i\_{j-1}, t\_{j-1}, \theta, \mathring{b}^i \right) \Delta^i\_l \right] \right)$$

$$- \frac{1}{2} \left( \mathring{b}^i - \nu \right)' \phi^{-1} \left( \mathring{b}^i - \nu \right) + \frac{q}{2} \log \left( 2\pi \right) - \frac{1}{2} \log \left( \det \left( -H \left( \mathring{b}^i | \theta, \Psi \right) \right) \right) \right]. \tag{12}$$

The MLEs of ð Þ *θ*, Ψ can be obtained using one of the optimization tools and numerical computation software, especially when it is complicated to compute the *Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential… DOI: http://dx.doi.org/10.5772/intechopen.90751*

gradients of the likelihood analytically. Here, we propose to use the genetic algorithm as an optimization tool to maximize the approximate likelihood function Eq. (12) using MATLAB software:

$$\left(\hat{\theta}, \hat{\Psi}\right) = \operatorname\*{argmin}\_{GA} \left(-\log\left(L^{(a)}(\theta, \Psi)\right)\right). \tag{13}$$

#### **2.3 Genetic algorithm**

system, where the fluctuations are important as for some cases in physics, e.g., the position and the speed of the Brownian motion of a small particle. However, it can be also used for the larger system where, in spite of their small fluctuations, the stochastic description remains necessary when the deterministic equations may not

process *Y* in the SDME model Eq. (1) is the solution of the following functional

*<sup>∂</sup><sup>t</sup>* <sup>¼</sup> *LFPPY <sup>y</sup><sup>i</sup>*

th element of the diffusion matrix. We assume that the diffusion matrix is positive

its solution methods with some of its applications especially for problems of Brownian motion. So, here, we propose to characterize the transition density of the

½ � *Det*<sup>Σ</sup> �<sup>1</sup>

*Yi j* � � *k* � *<sup>Y</sup><sup>i</sup> j*�1 � � *k*

<sup>2</sup> ∗ *exp* � � 1 4Δ*<sup>j</sup>* X*N l*¼1

To test the effectiveness of this approach, we will guide our statistical methodology by simulation studies in order to examine the flexibility of its application to deduce its

by solving the integral with respect to the random effects density and maximizing the

For a nonlinear SDME model with Gaussian random effects using Eqs. (6), (7),

Δ*i j*

obtained likelihood function, we get the approximated estimators ^*θ* and Ψ^.

and (11), we obtain the following approximated likelihood function:

<sup>2</sup> *log det* ð Þ� ð Þ *<sup>ϕ</sup> ni*

*<sup>j</sup>*�1, *<sup>t</sup>*, *<sup>θ</sup>*, <sup>~</sup> *b <sup>i</sup>* � �

The MLEs of ð Þ *θ*, Ψ can be obtained using one of the optimization tools and numerical computation software, especially when it is complicated to compute the

<sup>2</sup> *log* ð Þ� <sup>2</sup>*<sup>π</sup>* <sup>1</sup>

formal solution of Eq. (10) proposed in [34], p: 4.109, denoted by

Δ*i j*

In [34], H. Risken deals in this book with the derivation of the FP equation and

*j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�1, *<sup>b</sup><sup>i</sup>* , *θ*

*k*¼1 *∂ <sup>∂</sup>Yk <sup>μ</sup><sup>k</sup> <sup>Y</sup><sup>i</sup>*

*<sup>t</sup>* in Eq. (1) using the Risken approximation based on the

X*N k*¼1

� �

�X*ni j*¼1

*log* 2 ffiffiffiffiffiffiffiffi ΠΔ*<sup>j</sup>* � � p �*<sup>N</sup>* � �

h i�

*:*

*<sup>j</sup>*�1, *tj*�1, *<sup>θ</sup>*, <sup>~</sup> *b*

*<sup>i</sup>* � �

� *<sup>μ</sup><sup>k</sup> <sup>Y</sup><sup>i</sup>*

*<sup>Y</sup> Y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*Yi <sup>j</sup>*�1, *<sup>b</sup><sup>i</sup>* , *θ*

<sup>2</sup> *log det* ð Þþ ð Þ <sup>Σ</sup>

*b i* j*θ*, Ψ � � � � � � �

*Yi j* � � *k* � *<sup>Y</sup><sup>i</sup> j*�1 � � *k* � *<sup>μ</sup><sup>k</sup> <sup>Y</sup><sup>i</sup>*

<sup>2</sup> *log det* �*<sup>H</sup>* <sup>~</sup>

h i�

Σ�<sup>1</sup> � � *lk*

*<sup>j</sup>*�1, *tj*�1, *<sup>θ</sup>*, *bi* � �

Δ*i j* (11)

in Eq. (6), and

� 1 4Δ*<sup>j</sup>*

> Δ*i j*

(12)

*<sup>∂</sup>Yk∂Yl* Σ*kl* ½ � and *μ<sup>k</sup>* is the k-th element of the drift vector and Σ*kl* is the kl-

� �

*j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�1, *bi* , *θ*

, *<sup>t</sup>*, *<sup>θ</sup>*, *<sup>b</sup><sup>i</sup>* � � � � <sup>þ</sup>

� �

of the

, (10)

So, under some assumptions, the transition density *PY y<sup>i</sup>*

� �

where *LFP* is the FP operator and *LFP* ¼ �P*<sup>N</sup>*

be stable for this type of system.

*∂*2

N-dimensional process *Y<sup>i</sup>*

� �

� �

*2.2.1 Approximated estimators*

*i*¼1

*lk <sup>Y</sup><sup>i</sup> j* � � *l* � *<sup>Y</sup><sup>i</sup> j*�1 � � *l* � *<sup>μ</sup><sup>l</sup> <sup>Y</sup><sup>i</sup>*

� � *q*

*logL*ð Þ *<sup>a</sup>* ð Þ *<sup>θ</sup>*, <sup>Ψ</sup> <sup>≃</sup> <sup>X</sup>*<sup>M</sup>*

Σ�<sup>1</sup> � �

� 1 2 ~ *b i* � *ν* � �<sup>0</sup>

**96**

1 2 P*<sup>N</sup> k*¼1 P*<sup>N</sup> l*¼1

*P*ð Þ *<sup>a</sup> <sup>Y</sup> yi j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�<sup>1</sup>, *bi* , *θ*

*P*ð Þ *<sup>a</sup> <sup>Y</sup> Y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*Yi <sup>j</sup>*�1, *bi* , *θ*

*Yi j* � � *l* � *<sup>Y</sup><sup>i</sup> j*�1 � � *l* � *<sup>μ</sup><sup>l</sup> <sup>Y</sup><sup>i</sup>*

partial differential Equation [23, 34]:

*Numerical Modeling and Computer Simulation*

*∂PY y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�1, *bi* , *θ*

definite so that its inverse exists *Det*ð Þ Σ 6¼ 0.

:

h i

<sup>¼</sup> <sup>2</sup> ffiffiffiffiffiffiffiffi ΠΔ*<sup>j</sup>* � � p �*<sup>N</sup>*

*<sup>j</sup>*�1, *<sup>t</sup>*, *<sup>θ</sup>*, *bi* � �

advantages and disadvantages. So, we substitute *P*ð Þ *<sup>a</sup>*

<sup>2</sup> *log* ð Þ� ð Þ <sup>2</sup>*<sup>π</sup>* <sup>1</sup>

*ϕ*�<sup>1</sup> ~ *b i* � *ν* � �

h i

þ *q*

The genetic algorithm is a random search technique to look for an exact or approximated optimum points for optimization problems [36–38]. It is based on the concepts of natural genetic evolution which contains the following stages: the reproduction, the crossing, and the mutation of a constantly evolving population. It sets up the evolution of a random population of potential solutions of the N cardinal; then, the N simultaneous iterative trajectories interact with each other by following or imitating the biological evolution, for a convergence of certain elements of the population towards an optimal point of the fitness function.

The GA can search in multiple directions to explore all the search space by the possibility of jumping across them, so that the seeds spread uniformly over the whole search space. In this algorithm, we have a diversity of initial populations which gives the global optimum faster than other algorithms, where the initial value is very important and should be enough close to the global optimum. All of these features allows the GA to be regarded as a driving tool of evolution giving good results for optimization processes [37, 39]. In the literature, there were many works about the application of GA in optimizing problems specially for likelihood function [40, 41]. For the use of GA, we must first define some parameters of the algorithm: Population size N, EN, SR, CP, MP, fitness function, and convergence criteria. In the following we present the GA steps:

#### **Steps of GA:**

**1**. Generate initial population *β* ð Þ 0 <sup>1</sup> , *β* ð Þ 0 <sup>2</sup> , … , *<sup>β</sup>*ð Þ <sup>0</sup> *N* n o, *<sup>m</sup>* <sup>¼</sup> 0 via an initialization strategy (random generation), in our case *β* ¼ ð Þ *θ*, Ψ .

For *m* ¼ 0:

**<sup>2</sup>**. Evaluate the Fitness function *log* �*L*ð Þ *<sup>a</sup> <sup>θ</sup>*ð Þ *<sup>m</sup>* , <sup>Ψ</sup>ð Þ *<sup>m</sup>* � � � � .

**3**. While (convergence criteria are not satisfied):

**Do**:

**4**. Replacement step (by using SR and EN): At the SR rate, individuals with the worst results in step 2 of fitness function are replaced by new ones randomly generated, and a number EN of individuals is selected and accepted for the next step.

**5**. Selection operator by using roulette wheel method, based on the fact that the more the individual has a good result of fitness function, the more likely he will be selected.

**6**. Crossover operator by using CP and mutation operator by using MP: It is a mechanism of perturbation on the candidate individuals (parents) according to CP and MP to generate new groups of individuals, and we obtain a new ð Þ *m* þ 1 *nd* population *β*ð Þ *<sup>m</sup>*þ<sup>1</sup> <sup>1</sup> , *<sup>β</sup>*ð Þ *<sup>m</sup>*þ<sup>1</sup> <sup>2</sup> , … , *<sup>β</sup>*ð Þ *<sup>m</sup>*þ<sup>1</sup> n o.

*N* **Else**: **7**. Evolution stops; get GA output. **8**. *m* ¼ *m* þ 1. End For.

In the simulation studies paragraphs, the GA is implemented using the MATLAB software, where the function "ga," to generate the genetic algorithm, requires inputs that are chosen according to the constraints of each example (see the help window in MATLAB). Moreover, the algorithm parameters are chosen according to what is early used in the literature [39], *EN* ¼ 4, *MP* ¼ 0*:*2, *CP* ¼ 0*:*8 *and SR* ¼ 1*=*3, and the search spaces are around the confidence interval of the minimal model parameters (see [42] and references therein).

with mean vector *<sup>μ</sup>* <sup>¼</sup> *<sup>α</sup>* <sup>þ</sup> *<sup>Y</sup><sup>i</sup>*

*DOI: http://dx.doi.org/10.5772/intechopen.90751*

<sup>∣</sup> <sup>j</sup>*β:bi*

approximated transition density of Y:

*j* i <sup>2</sup> <sup>þ</sup> <sup>Σ</sup>�<sup>1</sup> <sup>22</sup> � � *<sup>Y</sup>*ð Þ<sup>2</sup> *<sup>i</sup>*

*j* � � *<sup>τ</sup> exp* � *<sup>β</sup>:bi* � �<sup>0</sup>

> ΠΔ*<sup>j</sup>* � � p �<sup>2</sup>

scheme [45] with the following set of parameters:

*<sup>ς</sup>* <sup>¼</sup> *<sup>τ</sup>* � *exp* � *<sup>β</sup>:b<sup>i</sup>* � �Δ*<sup>i</sup>*

2*tr β:bi* � �∣*β:b<sup>i</sup>*

eigenvalues of *β:bi*

<sup>12</sup> *<sup>Y</sup>*ð Þ<sup>2</sup> *<sup>i</sup> tj*�<sup>1</sup> � *α*<sup>2</sup> � ��Δ*<sup>i</sup>*

*<sup>τ</sup>* <sup>¼</sup> <sup>1</sup>

*P*ð Þ *<sup>a</sup> <sup>Y</sup> Y<sup>i</sup> tj* , Δ*<sup>i</sup> j* j*Yi tj*�<sup>1</sup> , *bi* , *θ* � � <sup>¼</sup> <sup>2</sup> ffiffiffiffiffiffiffiffi

<sup>þ</sup> *<sup>β</sup>*12*b<sup>i</sup>*

**Figure 1.**

**99**

*tj*�<sup>1</sup> � *<sup>α</sup>*

� � � � , where

� � *exp* � *<sup>β</sup>:bi* � �Δ*<sup>i</sup>*

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential…*

We assume that the matrices *β:b<sup>i</sup>* and Σ have full rank and the real parts of the

exists. Under these assumptions, we derive from Eqs. (14) and (11) the following

<sup>2</sup> <sup>∗</sup> *exp* �

*tj*�<sup>1</sup> <sup>þ</sup> *<sup>β</sup>*21*b<sup>i</sup>*

In **Figure 1**, we illustrate in (a) the simulation of the OU process using the Euler

*A sample path of the OU process in the third graph of (a) for the given parameters set with the initial condition:* Y0 ¼ *(3,3) and time interval [3,10]; and the transition density for a transition from* Yj *to* Yjþ<sup>1</sup> *in (b).*

� 1 4Δ*<sup>j</sup>* Σ�<sup>1</sup> <sup>11</sup> � �h *Y*ð Þ<sup>1</sup> *<sup>i</sup> tj* � *<sup>Y</sup>*ð Þ<sup>1</sup> *<sup>i</sup> tj*�<sup>1</sup> þ � *β*11*b<sup>i</sup>* <sup>11</sup> *<sup>Y</sup>*ð Þ<sup>1</sup> *<sup>i</sup> tj*�<sup>1</sup> � *α*<sup>1</sup> � �

<sup>21</sup> *<sup>Y</sup>*ð Þ<sup>1</sup> *<sup>i</sup> tj*�<sup>1</sup> � *α*<sup>1</sup> � � <sup>þ</sup> *<sup>β</sup>*22*b<sup>i</sup>*

h i<sup>2</sup>

ð Þ <sup>Σ</sup>11Σ<sup>22</sup> �<sup>1</sup>

*tj* � *<sup>Y</sup>*ð Þ<sup>2</sup> *<sup>i</sup>*

Δ*i j*

<sup>j</sup>ΣΣ<sup>0</sup> <sup>þ</sup> *<sup>β</sup>:bi* � *tr <sup>β</sup>:bi* � �*:<sup>I</sup>* � �ΣΣ<sup>0</sup> *<sup>β</sup>:b<sup>i</sup>* � *tr <sup>β</sup>:b<sup>i</sup>* � �*:<sup>I</sup>* � �<sup>0</sup> � �*:* (17)

are positive definite in order that a stationary solution to Eq. (15)

*j*

� � and covariance matrix

<sup>22</sup> *<sup>Y</sup>*ð Þ<sup>2</sup> *<sup>i</sup> tj*�<sup>1</sup> � *α*<sup>2</sup>

*j*

� *:*

(18)

� � � � <sup>Δ</sup>*<sup>i</sup>*
