% EE 438 Digital Signal Processing with Applications 
% Random Signals
% file: random_signals_1.m
% Estimating mean of a sequence of random variables
%---------------------------------------------------------------
% Initialize
clear all;
close all;
%---------------------------------------------------------------
% Parameters
N = 100		% sample size
stddev = 1	% standard deviation of data
mu = 1		% mean of data
%---------------------------------------------------------------
n = 1:1:N;
x = randn(1,N);
y = stddev.*x + mu;

sum = 0;
sumsq = 0;
for m = 1:N
  sum = sum + y(m);
  sumsq = sumsq + y(m)^2;
  muhat(m) = sum/m;
  varhat = sumsq/m - muhat(m)^2;      
  if m == 1
    stddevhat(m) = sqrt(varhat);    % (biased estimate)
  else
    stddevhat(m) = sqrt(m./(m-1).*varhat); % (unbiased estimate)
  end
end

figure(1)
plot(n,y,'-',n,muhat,':',n,muhat+stddevhat,'--r',n,muhat-stddevhat,'--r')
msg=sprintf('\\mu=%.1f, \\sigma^2=%.1f',mu,stddev^2);
title(['Estimating mean of a sequence of random variables ' msg])
xlabel('Sample Number')
ylabel('Data')
legend('Data','Sample Mean','Sample Mean +/- Sample Std. Dev.')
%gtext('Data is i.i.d. Gaussian with mean 1 and std. dev. 1')