CODEWITHSUNDEEP
pythonnumpymatrix
JAA
JAA

Reputation: 19

multiply matrices on python 3

anyone know why I got the wrong result? I know the there is a mistake in the last line in the code which is f=p0*pn*p

import numpy as np

def passage(n,i,j):
    # this function calculate the first time passage distribution after n 
    #steps starting at i and end at j.

    p=np.matrix([[0,1,0,0],[0.5,0,0.5,0],[0,.5,0,.5],[0,0,1,0]])

    p0=p[:]
    for k in range(len(p)): # let elements in column j be zeros
        p0[k,j]=0

    p1=p0[:]
    for k in range(len(p)): #let element in column j and row j be zeros 
        p1[j,k]=0

    pn=np.linalg.matrix_power(p1,n-2)


    f=p0*pn*p      # this line gives us wrong result. why


    return f

Upvotes: 0

Views: 84

Answers (2)

hpaulj
hpaulj

Reputation: 231355

One potential source of problems:

p0=p[:]

produce a view, not a copy. That means that after those 2 loops (which probably can be written without looping), p, p0 and p1 have the same values.

Did you actually test this code line by line, making sure each step was correct? When I write functions in Python, and especially numpy have test out all steps interactively.

Another thing - unless you really need it, don't use np.matrix. Stick with np.array when creating arrays, even 2d ones. And use np.dot (or @) if you need matrix multiplication.

A quick rewrite of your function:

def passage(n,i,j):
    p=np.array([[0,1,0,0],[0.5,0,0.5,0],[0,.5,0,.5],[0,0,1,0]])
    p0=p.copy()
    p0[:,j] = 0
    p1=p0.copy()
    p1[j,:] = 0
    pn=np.linalg.matrix_power(p1,n-2)
    f = p0@pn@p    # or p0.dot(pn.dot(p))
    return f

In [15]: passage(4,0,0)
Out[15]: 
array([[ 0.125 ,  0.    ,  0.375 ,  0.    ],
       [ 0.    ,  0.1875,  0.    ,  0.1875],
       [ 0.1875,  0.    ,  0.5625,  0.    ],
       [ 0.    ,  0.375 ,  0.    ,  0.375 ]])

Upvotes: 1

tommy.carstensen
tommy.carstensen

Reputation: 9622

The solution is to use np.copy() to copy your p matrix.

p0 = np.copy(p)
p1 = np.copy(p0)

If you still don't get the expected answer, then you might want to have a look at np.dot() and np.matmul().

Also note that PEP465 introduced @ as a matrix multiplication operator.

Upvotes: 0

Related Questions