A Dynamic Approach to Linear Statistical Calibration with an Application in Microwave Radiometry

Derick L. Rivers and Edward L. Boone Corresponding author. Email: riversdl@vcu.edu Department of Statistical Sciences and Operations Research, Virginia Commonwealth University

Abstract

The problem of statistical calibration of a measuring instrument can be framed both in a statistical context as well as in an engineering context. In the first, the problem is dealt with by distinguishing between the “classical" approach and the “inverse" regression approach. Both of these models are static models and are used to estimate “exact" measurements from measurements that are affected by error. In the engineering context, the variables of interest are considered to be taken at the time at which you observe it. The Bayesian time series analysis method of Dynamic Linear Models (DLM) can be used to monitor the evolution of the measures, thus introducing an dynamic approach to statistical calibration. The research presented employs the use of Bayesian methodology to perform statistical calibration. The DLM’s framework is used to capture the time-varying parameters that maybe changing or drifting over time. Two separate DLM based models are presented in this paper. A simulation study is conducted where the two models are compared to some well known ’static’ calibration approaches in the literature from both the frequentist and Bayesian perspectives. The focus of the study is to understand how well the dynamic statistical calibration methods performs under various signal-to-noise ratios, $r$ . The posterior distributions of the estimated calibration points as well as the $95\%$ coverage intervals are compared by statistical summaries. These dynamic methods are applied to a microwave radiometry dataset.

1. Introduction

Calibrating measurement instruments is a important problem that engineers frequently need to address. There exist several statistical methods that address this problem that are based on a simple linear regression approach. In tradition simple linear regression the goal is to relate a known value of X to a uncertain value of Y using a linear relationship. In contrast, the statistical calibration problem seeks to utilize a simple linear regression model to relate a known value of Y to an uncertain value of X. This is why statistical calibration is sometimes called inverse regression due to its relationship to simple linear regression (Osborne 1991; Ott and Longnecker 2009). Recall in linear regression the model is given as follows:

{\bf Y}={\bf X}{\bm{\beta}}+{\bm{\epsilon}}

(1)

where Y is a $(n\times 1)$ response vector, X is a $(n\times p)$ matrix of independent variables with $p=k+1$ total model parameters, ${\bm{\beta}}$ is a $(p\times 1)$ vector of unknown fixed parameters and ${\bm{\epsilon}}$ is a $(n\times 1)$ vector of uncorrelated error terms with zero mean (Myers 1990; Draper and Smith 1998; Montgomery et al. 2012). It is assumed that the value of the predictor variable X = x are nonrandom and observed with negligible error, while the $n$ error terms are random variables with mean zero and constant variance $\sigma^{2}$ (Myers 1990). Typically, in regression, of interest is the estimation of the parameter vector; ${\bm{\beta}}$ , and possibly the prediction of a future value $\hat{\bf Y}_{i|new}$ corresponding to a new ${\bf X}=x_{i|new}^{\prime}$ value. The prediction problem is relatively straightforward, due to the fact that a future ${\bf Y}_{i}$ value can be made directly by substituting $x_{i|new}^{\prime}$ into (1) with $E[\epsilon]=0$ .
For the statistical calibration problem let $y_{0}$ be the known observed value of the response and $x_{0}$ be the corresponding regressor, $x_{0}$ which is to be estimated. This problem is conducted in two stages: first measurement pairs $(x_{i},y_{i})$ of data is observed and a simple linear regression line is fit by estimating ${\bm{\beta}}$ ; secondly, $m$ observations of the response are observed, all corresponding to a single $x_{0}$ (Özyurt and Erar 2003). Since $y_{0}$ is fixed, inferences are different than those in a traditional regression (or prediction) problem (Osborne 1991; Eno 1999; Eno and Ye 2000).

1. Classical Calibration Methods

Eisenhart (1939) offered the first solution to the calibration problem, and is commonly known as the $``classical"$ estimator to the linear calibration problem. They assumed that the relationship between $x$ and $y$ was of a simple linear form:

E(Y|X=x)=\beta_{0}+\beta_{1}x.

The estimated regression line for the first stage of the experiment is given by

\hat{Y}=\hat{\beta_{0}}+\hat{\beta}_{1}X,

(2)

where $\hat{\beta_{0}}$ and $\hat{\beta_{1}}$ are the least squares estimate of $\beta_{0}$ and $\beta_{1}$ , respectively. Using the data collected at the first stage of experimentation, Eisenhart (1939) inverts Equation (2) to estimate the unknown regressor value $x_{0}$ for an observed response value $y_{0}$ , by:

\hat{x}_{0,c}=\frac{y_{0}-\hat{\beta}_{0}}{\hat{\beta}_{1}}

(3)

where $\hat{x}_{0,c}$ denotes the $``classical"$ estimator for $x_{0}$ . Since division by $\hat{\beta}_{1}$ is used there is an implicit assumption that $\lvert\hat{\beta}_{1}\rvert>0$ .
Assuming that $\lvert\hat{\beta}_{1}\rvert>0$ , Brown (1993) describes the following interval estimate corresponding to Eisenhart (1939):

\frac{y_{0}-\hat{\beta}_{0}}{\hat{\beta}_{1}}\left(1+\frac{\hat{\sigma}^{2}t^{2}}{\hat{\beta}_{1}^{2}S_{xx}}\right)\pm\frac{\hat{\sigma}t}{\hat{\beta}_{1}}\left(1+\frac{1}{2n}+\frac{(y_{0}-\hat{\beta}_{0})^{2}+\hat{\sigma}^{2}t^{2}}{2\hat{\beta}_{1}^{2}S_{xx}}\right),

where

\hat{\sigma}=\sqrt{\frac{\sum_{i=1}^{n}(y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1}x_{i})^{2}}{n-2}},

S_{xx}=\sum_{i=1}^{n}(x-\bar{x})^{2},

Krutchkoff (1967) proposed a competitive approach to Eisenhart’s (1939) classical linear calibration solution, which he called the $``inverse"$ regression calibration method and is written as:

X_{i}=\phi+\delta Y_{i}+\epsilon^{{}^{\prime}}_{i},

where $\phi$ and $\delta$ are the parameters in the linear relationship and $\epsilon^{{}^{\prime}}_{i}$ are independent identically distributed measurement errors with a zero mean and finite variance. Here $\phi$ and $\delta$ are estimated via least squares. The unknown $x_{0}$ can be estimated directly by substituting $y_{0}$ into the fitted equation:

\hat{x}_{0,I}=\hat{\phi}+\hat{\delta}y_{0}.

(4)

We let $\hat{x}_{0,I}$ denote the $``inverse"$ estimator of $x_{0}$ . The $100(1-\alpha)\%$ confidence interval for $E(x_{0,I}|y_{0})$ can be written as

x_{0,I}(y_{0})\pm t_{\nicefrac{{\alpha}}{{2}}}\hat{\sigma}\sqrt{\frac{1}{n}+\frac{(y_{0}-\bar{y})^{2}}{S_{yy}}}

where

S_{yy}=\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}.

Krutchkoff (1967) used a simulation study, where he found that the mean squared error of estimation for $x_{0}$ was uniformly less for this estimator versus the classical estimator. The inverse approach was later supported by Lwin and Maritz (1982). For criticisms of Krutchkoff’s (1967) approach such as bias see Osborne (1991).

2. Bayesian Calibration Methods

The first noted Bayesian solution to the calibration problem was presented by Hoadley (1970). His work was motivated by the unanswered question in the Frequentist community of whether $\beta_{1}$ is zero (or close to zero). Hoadley (1970) justified the use of the $``inverse"$ estimator (Krutchkoff, 1967) by considering the ususal $F$ -statistic to test the hypothesis that $\beta_{1}=0$ where $F=\hat{\beta}_{1}^{2}S_{xx}/\hat{\sigma}^{2}$ ,

\hat{\sigma}^{2}=\frac{\left\{\sum_{i=1}^{n}\left(y_{1i}-(\hat{\beta}_{0}+\hat{\beta}_{1}x_{i})\right)^{2}+\sum_{j=1}^{m}\left(y_{2j}-\bar{y}_{2}\right)^{2}\right\}}{(n+m-3)}.

The assumption made by Hoadley (1970) reflects that $x_{0}$ is random and a priori independent of $\pi(\beta_{0},\beta_{1},\sigma^{2})$ , so that the joint prior distribution of $\pi(\beta_{0},\beta_{1},\sigma^{2},x_{0})\propto\pi(\beta_{0},\beta_{1},\sigma^{2})\pi(x_{0})$ . Hoadley (1970) first assumed that $(\beta_{0},\beta_{1},\sigma^{2})$ had a uniform distribution,

\pi(\beta_{0},\beta_{1},\sigma^{2})\propto\sigma^{-2},

but the prior distribution for $x_{0}$ was not given.
Hoadley (1970) shows for $m=1$ (one observation at the prediction stage), that if $x_{0}$ has a prior density from a Student t distribution with $n-3$ degrees of freedom, a mean of 0, and a scale parameter

\sigma=\frac{n+1}{n-3},

the posterior distribution is

\pi(x_{0}|{\bf Data})=t_{n-2}\left(\hat{x}_{0,I},\left[\frac{n+1+(\hat{x}_{0,I})^{2}/R}{F+n-2}\right]\right),

(5)

where $\hat{x}_{0,I}$ is the inverse estimator given by (4) $,R=\frac{F}{F+n-2}$ and $F=\hat{\beta}_{1}^{2}S_{xx}/\hat{\sigma}^{2}$ .
Hunter and Lamboy (1981) also considered the calibration problem from a Bayesian point of view and is similar to that of Hoadley (1970) because both assume the prior distribution to be

\pi(\beta_{0},\beta_{1},\sigma^{2},\eta)\propto\sigma^{-2}

where $\eta=\beta_{0}+\beta_{1}x_{0}$ which is the predicted $y_{0}$ . The primary difference between their approach and the approach of Hoadley (1970) is that a priori they assume that $\eta$ and $(\beta_{0},\beta_{1},\sigma^{2})$ are independent while Hoadley (1970) assumed a priori that $x_{0}$ and $(\beta_{0},\beta_{1},\sigma^{2})$ are independent.
Hunter and Lamboy (1981) uses an approximation to the posterior distribution of the unknown regressor $x_{0}$ by

\pi(x_{0}|{\bf Data})=N\left(\hat{x}_{0,c},\frac{(s_{11}+s_{33})s_{22}-s_{12}^{2}}{s_{22}\hat{\beta}_{1}^{2}}\right),

(6)

where

{\bf S}=\{s_{i,j}\}=\left[\begin{array}[]{ccc}s_{11}&s_{12}&0\\ s_{12}&s_{22}&0\\ 0&0&s_{33}\end{array}\right]=\left[\begin{array}[]{cc}({\bf X}^{{}^{\prime}}{\bf X})^{-1}\hat{\sigma}^{2}&{\bf 0}\\ {\bf 0}&\nicefrac{{\hat{\sigma}^{2}}}{{m}}\end{array}\right],

with $\hat{x}_{0,c}$ being the classical estimator given in Equation (3), $s_{i,j}$ denote the element of the $i^{th}$ row and $j^{th}$ column from variance-covariance matrix of the joint posterior density of ( $\beta_{0}$ , $\beta_{1}$ , $\eta$ ).
The remainder of this paper is organized as follows. Section 2 presents the development of the dynamic approaches to the statistical calibration problem. In Section 3 the results from the simulation study where the dynamics methods are evaluated along with the static approaches are presented. In Section 4 the proposed methods are applied to microwave radiometer data. In Section 5 future work and other considerations are given.

2. Dynamic Calibration Approach

Traditional calibration methods assume the regression relationship is “static” in time. In many cases this is false, for example in microwave radiometry the static nature of the relationship is known to change across time. A dynamic approach can be created by letting the regression coefficients vary through time,

y_{t}=\beta_{0t}+\beta_{1t}x_{t}+\epsilon_{t},

where $\epsilon_{t}\stackrel{{\scriptstyle iid}}{{\sim}}N[0,\sigma^{2}_{t}]$ and is known as the $observational$ error.
The model may have different defining parameters at different times. One approach is to model $\beta_{0t}$ and $\beta_{1t}$ by using random walk type evolutions for the defining parameters, such as:

	$\displaystyle\beta_{0t}$	$\displaystyle=$	$\displaystyle\beta_{0(t-1)}+\omega_{\beta_{0t}},$
	$\displaystyle\beta_{1t}$	$\displaystyle=$	$\displaystyle\beta_{1(t-1)}+\omega_{\beta_{1t}},$

where $\omega_{\beta_{0t}}$ and $\omega_{\beta_{1t}}$ are independent zero-mean error terms with finite variances. At any time $t$ the calibration problem is given by:

y_{0t}=\beta_{0t}+\beta_{1t}x_{0t}+\epsilon_{t},\hskip 30.0ptt=1,2,\dots,T.

Bayesian Dynamic Linear Models (DLMs) approach of West et al. (1985); West and Harrison (1997) can be employed to achieve this goal. Recall the DLM framework is:

$\displaystyle\mbox{Observation equation}:\hskip 28.45274pt$	$\displaystyle{\bf Y}_{t}={\bf X}_{t}{\bf X}_{t}{\bm{\theta}_{t}}+{\bm{\epsilon}}_{t},\hskip 25.6073pt$	$\displaystyle{\bm{\epsilon}}_{t}\sim N_{r}[{\bf 0,E}]$
$\displaystyle\mbox{System equation}:\hskip 28.45274pt$	$\displaystyle{\bm{\theta}_{t}}={\bf G}_{t}{\bm{\theta}_{t-1}}+{\bm{\omega}_{t}},$	$\displaystyle{\bm{\omega}}_{t}\sim N_{d}[{\bf 0,W}]$
$\displaystyle\mbox{Initial information}:\hskip 28.45274pt$	$\displaystyle({\bm{\theta}_{0}}\|D_{0})\sim N_{d}[{\bf m_{0},C_{0}}],$

for some prior mean $\bf m_{0}$ and variance $\bf C_{0}$ with the vector of error terms, ${\bm{\epsilon}}_{t}$ and ${\bm{\omega}}_{t}$ independent across time and at any time.
To update the model through time West and Harrison (1997) give the following method:

(a)

Posterior distribution at $t-1$ : For some mean ${\bf m}_{t-1}$ and variance ${\bf C}_{t-1}$ ,
$({\bm{\theta}_{t-1}}|D_{t-1})\sim N_{d}[{\bf m}_{t-1},{\bf C}_{t-1}]$ .
(b)

Prior distribution at time $t$ : $({\bm{\theta}_{t}}|D_{t-1})\sim N_{d}[{\bf a}_{t},{\bf R}_{t}]$ , where
${\bf a}_{t}={\bf G}_{t}{\bf m}_{t-1}$ and ${\bf R}_{t}={\bf G}_{t}{\bf C}_{t-1}{\bf G}^{{}^{\prime}}_{t}+{\bf W}$ .
(c)

One-step forecast: $({\bf Y}_{t}|D_{t-1})\sim N_{r}[{\bf f}_{t},{\bf Q}_{t}]$ , where
${\bf f}_{t}={\bf X}_{t}{\bf a}_{t}$ and ${\bf Q}_{t}={\bf X}_{t}{\bf R}_{t}{\bf X}_{t}^{{}^{\prime}}+{\bf E}$ .
(d)

Posterior distribution at time $t$ : $({\bm{\theta}_{t}}|D_{t})\sim N_{d}[{\bf m}_{t},{\bf C}_{t}]$ , with
${\bf m}_{t}={\bf a}_{t}+{\bf A}_{t}{\bf e}_{t}$ and ${\bf C}_{t}={\bf R}_{t}-{\bf A}^{{}^{\prime}}_{t}{\bf Q}_{t}{\bf A}_{t},$
where
${\bf A}_{t}={\bf Q}^{-1}_{t}{\bf X}_{t}{\bf R}_{t}$ and ${\bf e}_{t}={\bf Y}_{t}-{\bf f}_{t}$ .

The DLM framework is used to establish the evolving relationship between the fixed design matrix ${\bf X}_{t}$ and ${\bf Y}_{t}$ by estimating ${\bm{\theta}}_{t}$ , which is a $(d\times n)$ matrix of time-varying regression coefficients $\beta_{0t}$ and $\beta_{1t}$ . For our calibration situation ${\bf Y}_{t}$ is a $(r\times n)$ matrix of responses and ${\bf G}_{t}$ is a known ( $d\times d$ ) system matrix. The error ${\bm{\epsilon}}_{t}$ and ${\bm{\omega}}_{t}$ are independent normally distributed random $(r\times n)$ matrices with zero mean and constant variance-covariance matrices E and W. For simplification ${\bf G}_{t}$ is set equal to ${\bf I}_{(d\times d)}$ , ${\bf E}$ is set equal to $\sigma^{2}_{E}{\bf I}_{(r\times r)}$ and ${\bf W}$ is $\sigma^{2}_{W}\left[{\bf X}_{t}^{{}^{\prime}}{\bf X}_{t}\right]^{-1}$ . The past information is contained in the set $D_{0}$ .
We specify a prior in the first stage of calibration for the unknown variances and derive an algorithm to draw from the posterior distribution of the unknown parameters,

\pi({\bm{\theta}}_{t},\sigma^{2}_{E},\sigma^{2}_{W}|{\bf Y}_{t})\propto\pi({\bm{\theta}}_{t}|\sigma^{2}_{E},\sigma^{2}_{W},{\bf Y}_{t})\pi(\sigma^{2}_{E},\sigma^{2}_{W}|{\bf Y}_{t}).

The second stage of the calibration experiment consists of using the joint posterior distribution $\pi({\bm{\theta}}_{t},\sigma^{2}_{E},\sigma^{2}_{W}|{\bf Y}_{t})$ to derive $x_{0t}|{\bm{\theta}}_{t},\sigma^{2}_{E},\sigma^{2}_{W}$ for each draw of $\pi(\sigma^{2}_{E},\sigma^{2}_{W}|{\bf Y}_{t})$ . The estimator for the parameter of interest, $x_{0t}$ , is defined in a manner akin to Eisenhart (1939); Hunter and Lamboy (1981); Eno (1999), where

x_{0t}=\frac{y_{0t}-\beta_{0t}}{\beta_{1t}}.

(7)

In the final stage of the calibration experiment, the posterior distribution summary statistics are gathered at each time point $t$ . The posterior median and credible intervals are taken for each $t$ across the draws of $x_{0t}|{\bm{\theta}}_{t},\sigma^{2}_{E},\sigma^{2}_{W}$ . The result of the dynamic calibration experiment is a time series of calibration distributions across time. We will be able to observe the distributional changes of the system with respect to the calibration reference.
The proposed calibration estimator is developed by first considering the joint posterior distribution $\pi(\sigma^{2}_{E},\sigma^{2}_{W}|{\bf Y}_{t})$ . We let ${\bm{\Gamma}}$ denote the vector of unknown DLM dispersion parameters where ${\bm{\Gamma}}^{\prime}=(\sigma^{2}_{E},\sigma^{2}_{W})$ . The prior information for the dispersion parameters is described by a prior density $\pi({\bm{\Gamma}})$ which summarizes what is known about the variance parameters before any data are observed. Using the Bayesian inferential approach, the prior information about the parameters must be combined with information contained in the data. The information provided by the data is captured by the likelihood functions, $f_{{\bf Y}}({\bf Y}_{t}|{\bm{\theta}}_{t},\sigma^{2}_{E},\sigma^{2}_{W})$ and $f_{\bm{\theta}}({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1},\sigma^{2}_{W})$ for the observation equation and the system equation, respectively. The combined information is described by the posterior density using the Bayes theorem (Bernardo and Smith 1994) as

\pi({\bm{\Gamma}}|{\bf Y}_{t})\propto f_{{\bf Y}}({\bf Y}_{t}|{\bm{\theta}}_{t},\sigma^{2}_{E},\sigma^{2}_{W})\cdot f_{\bm{\theta}}({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1},\sigma^{2}_{W})\cdot\pi({\bm{\Gamma}}).

For our calibration problem it is believe that $\sigma^{2}_{E}>\sigma^{2}_{W}$ . To deal with the variance relationship we specify the following prior distributions:

	$\displaystyle\sigma^{2}_{E}$	$\displaystyle\sim$	$\displaystyle Uniform(0,1)$		(8)
	$\displaystyle\sigma^{2}_{W}\|\sigma^{2}_{E}$	$\displaystyle\sim$	$\displaystyle Uniform(0,\sigma^{2}_{E}).$		(9)

Prior distributions (8) and (9) ensures the system variance to be less than the observation variance. Since these are proper prior distributions the resulting posterior distribution will also be proper.
In the first stage of calibration, the joint distribution of the observations, states, and unknown parameters is as follows:

$\displaystyle\pi({\bf Y}_{1:T},{\bm{\theta}}_{0:T},\sigma^{2}_{E},\sigma^{2}_{W})$	$\displaystyle=$	$\displaystyle f_{{\bf Y}}({\bf Y}_{1:T}\|{\bm{\theta}}_{0:T},\sigma^{2}_{E},\sigma^{2}_{W})\cdot f_{\bm{\theta}}({\bm{\theta}}_{0:T}\|\sigma^{2}_{W})\cdot\pi({\bm{\Gamma}})$
	$\displaystyle=$	$\displaystyle\prod_{t=1}^{T}f_{{\bf Y}}({\bf Y}_{t}\|{\bm{\theta}}_{t},\sigma^{2}_{E})\cdot\prod_{t=1}^{T}f_{\bm{\theta}}({\bm{\theta}}_{t}\|{\bm{\theta}}_{t-1},\sigma^{2}_{W})$
		$\displaystyle\hskip 10.0pt\cdot\pi({\bm{\theta}}_{0})\cdot\pi(\sigma^{2}_{E})\cdot\pi(\sigma^{2}_{W}\|\sigma^{2}_{E}).$

where the likelihood for the observation equation is

f_{{\bf Y}}({\bf Y}_{t}|{\bm{\theta}}_{t},\sigma^{2}_{E})\propto\sigma^{-T}_{E}\mbox{exp}\left\{-\frac{1}{2\sigma^{2}_{E}}\sum_{t=1}^{T}({\bf Y}_{t}-{\bf X}_{t}{\bm{\theta}}_{t})^{2}\right\}

and the likelihood for the system equation is

f_{{\bm{\theta}}}({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1},\sigma^{2}_{W})\propto\sigma^{-T}_{W}\mbox{exp}\left\{-\frac{1}{2\sigma^{2}_{W}}\sum_{t=1}^{T}({\bm{\theta}}_{t}-{\bm{\theta}}_{t-1})^{2}\right\}.

Given the joint distribution above, the posterior distribution is

\pi({\bf x}_{0t}|{\bm{\theta}}_{t},{\bm{\Gamma}},{\bf Y}_{t})

(10)

where

{\bf x}_{0t}=\frac{\bf{y}^{*}_{0t}}{{\bm{\theta}}_{t}}

(11)

and ${\bf y}^{*}_{0t}={\bf y}_{0t}-\bar{y}_{t}$ (i.e. $\bar{y}_{t}$ is the cumulative mean of the observations up to time $t$ ) and ${\bm{\theta}}_{t}=\hat{\bm{\beta}}_{1t}$ . Samples from the posterior distribution in Equation (10) are drawn by implementing the Sampling Importance Resampling (Albert 2007; Givens and Hoeting 2005) approach.
The development of the estimator in Equation (11) is deterministic in approach. We present a fully Bayesian approach to dynamic calibration that incorporates the uncertainty in estimation. The second dynamic calibration model is derived by Bayes’ theorem

\pi({\bf x_{0t}}|{\bf Y}_{t})\propto\pi({\bf x_{0t}})f({\bf Y}_{t}|{\bf x_{0t}}),

where $\pi({\bf x_{0t}}|{\bf Y}_{t})$ is the posterior distribution for ${\bf x}_{0t}$ . The prior belief for the calibration values is denoted as $\pi({\bf x_{0t}})$ with the $f({\bf Y}_{t}|{\bf x_{0t}})$ denoting the likelihood function.
The objective of any Bayesian approach is to obtain the posterior distribution from which inferences can be made. Here the desired posterior is

\pi({\bf x_{0t}}|{\bf Y}_{t})

(12)

which must be dynamic through time. We determine the posterior distribution (12) in a similiar manner as described above in Equations (10) and (11). In the first stage of the calibration experiment the data is scaled and centered, therefore setting the $y-$ intercept equal to zero and the reference measurements centered at zero. Centering of the data is used to reduce the parameter space. The posterior distribution can be thought of as:

\pi({\bf z_{0t}}|{\bf Y}^{*}_{t})\propto\pi({\bf z_{0t}})f({\bf Y}^{*}_{t}|{\bf z_{0t}}),

(13)

with ${\bf z_{0t}}$ representing the transformed calibrated value at time $t$ and ${\bf Y}^{*}_{t}={\bf Y}_{t}-\bar{Y}_{t}$ , where $\bar{Y}_{t}$ is the cumulative mean of the observations. Given this information $a~priori$ we define the prior distribution

\pi({\bf z_{0t}})=N(0,1).

The posterior density in Equation (13) is defined as

\pi({\bf z_{0t}}|{\bf Y}^{*}_{t})\propto\mbox{exp}\left\{-\frac{1}{2}\left[\sigma^{-2}_{Y_{t}}\sum_{t=1}^{T}({\bm{\xi}}_{t}-{\bf z}_{0t})^{2}+{\bf z}^{2}_{0t}\right]\right\}

(14)

where ${\bm{\xi}}_{t}=\nicefrac{{{\bf Y}^{*}_{0t}}}{{\bm{\theta}_{t}}}$ . Applying Bayes theorem and completing the square, the posterior distribution is

\pi({\bf z_{0t}}|{\bf Y}^{*}_{t})\sim N(\mu_{z_{0t}},\sigma^{2}_{z_{0t}}),

(15)

with

	$\displaystyle\mu_{z_{0t}}$	$\displaystyle=$	$\displaystyle\frac{{\bm{\xi}}_{t}}{1+\sigma^{2}_{Y_{t}}},$
	$\displaystyle\sigma^{2}_{z_{0t}}$	$\displaystyle=$	$\displaystyle\frac{1}{1+\sigma^{2}_{Y_{t}}}$

and

\sigma^{2}_{Y_{t}}=\mbox{tr}({\bf Q}_{t}).

where tr( . ) denotes trace of the one-step forecast variance-covariance matrix. We derive the posterior in Equation (12) by drawing from Equation (15) and transforming the data back to the original scale as so:

{\bf x}_{0t}=\bar{X}+{\bf z}_{0t}\sigma_{X},

(16)

where $\bar{X}$ is the mean of the reference measurements vector and $\sigma_{X}$ is the standard deviation of the reference measurements vector.
The dynamic calibration algorithm is developed for both of the approaches using R (R Development Core Team, 2013) and is conducted as below.

Algorithm: Dynamic Calibration 1. Generate

M

proposal samples for

(\sigma^{2}_{E},\sigma^{2}_{W})

from

\pi(\sigma^{2}_{E})

and

\pi(\sigma^{2}_{W}|\sigma^{2}_{E})

; 2. Calibration data are fit using the DLM framework for each of the

M

proposal samples

(\sigma^{2(m)}_{E},\sigma^{2(m)}_{W})

, with the prior moments for

({\bm{\theta}_{0}}|D_{0})

\bf m_{0}=1_{d}

and

\bf C_{0}=100I_{(d\times d)}

, where

{\bf 1_{d}}

is a

d-

dimensional vector of ones.

a.

Data are scaled and shifted such that $\sum^{r}_{i=1}x_{i}=0$ , $\frac{1}{n}\sum^{r}_{i=1}x^{2}_{i}=1$ and $y-$ intercept $=0$ , where $y^{*}_{t}=y_{t}-\bar{y}_{t}$ for all $t$ (i.e. $\bar{y}_{t}$ is the cumulative mean up to time $t$ );
b.

Estimate ${\bm{\theta}}^{(m)}_{t}|\sigma^{2(m)}_{E},\sigma^{2(m)}_{W}$ for the $m^{th}$ proposal sample is calculated $\mbox{for all}~t$ ;
c.

Estimate $x^{(m)}_{0t}|{\bm{\theta}}^{(m)}_{t},\sigma^{2(m)}_{E},\sigma^{2(m)}_{W}$ for the $m^{th}$ proposal sample is calculated $\mbox{for all}~t$ , using either Equation (11) or drawing from Equation (15);
d.

Calculate log-likelihood density weights, $log[f({\bm{\Gamma}}^{(m)})]$ , for each $(\sigma^{2(m)}_{E},\sigma^{2(m)}_{W})$ pair

Sampling Importance Resampling (SIR) is used to simulate samples of $x_{0t}|{\bm{\theta}}_{t},\sigma^{2}_{E},\sigma^{2}_{W}$ by accepting a subset of $N=1,000$ from the proposal density to be distributed according to the posterior density $\pi({\bm{\Gamma}}|{\bf Y}_{t})$ with candidate density $\pi({\bm{\Gamma}})$ .

a.

Calculate the standardized importance weights, $w({{\bm{\Gamma}}^{(1)}}),\dots,w({{\bm{\Gamma}}^{(M)}})$ , where $w({{\bm{\Gamma}^{(m)}}})=log[f({\bm{\Gamma}}^{(m)})]-log[g({\bm{\Gamma}}^{(m)})]$ for the $m^{th}$ proposal sample;
b.

Sample $N$ calibrated time series from the $M$ proposal values with replacement given probabilities $p({{\bm{\Gamma}}^{(m)}})$ where

$p({{\bm{\Gamma}}^{(m)}})=\frac{e^{w({{\bm{\Gamma}^{(m)}}})}}{\sum_{j=1}^{M}e^{w({{\bm{\Gamma}^{(j)}}})}}.$

Rescale calibrated time series to original scale by Equation (16) and take summary statistics (i.e. medians and credible sets) across each time $t$ .

Constant $g_{t}=0$
			$r=10$			$r=100$			$r=1000$
Ref.	Model	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW
2	$M_{D1}$	0.0008	0.995	2.519	0.0035	0.983	2.523	0.0307	0.939	2.517
	$M_{D2}$	0.0012	1.000	3.782	0.0038	1.000	3.782	0.0308	1.000	3.782
	$M_{F1}$	0.0001	1.000	1.224	0.0012	1.000	3.868	0.0123	1.000	12.229
	$M_{F2}$	0.0001	1.000	1.223	0.0016	1.000	3.863	0.0335	1.000	12.168
	$M_{B1}$	0.0002	0.997	1.182	0.0022	1.000	3.866	0.0386	1.000	12.177
	$M_{B2}$	0.0014	1.000	1.458	0.0139	1.000	4.606	0.1391	1.000	14.565
5	$M_{D1}$	0.0008	0.995	2.496	0.0035	0.983	2.509	0.0307	0.941	2.514
	$M_{D2}$	0.0013	1.000	3.983	0.0039	1.000	3.983	0.0307	1.000	3.983
	$M_{F1}$	0.0001	1.000	1.223	0.0012	1.000	3.865	0.0123	1.000	12.220
	$M_{F2}$	0.0001	1.000	1.222	0.0022	1.000	3.860	0.0813	1.000	12.113
	$M_{B1}$	0.0002	1.000	1.223	0.0023	1.000	3.861	0.0792	1.000	12.116
	$M_{B2}$	0.0014	1.000	1.457	0.0139	1.000	4.604	0.1069	1.000	10.748

Constant $g_{t}=0$
			$r=2$			$r=20$			$r=200$
Ref.	Model	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW
2	$M_{D1}$	0.0012	0.992	2.519	0.0041	0.981	2.520	0.0323	0.939	2.528
	$M_{D2}$	0.0015	1.000	3.782	0.0044	1.000	3.782	0.0325	1.000	3.782
	$M_{F1}$	0.0001	1.000	1.230	0.0010	1.000	3.871	0.0114	1.000	12.231
	$M_{F2}$	0.0001	1.000	1.229	0.0012	1.000	3.866	0.0314	1.000	12.170
	$M_{B1}$	0.0001	1.000	1.230	0.0019	1.000	3.869	0.0371	1.000	12.179
	$M_{B2}$	0.0190	1.000	1.155	0.0243	1.000	3.767	0.1381	1.000	14.567
5	$M_{D1}$	0.0011	0.992	2.508	0.0041	0.981	2.510	0.032	0.939	2.514
	$M_{D2}$	0.0017	1.000	3.983	0.0045	1.000	3.983	0.032	1.000	3.983
	$M_{F1}$	0.0001	1.000	1.228	0.0010	1.000	3.868	0.011	1.000	12.222
	$M_{F2}$	0.0001	1.000	1.227	0.0019	1.000	3.863	0.081	1.000	12.114
	$M_{B1}$	0.0001	1.000	1.227	0.0021	1.000	3.864	0.082	1.000	12.118
	$M_{B2}$	0.0013	1.000	1.462	0.0137	1.000	4.607	0.138	1.000	14.560

Stepped $g_{t}=a_{i}$
			$r=10$			$r=100$			$r=1000$
Ref.	Model	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW
2	$M_{D1}$	0.0191	0.961	2.509	0.0198	0.953	2.506	0.0406	0.926	2.543
	$M_{D2}$	0.0196	1.000	3.782	0.0201	1.000	3.782	0.0408	1.000	3.783
	$M_{F1}$	0.0001	1.000	9.094	0.0004	1.000	9.813	0.0094	1.000	15.209
	$M_{F2}$	0.0046	1.000	9.065	0.0073	1.000	9.779	0.0528	1.000	15.098
	$M_{B1}$	0.0859	1.000	9.072	0.0838	1.000	9.786	0.1866	1.000	15.109
	$M_{B2}$	0.1399	1.000	10.830	0.0823	1.000	11.687	0.1836	1.000	18.115
5	$M_{D1}$	0.0191	0.961	2.510	0.0197	0.954	2.511	0.0405	0.924	2.516
	$M_{D2}$	0.0196	1.000	3.983	0.0201	1.000	3.983	0.0405	1.000	3.983
	$M_{F1}$	0.0001	1.000	9.087	0.0004	1.000	9.806	0.0094	1.000	15.199
	$M_{F2}$	0.0184	1.000	9.041	0.0267	1.000	9.749	0.1620	1.000	14.995
	$M_{B1}$	0.0199	1.000	9.044	0.0267	1.000	9.752	0.1559	1.000	14.999
	$M_{B2}$	0.0706	1.000	10.826	0.0618	1.000	8.625	0.1742	1.000	15.091

Stepped $g_{t}=a_{i}$
			$r=2$			$r=20$			$r=200$
Ref.	Model	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW
2	$M_{D1}$	0.0209	0.957	2.520	0.0219	0.950	2.522	0.0436	0.921	2.526
	$M_{D2}$	0.0214	1.000	3.782	0.0222	1.000	3.782	0.0438	1.000	3.782
	$M_{F1}$	0.0001	1.000	9.103	0.0003	1.000	9.822	0.0086	1.000	15.216
	$M_{F2}$	0.0047	1.000	9.075	0.0073	1.000	9.788	0.0511	1.000	15.105
	$M_{B1}$	0.0084	1.000	9.081	0.0115	1.000	9.795	0.0601	1.000	15.116
	$M_{B2}$	0.0709	1.000	10.842	0.0826	1.000	11.698	0.2054	1.000	18.122
5	$M_{D1}$	0.0209	0.957	2.509	0.0218	0.949	2.511	0.0436	0.920	2.516
	$M_{D2}$	0.0214	1.000	3.983	0.0221	1.000	3.983	0.0435	1.000	3.983
	$M_{F1}$	0.0001	1.000	9.096	0.0003	1.000	9.815	0.0086	1.000	15.205
	$M_{F2}$	0.0185	1.000	9.050	0.0267	1.000	9.758	0.1616	1.000	15.002
	$M_{B1}$	0.0199	1.000	9.053	0.0281	1.000	9.761	0.1641	1.000	15.006
	$M_{B2}$	0.0708	1.000	10.836	0.0825	1.000	11.693	0.2052	1.000	18.114

Sinusoidal $g_{t}=0.1\mbox{sin}(0.025t)$
			$r=10$			$r=100$			$r=1000$
Ref.	Model	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW	AvMSE	AvCP	AvIW
2	$M_{D1}$	4.4088	0.628	2.657	4.4794	0.629	2.648	4.7214	0.638	2.681
	$M_{D2}$	4.4002	0.829	3.783	4.4708	0.825	3.783	4.7123	0.810	3.783
	$M_{F1}$	0.0001	1.000	21.980	0.0012	1.000	22.307	0.0123	1.000	25.206
	$M_{F2}$	0.1541	1.000	21.665	0.1670	1.000	21.978	0.2943	1.000	24.738
	$M_{B1}$	0.1689	0.975	20.933	0.1868	1.000	21.994	0.3174	1.000	24.757
	$M_{B2}$	0.4127	1.000	26.178	0.4258	1.000	26.567	0.5531	1.000	30.020
5	$M_{D1}$	4.4087	0.628	2.646	4.4793	0.630	2.648	4.7214	0.635	2.653
	$M_{D2}$	4.3906	0.845	3.984	4.4609	0.839	3.984	4.7023	0.824	3.984
	$M_{F1}$	0.0001	1.000	21.964	0.0012	1.000	22.291	0.0123	1.000	25.188
	$M_{F2}$	0.5810	1.000	21.371	0.6218	1.000	21.671	1.0152	1.000	24.306
	$M_{B1}$	0.5956	1.000	21.377	0.5909	1.000	21.678	0.9628	1.000	24.314
	$M_{B2}$	0.4123	1.000	26.166	0.3087	1.000	18.973	0.4658	1.000	25.009

A Dynamic Approach to Linear Statistical Calibration with an Application in Microwave Radiometry

Abstract

1. Introduction

1. Classical Calibration Methods

2. Bayesian Calibration Methods

2. Dynamic Calibration Approach

3. Simulation Study

1. Interpolation case

2. Extrapolation case

4. Application to Microwave Radiometer

5. Discussion

References

2 References-Sinusoidal Gain- w $/$ Burn In $=200$
	Mean Squared Error	Coverage Probability	Interval Width
$M_{D1}$	0.63553	0.72185	1.15722
$M_{D2}$	0.63333	0.96380	3.77191