%% Determine Dumbbell Chain Mixing Time
% Create a "dumbbell" Markov chain containing 10 states in each "weight" and three states in the "bar."
% Specify random transition probabilities between states within each weight.
% If the Markov chain reaches the state in a weight that is closest to the bar, then specify a high probability of transitioning to the bar.
% Specify uniform transitions between states in the bar.

rng(1); % For reproducibility
w = 10;                              % Dumbbell weights
DBar = [0 1 0; 1 0 1; 0 1 0];        % Dumbbell bar
DB = blkdiag(rand(w),DBar,rand(w));  % Transition matrix

% Connect dumbbell weights and bar
DB(w,w+1) = 1;                       
DB(w+1,w) = 1; 
DB(w+3,w+4) = 1; 
DB(w+4,w+3) = 1;

mc = dtmc(DB);

%
% Visualize the transition matrix using a heatmap.

figure;
imagesc(mc.P);
colormap(jet);
axis square;
colorbar;

%
% Plot a directed graph of the Markov chain. Suppress node labels.

figure;
h = graphplot(mc);
h.NodeLabel = {};

%
% Plot the eigenvalues of the dumbbell chain.

figure;
eigplot(mc);

% The thin, red disc in the plot shows the spectral gap (the difference between the two largest eigenvalue moduli). The spectral gap determines the mixing time of the Markov chain. Large gaps indicate faster mixing, whereas thin gaps indicate slower mixing. In this case, the spectral gap is thin, indicating a long mixing time.
%
% Estimate the mixing time of the dumbbell chain and determine whether the chain is ergodic.

[~,tMix] = asymptotics(mc)
tf = isergodic(mc)

% On average, the time it takes for the total variation distance between any initial distribution and the stationary distribution to decay by a factor of e^1 is about 85 steps.
% Copyright 2012 The MathWorks, Inc.