PSTAT 160A PYTHON HOMEWORK
PYTHON HOMEWORK 6 { DUE 5/22 AT 5:00 PM
PSTAT 160A { S19
Professor Hohn
Instructions: Please note that you must work by yourself ! You will submit two les on Gau-
choSpace: (1) clear and concise explanations, graphics (if any), and results in PDF format (worth
40 points in total) and (2) your Python code in .py format (10 points). Note that if the grader
nds identical copies or very similar les, the grader cannot and will not grade them.
For Question #1, be sure to show all of your work! You may either (a) write at least a couple
sentences explaining your reasoning or (b) annotate your math work with brief explanations. Please
label any random variables or events that your use.
1. (10 points) (Without Python) Let (Xm)1 m=0 be a stationary discrete time Markov chain with
state space S = f1; 2; 3; 4g and transition matrix
P =
0
BB@
1=3 1=2 1=6 0
1=2 1=8 1=4 1=8
1=4 1=2 0 1=4
0 0 1=2 1=2
1
CCA (a) Brie
y justify why an invariant distribution exists.
(b) Calculate the invariant distribution.
2. (30 points) (With Python) In this Python exercise, you will compute the invariant distribution,
expected return times, and the long run frequency of states in a SDTMC. Be sure to annotate
your code with short explanations of what you are doing (worth 10 points).
Suppose (as above) that we have a stationary discrete time Markov chain with transition matrix
P.
(a) Find the invariant distribution of P by using Python to raise P to the power 100. Use
the print function to print your answer. Be sure to label your results. For example,
print('The invariant distribution that can be found by raising P to the
100th power is (%.3f,%.3f,%.3f,%.3f).' %tuple(inv_dist_via_matrix_power))
(b) For each i 2 S, given that the process starts in state i, estimate the expected return time
to i. That is, nd Ei[i] for each state i. Use the print function to print your answer. For
example, for each starting state,
print("The mean return time from state %s is %.3f"
% (starting_state, np.mean(rho)) )
As in previous assignments, start by simulating 10,000 sample paths of the Markov chain
(Xm)1 m=0 starting from state i until returning to state i.
(c) Use the information found in part (b) to calculate the invariant distribution for P. Use
the print function to print your answer.
(d) Simulate a sample path with 100,000 steps of the Markov chain (Xm)1m
=0 starting from
state i 2 S and calculate the long term frequency of hitting each state i 2 S. Use the
print function to print your answer.
(e) Take a screenshot showing your code and your results together (e.g. side by side).
Python Code Hints
numpy.linalg.matrix_power(P,3)
will give you P3 for a matrix P.
numpy.linalg.solve(A,b)
will solve the matrix equation Ax = b.
numpy.identity(n)
will give you the n n identity matrix.
numpy.matrix
will return a matrix from an array-like object or from a string of data. For example,
a = numpy.matrix([[1, 2], [3, 4]])
will give you the matrix
1 2
3 4
. It has certain special operators, such as * (matrix
multiplication) and ** (matrix power).
Change a matrix to a list via
nameOfMyMatrix.tolist()[0]
One way to make a bar graph is to use the library matplotlib.
matplotlib.pyplot.bar( x, height = y, align='center', alpha=0.5, color='g')
produces a bar graph with x-axis described by x and y axis described by y. Both x and y
are lists here.
matplotlib.pyplot.ylabel('Label Me')
creates a label for the y axis called Label Me.
matplotlib.pyplot.title('Title Me')
creates a title for the graph called Title Me.
matplotlib.pyplot.xticks(x, listOfNames)
will label the x axis tick marks with a list of names (e.g. Port numbers).
matplotlib.pyplot.show()
shows the graph/plot.
numpy.random.choice( alist, p = alist_prob )
will pick one element from the set alist using the probability distribution alist prob.
When Python indexes a list, the index of the list starts at 0. That is, to access the rst
entry of a list like a=[7, 8, 9], we need to write a[0].
alist = [1, 4, 7]
alist.append( 3 )
will add 3 to the end of your list. So, alist = [1, 4, 7, 3].
Page 2
set()
creates a set { an unordered collections of unique elements.
aSet = set( ['Stark', 'Lannister', 'Greyjoy'] )
aSet.add( 'Baratheon' )
will add Baratheon to the set aSet.
alist = [1, 4, 7, 1, 4, 4]
alist.count( 4 )
returns 3, the number of times 4 occurs in the list alist.
len (alist)
return the length (the number of items) of alist. alist can be a sequence (such as a list)
or a collection (like a set).
You may need to use a for loop or while statement in you code. See
http://www.openbookproject.net/books/bpp4awd/ch04.html for examples.
If you are using Python 2, you'll need to import division from Python 3 so that it acts
like Python 3 when using the / symbol in your calculations. If you are using Python 2,
write at the top of your py le,
from __future__ import division
You're importing the future!