Example: Estimation of a Gaussian Vector
From Gampmatlab
Main > Users' Guide and Examples > Estimation of a Gaussian Vector.
Next: Estimating a sparse vector.
Contents 
Program
The code for the example is in gampmatlab/trunk/code/examples/basic/simpleAWGN.m
.
To run the code, simply cd
to that directory and run simpleAWGN
at the MATLAB prompty. You will see a plot of a sparse vector and its estimate. In addition, the mean squared error (MSE) for both the leastsquares and GAMP methods will be printed.
Problem Description
In this first problem, we will show how to use the gampEst
function to estimate a Gaussian random vector from random linear measurements of the form
where is a Gaussian random matrix and is Gaussian noise. The solution to this problem is, of course, standard and be performed easily without GAMP: For example, if the components of and are i.i.d., Gaussian zeromean with variances and , respectively, then the MMSE estimate is given by the least squares solution:
Nevertheless, we will show how to get an approximate solution to this problem using GAMP
simply to illustrate the syntax and usage of the MATLAB command.
MATLAB Program Description
Generation of the Data
At the top of the program simpleAWGN.m
, you will see the lines:
% Set path to the main directory addpath('../../main/'); % Parameters nx = 1000; % Number of input components (dimension of x) nz = 2000; % Number of output components (dimension of y) snr = 20; % SNR in dB.
The first line adds the necessary directories to the MATLAB path. You can also set the MATLAB path permanently if you intend to use GAMP significantly. The remaining lines set the parameters of the problem. Feel free to change these to experiment.
The next lines generate a random Gaussian vector :
% Create a random Gaussian vector xmean0 = 0; xvar0 = 1; x = normrnd(xmean0, sqrt(xvar0),nx,1);
The matrix and the noise vector are generated with Gaussian
i.i.d. components, with the noise scaled according to the parameter snr
:
% Create a random measurement matrix A = (1/sqrt(nx))*randn(nz,nx); % Compute the noise level based on the specified SNR. Since the components % of A have power 1/nx, the Ey(i)^2 = Ex(j)^2 = xvar. wvar = 10^(0.1*snr)*xvar0; % Generate the noise w = normrnd(0, sqrt(wvar), nz, 1); y = A*x + w;
Solution via GAMP
The main GAMP command to perform the estimation of from the measurement vector appears near the bottom of the file:
[xhat, xvar] = gampEst(inputEst, outputEst, A, opt);
The gampEst
function takes as input arguments:

inputEst
is a class of typeEstimIn
that describes the prior on the components of the unknown vector ; 
A
is a matrix for the transform for the relation . This can also be specified by aLinTrans
class as we will see in subsequent examples. 
outputEst
is class of typeEstimOut
describing the transition probabilities for the measurements ; and 
opt
is a set of options.
The outputs xhat
and xvar
are the estimates of the mean and variance
of the vector . You can, in fact, get other outputs to estimate the full marginal posteriors
of each the variables  we will illustrate this in a later example.
The idea of specifying the prior and measurement channel
via the MATLAB classes inputEst
and
outputEst
may be confusing at first.
However, as we will see in this sequence of examples, the use of classes
allows the algorithm to be very general, modular and capable of handling many different
types of distributions. If you are not familiar with MATLAB classes or objectoriented programming, you can read more on the
MATLAB objectoriented webpage.
However, the gampmatlab
package contains prebuilt
classes for many common input and output distributions.
For most problems, you can just use one of these
prebuilt classes and not have to build one from scratch.
In this particular example, since the input is a Gaussian random vector,
the class inputEst
is created using the prebuilt AwgnEstimIn
class with the line:
inputEst = AwgnEstimIn(xmean0, xvar0);
The output estimator class is also created via a prebuilt class. Specifically, the line
outputEst = AwgnEstimOut(y, wvar);
creates outputEst
class using the function AwgnEstimOut
for an AWGN measurement
channel. Note that the observation vector y
is passed to the estimator class  not the
gampEst
function.
With the estimator classes defined, the program then calls gampEst
which returns
an estimate vector xhat
.
The GAMP estimate is compared to the MMSE estimator. For this problem, the MMSE estimate is given by the solution in equation (1) above to the leastsquares problem. This estimate is computed in the MATLAB code with the lines:
xhatLS = xvar0*(A'*A*xvar0 + wvar*eye(nx))\(A'*y);
The program concludes by plotting x
, the GAMP estimate xhat
and LS estimate
xhatLS
. The MSE for the GAMP and LS estimates are also printed. You should see that GAMP estimate obtains an error close to the LS estimate.