ABSTRACT (original article): This paper investigates stochastic finite matrices and the corresponding finite Markov chains constructed using recurrence matrices for general families of orthogonal polynomials and multiple orthogonal polynomials. The paper explores the spectral theory of transition matrices, utilizing both orthogonal and multiple orthogonal polynomials.Several properties are derived, including classes, periodicity, recurrence, stationary states, ergodicity, expected recurrence times, time-reversed chains and reversibility. Furthermore, the paper uncovers factorization in terms of pure birth and pure death processes. The case study focuses on hypergeometric orthogonal polynomials, where all the computations can be carried out effectively. Particularly within the Askey scheme, all descendants under Hahn (excluding Bessel), such as Hahn, Jacobi, Meixner, Kravchuk, Laguerre,Charlier and Hermite, present interesting examples of recurrent reversible birth and death finite Markov chains. Additionally, the paper considers multiple orthogonal polynomials including multiple Hahn, Jacobi-Piñeiro, Laguerre of the first kind and Meixner of the second kind, along with their hypergeometric representations and derives the corresponding recurrent finite Markov chains and time-reversed chains. CITATION (original article): Branquinho, A., Díaz, J. E., Foulquié-Moreno, A., & Mañas, M. (2023). Finite Markov chains and multiple orthogonal polynomials. arXiv:2308.00182. https://doi.org/10.48550/arXiv.2308.00182
Introduction
Introduction
This code implements the whole process described in the following associated paper to get tridiagonal stochastic matrices from the recurrence matrices of seven families of orthogonal polynomials (Hahn, Jacobi-Piñeiro, Meixner, Kravchuk, Laguerre, Charlier and Hermite). As well as their corresponding stochastic bidiagonal factorizations and steady states.
In order to use it, evaluate the whole notebook and go directly to section INSTRUCTIONS to see how the different commands work.
Weight Functions
Weight Functions
Hahn: and Discrete orthogonality over
a,b>-1
N=0,1,2...
{0,1,...,N}
In[]:=
wHahn[x_,a_,b_,N_]:=Gamma[a+x+1]Gamma[b+N-x+1]/Gamma[a+1]/Gamma[x+1]/Gamma[b+1]/Gamma[N-x+1]
Jacobi--Piñeiro: . Continuous orthogonality over
a,b>-1
[0,1]
In[]:=
wJP[x_,a_,b_]:=x^a
Meixner: and . Discrete orthogonality over
b>0
0<c<1
{0,1,2...}
In[]:=
wMeixner[x_,b_,c_]:=Gamma[b+x]c^x/Gamma[b]/Gamma[x+1]
Kravchuk: ; Discrete orthogonality over
p<0<1
N=0,1,2...
{0,1,...,N}
In[]:=
wKravchuk[x_,p_,N_]:=p^x(1-p)^(N-x)Gamma[N+1]/Gamma[x+1]/Gamma[N-x+1]
Laguerre: . Continuous orthogonality over
a>-1
[0,∞)
In[]:=
wLaguerre[x_,a_]:=E^(-x)x^a
Charlier: . Discrete orthogonality over
b>0
{0,1,2...}
In[]:=
wCharlier[x_,b_]:=b^x/Gamma[x+1]
Hermite: Continuous orthogonality over
(-∞,∞)
In[]:=
wHermite[x_]:=E^(-x^2)
Recurrence Coefficients
Recurrence Coefficients
b𝑚 of the main diagonal and c𝑚 of the first subdiagonal of the Jacobi matrix
In[]:=
bHahn[m_,a_,b_,N_]:=If[m0,(N)(a+1)/(a+b+2),(N-m)(a+m+1)(a+b+m+1)/(a+b+2m+1)/(a+b+2m+2)+m(b+m)(a+b+N+m+1)/(a+b+2m)/(a+b+2m+1)]
In[]:=
cHahn[m_,a_,b_,N_]:=If[m1,((N)(a+1)/(a+b+2))((b+1)(a+b+N+1+1)/(a+b+2)/(a+b+2+1)),((N-(m-1))(a+(m-1)+1)(a+b+(m-1)+1)/(a+b+2(m-1)+1)/(a+b+2(m-1)+2))(m(b+m)(a+b+N+m+1)/(a+b+2m)/(a+b+2m+1))]
bJP[m_,a_,b_]:=If[m==0,1/2-(b-a)/2/(2m+a+b+2),1/2-(b^2-a^2)/2/(2m+a+b)/(2m+a+b+2)]
In[]:=
cJP[m_,a_,b_]:=m(m+a)(m+b)(m+a+b)/(2m+a+b-1)/(2m+a+b)/(2m+a+b)/(2m+a+b+1)
In[]:=
bMeixner[m_,b_,c_]:=(m+(b+m)c)/(1-c)
In[]:=
cMeixner[m_,b_,c_]:=mc(b+m-1)/(1-c)/(1-c)
In[]:=
bKravchuk[m_,p_,N_]:=(N-m)p+m(1-p)
In[]:=
cKravchuk[m_,p_,N_]:=m(1-p)(N-m+1)p
bLaguerre[m_,a_]:=2m+a+1
In[]:=
cLaguerre[m_,a_]:=m(m+a)
In[]:=
bCharlier[m_,b_]:=b+m
In[]:=
cCharlier[m_,b_]:=bm
In[]:=
bHermite[m_]:=0
In[]:=
cHermite[m_]:=m/2
Jacobi Matrix
Jacobi Matrix
ofsizem,J𝑚=
b0 | 1 | 0... |
c1 | b1 | 1... |
0 | c2 | b2... |
In[]:=
JHahn[m_,a_,b_,N_]:=Table[cHahn[i,a,b,N]KroneckerDelta[i,j+1]+bHahn[i,a,b,N]KroneckerDelta[i,j]+KroneckerDelta[i,j-1],{i,0,m-1},{j,0,m-1}]
In[]:=
JJP[m_,a_,b_]:=Table[cJP[i,a,b]KroneckerDelta[i,j+1]+bJP[i,a,b]KroneckerDelta[i,j]+KroneckerDelta[i,j-1],{i,0,m-1},{j,0,m-1}]
In[]:=
JMeixner[m_,b_,c_]:=Table[cMeixner[i,b,c]KroneckerDelta[i,j+1]+bMeixner[i,b,c]KroneckerDelta[i,j]+KroneckerDelta[i,j-1],{i,0,m-1},{j,0,m-1}]
In[]:=
JKravchuk[m_,p_,N_]:=Table[cKravchuk[i,p,N]KroneckerDelta[i,j+1]+bKravchuk[i,p,N]KroneckerDelta[i,j]+KroneckerDelta[i,j-1],{i,0,m-1},{j,0,m-1}]
In[]:=
JLaguerre[m_,a_]:=Table[cLaguerre[i,a]KroneckerDelta[i,j+1]+bLaguerre[i,a]KroneckerDelta[i,j]+KroneckerDelta[i,j-1],{i,0,m-1},{j,0,m-1}]
In[]:=
JCharlier[m_,b_]:=Table[cCharlier[i,b]KroneckerDelta[i,j+1]+bCharlier[i,b]KroneckerDelta[i,j]+KroneckerDelta[i,j-1],{i,0,m-1},{j,0,m-1}]
In[]:=
JHermite[m_]:=Table[cHermite[i]KroneckerDelta[i,j+1]+bHermite[i]KroneckerDelta[i,j]+KroneckerDelta[i,j-1],{i,0,m-1},{j,0,m-1}]
Orthogonal Polynomials
Orthogonal Polynomials
p𝑚(x) in their monic version
In[]:=
Hahn[m_,x_,a_,b_,N_]:=Pochhammer[a+1,m]Pochhammer[-N,m]Sum[Pochhammer[-m,l]Pochhammer[-x,l]Pochhammer[a+b+m+1,l]/Pochhammer[-N,l]/Pochhammer[a+1,l]/l!,{l,0,m}]/Pochhammer[a+b+m+1,m]
In[]:=
JP[m_,x_,a_,b_]:=(-1)^mPochhammer[a+1,m]Sum[Pochhammer[-m,l]Pochhammer[a+b+m+1,l]x^l/Pochhammer[a+1,l]/l!,{l,0,m}]/Pochhammer[a+b+m+1,m]
In[]:=
Meixner[m_,x_,b_,c_]:=(c/(c-1))^mPochhammer[b,m]Sum[Pochhammer[-m,l]Pochhammer[-x,l]((c-1)/c)^l/Pochhammer[b,l]/l!,{l,0,m}]
In[]:=
Kravchuk[m_,x_,p_,N_]:=p^mPochhammer[-N,m]Sum[Pochhammer[-m,l]Pochhammer[-x,l](1/p)^l/Pochhammer[-N,l]/l!,{l,0,m}]
In[]:=
Laguerre[m_,x_,a_]:=(-1)^mPochhammer[a+1,m]Sum[Pochhammer[-m,l]x^l/Pochhammer[a+1,l]/l!,{l,0,m}]
In[]:=
Charlier[m_,x_,b_]:=(-b)^mSum[Pochhammer[-m,l]Pochhammer[-x,l](-1/b)^l/l!,{l,0,m}]
In[]:=
Hermite[m_,x_]:=(-1)^mPi^(1/2)Sum[Pochhammer[-m,l]x^l/Gamma[(-m+1+l)/2]/l!,{l,0,m}]
Tridiagonal Stochastic Matrix
Tridiagonal Stochastic Matrix
First we compute the conjugation matrices σ𝑚
In[]:=
SigmaHahn[m_,x_,a_,b_,N_]:=Table[Hahn[i,x,a,b,N]KroneckerDelta[i,j],{i,0,m-1},{j,0,m-1}]
In[]:=
SigmaJP[m_,x_,a_,b_]:=Table[JP[i,x,a,b]KroneckerDelta[i,j],{i,0,m-1},{j,0,m-1}]
In[]:=
SigmaMeixner[m_,x_,b_,c_]:=Table[Meixner[i,x,b,c]KroneckerDelta[i,j],{i,0,m-1},{j,0,m-1}]
In[]:=
SigmaKravchuk[m_,x_,p_,N_]:=Table[Kravchuk[i,x,p,N]KroneckerDelta[i,j],{i,0,m-1},{j,0,m-1}]
In[]:=
SigmaLaguerre[m_,x_,a_]:=Table[Laguerre[i,x,a]KroneckerDelta[i,j],{i,0,m-1},{j,0,m-1}]
In[]:=
SigmaCharlier[m_,x_,b_]:=Table[Charlier[i,x,b]KroneckerDelta[i,j],{i,0,m-1},{j,0,m-1}]
In[]:=
SigmaHermite[m_,x_]:=Table[Hermite[i,x]KroneckerDelta[i,j],{i,0,m-1},{j,0,m-1}]
Now we compute the last zero of the m-th polynomial p𝑚(x) numerically. We denote it by λ
In[]:=
LambdaHahn[m_,a_,b_,N_]:=x/.Solve[Hahn[m,x,a,b,N]==0,x][[m]]
In[]:=
LambdaJP[m_,a_,b_]:=x/.Solve[JP[m,x,a,b]==0,x][[m]]
In[]:=
LambdaMeixner[m_,b_,c_]:=x/.Solve[Meixner[m,x,b,c]==0,x][[m]]
In[]:=
LambdaKravchuk[m_,p_,N_]:=x/.Solve[Kravchuk[m,x,p,N]==0,x][[m]]
In[]:=
LambdaLaguerre[m_,a_]:=x/.Solve[Laguerre[m,x,a]==0,x][[m]]
In[]:=
LambdaCharlier[m_,b_]:=x/.Solve[Charlier[m,x,b]==0,x][[m]]
In[]:=
LambdaHermite[m_]:=Max[x/.Solve[Hermite[m,x]==0,x]]
With this we can compute the stochastic matrices of size m, P𝑚
In[]:=
PHahn[m_,a_,b_,N_]:=Inverse[SigmaHahn[m,LambdaHahn[m,a,b,N],a,b,N]].JHahn[m,a,b,N].SigmaHahn[m,LambdaHahn[m,a,b,N],a,b,N]/LambdaHahn[m,a,b,N]
In[]:=
PJP[m_,a_,b_]:=Inverse[SigmaJP[m,LambdaJP[m,a,b],a,b]].JJP[m,a,b].SigmaJP[m,LambdaJP[m,a,b],a,b]/LambdaJP[m,a,b]
In[]:=
PMeixner[m_,b_,c_]:=Inverse[SigmaMeixner[m,LambdaMeixner[m,b,c],b,c]].JMeixner[m,b,c].SigmaMeixner[m,LambdaMeixner[m,b,c],b,c]/LambdaMeixner[m,b,c]
In[]:=
PKravchuk[m_,p_,N_]:=Inverse[SigmaKravchuk[m,LambdaKravchuk[m,p,N],p,N]].JKravchuk[m,p,N].SigmaKravchuk[m,LambdaKravchuk[m,p,N],p,N]/LambdaKravchuk[m,p,N]
In[]:=
PLaguerre[m_,a_]:=Inverse[SigmaLaguerre[m,LambdaLaguerre[m,a],a]].JLaguerre[m,a].SigmaLaguerre[m,LambdaLaguerre[m,a],a]/LambdaLaguerre[m,a]
In[]:=
PCharlier[m_,b_]:=Inverse[SigmaCharlier[m,LambdaCharlier[m,b],b]].JCharlier[m,b].SigmaCharlier[m,LambdaCharlier[m,b],b]/LambdaCharlier[m,b]
In[]:=
PHermite[m_]:=N[Inverse[SigmaHermite[m,LambdaHermite[m]]].JHermite[m].SigmaHermite[m,LambdaHermite[m]]/LambdaHermite[m]]
Steady States
Steady States
π𝑚
In[]:=
SteadyHahn[m_,a_,b_,N_]:=Table[Hahn[i,LambdaHahn[m,a,b,N],a,b,N]^2/Sum[Hahn[i,k,a,b,N]^2wHahn[k,a,b,N],{k,0,N}],{i,0,m-1}]/Sum[Hahn[i,LambdaHahn[m,a,b,N],a,b,N]^2/Sum[Hahn[i,k,a,b,N]^2wHahn[k,a,b,N],{k,0,N}],{i,0,m-1}]
In[]:=
SteadyJP[m_,a_,b_]:=Table[JP[i,LambdaJP[m,a,b],a,b]^2/Integrate[JP[i,x,a,b]^2wJP[x,a,b](1-x)^b,{x,0,1}],{i,0,m-1}]/Sum[JP[i,LambdaJP[m,a,b],a,b]^2/Integrate[JP[i,x,a,b]^2wJP[x,a,b](1-x)^b,{x,0,1}],{i,0,m-1}]
In[]:=
SteadyMeixner[m_,b_,c_]:=Table[Meixner[i,LambdaMeixner[m,b,c],b,c]^2/Sum[Meixner[i,k,b,c]^2wMeixner[k,b,c],{k,0,Infinity}],{i,0,m-1}]/Sum[Meixner[i,LambdaMeixner[m,b,c],b,c]^2/Sum[Meixner[i,k,b,c]^2wMeixner[k,b,c],{k,0,Infinity}],{i,0,m-1}]
In[]:=
SteadyKravchuk[m_,p_,N_]:=Table[Kravchuk[i,LambdaKravchuk[m,p,N],p,N]^2/Sum[Kravchuk[i,k,p,N]^2wKravchuk[k,p,N],{k,0,N}],{i,0,m-1}]/Sum[Kravchuk[i,LambdaKravchuk[m,p,N],p,N]^2/Sum[Kravchuk[i,k,p,N]^2wKravchuk[k,p,N],{k,0,N}],{i,0,m-1}]
Expected Return Times
Expected Return Times
LU Decomposition
LU Decomposition
This decomposition is not possible for the Hermite family since the main diagonal of the recurrence matrix is null
Stochastic Bidiagonal Factorization
Stochastic Bidiagonal Factorization
We start by defining the auxiliar matrices
Instructions
Instructions
Here comes a summary of the commands to get the stochastic matrices, their respective bidiagonal factorizations and the steady states for the seven families.
Let's remind first all these families and their possible choice of parameters
◼
HERMITE: There are no parameters
Stochastic Matrix of size m corresponding to the Hahn/Jacobi-Piñeiro/Meixner/Kravchuk/Laguerre/Charlier/Hermite polynomials
PHahn[m,a,b,N]/PJP[m,a,b]/PMeixner[m,b,c]/PKravchuk[m,p,N]/PLaguerre[m,a]/PCharlier[m,b]/PHermite[m]
PHahn[m,a,b,N]/PJP[m,a,b]/PMeixner[m,b,c]/PKravchuk[m,p,N]/PLaguerre[m,a]/PCharlier[m,b]/PHermite[m]
Examples
Examples
Lower Bidiagonal Stochastic Matrix of size m corresponding to the Hahn/Jacobi-Piñeiro/Meixner/Kravchuk/Laguerre/Charlier polynomials
PiHahn[m,a,b,N]/PiJP[m,a,b]/PiMeixner[m,b,c]/PiKravchuk[m,p,N]/PiLaguerre[m,a]/PiCharlier[m,b]
PiHahn[m,a,b,N]/PiJP[m,a,b]/PiMeixner[m,b,c]/PiKravchuk[m,p,N]/PiLaguerre[m,a]/PiCharlier[m,b]
Examples
Examples
Upper Bidiagonal Stochastic Matrix of size m corresponding to the Hahn/Jacobi-Piñeiro/Meixner/Kravchuk/Laguerre/Charlier polynomials
UpsilonHahn[m,a,b,N]/UpsilonJP[m,a,b]/UpsilonMeixner[m,b,c]/UpsilonKravchuk[m,p,N]/UpsilonLaguerre[m,a]/UpsilonCharlier[m,b]
UpsilonHahn[m,a,b,N]/UpsilonJP[m,a,b]/UpsilonMeixner[m,b,c]/UpsilonKravchuk[m,p,N]/UpsilonLaguerre[m,a]/UpsilonCharlier[m,b]
Example
Example
Steady State of the tridiagonal stochastic matrix of size m corresponding to Hahn/Jacobi-Piñeiro/Meixner/Kravchuk/Laguerre/Charlier/Hermite polynomials
SteadyHahn[m,a,b,N]/SteadyJP[m,a,b]/SteadyMeixner[m,b,c]/SteadyKravchuk[m,p,N]/SteadyLaguerre[m,a]/SteadyCharlier[m,b]/SteadyHermite[m]
SteadyHahn[m,a,b,N]/SteadyJP[m,a,b]/SteadyMeixner[m,b,c]/SteadyKravchuk[m,p,N]/SteadyLaguerre[m,a]/SteadyCharlier[m,b]/SteadyHermite[m]
Examples
Examples
Expected Return Times for the m states of the Markov chain associated to the tridiagonal stochastic matrix corresponding to the Hahn/Jacobi Piñeiro/Meixner/Kravchuk/Laguerre/Charlier/Hermite polynomials
ReturnTimesHahn[m,a,b,N]/ReturnTimesJP[m,a,b]/ReturnTimesMeixner[m,b,c]/ReturnTimesKravchuk[m,p,N]/ReturnTimesLaguerre[m,a]/ReturnTimesCharlier[m,b]/ReturnTimesHermite[m]
ReturnTimesHahn[m,a,b,N]/ReturnTimesJP[m,a,b]/ReturnTimesMeixner[m,b,c]/ReturnTimesKravchuk[m,p,N]/ReturnTimesLaguerre[m,a]/ReturnTimesCharlier[m,b]/ReturnTimesHermite[m]
Examples
Examples
To get a matrix format, implement the instruction //MatrixForm after the matrix.
IMPORTANT: Be careful since this instruction does NOT allow to do operations with the matrices.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Finite Markov chains and orthogonal polynomials
by Amílcar Branquinho, Juan E. F. Díaz, Ana Foulquié-Moreno & Manuel Mañas
Wolfram Community, STAFF PICKS, May 6, 2024
https://community.wolfram.com/groups/-/m/t/3190529
by Amílcar Branquinho, Juan E. F. Díaz, Ana Foulquié-Moreno & Manuel Mañas
Wolfram Community, STAFF PICKS, May 6, 2024
https://community.wolfram.com/groups/-/m/t/3190529