MATLAB Examples

# Compare Markov Chain Mixing Times

This example compares the estimated mixing times of several Markov chains having differing structures. Convergence theorems typically require ergodic unichains. Therefore, before comparing mixing time estimates, this example ensures that the Markov chains are ergodic unichains.

## Fast-Mixing Chain

Create a 23-state Markov chain from a random transition matrix containing 250 infeasible transitions of 529 total transitions. Plot a digraph of the Markov chain and identify classes using node color and markers.

```rng(1); % For reproducibility numStates = 23; Zeros1 = 250; mc1 = mcmix(numStates,'Zeros',Zeros1); figure; graphplot(mc1,'ColorNodes',true); ```

mc1 represents a unichain because it is a single recurrent, aperiodic class.

Determine whether the Markov chain is ergodic.

```tf1 = isergodic(mc1) ```
```tf1 = logical 1 ```

tf1 = 1 indicates that mc1 represents an ergodic unichain.

Plot the eigenvalues of the Markov chain on the complex plane.

```figure; eigplot(mc1); ```

The spectral gap is the length of the radius of the disc created by the circle having a radius of length equal to the second largest eigenvalue modulus (SLEM) and the circle of radius one. 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 gap is large indicating a fast-mixing chain.

Estimate the mixing time of the chain.

```[~,tMix1] = asymptotics(mc1) ```
```tMix1 = 0.8357 ```

On average, it takes 0.6121 steps for the total variation distance to decay by a factor of .

## Slow-Mixing Chain

Create another 23-state Markov chain from a random transition matrix containing 475 infeasible transitions. Because there are fewer feasible transitions, this chain should take longer to mix. Plot a digraph of the Markov chain and identify classes using node color and markers.

```Zeros2 = 475; mc2 = mcmix(numStates,'Zeros',Zeros2); figure; graphplot(mc2,'ColorNodes',true); ```

mc2 represents a unichain because it has a single recurrent, aperiodic class and several transient classes.

Determine whether the Markov chain is ergodic.

```tf2 = isergodic(mc2) ```
```tf2 = logical 0 ```

tf2 = 0 indicates that mc2 is not ergodic.

Extract the recurrent subchain from mc2. Determine whether the subchain is ergodic.

```[bins,~,ClassRecurrence] = classify(mc2); recurrentClass = find(ClassRecurrence,1); recurrentState = find((bins == recurrentClass),1); sc2 = subchain(mc2,recurrentState); tf2 = isergodic(sc2) ```
```tf2 = logical 1 ```

sc2 represents an ergodic unichain.

Plot the eigenvalues of the subchain on the complex plane.

```figure; eigplot(sc2); ```

The spectral gap in the subchain is much thinner than the gap in mc1, which indicates that the subchain mixes more slowly.

Estimate the mixing time of the subchain.

```[~,tMix2] = asymptotics(sc2) ```
```tMix2 = 5.0201 ```

On average, it takes 16.4417 steps for the total variation distance to decay by a factor of .

## 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.
```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; mc3 = dtmc(DB); ```

Plot a directed graph of the dumbbell chain. Suppress node labels.

```figure; h = graphplot(mc3,'ColorNodes',true); h.NodeLabel = {}; ```

mc3 represents a unichain because it has a single recurrent, aperiodic.

Determine whether the Markov chain is ergodic.

```tf3 = isergodic(mc3) ```
```tf3 = logical 1 ```

tf3 = 1 indicates that mc3 is ergodic.

Plot the eigenvalues of the dumbbell on the complex plane.

```figure; eigplot(mc3); ```

The spectral gap in the subchain is very thin, which indicates that the dumbbell chain mixes very slowly.

Estimate the mixing time of the dumbbell chain.

```[~,tMix3] = asymptotics(mc3) ```
```tMix3 = 90.4334 ```

On average, it takes 36.3541 steps for the total variation distance to decay by a factor of .