where you're checking that the transition matrix sums to one across the rows, I think you can have problems due to finite precision. The sum can be epsilon different than unity. A reasonable check would be that abs(1-sum) < n*eps where eps is the floating-point precision (say, 4e-16).
In fact, a strict check like this can sometimes be annoying because you might save or transmit the values in a JSON format with less than full precision.
Since the row has "n" entries, and errors can accumulate in each add, something like (n+2)*eps may be more appropriate. (I wrote a similar library in Matlab.)
That was a fast delete! I was trying to reply that the form in your code is detailed in the excellent book Simulation Modeling and Analysis by Law (previous editions by Law and Kelton. I'm sure there's a story there...)
In fact, a strict check like this can sometimes be annoying because you might save or transmit the values in a JSON format with less than full precision.