import numpy M = numpy.array([[0, 1/2, 1, 0], [1/3, 0, 0, 1/2], [1/3, 0, 0, 1/2], [1/3, 1/2, 0, 0]]) beta = 0.8 es1 = numpy.array([1, 0, 0, 0]) r0 = numpy.array([1/4, 1/4, 1/4, 1/4]) r1 = beta * M * r0 + (1-beta) * es1 print(r1)