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_{21}$$

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)

• What did you observe?

Endnotes

• Version: 16 October 2022
• Rmd Source