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())  # eigenvalue, multiplicity, eigenvector
## [(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)\).

Solutions

import numpy as np

a = sym.Symbol('a')
b = sym.Symbol('b')
P =  sym.Symbol('P')
P=sym.Matrix([[a, 1-a], [b,1-b]])
Pt =  sym.Symbol('Pt')
Pt= P.T
print("#############")
## #############
print(Pt.eigenvects())
## [(1, 1, [Matrix([
## [-b/(a - 1)],
## [         1]])]), (a - b, 1, [Matrix([
## [-1],
## [ 1]])])]
print("#############")
## #############
Q=P.subs(a,0.1).subs(b,0.4)  ## testing code for subs (will not be using these values)
print(Q)
## Matrix([[0.100000000000000, 0.900000000000000], [0.400000000000000, 0.600000000000000]])
print(Q[0])
## 0.100000000000000
print(Q[1])
## 0.900000000000000
print(Q[2])
## 0.400000000000000
print(Q[3])
## 0.600000000000000
v1= np.array(  [Q[0],Q[1]]  )
v2= np.array(  [Q[2],Q[3]]  )  
Q = np.array([v1,v2])
print(Q)
## [[0.100000000000000 0.900000000000000]
##  [0.400000000000000 0.600000000000000]]

So we read off that the eigenvector we are after is \([-b/(a-1), 1]\), and after normalization, we have \[\frac{1}{1-a+b} *[b, 1-a].\]

Note, we would get the nicer looker expression, if we had chosen: \(a \to 1-a\) at the outset.


def step(i,A,B):    # advancing the MC by one step, when you are at state i, with values for a and b
    Q = P.subs(a,A).subs(b,B)
    v1= np.array([Q[0],Q[1]])
    v2=np.array([Q[2],Q[3]])
    Q = np.array([v1,v2])
    q = Q[i-1]    # this the transition vector for state i
    x=-1
    u = np.random.uniform()   # imagine the interval [0,1] split up into smaller intervals with probabilities q_j summing to 1
    j=0
    cumq = np.cumsum(q)
    while(x==-1):            
        j = j+1
        if (u <= cumq[j-1]):        # if u lands in the interval of length q_j, then we jump to state j
            x = j
    return x

def steps(i,n,A,B):  # advances the MC by n steps, starting at state i, (with values A and B) keeping a complete history
    x = np.array([i])
    for j in range(n):
        x = np.append(x,  step(x[j],A,B)   )
    return x

print(steps(1,10, 0.3,0.2))
## [1 2 2 2 2 2 2 2 2 1 2]
x= [steps(1,100,0.3,0.2)[100] for _ in range(300)]
freq = np.unique(x,return_counts = True)
print(freq)
## (array([1, 2]), array([ 68, 232], dtype=int64))
stat = (1/300)* np.array( [freq[1][0], freq[1][1] ])
thstat = (1/(1-0.3+0.2))*np.array([0.2, 1-0.3])
print("#############")
## #############
print(stat)
## [0.22666667 0.77333333]
print(thstat)
## [0.22222222 0.77777778]
print("#############")
## #############

In the sample path average, we have


path = steps(1, 1000, 0.3,0.2)
freq = np.unique(path,return_counts = True)
print(freq)
## (array([1, 2]), array([185, 816], dtype=int64))
stat = (1/1001)* np.array( [freq[1][0], freq[1][1] ])
print(stat)
## [0.18481518 0.81518482]

In both cases, we recover the stationary distribution.

And in a final example, we count the occurrences of \((1,1,2)\) occurs


path = steps(1, 1000, 0.3,0.2)
oneonetwo = 0
for i in range(998):
  if (path[i]==1 and path[i+1]==1 and path[i+2]==2):
    oneonetwo = oneonetwo +1

print(oneonetwo/1000)
## 0.045
tans= (thstat[0])*(0.3) * (0.7)
print(tans)
## 0.04666666666666666

Endnotes