- Email: [email protected]
q-Gaussian Density Model and its Application to State Estimation of Nonlinear Systems
3rd IFAC International Conference on Intelligent Control and Automation Science. September 2-4, 2013. Chengdu, China
q-Gaussian Density Model and its Application to State Estimation of Nonlinear Systems Xifeng Li*, Yongle Xie.** * University of Electronic Science and Technology of China, Chengdu, Sichuan, 611731 China (Tel: +86-13458530206; e-mail: xfengl@ hotmail.com). ** University of Electronic Science and Technology of China, Chengdu, Sichuan, 611731 China (e-mail: [email protected]) Abstract: Probability density function (PDF) plays a vital role in system analysis involving stochastic factors. A good estimate of true PDF conditioned under certain performance criterion could help acquire more information of the system. With help of the new information, many features of the system that we are concerning can be revealed effectively, especially for nonlinear non-Gaussian stochastic systems. In this paper, based on the Tsallis entropy, we derive a class of PDFs with explicit form called q-Gaussian PDFs. These PDFs have a parameter that indicates the fractal feature of the system. Based on the explicit form of q-Gaussian PDFs, we propose an extension of Gaussian particle filter (GPF) called q-Gaussian particle filter (q-GPF). The experimental results show that the q-GPF is a more effective method to estimate the state of nonlinear stochastic system compared with the GPF. Keywords: Tsallis entropy, system uncertainty, Gaussian particle filter, Nonlinear system. performance and lower computation cost. It is noted that all results derived from GPF are based on Gaussian assumption. However, the true PDF is an approximation of Gaussian (or Gaussian-like) due to insufficient sampling points in actual situations (C. Tsallis, 1988).
1. INTRODUCTION The most important problems arising in system analysis is probably the representation and suppression of noise, or uncertainty. It means that the real state of a system, whether measured or estimated, is rarely known perfectly because (a) limit precision of measurement and processes, and/or (b) estimates of evolving systems are based on process models that fail to include all governing parameters. The system uncertainty associated with the state can be depicted via PDF which incorporates all knowledge about the state (S. Julier & J. K. Uhlmann, 1996). In practice, the PDF evolves due to the non-equilibrium of the system state. One effective way suited for representing this situation is the use of the so-called dynamic state space (DSS) model (P. J. Harrison & C. F. Stevens, 1976).
In this paper, based on Tsallis entropy (C. Tsallis, 1988), unlike the implicit form of q-Gaussian PDFs given by H. Suyari and M. Tsukada (2005), we develop a group of explicit form of PDFs under normalization and second moment constraints. Based on these easy-use PDFs with explicit form, we propose an extension form of Gaussian particle filter called q-Gaussian particle filter (q-GPF). Two nonlinear system models are presented to demonstrate that the q-GPF using q-Gaussian PDFs in placed of the Gaussian assumption improves the GPF obviously and effectively so that a more accurate estimate of the system state is obtained.
The DSS model represents the time-varying dynamics of an unobserved state variable xn of a system as the distribution p ( xn | xn −1 ) , where n indicates time (or any other physical
evolving parameter). It is known that the estimation of distribution is a standard filtering problem. Unfortunately, closed-form analytical and optimal solutions to this filtering problem exist in a small number of cases, and in many realworld scenarios, closed-form solutions cannot be obtained (M. Vemula, M. F. Bugallo & P. M. Djuric, 2007). Recently, one application of DSS model called Gaussian particle filtering (GPF) (J. H. Kotecha & P. M. Djuric, 2003) has been proved to be powerful in tackling with the non-closed-form problems. According to Y. Wu, X. Hu, D. Hu & M. Wu (2005), the Gaussian particle filter extends the conventional Gaussian filter using Monte Carlo integration and Bayesian update rule. Compared with other particle filters, the GPF has improved 978-3-902823-45-8/2013 © IFAC
602
The paper is organized as follows. We start by presenting Tsallis entropy in Section 2 so that we can develop our results with his later on. In section 3, we show that the PDFs we calculated are all Gaussian-like and explain why they would be proper candidates to replace Gaussian assumption in Gaussian particle filtering. In Section 4, we illustrate our analysis of q-GPF. Section 5 provides some simulation results of the q-GPF as compared with GPF. Finally, Section 6 offers some conclusions. 2. TSALLIS ENTROPY In 1988, Tsallis published a powerful generalization form of Shannon entropy. Partly owing to its corresponding probability density function having a close connection to Levy process, which is ubiquitous in nature (C. Tsallis, S. V. F. Levy, A. M. C. Souza & R. Maynard, 1995), Tsallis 10.3182/20130902-3-CN-3020.00157
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
entropy is widely used in physics and field of information theory effectively (G. Wilk & Z. Wlodarczyk, 2000), for example, a theoretical model of loss systems is proposed within the framework of maximum Tsallis entropy principle (Shachi Sharma & Karmeshu, 2008), a medical image algorithm based on Tsallis entropy measures is introduced (Nicholas J. Tustision, S.P. Awate, G. Song, T.S. Cook & J. C. Gee, 2011). With the help of Tsallis entropy, a novel method to quantify brain macro-states is presented (Charles Francois, V. Latchoumane & J. Jeong, 2011). Tsallis entropy with discrete form is defined as N
∑
1− Sq =
p iq
⎛ pq ( x) − pqq ( x) ⎞ 2 + λ ⋅ x ⋅ p x + µ ⋅ p x ( ) ( ) ⎜ ⎟⎟dx (5) q q ∫−∞ ⎜ q − 1 ⎝ ⎠ ∞
The condition for this is
1 ⋅ (1 − q ⋅ pqq −1 ( x) ) + λ ⋅ x 2 + µ = 0 (6) q −1 So we get
(
(1)
i =1
q −1
⎛ 1 ⎞ 2(1 − q ) ⎟ 1 ⎝ 1 − q ⎠ 3q − 1 ⋅ pq ( x ) = 1⎞ ⎛ 1 2πσ Γ⎜ − ⎟ ⎝ 1− q 2 ⎠
pi
is the probability of the
It is easily seen that Tsallis entropy reduces to Shannon entropy as the parameter q tends to unity. Thus, according to M. Vemula, M. F. Bugallo and P. M. Djuric (2007), Tsallis entropy is a proper tool to deal with non-extensive systems having multi-fractal features because its index q appears to characterize universal class of nonadditivity.
1
⎛ 1 2(1 − q) 2 ⎞ q −1 x ⎟ ⋅ ⎜1 + 2 ⎝ 2σ 3q − 1 ⎠ where Γ(⋅) denotes Γ function and
In practice, limited sampling points are not enough to guarantee of a real Gaussian distribution that only exists in the ideal situations where sampling points are infinity. Thus, in order to acquire better accuracy, a reasonable candidate for Gaussian distribution is needed and it should be somewhat Gaussian-like for the practical use due to the robustness of the stable system, that is, because almost all system available are all under stationary state so that the PDF associated with this state should be Gaussian according to the statistical physics.
(2)
where p ( x) denotes PDF of the current state of the system. Suppose that the PDF are to obey the following constraints: (3)
Using Lagrange parameters, the PDF that maximizes the Shannon entropy is as follows:
p ( x) =
⎛ x2 ⎞ 1 exp ⎜ − 2 ⎟ 2πσ ⎝ 2σ ⎠
< q < 1 . We will also
Here, we get a group of PDFs, while parameter q in (8) changes in the range (1/3,1), from Figure 1, it is easily seen that the corresponding PDFs are all Gaussian-like, and Gaussian density is the limit case when q is very close to one.
The continuous case of Shannon entropy (C. E. Shannon, 1948) associated with one dimension is given by
p( x)dx = σ 2
3
is, < X >= 0 , since this causes no loss of generality. The pq ( x) with some typical value of q is depicted in Fig.1.
The q-Gaussian distribution proposed by H. Suyari and M. Tsukada (2005), is inconvenient in practice use due to its implicit formula, thus here, we derive the explicit form of qGaussian distributions under certain constraints.
H ( p) = − ∫ p( x) log p ( x)dx
1
(8)
assume that random variables X i have zero mean value, that
3. EXPLICIT FORM OF Q-GAUSSIAN PDFS BASED ON TSALLIS ENTROPY
2
(7)
Γ⎜
system being in state i and N denotes the number of states.
∫ p( x)dx = 1, ∫ x
)
And consequently (adjusting to satisfying the constraints (3))
where positive Tsallis parameter q measures the degree of non-extensivity of a system, that is, q denotes the complex class to which the system belongs.
1/( q −1)
⎛1 ⎞ pq ( x) = ⎜ ⋅ 1 + ( λ x 2 + µ ) ( q − 1) ⎟ ⎝q ⎠
Actually, applying q-Gaussian density instead of Gaussian assumption to implement Monte Carlo integration could lead to a more accurate result compared with results given by the traditional Gaussian particle filter algorithm, and this enhanced GPF algorithm that we developed is named qGaussian particle filter algorithm. The effectiveness of the proposed algorithm has been proved by the simulations in Section 5.
(4)
which is Gaussian distribution. Now, we wish to optimize continuous case of (1) under the constraints given by (3). By introducing Lagrange parameters, we perform the above optimization, this requires, by the calculus of variations, maximizing
603
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
closed-form analytic expression for the posterior distribution does not exist in general. Due to high-dimensional integrations, numerical solutions are not practical to implement. Recently, importance sampling-based filters have been used to update the posterior distributions. There, a distribution is represented by a weighted set of samples (or particles) from the distributions, which are propagated through the dynamic system using importance sampling to sequentially update the posterior distributions. These methods are collectively called sequential importance sampling (SIS) filters, or particle filters, and provide optimal results asymptotically in the number of particles.
Fig. 1. q-Gaussian PDF. From appearance of the curve, the distributions are similar to Gaussian density, with parameter q changing from 1/3 to 1, the q-Gaussian distributions approaches Gaussian. Furthermore, a weighted sum (or product) of the q-Gaussian PDF can be used to arbitrarily approximate closely another density just like the weighted sum of Gaussian density does (H. Suyari & M. Tsukada, 2005). So, the q-Gaussian PDF that we derived with explicit form can be a building block either for more complex cases or for some special applications to analysis of nonlinear systems.
4.2 Particle Filter The particle filter consists of recursive propagation of the weights and support points as each measurement is received i
i
N
sequentially. Let {x0:k , wk }i =s1 denote a random measure that characterizes the posterior PDF
p ( x0:k | y1:k ) , where
4. Q-GAUSSIAN PARTICLE FILTER ALGORITHM
{x0:i k , i = 0,..., N s } is a set of support points with associated
At first, we outline the DSS model that is the basis of the particle filter, and then the particle filter algorithm for system analysis is illustrated. Based on the GPF, the q-GPF algorithm is represented in the last.
weights {w0:k , i = 0,..., N s } and x0:k is the set of all state up
to time k . Suppose p ( x) ∝ π ( x ) is a PDF from which it is
4.1 DSS Model
evaluated. In addition, let x ~ q ( x ), i = 1,...N s be samples
i
difficult to draw samples but for which
The DSS model can be written as xn = f ( xn −1 , un ) (Process Equation)
(9)
[{xki , wki }iN=s1 ] = SIS [{xki −1 , wki −1}iN=s1 , yk ] • FOR i = 1: N s o Draw xki ~ q ( xk | xki −1 , yk ) o Assign the particle a weight , wki ∝ wki −1
the prior knowledge of the initial state is given by the probability distribution p ( x0 ) . Our objective is to estimate
• END FOR
recursively in time: the filtering distribution or the marginal posterior of the state at time n p( xn | y0:n ) = Cn p( xn | y0:n −1 ) p( yn | xn ) (11)
4.3 Gaussian Particle Filter
n
| y0:n −1 ) p( yn | xn )dxn
)
p ( yk | xki ) p ( xki | xki −1 ) q ( xki | xki −1 , yk )
A major disadvantage of particle filters is the computational complexity, a large part of which comes from a procedure called re-sampling. J. H. Kotecha and P. M. Djuric (2003) proposed a new particle filter called Gaussian particle filter (GPF). The GPF approximates the posterior mean and covariance of the unknown state variable using importance sampling. Unlike the SIS filters (particle filter), re-sampling is not required in the GPF. This results in a reduced complexity as compared with the SIS with re-sampling, and this is a major advantage.
where Cn is the normalizing constant given by
( ∫ p( x
can be
that are easily generated from a proposal q (⋅) called an importance density. A pseudo-code description of this algorithm is given by algorithm I [15], which is as follows
yn = h( xn , vn ) (Observation Equation) (10) where f (⋅) and h(⋅) are some known functions, and un and vn are the random noise vectors of given distributions. The process equation represents a system evolving with time n , where the system is represented by the hidden state xn , and
Cn =
π ( x)
i
−1
And the predictive distribution
p( xn +1 | y0:n ) = ∫ p ( xn +1 | xn ) p( xn | y0:n )dxn (12) Once the DSS model is linear with Gaussian noise and the prior distribution p ( x0 ) is Gaussian, the posterior and
It is worth pointing out that the GPF totally depends on the assumption that the PDF of the DSS model is Gaussian. This premise, however, is not always right in all the cases.
predictive distribution are Gaussian, and the Kalman filter supplies the mean and covariance sequentially. However, for most nonlinear models and non-Gaussian noise problems, 604
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
4.4 q-Gaussian Particle Filter Due to the limitation of the number of the sampling points, the true distribution of the model in practice is approximation of Gaussian distribution. To improve the accuracy of the GPF without adding the computational cost significantly, here, we propose an enhanced GPF named q-Gaussian particle filter based on the q-Gaussian PDFs with explicit form shown in (8) of Section III. The q-GPF algorithm is shown as follows and it is easily noted that the choice of the density function is the only difference between q-GPF and GPF; in the former, it is qGaussian PDF while in the later, it is Gaussian distribution.
q-GPF- Measurement update algorithm 1. Draw sample from the important function π ( xn | y0:n ) and denote them as {xn( j ) }Mj =1 .
3.
2
π (x
( j) n
| y0:n )
parameter
wn( j ) M
∑w
µn = ∑ w x j =1
We apply the particle filter, GPF and the q-GPF using
n = 100 particles to re-estimate X . We then test how well
( j) ( j) n n
each
Σ n = ∑ wn( j ) ( µ n − xn( j ) )( µ n − xn( j ) ) H j =1
.
j = 1,L, M , sample ( j) p( xn +1 | xn = xn ) to obtain {xn( +j )1}Mj =1 .
2.
For
3.
Compute the mean
Σ n +1 =
1 M
from
µ n +1 and covariance Σ n+1 as
( j) n +1
∑ (µ j =1
n +1
simulated
(22)
gives a better result of estimation than the others. What we not easily seen in Figure 2, is a comparison with respect to both performance and efficiency. We refer to Table I to illustrate this comparison. The q-GPF with a proper parameter q provides an appreciable improvement in the error when compared with both particle filter and GPF; for example, for the case q = 0.98 , the mean squared error
∑x M
initial
From Figure 2, we can see that each method is able to do a reasonable job in re-estimating the state X 1:T , but q-GPF
( j) M
j =1
the
1/ 2
them as { xn } j =1 .
M
fits
⎛1 T 2⎞ RMSE = ⎜ ∑ ( xn − xˆn ) ⎟ ⎝ T n =1 ⎠
Draw samples from pq ( xn ; µ n , Σ n ) and denote
1 M
procedure
The Root Mean Square Error (RMSE) is employed to evaluate the estimation accuracy, which is defined
q-GPF- Time update algorithm
µn +1 =
estimation
unobserved time series X .
.
M
1.
are known quantities. We want
θ to obtain X and Y from the time update process and measurement systems.
. Estimate the mean and covariance by M
θ = [a, b, c, d ]
simulate T = 100 observations under
( j) n
j =1
2
to obtain estimates of the state process X given observation from the data vector Y . To compare the performance of the estimation routines of the particle filter, the Gaussian particle filter and the proposed q-Gaussian particle filter, we apply these procedures on a simple nonlinear model where θ = [0.6, 0.2, 2, 0.2] . Equations (20) and (21) are used to
.
Normalize the weights as
(21)
where W ~ N (0, b ), V ~ N (0, d ) .We assume that the
p( yn | xn( j ) ) pq ( xn = xn( j ) ; µn , Σ n )
wn ( j ) = 4.
yn = c ⋅ exp( xn ) + V
Obtain the respective weights by
wn( j ) =
5.1 Simulation I Here, we present the result for one model: the short-term variations and long-term dynamics model (E. S. Schartz & J. Smith, 2000), which is popularly used in econometrics. We choose this model since it is highly nonlinear. The DSS equation for this model can be written as xn = a ⋅ xn −1 + W (20)
q-GPF Algorithm
2.
model is chosen which comes from the economy. Here, the measurement noise is considered to be a Gaussian noise. In example II, another nonlinear model, which has been used extensively in the literature for benchmarking numerical filtering techniques (M. S. Arulampalam, S. Maskell, N. Gordon & T. Clapp, 2002), is presented under the Gaussian noise. The performance of the q-GPF is compared with a conventional SIS filter (particle filter) and GPF.
(MSE), where MSE = RMSE (Simon Godsill & Tim Clapp, 2000), of GPF is 0.0049 while that of the q-GPF is 0.0015. The q-GPF with a proper q also outperforms a little the GPF in time cost; for example, for the q-GPF with q=0.98, the time cost is 0.050s while that of the GPF is 0.053s. 2
− xn( +j )1 )( µn +1 − xn( +j )1 ) H
5. SIMULATIONS The q-GPF proposed in the Section 4 is now applied to two different systems. In simulation I, a nonlinear time series 605
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
100 time samples. The performance of the q-GPF algorithm has also been compared with the particle filter and GPF. From Figure 5, it is easily seen that each method is also able to estimate the state of the system correctly, while q-GPF offers a better results of estimation than the others. Figure 6 and 7 show the MSE with 100 particles and the average MSE in filtering while using different number of particles ranging from 50 to 500. Fig.2. Simulation I: This figure compares the plot of simulated state, represented by the solid black line, to that of the estimated state, represented by the dotted asterisk, circle and square lines. These colours correspond to the particle filter, GPF and q-GPF with q=0.98 respectively. The number of particles is 100.
It is clear that the estimation accuracy of q-GPF is the best all, especially while using smaller number of particles. The processing time of 1000 iterations when using the number of particle is 50 is given in Table II. It is also easily noted that q-GPF has cost nearly equal processing time compared with the GPF.
The MSE of each filter was subsequently computed. Fig. 3 shows the MSE of each filter with 100 particles for 50 Monte Carlo runs. The comparison of the average MSE over the 50 times samples among filters, in terms of the number of particles, is shown in Fig. 4. The average MSE was also compared for the number of particles 50, 100 and 500 shown in Fig. 4. From Fig. 4, it is easily noted that even for small number of particles the qGPF performs significantly better than GPF. An increase in the number of particles reduces the MSE further for the qGPF. Table I: Simulation I: Comparison of the particle filters Algorithm PF MSE Time/s GPF MSE Time/s q-GPF MSE Time/s
1
2
3
0.0087 0.039
0.0085 0.038
0.0091 0.039
0.0049 0.053 q=0.98 0.0015 0.050
0.0050 0.050 q=0.97 0.0011 0.048
0.0070 0.049 q=0.96 0.0027 0.048
Figure 3: Simulation I: MSE for 50 Monte Carlo runs (the number of particles is 100), with q=0.98 for q-GPF. Table II: Simulation II: Performance comparison (1000 iterations) Algorithm Processing time (s)
Particle filter 0.1567
GPF 0.5631
qGPF(q=0.85) 0.5624
The performance of the q-GPF becomes very close for the number of particles 500 and 100. In conclusion, the performance of the q-GPF appears to be better than the other ones taken here for comparison. 5.2 Simulation II In this section we apply the q-GPF to estimate the state of another nonlinear time series. The continuous time dynamics of the system are
xt = 0.5 xt −1 + 25
xt −1 + 8cos(1.2(t − 1)) + ut (23) 1 + xt2−1
yt = where
xt2 + vt 20
Figure 4: Simulation I: Average MSE over time samples with the number of particles 50, 100 and 500 (for 50 Monte Carlo runs). The q of q-GPF is 0.98.
(24)
Gaussian
distribution noise ut ~ N (0,1), vt ~ N (0,1) . The simulation carried out for 606
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
REFERENCES C. Tsallis (1988), Possible generalization of BoltzmannGibbs statistics, Journal of Statistical Physics, vol. 51 (1/2), pp. 479-487. C. Tsallis, S. V. F. Levy, A. M. C. Souza, and R. Maynard (1995), Statistical mechanical foundation of the ubiquity of levy distributions in nature, Physical Review Letters, vol.75 (20), pp.3589-3593. Charles Francois, V. Latchoumane, and J. Jeong (2011), Quantification of brain macrostates using dynamical nonstationary of physiological time series, IEEE Trans. Biomedical Engineering, vol. 58 (4), pp.1084-1093. C. Tsallis (2009), An introduction to nonextensive statistical mehcanics, Springer Science +Business Media LLC, NewYork. C. E. Shannon (1948), A mathematical theory of communication, The Bell System Technical Journal, vol. 27, pp.379-423, 623-656. E. S. Schartz, J. Smith (2000), Short-term variations and long-term dynamics in commodity prices, Management Science, vol.46 (7), pp.893-911. G. Wilk and Z. Wlodarczyk (2000), Interpretation of the nonextensity parameter q in some applications of Tsallis statistics and Levy distributions, Physical Review Letter, vol. 84 (13), 2770-2773 H. Suyari, M. Tsukada (2005), Law of error in Tsallis statistics, IEEE Trans. Information Theory, vol. 51 (2), pp. 753-757. J. H. Kotecha, P. M. Djuric (2003), Gaussian particle filtering, IEEE Trans. Signal Processing, vol. 51 (10), pp. 25922601 M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp (2002), A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking, IEEE Trans. Signal Processing, vol 50 (2), pp.174-188. M. Vemula, M. F. Bugallo and P. M. Djuric (2007), Performance comparison of Gaussian-based filters using information measures, IEEE Signal Processing Lett., vol. 14 (12), pp. 1020-1023. Nicholas J. Tustision, S.P. Awate, G. Song, T.S. Cook, and J. C. Gee (2011), Point set registration using HavrdaCharvat-Tsallis entropy measures, IEEE Trans. Medical Imaging, vol. 30 (2), pp.451-460. P. J. Harrison and C. F. Stevens (1976), Bayesian forecasting (with discussion), J. R. Dtatist. Soc. B, vol. 38, pp. 205247. Shachi Sharma and Karmeshu (2008), Bimodal packet distribution in loss systems using maximum Tsallis entropy principle, IEEE Trans. Communications, vol. 56 (9), pp.1530-1535. Simon Godsill and Tim Clapp (2000), Improvement strategies for Monte Carlo Practice Filters, pp. 139-158, Springer-Verlag. S. Julier and J. K. Uhlmann (1996), A general method for approximation nonlinear transformations of probalility distributions, http://www.Robots.ox .ac.uk/siju/work/work.html. Y. Wu, X. Hu, D. Hu and M. Wu (2005), Comment on Gaussian particle filtering, IEEE Trans. Signal Processing, vol. 53 (8), pp. 3350-3351.
Figure 5: Simulation II: This figure compares the plot of simulated state, represented by the solid black line, to that of the estimated state, represented by the dotted asterisk, circle and square lines. These colours correspond to the particle filter, GPF and q-GPF with q=0.85 respectively. The number of particles is 100. 6. CONCLUSIONS This paper develops a group of Gaussian-like PDFs with explicit forms called q-Gaussian PDFs based on Tsallis entropy. With the aid of these q-Gaussian PDFs, an effective application to nonlinear system state estimation using the extension of Gaussian particle filter, called the q-Gaussian particle filter, is discussed. Simulations on estimating the system state demonstrate that the q-Gaussian particle filter is superior to traditional Gaussian particle filter in accuracy.
Figure 6: Simulation II: MSE for 50 Monte Carlo runs (the number of particles is 100), with q=0.85 for q-GPF.
Figure 7: Simulation II: Average MSE over time samples with the number of particles 50, 100 and 500 (for 50 Monte Carlo runs). The q of q-GPF is 0.85.
607
q-Gaussian Density Model and its Application to State Estimation of Nonlinear Systems Xifeng Li*, Yongle Xie.** * University of Electronic Science and Technology of China, Chengdu, Sichuan, 611731 China (Tel: +86-13458530206; e-mail: xfengl@ hotmail.com). ** University of Electronic Science and Technology of China, Chengdu, Sichuan, 611731 China (e-mail: [email protected]) Abstract: Probability density function (PDF) plays a vital role in system analysis involving stochastic factors. A good estimate of true PDF conditioned under certain performance criterion could help acquire more information of the system. With help of the new information, many features of the system that we are concerning can be revealed effectively, especially for nonlinear non-Gaussian stochastic systems. In this paper, based on the Tsallis entropy, we derive a class of PDFs with explicit form called q-Gaussian PDFs. These PDFs have a parameter that indicates the fractal feature of the system. Based on the explicit form of q-Gaussian PDFs, we propose an extension of Gaussian particle filter (GPF) called q-Gaussian particle filter (q-GPF). The experimental results show that the q-GPF is a more effective method to estimate the state of nonlinear stochastic system compared with the GPF. Keywords: Tsallis entropy, system uncertainty, Gaussian particle filter, Nonlinear system. performance and lower computation cost. It is noted that all results derived from GPF are based on Gaussian assumption. However, the true PDF is an approximation of Gaussian (or Gaussian-like) due to insufficient sampling points in actual situations (C. Tsallis, 1988).
1. INTRODUCTION The most important problems arising in system analysis is probably the representation and suppression of noise, or uncertainty. It means that the real state of a system, whether measured or estimated, is rarely known perfectly because (a) limit precision of measurement and processes, and/or (b) estimates of evolving systems are based on process models that fail to include all governing parameters. The system uncertainty associated with the state can be depicted via PDF which incorporates all knowledge about the state (S. Julier & J. K. Uhlmann, 1996). In practice, the PDF evolves due to the non-equilibrium of the system state. One effective way suited for representing this situation is the use of the so-called dynamic state space (DSS) model (P. J. Harrison & C. F. Stevens, 1976).
In this paper, based on Tsallis entropy (C. Tsallis, 1988), unlike the implicit form of q-Gaussian PDFs given by H. Suyari and M. Tsukada (2005), we develop a group of explicit form of PDFs under normalization and second moment constraints. Based on these easy-use PDFs with explicit form, we propose an extension form of Gaussian particle filter called q-Gaussian particle filter (q-GPF). Two nonlinear system models are presented to demonstrate that the q-GPF using q-Gaussian PDFs in placed of the Gaussian assumption improves the GPF obviously and effectively so that a more accurate estimate of the system state is obtained.
The DSS model represents the time-varying dynamics of an unobserved state variable xn of a system as the distribution p ( xn | xn −1 ) , where n indicates time (or any other physical
evolving parameter). It is known that the estimation of distribution is a standard filtering problem. Unfortunately, closed-form analytical and optimal solutions to this filtering problem exist in a small number of cases, and in many realworld scenarios, closed-form solutions cannot be obtained (M. Vemula, M. F. Bugallo & P. M. Djuric, 2007). Recently, one application of DSS model called Gaussian particle filtering (GPF) (J. H. Kotecha & P. M. Djuric, 2003) has been proved to be powerful in tackling with the non-closed-form problems. According to Y. Wu, X. Hu, D. Hu & M. Wu (2005), the Gaussian particle filter extends the conventional Gaussian filter using Monte Carlo integration and Bayesian update rule. Compared with other particle filters, the GPF has improved 978-3-902823-45-8/2013 © IFAC
602
The paper is organized as follows. We start by presenting Tsallis entropy in Section 2 so that we can develop our results with his later on. In section 3, we show that the PDFs we calculated are all Gaussian-like and explain why they would be proper candidates to replace Gaussian assumption in Gaussian particle filtering. In Section 4, we illustrate our analysis of q-GPF. Section 5 provides some simulation results of the q-GPF as compared with GPF. Finally, Section 6 offers some conclusions. 2. TSALLIS ENTROPY In 1988, Tsallis published a powerful generalization form of Shannon entropy. Partly owing to its corresponding probability density function having a close connection to Levy process, which is ubiquitous in nature (C. Tsallis, S. V. F. Levy, A. M. C. Souza & R. Maynard, 1995), Tsallis 10.3182/20130902-3-CN-3020.00157
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
entropy is widely used in physics and field of information theory effectively (G. Wilk & Z. Wlodarczyk, 2000), for example, a theoretical model of loss systems is proposed within the framework of maximum Tsallis entropy principle (Shachi Sharma & Karmeshu, 2008), a medical image algorithm based on Tsallis entropy measures is introduced (Nicholas J. Tustision, S.P. Awate, G. Song, T.S. Cook & J. C. Gee, 2011). With the help of Tsallis entropy, a novel method to quantify brain macro-states is presented (Charles Francois, V. Latchoumane & J. Jeong, 2011). Tsallis entropy with discrete form is defined as N
∑
1− Sq =
p iq
⎛ pq ( x) − pqq ( x) ⎞ 2 + λ ⋅ x ⋅ p x + µ ⋅ p x ( ) ( ) ⎜ ⎟⎟dx (5) q q ∫−∞ ⎜ q − 1 ⎝ ⎠ ∞
The condition for this is
1 ⋅ (1 − q ⋅ pqq −1 ( x) ) + λ ⋅ x 2 + µ = 0 (6) q −1 So we get
(
(1)
i =1
q −1
⎛ 1 ⎞ 2(1 − q ) ⎟ 1 ⎝ 1 − q ⎠ 3q − 1 ⋅ pq ( x ) = 1⎞ ⎛ 1 2πσ Γ⎜ − ⎟ ⎝ 1− q 2 ⎠
pi
is the probability of the
It is easily seen that Tsallis entropy reduces to Shannon entropy as the parameter q tends to unity. Thus, according to M. Vemula, M. F. Bugallo and P. M. Djuric (2007), Tsallis entropy is a proper tool to deal with non-extensive systems having multi-fractal features because its index q appears to characterize universal class of nonadditivity.
1
⎛ 1 2(1 − q) 2 ⎞ q −1 x ⎟ ⋅ ⎜1 + 2 ⎝ 2σ 3q − 1 ⎠ where Γ(⋅) denotes Γ function and
In practice, limited sampling points are not enough to guarantee of a real Gaussian distribution that only exists in the ideal situations where sampling points are infinity. Thus, in order to acquire better accuracy, a reasonable candidate for Gaussian distribution is needed and it should be somewhat Gaussian-like for the practical use due to the robustness of the stable system, that is, because almost all system available are all under stationary state so that the PDF associated with this state should be Gaussian according to the statistical physics.
(2)
where p ( x) denotes PDF of the current state of the system. Suppose that the PDF are to obey the following constraints: (3)
Using Lagrange parameters, the PDF that maximizes the Shannon entropy is as follows:
p ( x) =
⎛ x2 ⎞ 1 exp ⎜ − 2 ⎟ 2πσ ⎝ 2σ ⎠
< q < 1 . We will also
Here, we get a group of PDFs, while parameter q in (8) changes in the range (1/3,1), from Figure 1, it is easily seen that the corresponding PDFs are all Gaussian-like, and Gaussian density is the limit case when q is very close to one.
The continuous case of Shannon entropy (C. E. Shannon, 1948) associated with one dimension is given by
p( x)dx = σ 2
3
is, < X >= 0 , since this causes no loss of generality. The pq ( x) with some typical value of q is depicted in Fig.1.
The q-Gaussian distribution proposed by H. Suyari and M. Tsukada (2005), is inconvenient in practice use due to its implicit formula, thus here, we derive the explicit form of qGaussian distributions under certain constraints.
H ( p) = − ∫ p( x) log p ( x)dx
1
(8)
assume that random variables X i have zero mean value, that
3. EXPLICIT FORM OF Q-GAUSSIAN PDFS BASED ON TSALLIS ENTROPY
2
(7)
Γ⎜
system being in state i and N denotes the number of states.
∫ p( x)dx = 1, ∫ x
)
And consequently (adjusting to satisfying the constraints (3))
where positive Tsallis parameter q measures the degree of non-extensivity of a system, that is, q denotes the complex class to which the system belongs.
1/( q −1)
⎛1 ⎞ pq ( x) = ⎜ ⋅ 1 + ( λ x 2 + µ ) ( q − 1) ⎟ ⎝q ⎠
Actually, applying q-Gaussian density instead of Gaussian assumption to implement Monte Carlo integration could lead to a more accurate result compared with results given by the traditional Gaussian particle filter algorithm, and this enhanced GPF algorithm that we developed is named qGaussian particle filter algorithm. The effectiveness of the proposed algorithm has been proved by the simulations in Section 5.
(4)
which is Gaussian distribution. Now, we wish to optimize continuous case of (1) under the constraints given by (3). By introducing Lagrange parameters, we perform the above optimization, this requires, by the calculus of variations, maximizing
603
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
closed-form analytic expression for the posterior distribution does not exist in general. Due to high-dimensional integrations, numerical solutions are not practical to implement. Recently, importance sampling-based filters have been used to update the posterior distributions. There, a distribution is represented by a weighted set of samples (or particles) from the distributions, which are propagated through the dynamic system using importance sampling to sequentially update the posterior distributions. These methods are collectively called sequential importance sampling (SIS) filters, or particle filters, and provide optimal results asymptotically in the number of particles.
Fig. 1. q-Gaussian PDF. From appearance of the curve, the distributions are similar to Gaussian density, with parameter q changing from 1/3 to 1, the q-Gaussian distributions approaches Gaussian. Furthermore, a weighted sum (or product) of the q-Gaussian PDF can be used to arbitrarily approximate closely another density just like the weighted sum of Gaussian density does (H. Suyari & M. Tsukada, 2005). So, the q-Gaussian PDF that we derived with explicit form can be a building block either for more complex cases or for some special applications to analysis of nonlinear systems.
4.2 Particle Filter The particle filter consists of recursive propagation of the weights and support points as each measurement is received i
i
N
sequentially. Let {x0:k , wk }i =s1 denote a random measure that characterizes the posterior PDF
p ( x0:k | y1:k ) , where
4. Q-GAUSSIAN PARTICLE FILTER ALGORITHM
{x0:i k , i = 0,..., N s } is a set of support points with associated
At first, we outline the DSS model that is the basis of the particle filter, and then the particle filter algorithm for system analysis is illustrated. Based on the GPF, the q-GPF algorithm is represented in the last.
weights {w0:k , i = 0,..., N s } and x0:k is the set of all state up
to time k . Suppose p ( x) ∝ π ( x ) is a PDF from which it is
4.1 DSS Model
evaluated. In addition, let x ~ q ( x ), i = 1,...N s be samples
i
difficult to draw samples but for which
The DSS model can be written as xn = f ( xn −1 , un ) (Process Equation)
(9)
[{xki , wki }iN=s1 ] = SIS [{xki −1 , wki −1}iN=s1 , yk ] • FOR i = 1: N s o Draw xki ~ q ( xk | xki −1 , yk ) o Assign the particle a weight , wki ∝ wki −1
the prior knowledge of the initial state is given by the probability distribution p ( x0 ) . Our objective is to estimate
• END FOR
recursively in time: the filtering distribution or the marginal posterior of the state at time n p( xn | y0:n ) = Cn p( xn | y0:n −1 ) p( yn | xn ) (11)
4.3 Gaussian Particle Filter
n
| y0:n −1 ) p( yn | xn )dxn
)
p ( yk | xki ) p ( xki | xki −1 ) q ( xki | xki −1 , yk )
A major disadvantage of particle filters is the computational complexity, a large part of which comes from a procedure called re-sampling. J. H. Kotecha and P. M. Djuric (2003) proposed a new particle filter called Gaussian particle filter (GPF). The GPF approximates the posterior mean and covariance of the unknown state variable using importance sampling. Unlike the SIS filters (particle filter), re-sampling is not required in the GPF. This results in a reduced complexity as compared with the SIS with re-sampling, and this is a major advantage.
where Cn is the normalizing constant given by
( ∫ p( x
can be
that are easily generated from a proposal q (⋅) called an importance density. A pseudo-code description of this algorithm is given by algorithm I [15], which is as follows
yn = h( xn , vn ) (Observation Equation) (10) where f (⋅) and h(⋅) are some known functions, and un and vn are the random noise vectors of given distributions. The process equation represents a system evolving with time n , where the system is represented by the hidden state xn , and
Cn =
π ( x)
i
−1
And the predictive distribution
p( xn +1 | y0:n ) = ∫ p ( xn +1 | xn ) p( xn | y0:n )dxn (12) Once the DSS model is linear with Gaussian noise and the prior distribution p ( x0 ) is Gaussian, the posterior and
It is worth pointing out that the GPF totally depends on the assumption that the PDF of the DSS model is Gaussian. This premise, however, is not always right in all the cases.
predictive distribution are Gaussian, and the Kalman filter supplies the mean and covariance sequentially. However, for most nonlinear models and non-Gaussian noise problems, 604
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
4.4 q-Gaussian Particle Filter Due to the limitation of the number of the sampling points, the true distribution of the model in practice is approximation of Gaussian distribution. To improve the accuracy of the GPF without adding the computational cost significantly, here, we propose an enhanced GPF named q-Gaussian particle filter based on the q-Gaussian PDFs with explicit form shown in (8) of Section III. The q-GPF algorithm is shown as follows and it is easily noted that the choice of the density function is the only difference between q-GPF and GPF; in the former, it is qGaussian PDF while in the later, it is Gaussian distribution.
q-GPF- Measurement update algorithm 1. Draw sample from the important function π ( xn | y0:n ) and denote them as {xn( j ) }Mj =1 .
3.
2
π (x
( j) n
| y0:n )
parameter
wn( j ) M
∑w
µn = ∑ w x j =1
We apply the particle filter, GPF and the q-GPF using
n = 100 particles to re-estimate X . We then test how well
( j) ( j) n n
each
Σ n = ∑ wn( j ) ( µ n − xn( j ) )( µ n − xn( j ) ) H j =1
.
j = 1,L, M , sample ( j) p( xn +1 | xn = xn ) to obtain {xn( +j )1}Mj =1 .
2.
For
3.
Compute the mean
Σ n +1 =
1 M
from
µ n +1 and covariance Σ n+1 as
( j) n +1
∑ (µ j =1
n +1
simulated
(22)
gives a better result of estimation than the others. What we not easily seen in Figure 2, is a comparison with respect to both performance and efficiency. We refer to Table I to illustrate this comparison. The q-GPF with a proper parameter q provides an appreciable improvement in the error when compared with both particle filter and GPF; for example, for the case q = 0.98 , the mean squared error
∑x M
initial
From Figure 2, we can see that each method is able to do a reasonable job in re-estimating the state X 1:T , but q-GPF
( j) M
j =1
the
1/ 2
them as { xn } j =1 .
M
fits
⎛1 T 2⎞ RMSE = ⎜ ∑ ( xn − xˆn ) ⎟ ⎝ T n =1 ⎠
Draw samples from pq ( xn ; µ n , Σ n ) and denote
1 M
procedure
The Root Mean Square Error (RMSE) is employed to evaluate the estimation accuracy, which is defined
q-GPF- Time update algorithm
µn +1 =
estimation
unobserved time series X .
.
M
1.
are known quantities. We want
θ to obtain X and Y from the time update process and measurement systems.
. Estimate the mean and covariance by M
θ = [a, b, c, d ]
simulate T = 100 observations under
( j) n
j =1
2
to obtain estimates of the state process X given observation from the data vector Y . To compare the performance of the estimation routines of the particle filter, the Gaussian particle filter and the proposed q-Gaussian particle filter, we apply these procedures on a simple nonlinear model where θ = [0.6, 0.2, 2, 0.2] . Equations (20) and (21) are used to
.
Normalize the weights as
(21)
where W ~ N (0, b ), V ~ N (0, d ) .We assume that the
p( yn | xn( j ) ) pq ( xn = xn( j ) ; µn , Σ n )
wn ( j ) = 4.
yn = c ⋅ exp( xn ) + V
Obtain the respective weights by
wn( j ) =
5.1 Simulation I Here, we present the result for one model: the short-term variations and long-term dynamics model (E. S. Schartz & J. Smith, 2000), which is popularly used in econometrics. We choose this model since it is highly nonlinear. The DSS equation for this model can be written as xn = a ⋅ xn −1 + W (20)
q-GPF Algorithm
2.
model is chosen which comes from the economy. Here, the measurement noise is considered to be a Gaussian noise. In example II, another nonlinear model, which has been used extensively in the literature for benchmarking numerical filtering techniques (M. S. Arulampalam, S. Maskell, N. Gordon & T. Clapp, 2002), is presented under the Gaussian noise. The performance of the q-GPF is compared with a conventional SIS filter (particle filter) and GPF.
(MSE), where MSE = RMSE (Simon Godsill & Tim Clapp, 2000), of GPF is 0.0049 while that of the q-GPF is 0.0015. The q-GPF with a proper q also outperforms a little the GPF in time cost; for example, for the q-GPF with q=0.98, the time cost is 0.050s while that of the GPF is 0.053s. 2
− xn( +j )1 )( µn +1 − xn( +j )1 ) H
5. SIMULATIONS The q-GPF proposed in the Section 4 is now applied to two different systems. In simulation I, a nonlinear time series 605
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
100 time samples. The performance of the q-GPF algorithm has also been compared with the particle filter and GPF. From Figure 5, it is easily seen that each method is also able to estimate the state of the system correctly, while q-GPF offers a better results of estimation than the others. Figure 6 and 7 show the MSE with 100 particles and the average MSE in filtering while using different number of particles ranging from 50 to 500. Fig.2. Simulation I: This figure compares the plot of simulated state, represented by the solid black line, to that of the estimated state, represented by the dotted asterisk, circle and square lines. These colours correspond to the particle filter, GPF and q-GPF with q=0.98 respectively. The number of particles is 100.
It is clear that the estimation accuracy of q-GPF is the best all, especially while using smaller number of particles. The processing time of 1000 iterations when using the number of particle is 50 is given in Table II. It is also easily noted that q-GPF has cost nearly equal processing time compared with the GPF.
The MSE of each filter was subsequently computed. Fig. 3 shows the MSE of each filter with 100 particles for 50 Monte Carlo runs. The comparison of the average MSE over the 50 times samples among filters, in terms of the number of particles, is shown in Fig. 4. The average MSE was also compared for the number of particles 50, 100 and 500 shown in Fig. 4. From Fig. 4, it is easily noted that even for small number of particles the qGPF performs significantly better than GPF. An increase in the number of particles reduces the MSE further for the qGPF. Table I: Simulation I: Comparison of the particle filters Algorithm PF MSE Time/s GPF MSE Time/s q-GPF MSE Time/s
1
2
3
0.0087 0.039
0.0085 0.038
0.0091 0.039
0.0049 0.053 q=0.98 0.0015 0.050
0.0050 0.050 q=0.97 0.0011 0.048
0.0070 0.049 q=0.96 0.0027 0.048
Figure 3: Simulation I: MSE for 50 Monte Carlo runs (the number of particles is 100), with q=0.98 for q-GPF. Table II: Simulation II: Performance comparison (1000 iterations) Algorithm Processing time (s)
Particle filter 0.1567
GPF 0.5631
qGPF(q=0.85) 0.5624
The performance of the q-GPF becomes very close for the number of particles 500 and 100. In conclusion, the performance of the q-GPF appears to be better than the other ones taken here for comparison. 5.2 Simulation II In this section we apply the q-GPF to estimate the state of another nonlinear time series. The continuous time dynamics of the system are
xt = 0.5 xt −1 + 25
xt −1 + 8cos(1.2(t − 1)) + ut (23) 1 + xt2−1
yt = where
xt2 + vt 20
Figure 4: Simulation I: Average MSE over time samples with the number of particles 50, 100 and 500 (for 50 Monte Carlo runs). The q of q-GPF is 0.98.
(24)
Gaussian
distribution noise ut ~ N (0,1), vt ~ N (0,1) . The simulation carried out for 606
IFAC ICONS 2013 September 2-4, 2013. Chengdu, China
REFERENCES C. Tsallis (1988), Possible generalization of BoltzmannGibbs statistics, Journal of Statistical Physics, vol. 51 (1/2), pp. 479-487. C. Tsallis, S. V. F. Levy, A. M. C. Souza, and R. Maynard (1995), Statistical mechanical foundation of the ubiquity of levy distributions in nature, Physical Review Letters, vol.75 (20), pp.3589-3593. Charles Francois, V. Latchoumane, and J. Jeong (2011), Quantification of brain macrostates using dynamical nonstationary of physiological time series, IEEE Trans. Biomedical Engineering, vol. 58 (4), pp.1084-1093. C. Tsallis (2009), An introduction to nonextensive statistical mehcanics, Springer Science +Business Media LLC, NewYork. C. E. Shannon (1948), A mathematical theory of communication, The Bell System Technical Journal, vol. 27, pp.379-423, 623-656. E. S. Schartz, J. Smith (2000), Short-term variations and long-term dynamics in commodity prices, Management Science, vol.46 (7), pp.893-911. G. Wilk and Z. Wlodarczyk (2000), Interpretation of the nonextensity parameter q in some applications of Tsallis statistics and Levy distributions, Physical Review Letter, vol. 84 (13), 2770-2773 H. Suyari, M. Tsukada (2005), Law of error in Tsallis statistics, IEEE Trans. Information Theory, vol. 51 (2), pp. 753-757. J. H. Kotecha, P. M. Djuric (2003), Gaussian particle filtering, IEEE Trans. Signal Processing, vol. 51 (10), pp. 25922601 M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp (2002), A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking, IEEE Trans. Signal Processing, vol 50 (2), pp.174-188. M. Vemula, M. F. Bugallo and P. M. Djuric (2007), Performance comparison of Gaussian-based filters using information measures, IEEE Signal Processing Lett., vol. 14 (12), pp. 1020-1023. Nicholas J. Tustision, S.P. Awate, G. Song, T.S. Cook, and J. C. Gee (2011), Point set registration using HavrdaCharvat-Tsallis entropy measures, IEEE Trans. Medical Imaging, vol. 30 (2), pp.451-460. P. J. Harrison and C. F. Stevens (1976), Bayesian forecasting (with discussion), J. R. Dtatist. Soc. B, vol. 38, pp. 205247. Shachi Sharma and Karmeshu (2008), Bimodal packet distribution in loss systems using maximum Tsallis entropy principle, IEEE Trans. Communications, vol. 56 (9), pp.1530-1535. Simon Godsill and Tim Clapp (2000), Improvement strategies for Monte Carlo Practice Filters, pp. 139-158, Springer-Verlag. S. Julier and J. K. Uhlmann (1996), A general method for approximation nonlinear transformations of probalility distributions, http://www.Robots.ox .ac.uk/siju/work/work.html. Y. Wu, X. Hu, D. Hu and M. Wu (2005), Comment on Gaussian particle filtering, IEEE Trans. Signal Processing, vol. 53 (8), pp. 3350-3351.
Figure 5: Simulation II: This figure compares the plot of simulated state, represented by the solid black line, to that of the estimated state, represented by the dotted asterisk, circle and square lines. These colours correspond to the particle filter, GPF and q-GPF with q=0.85 respectively. The number of particles is 100. 6. CONCLUSIONS This paper develops a group of Gaussian-like PDFs with explicit forms called q-Gaussian PDFs based on Tsallis entropy. With the aid of these q-Gaussian PDFs, an effective application to nonlinear system state estimation using the extension of Gaussian particle filter, called the q-Gaussian particle filter, is discussed. Simulations on estimating the system state demonstrate that the q-Gaussian particle filter is superior to traditional Gaussian particle filter in accuracy.
Figure 6: Simulation II: MSE for 50 Monte Carlo runs (the number of particles is 100), with q=0.85 for q-GPF.
Figure 7: Simulation II: Average MSE over time samples with the number of particles 50, 100 and 500 (for 50 Monte Carlo runs). The q of q-GPF is 0.85.
607