Variance Quantification

This project addresses quantification of the measurement noise induced error in the estimation of transfer function type model structures of the form

 y_t = G(q,theta)u_t + H(q,theta)e_t

where y_t is an oberved output, u_t is an observed input with power spectral density Phi_u(omega) and e_t is a zero mean white noise corruption with variance E{e_t^2}=sigma^2. The above transfer functions are assumed of the form

 G(q,theta) = q^{-k}frac{B(q,theta)}{A(q,theta)},quad H(q,theta) = frac{C(q,theta)}{D(q,theta)}


 A(q,theta) = 1+a_1q^{-1}+a_2q^{-2}+cdots+a_{n_a}q^{-n_a},qquad B(q,theta) = b_0+b_1q^{-1}+b_2q^{-2}+cdots+b_{n_b}q^{-n_b}
 D(q,theta) = 1+d_1q^{-1}+d_2q^{-2}+cdots+d_{n_d}q^{-n_d},qquad C(q,theta) = 1+c_1q^{-1}+c_2q^{-2}+cdots+c_{n_c}q^{-n_c}.

The parameter estimate widehattheta_N is assumed to be found by least squares estimation:

 widehattheta_N = mbox{arg}min_theta sum_{t=1}^N varepsilon_t^2(theta),qquad varepsilon_t(theta) = y_t-widehat y_{tmid t-1}(theta) = H^{-1}(q,theta)left[y_t - G(q,theta)u_tright].

It has been widely accepted that the variability of the ensuing estimate G(e^{jomega},widehattheta_N) of the system frequency response can be approximated by the formula

Existing Variance Quantification
 mbox{Var}{G(e^{jomega},widehattheta_N)}approx frac{m}{N}frac{sigma^2 |H(e^{jomega})|^2}{Phi_u(omega)}.

In the Ouput Error or Box–Jenkins model structure cases, if m_a=m_b then m=m_a.

This approximation depends on theoretical analysis that calculates the variance when m=infty, and then uses that quantity for the finite m in question. The quality of the approximation therefore depends on the speed of convergence with increasing m.

This project examines new analysis methods based on the principles of reproducing kernels and orthonormal parameterizations of model structures to derive a new variance approximation

New Variance Quantification
 mbox{Var}{G(e^{iomega},widehattheta_N)}approx frac{(m_a+m_b)}{N}frac{sigma^2 |H(e^{jomega})|^2}{Phi_u(omega)} sum_{k=1}^mfrac{1-|xi_k|^2}{|e^{iomega}-xi_k|^2}.

where the xi_k are the zeros of the true underlying denominator A(q)

 A(q) = prod_{k=1}^{m_a}(1-q^{-1}xi_k).

In some cases, such as white input excitation Phi_u(omega), this quantification is exact for finite m_a,m_b. In the non-white case it has been observed to offer quicker convergence with increasing m_a,m_b, and hence better approximation for finite m_a,m_b.

For example, consider the system

 G(s) = frac{0.0012(1-3.33s)^3}{(s+0.9163)^2(s+0.3567)^2(s+0.2231)^3 }.

At 1 second sample period, this has zero order hold equivalent G(q) representation with poles as

 A(q) = {(q-0.8)^3(q-0.7)^2(q-0.4)^2}

We conduct a simulation experiment where a 7'th order output error model structure is estimated on the basis of observing a length N=10000 sample input-output record from this system, and for which the output is corrupted by white Gaussian noise of variance sigma^2 = 0.01, and with input which is a realisation of a stationary Gaussian process with spectral density

 Phi_u(omega) = frac{1}{1.25-cosomega}.

The sample mean square error over 10000 estimation experiments with different input and noise realisations is used as an estimate of mbox{Var}{G(e^{jomega},widehattheta^n_N)}.

This is plotted as the green solid line in the figure below. The pre-exiting approximation is shown as the magenta dash-dot line. Our new approximation is shown as the blue dashed line.

Clearly, in this case it is a much better representation of the true variability (green line) than the pre-existing approximation.

Of course, this is just one example. If you are interesting exploring further, you may wish to download the following Matlab script which profiles both variance expressions together with the true variance (estimated via Monte-Carlo analysis) for a range of models.

Matlab script profiling variance expressions

The matlab script avar.m requires either the matlab system identification toolbox, or the alternative toolbox to run.

For more details on the theoretical analysis and other aspects underlying this approximation, please refer to the following publications. As you will see from the author list, this work is the results of a collaboration with Hakan Hjalmarsson and Fredrik Gustafsson.

Journal Papers

Conference Papers

Maintained by Prof. Brett Ninness
University of Newcastle
19 May 2010, © Copyright