To to complete the specification of the statistical model (1)–(3) we need to model the additive failure componentssk. Let
hk(·) =
hk,1(·) ...
hk,m(·)
, (6)
wherehk,1(·),...,hk,m(·)is a partition of the observations such that the additive errors have a compatible partition
sk =
sk,1
...
sk,m
, (7)
where sk,1,...,sk,m are mutually independent. This kind of model would be natural for a situation in which we have pseudorange obser-vations from several satellites, and depending whether the individual signals between the satellites and the receiver is unobstructed or not, there may be an additional error component present. Of course, we would have to assume that the visibility to one satellite is inde-pendent of the visibility to any other satellite. To model the presence of the sensor errors in the observations at timek, we use indicator variablesλk,l ∈ {0,1}so that
sk,l =λk,lrk,l, (8) whererk,l is the magnitude of the error. We use two models for the indicator variables throughout the thesis. The first model for the indicator variable is the Bernoulli-distribution
P(λk,l =1) =θk =1−P(λk,l =0) (9) whereθk is the probability of an additive sensor error being present inyk,l. The notationP(·)is used for the probability mass functions of discrete random variables.
In the case of dynamic systems it is sometimes reasonable to model the indicator variable as a Markov chain[68]with transition probab-ilities
P(λk,l =j |λk−1,l =i) =θj i, (10) P(λ0,l =1) =θ0. (11)
Note that in both models, we assume for simplicity that the probabil-ities of the indicator variable are the same for each of the elements of the observation vector.
We use a linear state transition model for the sensor error size
rk+1,l =φk,lrk,l +εk+1,l (12) r0,l �N�r0|0,l,P0r|0,l�
, (13)
where εk+1,l �N�0,Σk+1,l�
is a Gaussian white noise process. The choiceφk,l =0 results in the Gaussian white noise process for the er-ror size that could be used to model outliers in the observation noise [55],[13]. The state transition model with the coefficientφk,l =1 is the Gaussian random walk that can be used to model the evolution of the multipath bias in GPS pseudorange measurements[23],[24], [P5].
The system (1)–(3) can be expressed using the sensor error models as
�xk+1
rk+1
�
=
�Fk 0 0 Φk
��xk
rk
� +
�wk
εk+1
�
(14) yk =hk(xk) + Λkrk+vk (15)
�x0
r0
�
�N
��x0|0
r0|0
� ,
�P0|0 0 0 P0r|0
��
(16)
where we have defined
Φk :=
φk,1
...
φk,m
,
Λk :=
λk,1Iny1
...
λk,mInym
,
P0r|0=
P0r|0,1
...
P0r|0,m
, r0|0:=�
r0T|0,1 ··· r0T|0,m�T
, εk :=�
εTk,1 ··· εTk,m�T
.
The augmented system (14)–(16) is a standard state-space system if the indicator variables would be known. Naturally one could consider the indicator variable as a part of the state and formulate a larger non-linear and non-Gaussian system, but we handleΛk separately due to the approximation methods discussed later on.
Furthermore, in the caseφi,j =0, ∀i,j we can write the system as xk+1=Fkxk+wk (17)
yk =hk(xk) +vk(Λk) (18) x0�N�x0|0,P0|0�
, (19)
where the observation noise givenΛk is a Gaussian white noise pro-cess
vk(Λk):= Λkrk+vk �N(0,Rk(Λk)) (20) Rk(Λk):=Rk+ ΛkΣkΛTk (21) The system (17)–(19) is convenient if we are considering the additive sensor errors simply as nuisance parameters causing the deteriora-tion of the observadeteriora-tion quality.
Heavy-tailed distributions
Another approach for modeling faulty data would be to use non-Gaussian distributions for the error components. In the robust fil-ter design similar to M-estimation [31, 27], very large outliers are modeled with the noise process
sk+vk �� (22)
where� is a member of a family of distributions with fat tails. The robust filter design is then based on finding a state estimator that performs the best should the noise have any of the distributions in the family� [50],[51],[56],[P1].
Outliers in the observation noise process can also be modeled simply using a heavy-tailed distribution such as Student-t distribution for the observation noise instead of the Gaussian noise[69, 1]. A sample drawn from a heavy-tailed distribution would result in more realiza-tions that are far away from the bulk of the data.
Hierarchical modeling of varying environment
A Bayesian approach to the state-space estimation problem enables us to use hierarchical modeling to describe more complex real world phenomena. One application of hierarchical modeling is to describe the probability of the presence of the additive sensor error (9) as a time-evolving parameter. Given the model we can solve it jointly with the state[P5]. The approach is reasonable because the probability of the faulty observation is often dependent on the surrounding environments that change gradually when the MS is moving with a reasonably low velocity.
Now the probability mass function ofλk,l is defined with a hierarch-ical model depending on the time-varying unknown variableθk as follows
P(λk,l =0|θk) =1−P(λk,l =1|θk) =1−θk. (23) In the state-transition model for the uncertainty parameter θk we have to take into account thatθk∈[0,1]. The parameter is modeled as
a Markov process, and we take the densityp(θk+1|θk)to be unimodal, with the mode near to the value ofθk. A probability density fulfilling these criteria would be
beta(ξ|α,β) = Γ(α+β)
Γ(α)Γ(β)ξα−1(1−ξ)β−1, (24) that is a beta density with parametersαandβ evaluated atξ. The mode and variance of a beta distributed random variable are
mode(ξ) = α−1
α+β −2, V(ξ) = αβ
(α+β)2(α+β−1). (25) The beta density function is unimodal whenα,β >1, and the vari-ance of a beta distributed random variable depends onαandβ, with variance→0 asα,β → ∞.
The model probability is modeled as having the state-transition dens-ity
p(θk+1|θk,S) =beta(θk+1|θk(S−2) +1,(1−θk)S+2θk−1), (26) whereS is a tuning parameter. The state transition density is illus-trated in Figure 1. In Figure 2 we have drawn a few example sample paths of the process with a corresponding sample path ofλk |θk.
0 0.2 0.4 0.6 0.8 1
θk+1
p(θk+1|0.3,10)
p(θk+1|0.5,100)
p(θk+1|0.7,500)
❅❅❘
❅❅❘
��
✠
Figure 1:State-transition densities with different parameter values.
The mode and variance of (26) are
mode(θk+1|θk,S) =θk, V(θk+1|θk,S) = (1−θk)θk
S−1 . (27)
0 10 20 30 40 50 60 70 80 90 100 0
0.2 0.4 0.6 0.8 1
0 10 20 30 40 50 60 70 80 90 100
0 0.2 0.4 0.6 0.8 1 θk+1|θk,100
λk|θk
k
Figure 2:Five sample paths of the processθk+1|θk,100 and a sample path ofλk |θk corresponding to the red coloredθk+1|θk,100.
The most likely value of the model uncertaintyθk+1corresponds to the model uncertainty at the previous time stepθk, and increasing the value of the tuning parameterSreduces the probability ofθk+1
deviating significantly fromθk. The hierarchical state-space model for the system (17)–(19) is illustrated by the directed acyclic graph (DAG) in Figure 3.
3 Estimation methods
From the whiteness and the mutual independence assumptions of the noise processes it follows thatxk,rk andΛk are Markov processes
p(xk,rk,Λk |x0:k−1,r0:k−1,Λ0:k−1)
=p(xk |xk−1)p(rk |rk−1)P(Λk|Λk−1), (28)
x0
θ0
x1
y1
Λ1
θ1
x2
y2
Λ2
θ2
xk
yk
Λk
θk
Figure 3:DAG of the hierarchical state-space model
and therefore
p(x0:k,r0:k,Λ0:k)
=p(x0)p(r0)P(Λ0)
�k i=1
p(xi |xi−1)p(ri |ri−1)P(Λi |Λi−1), (29) where
p(xk |xk−1) =pwk(xk−Fk−1xk−1) =N(xk |Fk−1xk−1,Qk−1) (30) p(rk |rk−1) =pεk(rk−Φk−1rk−1) =N(rk |Φk−1rk−1,Σk), (31) when the state transition models are defined with additive white Gaussian process noises. The subscripted expressionpx(·)is occa-sionally used to emphasize that we are considering the probability density of the random variablex, but we omit the subscript whenever it is clear from the context what random variable we are considering.
Also, the observations are conditionally independent p(y1:k |x0:k,r0:k,Λ0:k) =
�k i=1
p(yi |xi,ri,Λi). (32) Thelikelihood function(32) can be expressed as
�k i=1
pvi(yi −hi(xi)−Λiri) =
�k i=1
N�yi |hi(xi) + Λiri,Ri�
, (33)
in the case of additive white Gaussian observation noise (2). If the noise is non-additive, then the form of the likelihood function can be significantly more complex.