Understanding Ergodicity and Ensemble Averaging

Question

Literature says that a stationary signal is ergodic, if its ensemble average = time averages. Should it be the statistics computed by time averaging = statistics computed by ensemble averaging?The way I understood from the book in the link is as follows :

Case1: There are N samples in a time series I compute the mean1,skewness1,kurtosis1 .

Case2: Next, with the N samples, T experiments are conducted and in each of the experiment is calculate the statistics. After, T experiments are done, we calculate the average of these statistics - mean_avg,skewness_avg, kurtosis_avg. This is called as Ensemble averaging or Probabilistic averaging.

If mean1 = mean_avg and so on, the the process/signal is ergodic. Is my understanding correct?

Answer 1

I will try to explain this practically using MATLAB notation.
Yet before that I must say the ergodic property sometime is limited to a level of moment, namely ergodic in the 1st , 2nd, 3rd moment, etc...
If the process is IID it is promised to be ergodic.

Now, assume we have a function myRandomProcess in MATLAB.
It returns a row vector of length n with the random properties of our random process.

Now, assume we build a matrix A which is $ m $ by $ n $ matrix which each row is a different vector created using the above function.

Now, using the mean and var functions in MATLAB we can compute the properties of the process.

If we apply mean a long the columns (Namely, on each row) we will get $ m $ results.
If those are close enough to the $ n $ results yielded by applying the mean function along the rows (On each column by itself) we can assume the process is ergodic in its first moment.

Now, if $ m $ and $ n $ approaches infinity "close enough" can be small as we want.

Most properties of random process related to its ergodicity can be inferred by looking on its auto correlation function.

Update

Simple MATLAB code to convey the idea:

% Ergodic Random Process

% The Prbability dimension is a long rows (Vertical) and the time dimension
% is a long the columns (Horizontal).

numSamples  = 1000;
numTrials   = 1000;

% Ergodic Process - IID Random Samples

mRandomProcess = randn(numTrials, numSamples);

vTimeAverage = mean(mRandomProcess, 2);
vSampleAverage = mean(mRandomProcess, 1);

vTimeVariance = var(mRandomProcess, 1, 2);
vSampleVariance = var(mRandomProcess, 1, 1);

disp(['Ergodic Process Statistics']);
disp(['1st Moment Ergodicity - ', num2str(mean([vTimeAverage(:) - vSampleAverage(:)] .^ 2))]);
disp(['2nd Moment Ergodicity - ', num2str(mean([vTimeVariance(:) - vSampleVariance(:)] .^ 2))]);

% Un Ergodic Process

% The DC Level is drawed once per random process (Realization)
vDcLevel = (10 * rand(numTrials, 1)) - 5;

mRandomProcess = randn(numTrials, numSamples) + repmat(vDcLevel, [1, numSamples]);

vTimeAverage = mean(mRandomProcess, 2);
vSampleAverage = mean(mRandomProcess, 1);

vTimeVariance = var(mRandomProcess, 1, 2);
vSampleVariance = var(mRandomProcess, 1, 1);

disp(['Un Ergodic Process Statistics']);
disp(['1st Moment Ergodicity - ', num2str(mean([vTimeAverage(:) - vSampleAverage(:)] .^ 2))]);
disp(['2nd Moment Ergodicity - ', num2str(mean([vTimeVariance(:) - vSampleVariance(:)] .^ 2))]);

As you'll be able to see, the first process is Ergodic (As any IID process), as it is equivalent to calculate its moment along the time dimension or the sample dimension.

On the other hand, the second one isn't Ergodic.
Since its mean along the rows is zero (Drawn from a uniformly random variable [-5, 5]) yet along its columns its mean is the value of the specific realization of the uniform variable.

I hope this clarifies it.