Markov Chains in Maple

The following is a procedure for Maple that simulates a Markov chain by executing a specified number of trials:

MarkovSim := proc(Mat, startstage, ntrials)
  local count, nvisit, a, i, j, firstAbsorbingStage, run, stage, sumQuant;
  count:=Vector(ColumnDimension(Mat));
  nvisit:=Vector(ColumnDimension(Mat));
  firstAbsorbingStage:=0;
  for i from 1 to ColumnDimension(Mat) do
    for j from 1 to RowDimension(Mat) do
      if Mat[j,i]>=1.0 then firstAbsorbingStage:=i; break; fi;
    od;
    if firstAbsorbingStage>0 then break; fi;
  od;
  for i from 1 to ColumnDimension(Mat) do count[i]:=0; od;
  for run from 1 to ntrials do
    for i from 1 to ColumnDimension(Mat) do nvisit[i]:=0; od;
    nvisit[startstage]:=1;
    stage:=startstage;
    while stage
      a:=rand()/10.^12;
      sumQuant:=0;
        for i from 1 to RowDimension(Mat) do
        if a<(sumQuant+Mat[i,stage]) and a>sumQuant then
          stage:=i; nvisit[i]:=nvisit[i]+1; break;
        fi;
        sumQuant:=sumQuant+Mat[i,stage];
      od;
    od;
    count:=count+nvisit;
  od;
  printf("Average number of steps in different stages before absorption:\n");
  for i from 1 to firstAbsorbingStage-1 do printf("%g ", count[i]/ntrials); od;
  j:=0; printf("\n");
  for i from 1 to ColumnDimension(Mat) do j:=j+count[i]; od;
  printf("Average number of steps until absorption:\n %g", j/ntrials);
end proc: