## DESCRIPTION
## Determine which states are accessible from other states in a Discrete Time Markov Chain with a finite state space.
## ENDDESCRIPTION

## KEYWORDS('transition matrix')

## DBsubject('Probability')
## DBchapter('Stochastic process')
## DBsection('Markov chain')
## Date('20250330')
## Author('Joe Pedersen')
## Institution('USMA')

DOCUMENT();

loadMacros('PGstandard.pl', 'PGML.pl',"MathObjects.pl",
"parserFunction.pl", 'PGtikz.pl', 'niceTables.pl', 'PGcourse.pl');

Context('Numeric')->flags->set(
    tolerance => 0.001,
    tolType   => 'absolute',
);

$showPartialCorrectAnswers = 1;

@rr = map { $_ / 100} random_subset(13, 20..80);

@am = (
    [  0, $rr[0],   0, 1-$rr[0],   0,   0,   0,   0,   0,   0,   0,   0,   0], 
    [  0,   0, $rr[1],   0, 1-$rr[1],   0,   0,   0,   0,   0,   0,   0,   0], 
    [$rr[2],   0,   0,   0,   0, 1-$rr[2],   0,   0,   0,   0,   0,   0,   0], 
    [  0,   0, $rr[3],   0, 1-$rr[3],   0,   0,   0,   0,   0,   0,   0,   0], 
    [$rr[4],   0,   0,   0,   0, 1-$rr[4],   0,   0,   0,   0,   0,   0,   0], 
    [  0, $rr[5],   0, 1-$rr[5],   0,   0,   0,   0,   0,   0,   0,   0,   0], 
    [  0,   0,   0,  0,  0, $rr[6]/2,  0, $rr[6]/2, 1-$rr[6],  0,  0,  0,  0], 
    [  0,   0,   0,   0,   0,   0, $rr[7],   0, 1-$rr[7],   0,   0,   0,   0], 
    [0, 0, $rr[8]/4, (1-$rr[8])/4, 0, 0, $rr[9]/2, 0, (1-$rr[9])/2, 0, 0, 0, 1/4], 
    [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   1,   0,   0], 
    [  0,   0,   0,   0,   0,   0,   0,   0,   0, $rr[10],   0, 1-$rr[10],   0], 
    [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   1,   0], 
    [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   1,   0], 
);

@am2 = (
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
);

@i13 = (
    [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
);

@this_order = random_subset(13, (0..12));
@rev_order = ();
for $i (0..12) {
    $rev_order[$this_order[$i]] = $i;
}

for $i (0..12) {
    for $j (0..12) {
        $am2[$i][$j] = $am[$this_order[$i]][$this_order[$j]];
    }
}

Context("Matrix");
$M = Matrix(@am2);
$MN = Matrix(@am2);
$MT = Matrix(@i13);


for $i (0..12) {
    $MT = $MT + $MN; # Expected number of visits from i to j in N steps
    $MN = $MN * $M;  # Powers of the transition matrix
}

# answer order
@ao = random_subset(4, ($rev_order[3], $rev_order[6], $rev_order[9], $rev_order[12]));

@from0 = ();
for $j (0..12) {
    if ( $MT->element($ao[0]+1,$j+1) > 0 ) {
        push(@from0, $j);
    }
}
if ((scalar @from0) > 1) {
    $ans0 = Set(@from0);
} else {
    $ans0 = List(@from0);
}

@from1 = ();
for $j (0..12) {
    if ( $MT->element($ao[1]+1,$j+1) > 0 ) {
        push(@from1, $j);
    }
}
if ((scalar @from1) > 1) {
    $ans1 = Set(@from1);
} else {
    $ans1 = List(@from1);
}

@from2 = ();
for $j (0..12) {
    if ( $MT->element($ao[2]+1,$j+1) > 0 ) {
        push(@from2, $j);
    }
}
if ((scalar @from2) > 1) {
    $ans2 = Set(@from2);
} else {
    $ans2 = List(@from2);
}

@from3 = ();
for $j (0..12) {
    if ( $MT->element($ao[3]+1,$j+1) > 0 ) {
        push(@from3, $j);
    }
}
if ((scalar @from3) > 1) {
    $ans3 = Set(@from3);
} else {
    $ans3 = List(@from3);
}

BEGIN_PGML
Given the following probability transition matrix for a Markov chain:

Note: assume rows/columns are zero-indexed, with the top left entry containing:  [`p_{0,0} = P(X_{i+1} = 0|X_i = 0)`]

[#
  [. [$am2[0][0]] .]
  [. [$am2[0][1]] .]
  [. [$am2[0][2]] .]
  [. [$am2[0][3]] .]
  [. [$am2[0][4]] .]
  [. [$am2[0][5]] .]
  [. [$am2[0][6]] .]
  [. [$am2[0][7]] .]
  [. [$am2[0][8]] .]
  [. [$am2[0][9]] .]
  [. [$am2[0][10]] .]
  [. [$am2[0][11]] .]
  [. [$am2[0][12]] .]*

  [. [$am2[1][0]] .]
  [. [$am2[1][1]] .]
  [. [$am2[1][2]] .]
  [. [$am2[1][3]] .]
  [. [$am2[1][4]] .]
  [. [$am2[1][5]] .]
  [. [$am2[1][6]] .]
  [. [$am2[1][7]] .]
  [. [$am2[1][8]] .]
  [. [$am2[1][9]] .]
  [. [$am2[1][10]] .]
  [. [$am2[1][11]] .]
  [. [$am2[1][12]] .]*

  [. [$am2[2][0]] .]
  [. [$am2[2][1]] .]
  [. [$am2[2][2]] .]
  [. [$am2[2][3]] .]
  [. [$am2[2][4]] .]
  [. [$am2[2][5]] .]
  [. [$am2[2][6]] .]
  [. [$am2[2][7]] .]
  [. [$am2[2][8]] .]
  [. [$am2[2][9]] .]
  [. [$am2[2][10]] .]
  [. [$am2[2][11]] .]
  [. [$am2[2][12]] .]*

  [. [$am2[3][0]] .]
  [. [$am2[3][1]] .]
  [. [$am2[3][2]] .]
  [. [$am2[3][3]] .]
  [. [$am2[3][4]] .]
  [. [$am2[3][5]] .]
  [. [$am2[3][6]] .]
  [. [$am2[3][7]] .]
  [. [$am2[3][8]] .]
  [. [$am2[3][9]] .]
  [. [$am2[3][10]] .]
  [. [$am2[3][11]] .]
  [. [$am2[3][12]] .]*

  [. [$am2[4][0]] .]
  [. [$am2[4][1]] .]
  [. [$am2[4][2]] .]
  [. [$am2[4][3]] .]
  [. [$am2[4][4]] .]
  [. [$am2[4][5]] .]
  [. [$am2[4][6]] .]
  [. [$am2[4][7]] .]
  [. [$am2[4][8]] .]
  [. [$am2[4][9]] .]
  [. [$am2[4][10]] .]
  [. [$am2[4][11]] .]
  [. [$am2[4][12]] .]*

  [. [$am2[5][0]] .]
  [. [$am2[5][1]] .]
  [. [$am2[5][2]] .]
  [. [$am2[5][3]] .]
  [. [$am2[5][4]] .]
  [. [$am2[5][5]] .]
  [. [$am2[5][6]] .]
  [. [$am2[5][7]] .]
  [. [$am2[5][8]] .]
  [. [$am2[5][9]] .]
  [. [$am2[5][10]] .]
  [. [$am2[5][11]] .]
  [. [$am2[5][12]] .]*

  [. [$am2[6][0]] .]
  [. [$am2[6][1]] .]
  [. [$am2[6][2]] .]
  [. [$am2[6][3]] .]
  [. [$am2[6][4]] .]
  [. [$am2[6][5]] .]
  [. [$am2[6][6]] .]
  [. [$am2[6][7]] .]
  [. [$am2[6][8]] .]
  [. [$am2[6][9]] .]
  [. [$am2[6][10]] .]
  [. [$am2[6][11]] .]
  [. [$am2[6][12]] .]*

  [. [$am2[7][0]] .]
  [. [$am2[7][1]] .]
  [. [$am2[7][2]] .]
  [. [$am2[7][3]] .]
  [. [$am2[7][4]] .]
  [. [$am2[7][5]] .]
  [. [$am2[7][6]] .]
  [. [$am2[7][7]] .]
  [. [$am2[7][8]] .]
  [. [$am2[7][9]] .]
  [. [$am2[7][10]] .]
  [. [$am2[7][11]] .]
  [. [$am2[7][12]] .]*

  [. [$am2[8][0]] .]
  [. [$am2[8][1]] .]
  [. [$am2[8][2]] .]
  [. [$am2[8][3]] .]
  [. [$am2[8][4]] .]
  [. [$am2[8][5]] .]
  [. [$am2[8][6]] .]
  [. [$am2[8][7]] .]
  [. [$am2[8][8]] .]
  [. [$am2[8][9]] .]
  [. [$am2[8][10]] .]
  [. [$am2[8][11]] .]
  [. [$am2[8][12]] .]*

  [. [$am2[9][0]] .]
  [. [$am2[9][1]] .]
  [. [$am2[9][2]] .]
  [. [$am2[9][3]] .]
  [. [$am2[9][4]] .]
  [. [$am2[9][5]] .]
  [. [$am2[9][6]] .]
  [. [$am2[9][7]] .]
  [. [$am2[9][8]] .]
  [. [$am2[9][9]] .]
  [. [$am2[9][10]] .]
  [. [$am2[9][11]] .]
  [. [$am2[9][12]] .]*

  [. [$am2[10][0]] .]
  [. [$am2[10][1]] .]
  [. [$am2[10][2]] .]
  [. [$am2[10][3]] .]
  [. [$am2[10][4]] .]
  [. [$am2[10][5]] .]
  [. [$am2[10][6]] .]
  [. [$am2[10][7]] .]
  [. [$am2[10][8]] .]
  [. [$am2[10][9]] .]
  [. [$am2[10][10]] .]
  [. [$am2[10][11]] .]
  [. [$am2[10][12]] .]*

  [. [$am2[11][0]] .]
  [. [$am2[11][1]] .]
  [. [$am2[11][2]] .]
  [. [$am2[11][3]] .]
  [. [$am2[11][4]] .]
  [. [$am2[11][5]] .]
  [. [$am2[11][6]] .]
  [. [$am2[11][7]] .]
  [. [$am2[11][8]] .]
  [. [$am2[11][9]] .]
  [. [$am2[11][10]] .]
  [. [$am2[11][11]] .]
  [. [$am2[11][12]] .]*

  [. [$am2[12][0]] .]
  [. [$am2[12][1]] .]
  [. [$am2[12][2]] .]
  [. [$am2[12][3]] .]
  [. [$am2[12][4]] .]
  [. [$am2[12][5]] .]
  [. [$am2[12][6]] .]
  [. [$am2[12][7]] .]
  [. [$am2[12][8]] .]
  [. [$am2[12][9]] .]
  [. [$am2[12][10]] .]
  [. [$am2[12][11]] .]
  [. [$am2[12][12]] .]*
#]{align => '|c|c|c|c|c|c|c|c|c|c|c|c|c|', horizontalrules => 1}

Determine the quantities below. 

TIP:  do *NOT* type the numbers by hand, because you'll make an error somewhere and get the answers wrong. Instead, you can copy and paste the matrix above  
1) into MS Excel, save it as a CSV, and read that with pandas read_csv, or
2) into python as a triple-quoted string, and read that with pandas read_table


What is the set of states that are accessible from each of the following states? (Remember that the states are numbered 0 to 12)

Surround your answers with curly brackets, as in:  [` \{0,1,4,9\} `]

The states accessible from [`[$ao[0]] `] are: [_]{$ans0->cmp(requireParenMatch => 0, implicitList => 0)}  
The states accessible from [`[$ao[1]] `] are: [_]{$ans1->cmp(requireParenMatch => 0, implicitList => 0)}  
The states accessible from [`[$ao[2]] `] are: [_]{$ans2->cmp(requireParenMatch => 0, implicitList => 0)}  
The states accessible from [`[$ao[3]] `] are: [_]{$ans3->cmp(requireParenMatch => 0, implicitList => 0)}  
END_PGML

ENDDOCUMENT();