Bayesian approach to linear state space model with unknown measurement matrix

Abstract

A linear state space model is of two equations, one of state space and the other of measurement. Estimation methods for the parameters of the model have been developed from the historic Kalman filter method. The Bayes estimation has also been used under a variety of conditions on the parameter space. We explored the availability of the Bayes method for parameter estimation with no constraints on the parameter space and found that the estimation for the state space is acceptable as long as the priors are not vague on both the state and the parameter space. We also investigated the model where the measurement matrix is contaminated with noise and found that the estimates for the state space were more accurate than those by the methods in literature. We made remarks on extended applications of the Bayes method for the linear state space model where a variety of constraints are imposed on the parameter space.

Keywords

Bayes estimation gibbs sampling kalman filter mean distance error measurement error

1. Introduction

Consider a linear state space model consisting of two equations as

$\displaystyle s_{t}=F_{t-1}s_{t-1}+w_{t}$ (1) $\displaystyle y_{t}=\Phi_{t}s_{t}+e_{t},$

where $s_{t}$ is a $d\times 1$ state vector, $y_{t}$ is a $p\times 1$ measurement vector, $F_{t}$ and $\Phi_{t}$ are $d\times d$ and $p\times d$ coefficient matrices, respectively, and $w_{t}$ and $e_{t}$ are noises with $w_{t}\sim N_{d}(0,Q_{t})$ and $e_{t}\sim N_{p}(0,R_{t})$ where $N_{d}$ stands for ‘ $d$ -variate Normal distribution’. The state changes with time and the measurement is made on the state at every time point. The first equation is called a state equation and the second a measurement equation. $F_{t}$ and $\Phi_{t}$ are called system and measurement matrices, respectively.

This model is also called a dynamic linear model which was introduced by Kalman [12] and Kalman and Bucy [13]. The model was introduced as a model for aerospace-related research and its application has been extended to a wide range of research fields for analyzing data from economics, medicine, earth sciences, and signal processing among others [4, 5, 6, 8, 11].

The Kalman filtering method is one of the most prevalent methods for the state space model and related research is abundant in literature. In essence, the model is to be analyzed in a two-step recursive process. We first make a one-time ahead prediction of the state and then the state prediction is updated based on the current data. The process we mentioned can be regarded as the Bayesian method and this viewpoint has been reflected in literature by many researchers. [16, 19] showed that the Kalman filter is derived within the Bayesian framework. Chen [20] derived the Kalman filter within a maximum a posterior (MAP) method.

One of the key assumptions in Kalman filtering is that all parameters are known except the state. However, in many real applications, the measurement matrix is not deterministic for the linear state space model [2, 7, 15, 17, 18]. Some stochastic assumptions were made on the measurement matrix in this case as described in the next section.

A main concern in the estimation of the measurement matrix and the state is that the two unknowns form one term in the measurement equation and that this may make the state $s_{t}$ unidentifiable. We found in this work that it is estimated by a Bayes method with good quality as long as at least one of the priors on the state and the measurement matrix is informative. For the case that the measurement matrix is contaminated with measurement noise, the Bayes method is shown to produce at least as good results as the existing methods in the context of estimation accuracy.

This paper is of 6 sections. Some previous works that are closely related to our work are briefly described in Section 2. In Section 3, we derived implicitly the formula of the posterior distribution of the state and the measurement matrix with Gaussian priors assigned on them individually. Numerical experiments were carried out in Section 4 by applying Gibbs sampling for parameter estimation and we checked on the sensitivity of the estimates on the priors. We also discussed the non-informative priors for the state. In Section 5, we considered the case where the measurement matrix is observable with error and compared our results with those in [18]. The experiments strongly suggest that the Bayesian approach is also available for the case where the measurement matrix is observable. The paper concludes in Section 6 with some discussions and remarks.

2. Related works

Yue and Shi [7] assumed that the measurement matrix is unknown but is limited to a special form. For instance, only the last row of the matrix is unknown and the entries for the remaining part are fixed as either 0 or 1 in accordance with a specified model.

Bayesian approaches have been applied for estimating the measurement matrix in literature. For instance, Wills et al. [1] developed a fully parametrized linear state space model by the Gibbs sampling algorithm. However, this was under the limitation that the parameters are time-invariant. We extended their work to the model where the measurement matrix is variant over time.

We also considered a linear state space model where the measurement matrix is contaminated with measurement noise as considered in Ra et al. [18]. This contamination is possible under a certain situation of data collection, which is not rare in practice and will be described in detail in Section 5.

3. Bayesian approach with measurement matrix unknown

Under the assumption that all the parameters are known, the Kalman filter method is used to estimate the state. However, if some parameters are unknown, the Kalman filter method per se is not available. We consider the case where the measurement matrix is unknown. For this model, we apply a Bayesian method to estimate both $s_{t}$ and $\Phi_{t}$ . Model (1) can be expressed as

$\displaystyle s_{t}\sim N_{d}(F_{t-1}s_{t-1},Q_{t})$ $\displaystyle y_{t}\sim N_{p}(\Phi_{t}s_{t},R_{t}).$

For the initial state $s_{0}$ , we assign a multivariate normal prior,

$s_{0}\sim N_{d}(\delta_{0},\Lambda_{0}).$

As for the states $s_{t},∼{}t=1,2\ldots$ , the posterior prediction of the mean and covariance, $\hat{s}_{t|t-1}^{B}$ and $\widehat{\Sigma}_{t|t-1}^{B}$ , of $s_{t}$ given the data $y_{1},\ldots,y_{t-1}$ are used as the hyper-parameters of the normal prior of $s_{t}$ as

$s_{t}\sim N_{d}\left(\hat{s}_{t|t-1}^{B},\widehat{\Sigma}_{t|t-1}^{B}\right)$

Note that, the posterior prediction of the mean and covariance of $s_{t}$ are computed as

$\displaystyle\hat{s}_{t|t-1}^{B}=F_{t-1}\hat{s}_{t-1|t-1}^{B}.$ $\displaystyle\widehat{\Sigma}_{t|t-1}^{B}=F_{t-1}\widehat{\Sigma}_{t-1|t-1}^{B% }F_{t-1}^{T}+Q_{t}.$

Here, $\hat{s}_{t_{1}|t_{2}}^{B}$ and $\widehat{\Sigma}_{t_{1}|t_{2}}^{B}$ are respectively the posterior mean and covariance of the state at time $t_{1}$ based on data $y_{1},\ldots,y_{t_{2}}$ . This Bayesian procedure is the same in spirit as for the ordinary Kalman filter [20].

Before assigning the prior distributions on $\Phi_{t}$ , we transform $\Phi_{t}$ into a vectorized form $\Phi_{t}^{*}$ by concatenating the columns of $\Phi_{t}$ . That is,

$\displaystyle\Phi_{t}^{*}=\left[\begin{tabular}[]{c}$\Phi_{1t}$\\ $\Phi_{2t}$\\ $\vdots$\\ $\Phi_{dt}$\end{tabular}\right]$

where $\Phi_{it}$ is the $i$ -th column of $\Phi_{t}$ . Then the measurement Eq. (1) can be expressed as

$\displaystyle y_{t}\sim N_{p}\left(\left(s_{t}^{T}\otimes I_{p}\right)\Phi_{t}% ^{*},R_{t}\right)$

where $I_{p}$ is an identity matrix of rank $p$ and $\otimes$ denotes the Kronecker product.

We assign the multivariate normal prior on $\Phi_{t}^{*}$ as

$\displaystyle\Phi_{t}^{*}\sim N_{pd}(m_{t},P_{t})\text{ for }t=0,1,\ldots,$

where $m_{t}$ and $P_{t}$ are the prior mean and covariance matrix at time $t$ , respectively.

As for the posterior distribution for $s_{t}$ and $\Phi_{t}^{*}$ at time $t$ ,

$\displaystyle p\left(s_{t},\Phi_{t}^{*}|y_{t}\right)\propto p\left(y_{t}|s_{t}% ,\Phi_{t}^{*}\right)p\left(s_{t}|\hat{s}_{t|t-1}^{B},\widehat{\Sigma}_{t|t-1}^% {B}\right)p\left(\Phi_{t}^{*}|m_{t},P_{t}\right).$ (2)

Since the posterior is not obtained in an explicit form, we used the Markov chain Monte Carlo (MCMC) method, especially, the Gibbs sampling algorithm [3] to compute the posterior distribution.

To obtain samples by the Gibbs sampling algorithm, we need to derive the full conditional distribution for $s_{t}$ and $\Phi_{t}^{*}$ from the posterior distribution. The full conditional distribution of $s_{t}$ is derived from $p(y_{t}|s_{t},\Phi_{t})p(s_{t}|\hat{s}_{t|t-1}^{B},\widehat{\Sigma}_{t|t-1}^{B})$ , which is the same as the posterior obtained by the maximum a posterior(MAP) method used for the ordinary Kalman filter [20]. We derive the full conditional distribution for $s_{t}$ from Eq. (2) as

$\displaystyle s_{t}|y_{t},\Phi_{t}\sim N_{d}\left(\hat{s}_{t|t}^{B},\widehat{% \Sigma}_{t|t}^{B}\right)$ (3)

where

$\displaystyle\hat{s}_{t|t}^{B}=\hat{s}_{t|t-1}^{B}+\widehat{\Sigma}_{t|t-1}^{B% }\Phi_{t}^{T}\left(\Phi_{t}\widehat{\Sigma}_{t|t-1}^{B}\Phi_{t}^{T}+R_{t}% \right)^{-1}\left(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{B}\right)$ (4) $\displaystyle\widehat{\Sigma}_{t|t}^{B}=\widehat{\Sigma}_{t|t-1}^{B}-\widehat{% \Sigma}_{t|t-1}^{B}\Phi_{t}^{T}\left(\Phi_{t}\widehat{\Sigma}_{t|t-1}^{B}\Phi_% {t}^{T}+R_{t}\right)^{-1}\Phi_{t}\widehat{\Sigma}_{t|t-1}^{B}.$

Similarly, the full conditional distribution for $\Phi_{t}^{*}$ is derived in the form of a multivariate normal distribution as

$\displaystyle\Phi_{t}^{*}|s_{t},y_{t}\sim N_{p\times d}\left(\left(\left(s_{t}% ^{T}\otimes I_{p}\right)^{T}R_{t}^{-1}\left(s_{t}^{T}\otimes I_{p}\right)+P_{t% }^{-1}\right)^{-1}\left(\left(s_{t}^{T}\otimes I_{p}\right)^{T}R_{t}^{-1}y_{t}% +P_{t}^{-1}m_{t}\right),\right.\left.∼{}∼{}∼{}\left(\left(s_{t}^{T}\otimes I_{% p}\right)^{T}R_{t}^{-1}\left(s_{t}^{T}\otimes I_{p}\right)+P_{t}^{-1}\right)^{% -1}\right).$ (5)

In the Gibbs sampling, we sample $s_{t}$ as $s_{t}^{(k)},∼{}k=1,\ldots,K$ where $s_{t}^{(k)}$ is from the full conditional distribution of $s_{t}$ , ${\cal L}(s_{t}|{\Phi_{t}^{*}}^{(k-1)},y_{t})$ . The sample of $\Phi_{t}^{*}$ is obtained as ${\Phi_{t}^{*}}^{(k)},∼{}k=1,\ldots,K$ where ${\Phi_{t}^{*}}^{(k)}$ is from the full conditional distribution of $\Phi_{t}^{*}$ , ${\cal L}(\Phi_{t}^{*}|s_{t}^{(k)},y_{t})$ . As for the initial values, $s_{t}^{(0)}$ and ${\Phi_{t}^{*}}^{(0)}$ , they are obtained as

$\displaystyle s_{t}^{(0)}\sim N_{d}\left(\hat{s}_{t|t-1}^{B},∼{}\hat{\Sigma}_{% t|t-1}^{B}\right).$ $\displaystyle{\Phi_{t}^{*}}^{(0)}\sim N_{pd}\left(\widehat{\Phi}^{*}_{t-1|t-1}% ,∼{}\widehat{P}_{t-1|t-1}\right).$

[t!] Gibbs sampling algorithm for Bayes estimation of

s_{t}

and

\Phi_{t}

Known parameters:

F, Q, R

for

t=1

to Tmax,

Given

y_{t}

\hat{s}_{t-1|t-1},∼{}\hat{\Sigma}_{t-1|t-1},∼{}(\hat{\Phi}^{*}_{t-1|t-1},∼{}% \hat{P}_{t-1|t-1})

or prior information for

\Phi_{t}^{*}

Predict

\hat{s}_{t|t-1}^{B},\hat{\Sigma}_{t|t-1}^{B},m_{t},P_{t}

Set starting point

s_{t}^{(0)}\sim N_{d}(\hat{s}_{t|t-1}^{B},∼{}\hat{\Sigma}_{t|t-1}^{B})

{\Phi_{t}^{*}}^{(0)}\sim N_{pd}(m_{t},P_{t})

for

k=1

K

Sample

s_{t}^{(k)}

from

{\cal L}(s_{t}|{\Phi_{t}^{*}}^{(k-1)},y_{t})

in Eq. (3)

Sample

{\Phi_{t}^{*}}^{(k)}

from

{\cal L}(\Phi_{t}^{*}|s_{t}^{(k)},y_{t})

in Eq. (5)

end

Burn in the first

K_{0}

samples and calculate posterior estimates

\hat{s}_{t|t}=\frac{1}{K-K_{0}}\sum_{n=K_{0}+1}^{K}s_{t}^{(n)},∼{}\hat{\Sigma}% _{t|t}=\text{cov}(s_{t}^{(n)}),∼{}n\ \text{is from }K_{0}+1\text{ to }K

\hat{\Phi}^{*}_{t|t}=\frac{1}{K-K_{0}}\sum_{n=K_{0}+1}^{K}{\Phi_{t}^{*}}^{(n)}% ,∼{}\hat{P}_{t|t}=\text{cov}({\Phi_{t}^{*}}^{(n)}),∼{}n\text{ is from }K_{0}+1% \text{ to }K

end

Before we estimate $s_{t}$ and $\Phi_{t}^{*}$ , we need to check the convergence of the samples. An initial part of the sample may not be appropriate for the desired posterior distribution. Therefore, we remove an initial part of the sample which is referred to as the ‘burn in’ part of the sample. A simple way to determine the burn-in part is to draw a trace plot of the sample and burn in an initial part of the sample before the stationarity of the sample is observed. Several burn-in diagnostic methods were developed including that by Geweke, Gelman-Rubin and Heidelberger-Welch [9]. After burning in the first $K_{0}$ samples, the posterior mean and covariances of $s_{t}$ and $\Phi_{t}^{*}$ are computed. This procedure is summarized in Algorithm 3.

4. Numerical experiments: Effect of priors

In this section, we will investigate the effect of the priors upon the estimates of the state $s_{t}$ and the measurement coefficient $\Phi_{t}$ by numerical experiments and suggest a guideline for priors for a safe estimation process.

We set values for the parameters of model Eq. (1) and the initial value of the state as in Table 1 to generate data for our experiment. We denote by $s_{0}$ the true initial value of the state and let $s_{t}=(x_{t},y_{t},v_{t})^{T}$ where the first two coordinates are the $x$ and $y$ coordinates of the target on a plane at time $t$ , and the target moves with the velocity $v_{t}(-0.2,-0.2)$ from the initial position $(40,40)$ until $t=M$ .

Table 1
The parameter values set for data generation in model Eq. (1)

Parameter values	$F_{t}=\left[\begin{array}[]{ccc}1&0&-0.2\\ 0&1&-0.2\\ 0&0&1\end{array}\right]$ , $Q_{t}=\left[\begin{array}[]{ccc}0.1^{2}&0&0\\ 0&0.1^{2}&0\\ 0&0&0.1^{8}\end{array}\right]$
	$\Phi_{t}=\left[\begin{array}[]{ccc}1&0.5&0\\ 0.5&0.7&0\end{array}\right]$ , $R_{t}=\left[\begin{array}[]{ccc}1&0\\ 0&1\end{array}\right]$
Initial state value	$s_{0}=(x_{0},y_{0},v_{0})^{T}=(40,40,1)^{T}$
Data size	$M=150$
Number of data set	$n=30$

Table 2

Prior distributions on the state and the measurement matrix. In the table, $t=1,2,\ldots$ , ‘diag $(b_{n\times 1}$ )’ denotes the size $n$ diagonal matrix with components $b$ and ‘diag ${}_{n}(a$ )’ denotes the size $n$ diagonal matrix with components $a$

Case		1	2	3
$s$		Informative	Non-informative	Informative
$\Phi$		Informative	Informative	Non-informative
$s_{0}$	$\delta_{0}$	$s_{0}$	$\left(0,0,0\right)^{T}$	$s_{0}$
	$\Lambda_{0}$	diag $((1,1,0.1^{2})^{T})$	diag $((20^{2},20^{2},1)^{T})$	diag $((1,1,0.1^{2})^{T})$
$\Phi_{0}^{*}$	$m_{0}$	$\Phi_{\textit{0,true}}^{*}$	$\Phi_{\textit{0,true}}^{*}$	$\left(0,0,0,0,0,0\right)^{T}$
	$P_{0}$	$\text{diag}_{6}(0.1^{2})$	$\text{diag}_{6}(0.1^{2})$	$\text{diag}_{6}(1)$
$s_{t}$	$\hat{s}_{t\|t-1}^{B}$	$F_{t-1}\hat{s}_{t-1\|t-1}^{B}$
	$\hat{\Sigma}_{t\|t-1}$	$F_{t-1}\hat{\Sigma}_{t-1\|t-1}^{B}F_{t-1}^{T}+Q_{t}$
$\Phi_{t}^{*}$	$m_{t}$	$\Phi_{\textit{0,true}}^{*}$	$\Phi_{\textit{0,true}}^{*}$	$\hat{\Phi}_{t-1}^{*}$
	$P_{t}$	$\text{diag}_{6}(0.1^{2})$	$\text{diag}_{6}(0.1^{2})$	$\hat{P}_{t-1}$

In addition to the specification on the parameters and the state, we impose priors on the state and the measurement matrix with different levels of uncertainty as in Table 2. We consider three cases of the uncertainty for the state $s_{t}$ and $\Phi_{t}$ , both informative and only one part informative. The priors are a reflection of the domain knowledge, represented by more informative distributions as we get more domain knowledge.

We will look into the case that both are non-informative separately. In the experiment, we take $\Phi_{t}$ as constant over time. In estimating $s_{t}$ and $\Phi_{t}$ , we used an iterative process applying an MCMC method with a burn-in of 1000 iterations.

For performance evaluation, we used the mean distance error (MDE) defined for $s_{t}$ as

$\displaystyle\textit{MDE}(\hat{s}_{t})=\frac{1}{n}\sum_{i=1}^{n}\sqrt{(\hat{s}% _{1,t}^{i}-s_{1,t})^{2}+(\hat{s}_{2,t}^{i}-s_{2,t})^{2}}$

where $\hat{s}_{j,t}^{i}$ is the $j$ th component of $\hat{s}_{t}$ based on the $i$ th data set. This MDE is the mean of the estimation error of $s_{t}$ . As for $\Phi_{t}$ , its estimation is evaluated by the MDE given by

$\displaystyle\textit{MDE}(\widehat{\Phi}_{t})=\frac{1}{n}\sum_{i=1}^{n}\sqrt{% \left(\sum_{a=1}^{p}\sum_{b=1}^{d}(\widehat{\Phi}_{a,b,t}^{i}-\Phi_{a,b,t})^{2% }\right)}$

where $\hat{\Phi}_{a,b,t}^{i}$ is the $(a,b)$ entry of the matrix $\widehat{\Phi}_{t}$ based on the $i$ th data set. Note that $\textit{MDE}(\widehat{\Phi}_{t})$ is the average of the Frobenius norms of $\widehat{\Phi}_{t}-\Phi_{t}$ .

In the evaluation, the result by the ordinary Kalman filter is included as a reference since $\Phi_{t}$ is assumed known in the Kalman filter method. The MDE’s are compared among the three cases of uncertainty levels of the priors, and the comparison is summarized in Fig. 1 and 2 regarding $s_{t}$ and $\Phi_{t}$ respectively.

Figure 1.

MDE( $\hat{s}_{t}$ ) for 3 cases of priors. The three cases of prior assignment on $s_{t}$ and $\Phi_{t}$ are: informative priors for both, non-informative and informative priors respectively, and informative and non-informative priors respectively. The MDE by the ordinary Kalman filter is given for reference.

In Fig. 1, it is strongly indicated that a high uncertainty (case 2 in Table 2) about the initial state $s_{0}$ is reflected in unstable estimates of $s_{t}$ at an early stage of the estimation process which ends up with an estimate as good as those in the other cases of uncertainty. This is as expected with non-informative priors. The estimation error decreases as data accumulates enough so that the initial uncertainty due to a vague prior may be dominated by the information in data. The instability, if any, at an early stage takes place in various forms in size and time depending on the information amounts in the prior and the data. The shape of instability in the figure is typical when the prior is of no or little information. It may also have some fluctuations at the beginning with the size decaying down at a rate as data accumulates.

Figure 2.

MDE( $\widehat{\Phi}_{t}$ ) for 3 cases of priors. See the caption of Fig. 1 for the prior assignment.

Figure 3.

Mean distance errors of $s_{t}$ and $\Phi_{t}$ with the non-informative priors assigned on $s_{t}$ and $\Phi_{t}$ respectively in Table 2. The upper graph is of MDE( $\hat{s}_{t}$ ) and MDE( $\hat{\Phi}_{t}$ ) in the other graph.

It is noteworthy in Fig. 2 that MDE( $\Phi_{t}$ ) is large and unstable when the prior of $\Phi_{0}$ is non-informative (case 3 in Table 2) while the estimate of $s_{t}$ is stable and as good as those under the other cases of the priors. This phenomenon is possible since many different values of $\Phi_{t}$ are possible for given values of $y_{t}$ and $s_{t}$ . One concern though is that an accurate estimate of $s_{t}$ with a poor estimate of $\Phi_{t}$ may not be welcome in the context of statistical modeling.

When both of $s_{t}$ and $\Phi_{t}$ are with non-informative priors, we could anticipate a bizarre estimation result as shown in Fig. 3. The scales are quite different between the graph in Fig. 1 and the upper one in Fig. 3. It is about 60 for the latter graph while it is about 10 for the former graph. If we ignore the case (case 2 in Table 2) of the non-informative prior on $s_{t}$ , the scale difference is much larger. We also see an analogous phenomenon with respect to MDE( $\Phi_{t}$ ).

We further investigated the case where the priors on $s_{t}$ are non-informative with a wide range of uncertainty levels on $s_{0}$ as listed in Table 3. The prior on $\Phi_{t}$ is assigned as informative. The estimation results are displayed in Fig. 4. The graph in panel (b) is based on one set of data rather than the average of the results based on 30 data sets in order to show a detailed description of the estimation process. The results for the other data sets were more or less the same as that in the figure.

Table 3

Prior distributions on $s_{0}$ . In the table, ‘diag $(b_{n\times 1}$ )’ denotes the size $n$ diagonal matrix with components $b$

$\delta_{0}$	$(100,100,5)^{T}$
$\Lambda_{0}$	diag $((10^{2},10^{2},0.5^{2})^{T})$
	diag $((50^{2},50^{2},2.5^{2})^{T})$
	diag $((100^{2},100^{2},5^{2})^{T})$
	diag $((200^{2},200^{2},10^{2})^{T})$

Figure 4.

The results of simulation under non-informative priors on $s$ in Table 3. (a) The MDE of $s$ under each priors. (b) The estimates of $s$ for a data set. The solid line in (b) indicates the true state of $s$ and the circled points indicate the estimates at the beginning. The priors on $\Phi$ are assigned informative as in Table 2.

We could see in this figure that the estimation was bizarre when $\Lambda_{0}$ is given as the first one (denote it by $\Lambda_{0,1}$ ) in Table 3. Panel (b) in the figure displays the estimates of the ( $x, y$ ) coordinate of $s_{t}$ . Note that the true initial target location and velocity $s_{0}$ is (40,40,1) while they are set as (100,100,5) in the prior. As for the prior $\Lambda_{0,1}$ , the estimate of $s_{t}$ starts at around (50,60), then moves smoothly towards the origin until the x coordinate reaches about 25, then the trace changes its direction drastically wandering for a short while (see the sharp up-and-down in panel (a)), then the trace moves slowly towards the origin along the actual trace line of the target. This bizarre movement of the estimate was not shown for the other values of $\Lambda_{0}$ in Table 3 that are of a much higher uncertainty level than $\Lambda_{0,1}$ . This is an apparent indication that vagueness in prior helps for a good level of estimation accuracy when imposing prior distributions with little or no prior information. Although the initial target location was set far away from the true location in the prior, the vagueness, well incorporated with data, took us near the target location a few iterations after the beginning of the estimation process.

In target localization problems, we may impose priors on $\Phi_{0}$ subject to the signal reception system and its function reliability. In this context, the priors on $\Phi_{0}$ could be informative in real world. Furthermore, $\Phi_{t}$ may be given as a function of data and parameters as considered below. In this case, the prior would be relatively easy to determine. We will address this issue in comparison with a traditional estimation method as proposed in [18].

5. When

\Phi

is contaminated with measurement noise

In target localization problems, we use signal data that are received from the target and there are several signal reception methods. The data may be collected as angles of signal arrival to sensors from the target, as time difference of signal arrival, or as frequency difference of signal arrival which utilizes the frequency Doppler shift phenomenon, among others[14]. In case we use the method of time difference of signal arrival for data collection, we can form the measurement Eq. (1) as consisting of the target location as an unknown parameter and the observed values of the distances between sensors and the time differences of signal arrival. We can easily imagine this phenomenon if we consider a triangle of the three points of two sensors and one target. An explicit expression of the equation is given in Eq. (11).

Ra et al. [18] applied an optimization approach by using a cost function for target localization when $\Phi_{t}$ is contaminated with measurement error. This method is called the non-conservative robust Kalman filter (NCRKF). They considered the model in Eq. (1) with $\Phi_{t}$ replaced with $\widetilde{\Phi}_{t}-\Delta\Phi_{t}$ as

$\displaystyle s_{t}=F_{t-1}s_{t-1}+w_{t}$ $\displaystyle y_{t}=[\tilde{\Phi}_{t}-\Delta\Phi_{t}]s_{t}+e_{t}$

with additional conditions that

$\displaystyle E[\Delta\Phi_{t}]=0.E[\Delta\Phi_{t}w_{t}]=0.$

(6) $\displaystyle E[\Delta\Phi_{t}^{T}R_{t}^{-1}\Delta\Phi_{t}]=W_{t}.E[\Delta\Phi% _{t}^{T}R_{t}^{-1}e_{t}]=V_{t}.$

They derived the formulae for update and prediction of $s_{t}$ as

( $i$ ) Update:

$\displaystyle\hat{s}_{t|t}^{N}=\left(I+\Sigma_{t|t}^{N}W_{t}\right)\hat{s}_{t|% t-1}^{N}+\Sigma_{t|t}^{N}\tilde{\Phi}_{t}^{T}R_{t}^{-1}\left(y_{t}-\tilde{\Phi% }\hat{s}_{t|t-1}^{N}\right)-\Sigma_{t|t}^{N}V_{t}.$ (7) $\displaystyle\left(\Sigma_{t|t}^{N}\right)^{-1}=\left(\Sigma_{t|t-1}^{N}\right% )^{-1}+\tilde{\Phi}_{t}^{T}R_{t}^{-1}\tilde{\Phi}_{t}-W_{t}.$

( $i i$ ) Prediction:

$\displaystyle\hat{s}_{t+1|t}^{N}=F_{t}\hat{s}_{t|t}^{N}.$ $\displaystyle\Sigma_{t+1|t}^{N}=F_{t}\Sigma_{t|t}^{N}F_{t}^{T}+Q_{t+1}.$

Provided the measurement matrix in model Eq. (1) is observed with error, we may take the observed $\Phi_{t}$ , $\widetilde{\Phi}_{t}$ , as the hyper-parameter as

$\displaystyle\Phi_{t}^{}\sim N_{p\times d}(\tilde{\Phi}_{t}^{},P_{t})\text{ % for }t=0,1,\ldots,.$

We will call this prior an empirical prior.

Note that $\hat{s}_{t|t}^{B}$ in Eq. (4) can be expressed as a function of $\hat{\Sigma}_{t|t}^{B}$ as in Eq. (9) below:

$\displaystyle\hat{s}_{t|t}^{B}=\hat{s}_{t|t-1}^{B}+\hat{\Sigma}_{t|t-1}^{B}% \Phi_{t}^{T}(\Phi_{t}\hat{\Sigma}_{t|t-1}^{B}\Phi_{t}^{T}+R_{t})^{-1}\left(y_{% t}-\Phi_{t}\hat{s}_{t|t-1}^{B}\right)=\hat{s}_{t|t-1}^{B}-\hat{\Sigma}_{t|t-1}% ^{B}\Phi_{t}^{T}\left[(\Phi_{t}\hat{\Sigma}_{t|t-1}^{B}\Phi_{t}^{T}+R_{t})^{-1% }\left(\Phi_{t}\hat{\Sigma}_{t|t-1}^{B}\Phi_{t}^{T}R_{t}^{-1}-\left(\Phi_{t}% \hat{\Sigma}_{t|t-1}^{B}\Phi_{t}^{T}R_{t}^{-1}+I\right)\right)\right]\cdot% \left(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{B}\right)=\hat{s}_{t|t-1}^{B}-\hat{\Sigma% }_{t|t-1}^{B}\Phi_{t}^{T}\left[\left(\Phi_{t}\hat{\Sigma}_{t|t-1}^{B}\Phi_{t}^% {T}+R_{t}\right)^{-1}\Phi_{t}\hat{\Sigma}_{t|t-1}^{B}\Phi_{t}^{T}R_{t}^{-1}-R_% {t}^{-1}\right]\left(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{B}\right)=\hat{s}_{t|t-1}^% {B}+\hat{\Sigma}_{t|t-1}^{B}\left(I-\Phi_{t}^{T}\left(\Phi_{t}\hat{\Sigma}_{t|% t-1}^{B}\Phi_{t}^{T}+R_{t}\right)^{-1}\Phi_{t}\hat{\Sigma}_{t|t-1}^{B}\right)% \Phi_{t}^{T}R_{t}^{-1}\left(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{B}\right)$ (8) $\displaystyle=\hat{s}_{t|t-1}^{B}+\hat{\Sigma}_{t|t}^{B}\Phi_{t}^{T}R_{t}^{-1}% \left(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{B}\right).$ (9)

Equation (8) is obtained by the matrix inversion lemma[10]. It is noteworthy that Eq. (9) is the same as the state estimate at time $t$ in the ordinary Kalman filter(OKF). However, $\Phi_{t}$ is unknown unlike for the OKF. Therefore, the estimation error by the OKF method is generally smaller than that by the Bayes method with the posterior distribution in Eq. (2) and we have checked this through numerical experiments in Section 4.

The $\hat{s}_{t|t}^{N}$ in Eq. (5) can be expressed, under the conditions in Eq. (5), as

$\displaystyle\hat{s}_{t|t}^{N}=\hat{s}_{t|t-1}^{N}+\Sigma_{t|t}^{N}\left(W_{t}% \hat{s}_{t|t-1}^{N}+\tilde{\Phi}_{t}^{T}R_{t}^{-1}\left(y_{t}-\tilde{\Phi}\hat% {s}_{t|t-1}^{N}\right)-V_{t}\right)=\hat{s}_{t|t-1}^{N}+\Sigma_{t|t}^{N}\Phi_{% t}^{T}R_{t}^{-1}\left(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{N}\right)+\Sigma_{t|t}^{N% }\left(W_{t}\hat{s}_{t|t-1}^{N}+(\Phi_{t}^{T}+\Delta\Phi_{t}^{T})R_{t}^{-1}% \left(y_{t}-(\Phi_{t}+\Delta\Phi_{t})\hat{s}_{t|t-1}^{N}\right)-V_{t}\right)=% \hat{s}_{t|t-1}^{N}+\Sigma_{t|t}^{N}\Phi_{t}^{T}R_{t}^{-1}\left(y_{t}-\Phi_{t}% \hat{s}_{t|t-1}^{N}\right)+\Sigma_{t|t}^{N}\left(\left(W_{t}-\Delta\Phi_{t}^{T% }R_{t}^{-1}\Delta\Phi_{t}\right)\hat{s}_{t|t-1}^{N}+\left(\Delta\Phi_{t}R_{t}^% {-1}\left(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{N}\right)-V_{t}\right)\right)=\hat{s}% _{t|t-1}^{N}+\Sigma_{t|t}^{N}\Phi_{t}^{T}R_{t}^{-1}\left(y_{t}-\Phi_{t}\hat{s}% _{t|t-1}^{N}\right)+\Sigma_{t|t}^{N}\left(\left(W_{t}-\widehat{W}_{t}\right)% \hat{s}_{t|t-1}^{N}+\left(\widehat{V}_{t}-V_{t}\right)\right)$ (10)

where $\widehat{W}_{t}=\Delta\Phi_{t}^{T}R_{t}^{-1}\Delta\Phi_{t}$ and $\widehat{V}_{t}=\Delta\Phi_{t}R_{t}^{-1}(y_{t}-\Phi_{t}\hat{s}_{t|t-1}^{N})$ . Note that $E[\hat{W}_{t}]=W_{t}$ and $E[\hat{V}_{t}]=V_{t}$ .

The update equations for $\hat{s}_{t|t}$ are different between the NCRKF and the Bayesian methods as shown in Eqs (9) and (10), respectively. The first two terms on the right hand side of Eq. (10) are analogous to those for $\hat{s}_{t|t}^{B}$ in Eq. (9), but $\hat{s}_{t|t}^{N}$ has another term,

$\displaystyle\Sigma_{t|t}^{N}\left((W_{t}-\widehat{W}_{t})\hat{s}_{t|t-1}^{N}+% (\widehat{V}_{t}-V_{t})\right),$

which may be regarded as an error term due to the measurement error and its related function $\Delta\Phi_{t}$ . What’s worse is that this error term is influenced by $\hat{s}_{t|t-1}^{N}$ . In case of target localization, when the target is far off from the sensors, $\hat{s}_{t|t-1}^{N}$ would be a main source of instability of estimation at the initial stage of the estimation process. This will be demonstrated in numerical experiments below.
5.1 Numerical experiment for performance comparison

In this subsection, we compare the performance between our Bayesian approach and the estimation approach as proposed in [18] in target localization. Han et al. [15] considered a target localization problem where the signal is received from the target by observing the time difference of signal arrivals on different sensors. Suppose we have 5 sensors with their positions denoted by $(\eta_{x_{j}},\eta_{y_{j}})$ for $j=$ 0, 1, 2, 3, 4 where the main sensor is labeled 0. Denote by $\tau_{j}$ the distance between the target and sensor $j$ and by $d_{j0}$ the distance between sensors 0 and $j$ . Let

$r_{j}=\tau_{j}-\tau_{0}$

and denote the target position by $(\theta_{x},\theta_{y})$ on a plane. For notational simplicity, we ignore the time index in these notations.

Then we get the following equation by applying a cosine law to the triangle formed by the positions of the target and the sensors $0$ and $j$ on a plane:

$\displaystyle r_{j}^{2}-d_{j0}^{2}=-2[\eta_{x_{j}}-\eta_{x_{0}},∼{}\eta_{y_{j}% }-\eta_{y_{0}},∼{}r_{j}]\left[\begin{array}[]{ccc}\theta_{x}-\eta_{x_{0}}\\ \theta_{y}-\eta_{y_{0}}\\ \tau_{0}\end{array}\right].$

If we denote the left-hand side of this equation plus a measurement error by $y_{t}$ and let $\phi_{j}=-2[\eta_{x_{j}}-\eta_{x_{0}},∼{}\eta_{y_{j}}-\eta_{y_{0}},∼{}r_{j}]$ , then we have a measurement equation as

$\displaystyle y_{t}=\Phi_{t}s_{t}+e_{t}=[\tilde{\Phi}-\Delta\Phi_{t}]s_{t}+e_{% t},$ (11)

where

$\displaystyle s_{t}=\left[\begin{array}[]{ccccccccc}\theta_{x}-\eta_{x_{0}}\\ \theta_{y}-\eta_{y_{0}}\\ d_{{t},0}\\ \dot{\theta}_{x}-\dot{\eta}_{x_{0}}\\ \dot{\theta}_{y}-\dot{\eta}_{y_{0}}\\ \dot{d}_{{t},0}\\ \ddot{\theta}_{x}-\ddot{\eta}_{x_{0}}\\ \ddot{\theta}_{x}-\ddot{\eta}_{x_{0}}\\ \ddot{d}_{{t},0}\end{array}\right],\Phi_{t}=\left[\begin{array}[]{ccc}\phi_{1}% &\\ \vdots&0^{4\times 6}\\ \phi_{4}&\end{array}\right],\tilde{\Phi}_{t}=\left[\begin{array}[]{ccc}\tilde{% \phi}_{1}&\\ \vdots&0^{4\times 6}\\ \tilde{\phi}_{4}&\end{array}\right],\Delta\Phi_{t}=\left[\begin{array}[]{ccc}% \Delta\phi_{1}&\\ \vdots&0^{4\times 6}\\ \Delta\phi_{4}&\end{array}\right].$

Here, $\dot{\theta}$ and $\dot{\eta}$ represent the velocities of the target and sensors, and $\ddot{\theta}$ and $\ddot{\eta}$ their accelerations, respectively. $\tilde{\phi}_{j}$ is the observation of $\phi_{j}$ and $\Delta\phi_{j}$ is the zero mean error which satisfies Eq. (5).

With the time interval $T$ , the state equation is given by

$\displaystyle s_{t}=F_{t-1}s_{t-1}+w_{t}$ (12)

with

$F_{t-1}=\left[\begin{array}[]{ccccc}I^{3\times 3}&T\cdot I^{3\times 3}&\frac{1% }{2}T^{2}\cdot I^{3\times 3}\\ 0^{3\times 3}&I^{3\times 3}&T\cdot I^{3\times 3}\\ 0^{3\times 3}&0^{3\times 3}&I^{3\times 3}\end{array}\right]$

and

$\text{Cov}(w_{t})=Q_{t}=\left[\begin{array}[]{ccccc}\frac{1}{2}T^{2}\cdot I^{3% \times 3}\\ T\cdot I^{3\times 3}\\ I^{3\times 3}\end{array}\right]\left[\begin{array}[]{ccccc}0.1^{2}&0&0\\ 0&0.1^{2}&0\\ 0&0&2\cdot 0.1^{2}\end{array}\right]\left[\begin{array}[]{ccccc}\frac{1}{2}T^{% 2}\cdot I^{3\times 3}\\ T\cdot I^{3\times 3}\\ I^{3\times 3}\end{array}\right]^{T}.$

A detailed description of the parameters are explained in [15].

For numerical experiments, the parameter values for the state and measurement equations were specified as in Table 4 in [15]. We generated data under the same specification for our experiment for performance comparison between the NCRKF and Bayes methods.

Table 4

Specification for data generation from the model consisting of Eqs (11) and (12)

Target	Velocity $=$ 20(m/s)
	Direction $=$ $135^{\circ}$
	Initial position $=$ (5000,0)(m)
Sensors	Initial velocity $=$ 270 (m/s)
	Direction $=$ $30^{\circ}$
	Initial position of leader sensor $=$ (0,0) (m)
	$(\eta_{x_{i}},\eta_{y_{i}})=(\eta_{x_{0}},\eta_{y_{0}})+d_{i0}(\cos((i-1)\pi/2% +\pi/6),∼{}\sin((i-1)\pi/2+\pi/6))$
	$d_{i0}=$ 300 m for $i=$ 1, 2, 3, 4
Noise variance	Range difference $\sigma_{r}^{2}=0.5^{2}m^{2}$
	Sensor position $\sigma_{x}^{2}=0.5^{2}m^{2}$ , $\sigma_{y}^{2}=0.5^{2}m^{2}$
Data	Time interval $(T)=$ 0.04
	Number of time points $(M)=$ 450
	Number of data sets $(n)=$ 30

Figure 5.

Mean absolute error of the four components of $\hat{s}_{t|t}$ , the $x$ and $y$ coordinates of the target position and the target velocity. The results by NCRKF are given in black as 95% confidence intervals and those by the Bayes method in red as 95% credible intervals.

The priors for $s_{0}$ and $\Phi_{0}$ are given as

$\displaystyle s_{0}\sim N_{9}(\mu_{0},\Sigma_{0})$

and

$\displaystyle\Phi_{0}^{*}\sim N_{36}(\widetilde{\Phi}_{t},P_{t})∼{}\text{for}t% =0,1,\ldots,M$

where

$\displaystyle\mu_{0}=[5000,0,5000,-247.97,-120.86,-247.97,0,17.98,-12.30]^{T},$ $\displaystyle\Sigma_{0}=\textit{diag}([60^{2},60^{2},2\cdot 60^{2},5^{2},5^{2}% ,2\cdot 5^{2},3^{2},3^{2},2\cdot 3^{2}],\text{ and }$ (13) $\displaystyle P_{0}=\left[\begin{array}[]{cc}0.5I_{12}&\mathbb{0}\\ \mathbb{0}&0.01^{2}I_{24}\end{array}\right].$

The values of $\mu_{0}$ , $\Sigma_{0}$ , and $P_{0}$ are taken in accordance with the specification in Table 4.

The initial value $\hat{s}_{0}$ of $s_{0}$ was generated for the NCRKF method from the distribution $N_{9}(\mu_{0},\Sigma_{0})$ . Each of the 30 data sets generated was analyzed for the estimation of $s_{t}$ and the results are summarized in Fig. 5 as 95% confidence intervals for the NCRKF and as 95% credible intervals for the Bayes method. The interval types are different due to the difference in the estimation methods.

We can see in the figure that the Bayes method performs better than the NCRKF. In particular, the estimates by the NCRKF are quite unstable at an early stage of the process while those by the Bayes method look well under control. If we look at the $x$ coordinates of $\hat{s}_{t}$ , we can see in Table 5 that the mean absolute error by the Bayes method is about 0.84 times as large as that by the NCRKF on average. We also compared the MDE values of $\widehat{\Phi}_{t}$ between the two methods and the results are in Fig. 6. The Bayes estimates of $\Phi_{t}$ are closer to the true value than the estimates by the NCRKF where the estimates of $\Phi_{t}$ are the same as $\widetilde{\Phi}_{t}$ .

Table 5

Comparison of the mean absolute error of the $x$ -coordinate of the target position with the hyper-parameter $P_{0}$ in Eq. (5.1)

	Time( $t$ )
	2	4	6	8	10	12	14	16	18
Bayes method	54.66	48.61	45.49	34.85	29.39	21.79	12.23	7.89	1.83
NCRKF	76.86	76.66	57.51	37.65	35.36	24.77	12.65	7.70	2.21
Bayes/NCRKF	0.71	0.63	0.79	0.93	0.83	0.88	0.97	1.02	0.83

Table 6

Comparison of the mean absolute error of the $x$ -coordinate of the target position for $P_{0}^{a},P_{0}^{b},P_{0}^{c}$ and $P_{0}^{d}$

$P_{0}$	Time( $t$ )
	2	4	6	8	10	12	14	16	18
$P_{0}^{a}$	69.44	72.22	62.49	35.13	36.50	21.19	11.86	7.98	2.29
$P_{0}^{b}$	59.50	51.13	55.00	31.29	31.15	23.03	11.24	7.31	1.97
$P_{0}^{c}$	54.66	48.61	45.49	34.85	29.39	21.79	12.23	7.89	1.83
$P_{0}^{d}$	51.66	79.14	84.78	120.41	88.44	48.41	18.48	22.15	5.90

Figure 6.

MDE values for $\Phi_{t}$ . MDE( $\widehat{\Phi}_{t}$ ) is in red and MDE( $\widetilde{\Phi}$ ) in black.

Figure 7.

Mean absolute error of the four components of $\hat{s}_{t|t}$ as in Fig. 5 with the hyper-parameter $P_{0}=P_{0}^{a}$ in panel (a), $P_{0}=P_{0}^{b}$ in panel (b) and $P_{0}=P_{0}^{d}$ in panel (c).

The comparison above has been made with the prior on $\Phi_{0}$ with the hyper-parameter $P_{0}$ as in Eq. (5.1). The comparison is affected by the covariance structure $P_{0}$ . We considered three other structures for $P_{0}$ as given by

$P_{0}^{a}=\left[\begin{array}[]{cc}0.01I_{12}&\mathbb{0}\\ \mathbb{0}&0.01^{2}I_{24}\end{array}\right],P_{0}^{b}=\left[\begin{array}[]{cc% }0.1I_{12}&\mathbb{0}\\ \mathbb{0}&0.01^{2}I_{24}\end{array}\right],P_{0}^{d}=\left[\begin{array}[]{cc% }2I_{12}&\mathbb{0}\\ \mathbb{0}&0.01^{2}I_{24}\end{array}\right].$

The comparison results are summarized in Fig. 7. If we denote the $P_{0}$ in Eq. (5.1) by $P_{0}^{c}$ , the performance of estimating $s_{t}$ was the best with the hyper-parameter $P_{0}^{c}$ among the four $P_{0}$ values and the worst was with $P_{0}^{d}$ . We can also see in Table 6 that performance is the best with $P_{0}=P_{0}^{c}$ and worst with $P_{0}=P_{0}^{d}$ . This result indicates that the prior on $\Phi_{t}$ with a high uncertainty may yield estimates of $s_{t}$ with a lower accuracy than the estimates by a traditional method such as the NCRKF. This result is quite intuitive in the sense that a high uncertainty on the observed $\Phi_{t}$ would make the observed value less useful. In the same context, a lower uncertainty on $\widetilde{\Phi}_{t}$ would make the estimates of $s_{t}$ more influenced by $\widetilde{\Phi}_{t}$ . This result is observed in panel (a) where the estimation accuracy is more or less the same between the NCRKF and the Bayes method. So as far as the empirical prior is concerned, it seems desirable that the uncertainty level is low but not too low.

In the Gibbs sampling used for the Bayes estimation of $s_{t}$ and $\Phi_{t}$ , we set $K=$ 1500 and $K_{0}=$ 500 at an early stage of data accumulation and then we switched to $K=$ 1100 and $K_{0}=$ 100 to reduce the computing time with as good results as those without the switching.

6. Concluding remarks

In Bayes estimation for the state and the measurement matrix, we found that the state is estimated reasonably well as long as the priors are informative for at least one of the state and the measurement matrix. In reality, it is often the case for target localization problems that we have little or no information about the initial state. In this case, it is imperative that the prior on the measurement matrix is informative. If some of the entries of the measurement matrix are known, we may apply the Bayes method for estimation with highly informative priors imposed on the corresponding entries.

There are situations where the measurement matrix is a function of observed data. In the measurement equation, $y_{t}$ is a function of $\Phi_{t}s_{t}$ . If we express $\Phi_{t}$ , under a Bayesian framework, as a sum of $\widetilde{\Phi}_{t}$ and a matrix of white Gaussian noises, then we can confine the estimate of $\Phi_{t}$ within a range of $\widetilde{\Phi}_{t}$ in probability. This confinement and $y_{t}$ values will produce the full conditional posterior distribution of $s_{t}$ from which the conditional posterior mean of $s_{t}$ is obtained. Given this mean and the data $y_{t}$ , we can find the full conditional posterior distribution of $\Phi_{t}$ from which the conditional posterior mean of $\Phi_{t}$ is obtained. This iterative process ultimately converges since this process leads to the maximum of the posterior distribution of $s_{t}$ and $\Phi_{t}$ . In this process, we begin with a randomly generated value for $s_{0}$ and this could yield subsequently a conditional posterior mean of $\Phi_{t}$ which is far off the truth. But, if we have any information from data concerning $\Phi_{t}$ as is the case considered above, the conditional posterior mean of $\Phi_{t}$ could be controlled within a range of the truth. Another point to consider concerning the noise contaminated measurement matrix is that the measurement errors of $y_{t}$ and $\Phi_{t}$ may be correlated. This correlation may well be reflected in the full conditional posterior distributions of $s_{t}$ and $\Phi_{t}$ respectively. These are merits of the Bayes approach with empirical priors. A precaution is that an empirical prior with a too low uncertainty leads to data sensitive estimates.

We could apply the Bayes method to a larger set of state space models including nonlinear models as

$\displaystyle s_{t}=f_{t}(s_{t-1})+w_{t}.$ $\displaystyle y_{t}=h_{t}(s_{t},\Phi_{t})+e_{t}.$

where $f_{t}$ and $h_{t}$ are known nonlinear functions and $\Phi_{t}$ unknown. We may consider non-normal errors for $w_{t}$ and $e_{t}$ . In this case, the full conditional posterior distributions of $s_{t}$ and $\Phi_{t}$ may be given in complicated expressions. Stochastic computing algorithms such as Metropolis Hastings algorithm would be a useful tool for handling the full conditional posterior distributions.

Footnotes

Acknowledgments

This work was supported by the NRF grant (No. 2016R1D1A1B03936155) of the Republic of Korea.

References

Wills

Schön

T.B.

Lindsten

and Ninness

, Estimation of linear systems using a Gibbs sampler, in: IFAC Symposium on System Identification, 2012.

Zheng

, A robust H2 filtering approach and its application to equalizer design for communication systems, IEEE Transactions on Signal Processing 53(8) (2005), 2735–2747.

Casella

and George

E.I.

, Explaining the Gibbs sampler, The American Statistician 46(3) (1992), 167–174.

Hamilton

J.D.

, A new approach to the economic analysis of nonstationary time series and the business cycle, Econometrica 57 (1989), 357–384.

Hamilton

J.D.

, Time series analysis, Princeton: Princeton University Press, 1994.

Durbin

and Koopman

S.J.

, Time Series Analysis by State Space Methods 2nd ed., Oxford: Oxford University Press, 2012.

Yue

and Shi

, Multi-innovation self-tuning Kalman filter with unknown parameters systems, in: International Power, Electronics and Materials Engineering Conference, 2015.

Gordon

and Smith

A.F.M.

, Monitoring and modeling biomedical time series, Journal of the American Statistical Association 85 (1990), 328–337.

Sahlin

, Estimating convergence of Markov Chain Monte Carlo simulations, Master Thesis, Stockholm University, Stockholm, 2011.

10.

Woodbury

M.A.

, Inverting modified matrices, Princeton, NJ, 1950.

11.

Shumway

R.H.

and Stoffer

D.S.

, Dynamic linear models with switching, Journal of the American Statistical Association 86 (1991), 763–769.

12.

Kalman

R.E.

, A new approach to linear filtering and prediction problems, Journal of Basic Engineering 82 (1960), 35–45.

13.

Kalman

R.E.

and Bucy

, New results in linear filtering and prediction theory, The American Society of Mechanical Engineers 83 (1961), 95–107.

14.

Kim

S.H.

Park

J.H.

Yoon

W.S.

and Whang

I.H.

, A note on sensor arrangement for long-distance target localization, Signal Processing 133 (2017), 18–31.

15.

Han

S.K.

W.S.

Whang

I.H.

and Park

J.B

, Linear recursive passive target tracking filter for cooperative sea-skimming anti-ship missiles, IET Radar, Sonar and Navigation 8(7) (2014), 805–814.

16.

Peterka

, Bayesian approach to system identification, in: Trends and Progress in System Identification, Pergamon Press, 1981, pp. 239–304.

17.

W.S.

Whang

I.H.

and Ahn

J.Y.

, Robust horizontal line-of-sight rate estimator for sea skimming anti-ship missile with two-axis gimballed seeker, IEE Proceedings, Radar Sonar Navigations 152(1) (2005), 9–15.

18.

W.S.

Whang

I.H.

and Park

J.B.

, Non-conservative robust Kalman filtering using a noise corrupted measurement matrix, IET Control Theory and Applications 3(9) (2009), 1226–1236.

19.

Y.C.

and Lee

R.C.K.

, A Bayesian approach to problems in stochastic estimation and control, IEEE Transaction on Automatic Control 9 (1964), 333–339.

20.

Chen

, Bayesian Filtering: From Kalman Filters to Particle Filters, and Beyond, Adaptive Systems Laboratory, McMaster University, Hamilton, ON, Canada, 2003.

Bayesian approach to linear state space model with unknown measurement matrix

Abstract

Keywords

1. Introduction

3. Bayesian approach with measurement matrix unknown

Table 1 The parameter values set for data generation in model Eq. (1)

Footnotes

Acknowledgments

References

Table 1
The parameter values set for data generation in model Eq. (1)