2018-07-09 15:46:34 +00:00
|
|
|
/*!
|
|
|
|
* \file bayesian_estimation.cc
|
|
|
|
* \brief Interface of a library with Bayesian noise statistic estimation
|
|
|
|
*
|
|
|
|
* Bayesian_estimator is a Bayesian estimator which attempts to estimate
|
|
|
|
* the properties of a stochastic process based on a sequence of
|
|
|
|
* discrete samples of the sequence.
|
|
|
|
*
|
2019-07-21 17:32:52 +00:00
|
|
|
* [1]: LaMountain, Gerald, Vilà-Valls, Jordi, Closas, Pau, "Bayesian
|
|
|
|
* Covariance Estimation for Kalman Filter based Digital Carrier
|
2019-05-20 16:42:43 +00:00
|
|
|
* Synchronization," Proceedings of the 31st International Technical Meeting
|
2019-07-21 17:32:52 +00:00
|
|
|
* of the Satellite Division of The Institute of Navigation
|
2019-05-20 16:42:43 +00:00
|
|
|
* (ION GNSS+ 2018), Miami, Florida, September 2018, pp. 3575-3586.
|
|
|
|
* https://doi.org/10.33012/2018.15911
|
2018-07-09 15:46:34 +00:00
|
|
|
*
|
|
|
|
* \authors <ul>
|
|
|
|
* <li> Gerald LaMountain, 2018. gerald(at)ece.neu.edu
|
|
|
|
* <li> Jordi Vila-Valls 2018. jvila(at)cttc.es
|
|
|
|
* </ul>
|
2020-07-28 14:57:15 +00:00
|
|
|
* -----------------------------------------------------------------------------
|
2018-07-09 15:46:34 +00:00
|
|
|
*
|
2020-12-30 12:35:06 +00:00
|
|
|
* GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
|
2018-07-09 15:46:34 +00:00
|
|
|
* This file is part of GNSS-SDR.
|
|
|
|
*
|
2020-12-30 12:35:06 +00:00
|
|
|
* Copyright (C) 2010-2020 (see AUTHORS file for a list of contributors)
|
2020-02-08 00:20:02 +00:00
|
|
|
* SPDX-License-Identifier: GPL-3.0-or-later
|
2018-07-09 15:46:34 +00:00
|
|
|
*
|
2020-07-28 14:57:15 +00:00
|
|
|
* -----------------------------------------------------------------------------
|
2018-07-09 15:46:34 +00:00
|
|
|
*/
|
|
|
|
|
|
|
|
#include "bayesian_estimation.h"
|
2018-08-21 13:20:48 +00:00
|
|
|
|
2018-07-09 15:46:34 +00:00
|
|
|
|
|
|
|
Bayesian_estimator::Bayesian_estimator()
|
2021-12-15 01:48:18 +00:00
|
|
|
: kappa_prior(0),
|
|
|
|
nu_prior(0)
|
2018-07-09 15:46:34 +00:00
|
|
|
{
|
2021-12-15 01:48:18 +00:00
|
|
|
mu_prior = arma::zeros(1, 1);
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_est = mu_prior;
|
2021-12-15 01:48:18 +00:00
|
|
|
Psi_prior = arma::eye(1, 1) * (nu_prior + 1 + 1);
|
2018-08-21 12:57:22 +00:00
|
|
|
Psi_est = Psi_prior;
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
2019-07-21 17:32:52 +00:00
|
|
|
|
2018-07-09 15:46:34 +00:00
|
|
|
Bayesian_estimator::Bayesian_estimator(int ny)
|
2021-12-15 01:48:18 +00:00
|
|
|
: kappa_prior(0),
|
|
|
|
nu_prior(0)
|
2018-07-09 15:46:34 +00:00
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_prior = arma::zeros(ny, 1);
|
|
|
|
mu_est = mu_prior;
|
2021-12-15 01:48:18 +00:00
|
|
|
|
|
|
|
Psi_prior = arma::eye(ny, ny) * (nu_prior + ny + 1);
|
2018-08-21 12:57:22 +00:00
|
|
|
Psi_est = Psi_prior;
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
2019-07-21 17:32:52 +00:00
|
|
|
|
2021-12-15 01:48:18 +00:00
|
|
|
Bayesian_estimator::Bayesian_estimator(const arma::vec& mu_prior_0,
|
|
|
|
int kappa_prior_0,
|
|
|
|
int nu_prior_0,
|
|
|
|
const arma::mat& Psi_prior_0)
|
|
|
|
: mu_prior(mu_prior_0),
|
|
|
|
Psi_prior(Psi_prior_0),
|
|
|
|
kappa_prior(kappa_prior_0),
|
|
|
|
nu_prior(nu_prior_0)
|
2018-07-09 15:46:34 +00:00
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_est = mu_prior;
|
|
|
|
Psi_est = Psi_prior;
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
|
2018-08-21 12:57:22 +00:00
|
|
|
void Bayesian_estimator::init(const arma::mat& mu_prior_0, int kappa_prior_0, int nu_prior_0, const arma::mat& Psi_prior_0)
|
2018-08-14 01:18:01 +00:00
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_prior = mu_prior_0;
|
2018-08-14 01:18:01 +00:00
|
|
|
kappa_prior = kappa_prior_0;
|
2018-08-21 12:57:22 +00:00
|
|
|
nu_prior = nu_prior_0;
|
|
|
|
Psi_prior = Psi_prior_0;
|
2018-08-14 01:18:01 +00:00
|
|
|
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_est = mu_prior;
|
|
|
|
Psi_est = Psi_prior;
|
2018-08-14 01:18:01 +00:00
|
|
|
}
|
|
|
|
|
2019-07-21 17:32:52 +00:00
|
|
|
|
2018-07-09 15:46:34 +00:00
|
|
|
/*
|
|
|
|
* Perform Bayesian noise estimation using the normal-inverse-Wishart priors stored in
|
|
|
|
* the class structure, and update the priors according to the computed posteriors
|
|
|
|
*/
|
2018-08-21 12:57:22 +00:00
|
|
|
void Bayesian_estimator::update_sequential(const arma::vec& data)
|
2018-07-09 15:46:34 +00:00
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
int K = data.n_cols;
|
2018-07-09 15:46:34 +00:00
|
|
|
int ny = data.n_rows;
|
|
|
|
|
|
|
|
if (mu_prior.is_empty())
|
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_prior = arma::zeros(ny, 1);
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
if (Psi_prior.is_empty())
|
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
Psi_prior = arma::zeros(ny, ny);
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
arma::vec y_mean = arma::mean(data, 1);
|
2018-08-21 12:57:22 +00:00
|
|
|
arma::mat Psi_N = arma::zeros(ny, ny);
|
2018-07-09 15:46:34 +00:00
|
|
|
|
|
|
|
for (int kk = 0; kk < K; kk++)
|
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
Psi_N = Psi_N + (data.col(kk) - y_mean) * ((data.col(kk) - y_mean).t());
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
2018-08-21 12:57:22 +00:00
|
|
|
arma::vec mu_posterior = (kappa_prior * mu_prior + K * y_mean) / (kappa_prior + K);
|
|
|
|
int kappa_posterior = kappa_prior + K;
|
|
|
|
int nu_posterior = nu_prior + K;
|
2020-07-10 17:52:28 +00:00
|
|
|
arma::mat Psi_posterior = Psi_prior + Psi_N + (static_cast<float>(kappa_prior) * static_cast<float>(K)) / (static_cast<float>(kappa_prior) + static_cast<float>(K)) * (y_mean - mu_prior) * ((y_mean - mu_prior).t());
|
2018-07-09 15:46:34 +00:00
|
|
|
|
|
|
|
mu_est = mu_posterior;
|
|
|
|
if ((nu_posterior - ny - 1) > 0)
|
|
|
|
{
|
|
|
|
Psi_est = Psi_posterior / (nu_posterior - ny - 1);
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
Psi_est = Psi_posterior / (nu_posterior + ny + 1);
|
|
|
|
}
|
|
|
|
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_prior = mu_posterior;
|
2018-07-09 15:46:34 +00:00
|
|
|
kappa_prior = kappa_posterior;
|
2018-08-21 12:57:22 +00:00
|
|
|
nu_prior = nu_posterior;
|
|
|
|
Psi_prior = Psi_posterior;
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
/*
|
|
|
|
* Perform Bayesian noise estimation using a new set of normal-inverse-Wishart priors
|
|
|
|
* and update the priors according to the computed posteriors
|
|
|
|
*/
|
2018-08-21 12:57:22 +00:00
|
|
|
void Bayesian_estimator::update_sequential(const arma::vec& data, const arma::vec& mu_prior_0, int kappa_prior_0, int nu_prior_0, const arma::mat& Psi_prior_0)
|
2018-07-09 15:46:34 +00:00
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
int K = data.n_cols;
|
2018-07-09 15:46:34 +00:00
|
|
|
int ny = data.n_rows;
|
|
|
|
|
|
|
|
arma::vec y_mean = arma::mean(data, 1);
|
|
|
|
arma::mat Psi_N = arma::zeros(ny, ny);
|
|
|
|
|
|
|
|
for (int kk = 0; kk < K; kk++)
|
|
|
|
{
|
2018-08-21 12:57:22 +00:00
|
|
|
Psi_N = Psi_N + (data.col(kk) - y_mean) * ((data.col(kk) - y_mean).t());
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
2018-08-21 12:57:22 +00:00
|
|
|
arma::vec mu_posterior = (kappa_prior_0 * mu_prior_0 + K * y_mean) / (kappa_prior_0 + K);
|
|
|
|
int kappa_posterior = kappa_prior_0 + K;
|
|
|
|
int nu_posterior = nu_prior_0 + K;
|
2020-07-10 17:52:28 +00:00
|
|
|
arma::mat Psi_posterior = Psi_prior_0 + Psi_N + (Psi_prior_0 * static_cast<double>(K)) / (static_cast<double>(kappa_prior_0) + static_cast<double>(K)) * (y_mean - mu_prior_0) * ((y_mean - mu_prior_0).t());
|
2018-07-09 15:46:34 +00:00
|
|
|
|
|
|
|
mu_est = mu_posterior;
|
|
|
|
if ((nu_posterior - ny - 1) > 0)
|
|
|
|
{
|
|
|
|
Psi_est = Psi_posterior / (nu_posterior - ny - 1);
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
Psi_est = Psi_posterior / (nu_posterior + ny + 1);
|
|
|
|
}
|
|
|
|
|
2018-08-21 12:57:22 +00:00
|
|
|
mu_prior = mu_posterior;
|
2018-07-09 15:46:34 +00:00
|
|
|
kappa_prior = kappa_posterior;
|
2018-08-21 12:57:22 +00:00
|
|
|
nu_prior = nu_posterior;
|
|
|
|
Psi_prior = Psi_posterior;
|
2018-07-09 15:46:34 +00:00
|
|
|
}
|
|
|
|
|
2019-07-21 17:32:52 +00:00
|
|
|
|
2018-08-21 13:20:48 +00:00
|
|
|
arma::mat Bayesian_estimator::get_mu_est() const
|
2018-07-09 15:46:34 +00:00
|
|
|
{
|
|
|
|
return mu_est;
|
|
|
|
}
|
|
|
|
|
2019-07-21 17:32:52 +00:00
|
|
|
|
2018-08-21 13:20:48 +00:00
|
|
|
arma::mat Bayesian_estimator::get_Psi_est() const
|
2018-07-09 15:46:34 +00:00
|
|
|
{
|
|
|
|
return Psi_est;
|
|
|
|
}
|