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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 |
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: |