In class exercise: Two-by-two and symbolic computing

Symbolic computing inside Python (inside R)

Will we explore basic symbolic computing inside Python. In symbolic computing, the software is able to compute \(2x = x+x\) or \(\int_0^x x^2 dx = \tfrac{x^3}{3}\) without having a numerical value for \(x\). Popular software includes, Maple, Mathematica, and Sage; Python’s sympy is a basic in-house substitute, which we will explore.

import sympy as sym

x = sym.pi   # the digit pi
print(x)

## pi

print("#############")

## #############

print(x.evalf())

## 3.14159265358979

print("#############")

## #############

print(x.evalf(100))

## 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

print("#############")

## #############

x = sym.Symbol('x')   # you got to tell Python that x is symbol
print(sym.simplify((x + x)))

## 2*x

y = sym.Symbol('y')
y=sym.integrate(3 * x ** 2, x)  # integrate 3x^2
print(y)

## x**3

print("#############")

## #############

We can also do linear algebra

a = sym.Symbol('a')
b = sym.Symbol('b')
c = sym.Symbol('c')
d = sym.Symbol('d')
M =  sym.Symbol('M')
M=sym.Matrix([[a, b], [c, d]])
print(M)

## Matrix([[a, b], [c, d]])

print("#############")

## #############

print(M.det())

## a*d - b*c

print("#############")

## #############

print(M.T)

## Matrix([[a, c], [b, d]])

print("#############")

## #############

print(M.eigenvects())

## [(a/2 + d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2, 1, [Matrix([
## [-d/c + (a/2 + d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2)/c],
## [                                                         1]])]), (a/2 + d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2, 1, [Matrix([
## [-d/c + (a/2 + d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2)/c],
## [                                                         1]])])]

print("#############")

## #############

Msub = M.subs(a,1).subs(b,1).subs(c,1).subs(d,0)
values = Msub.eigenvects()
print("#############")

## #############

print(values)

## [(1/2 - sqrt(5)/2, 1, [Matrix([
## [1/2 - sqrt(5)/2],
## [              1]])]), (1/2 + sqrt(5)/2, 1, [Matrix([
## [1/2 + sqrt(5)/2],
## [              1]])])]

print("#############")

## #############

A Markov chain

Consider the two state Markov chain on \(\{1,2\}\) where \[p_{11} = a = 1-p_{12}\] and

\[p_{21} = b=1- p_{22}\]

We will consider the case where \(a,b \in (0,1)\).

• Write down the transition matrix.

• Find the stationary distribution of the transition matrix.

• Simulate the chain, (starting at \(1\)) and leave \(a\) and \(b\) as variables that you can choose.

• Choose \(a=0.3\) and \(b=0.2\).

• At \(t=100\) of your simulation, for this snapshot, record the value of the Markov chain; repeat this for \(N=300\) times, and find the proportion times that the chain is at \(1\) and \(2\).

• For a given simulation, of length \(t=1000\), also record the total number of times the Markov chain takes the value \(1\) and \(2\); take these numbers and divide it by \(1000\) (or a thousand and one depending how you coded)

  • Also record the number of times the sequence \((1,1,2)\) appears in the chain, and takes this number and divide by \(1000\); what do you expect to see?

• What did you observe?

Endnotes

• Version: 26 October 2023
• Rmd Source