evar

Variance estimation

Download
Syntax
Description
1-D examples
2-D examples
Information & Reference
See also
About the author

Download

evar.m

Syntax

   v = evar(y)

Description

Assuming that the deterministic function y has additive Gaussian noise, evar(y) returns an estimated variance of this noise.

A thin-plate smoothing spline model is used to smooth y. It is assumed that the model whose generalized cross-validation score is minimum can provide the variance of the additive noise. A few tests showed that evar works very well with "not too irregular" functions.

1-D examples

Create a deterministic signal with 1000 data points and add Gaussian noise whose variance ranges from 0 to 0.5

n = 1e3;
x = linspace(0,25,n);
y = round(sin(x));
sig2 = linspace(0,0.5,50);

As an example, this figure shows the effect of an additive noise with a variance of 0.2 (original signal in black, noisy signal in red).

yn = y + sqrt(0.2)*randn(size(x));
plot(x,yn,'r',x,y,'k')

Estimate the variance of the additive noise by using evar

sig2_est = zeros(1,50);
for i = 1:50
    yn = y + sqrt(sig2(i))*randn(1,n);
    sig2_est(i) = evar(yn);
end

Compare the estimated variance (sig2_est) with the true variance (sig2)

plot(sig2_est,sig2,'o',[0 0.5],[0 0.5],'k')
axis([-0.05 0.55 -0.05 0.55])
axis square
xlabel('Estimated variance')
ylabel('True variance')

2-D examples

Create a parametric surface and add Gaussian noise whose variance ranges from 0 to 0.5

[xp,yp] = deal(0:.02:1);
[x,y] = meshgrid(xp,yp);
f = exp(x+y) + sin((x-2*y)*3);
sig2 = linspace(0,0.5,50);

The following figure compares the original parametric surface with the same surface added with a Gaussian noise (variance = 0.1)

fn = f + sqrt(0.1)*randn(size(f));
figure, surf(xp,yp,f), zlim([0 8]), axis square
title('Original surface')
figure, surf(xp,yp,fn), zlim([0 8]), axis square
title('Noisy surface')

Estimate the variance of the additive noise by using evar...

sig2_est = zeros(1,50);
for i = 1:50
    fn = f + sqrt(sig2(i))*randn(size(f));
    sig2_est(i) = evar(fn);
end

... and compare the estimated variance (sig2_est) with the true variance (sig2)

plot(sig2_est,sig2,'o',[0 0.5],[0 0.5],'k')
axis([-0.05 0.55 -0.05 0.55])
axis square
xlabel('Estimated variance')
ylabel('True variance')

Information & Reference

Generalized cross-validation

Buckley MJ, Fast computation of a discretized thin-plate smoothing spline for image data, Biometrika, 81, 2, 1994, pp 247-258.

About the author

Damien Garcia, Eng., Ph.D.
Assistant professor, Department of radiology
CRCHUM, University of Montreal Hospital
Montreal, QC, Canada
Damien.Garcia.REMOVE-THIS@BiomeCardio.com