TF[xs_]:=TableForm[xs,TableDepth->2]

MF[M_]:=MatrixForm[M]

Binomial transform

BM[n_]:=Table[Binomial[x,y],{x,0,n-1},{y,0,n-1}]

DM[xs_]:=DiagonalMatrix[xs]

For weights γ(a,b) and γ(c) have the eigenvalues below

λ[a_,b_,d_]:=Binomial[a+d,d]/Binomial[a+b+d+1,d]

λs[n_,a_,b_]:=Table[λ[a,b,d],{d,0,n-1}]

λ[c_,d_]:=1/(c+1)^d

λs[n_,c_]:=Table[λ[c,d],{d,0,n-1}]

Also define a generic version

λs[n_]:=Table[λ[d],{d,0,n-1}]

Involution J

JEntry[n_,x_,y_]:=If[x+y==n-1,1,0]

J[n_]:=Table[JEntry[n,x,y],{x,0,n-1},{y,0,n-1}]

Binomial transform, defined so that the columns of binomial matrix are eigenvectors

BinomialTransform[X_]:=Module[{n=Length[X]},BM[n].X.Inverse[BM[n]]]

H[λs_List]:=BinomialTransform[DM[λs]];P[λs_List]:= H[λs].J[λs//Length];

H[n_,a_,b_]:=H[λs[n,a,b]]; P[n_,a_,b_]:=H[n,a,b].J[n]

H[n_,c_]:=H[λs[n,c]]; P[n_,c_]:=H[n,c].J[n]

H[n_]:=H[λs[n]];P[n_]:=H[n].J[n]

Hδ[n_,aP_,bP_]:=H[n,-aP,-bP];Pδ[n_,aP_,bP_]:=Hδ[n,aP,bP].J[n]

Example in introduction

MF/@{P[4,0,0],P[4,1,0],P[4,0,1],P[4,1/2],P[4,2],P[4,-4,-2]}

Number of anti-diagonal bands and (1/1-μ,2) weight for δ

MF/@{P[4,-4,-1],P[4,-5,-3],P[3,-1/(1-μ),-2]//Simplify}

Domain of δ: we need n < Ceiling(a’) so that extreme for Binomial(a’-1,y) is Binomial(a’-1,n-1) = (a’-1)...(a’-n+1)/n!, with no minus signs or zeros (when a’ is natural) provided that a’ > n-1 <=> Ceiling(a’) > n-1 <=> Ceiling(a’) >= n (comparing two natural numbers).

MF/@{Pδ[4,7/2,10],Pδ[5,7/2,10],Pδ[5,4,10],Pδ[5,5,10]}

Stochastic matrices, Corollary 6.3. Can have a more complicated zero region (when not reversible; in this case even have γ weight) than a triangle. Might be able to characterise irreducibility, but let’s not.

P[{1,λ[1],λ[2],λ[3],λ[3],λ[3]}]//MF

MF/@{P[{1,4/7,3/7,2/7,1/7}],P[{1,5/7+ϵ,4/7,3/7,2/7,1/7}]/.{ϵ->1/35},P[{1,9/10,9/10,9/10,8/10,8/10,8/10,8/10}]}//Simplify

Matrix in basis of columns of binomial matrix: for Ps is upper triangular with alternating sign eigenvalues

IBT[M_]:=Module[{n=Length[M]},Inverse[BM[n]].M.BM[n]]

MF/@IBT/@{P[4,0,0],P[4,1,0],P[4,0,1],P[4,1],P[4,-4,-2]}

MF/@IBT/@Table[P[n,0,1],{n,1,5}]

Bottom row entries: had a proposition finding them from eigenvalues, but this is unnecessarily complicated since given directly by the weight

HBottom[n_,a_,b_,y_]:=H[n,a,b][[n]][[y+1]]

HBottom[n_,a_,b_]:=Table[HBottom[n,a,b,y],{y,0,n-1}]

HBottomByWeight[n_,a_,b_,y_]:=Binomial[y+a,y]*Binomial[b+n-1-y,n-1-y]/Binomial[n+a+b,n-1]

HBottomByWeight[n_,a_,b_]:=Table[HBottomByWeight[n,a,b,y],{y,0,n-1}]

HBottomByAltFormula[n_,a_,b_,y_]:=Binomial[a+b+1,a]*(b+1)/Binomial[a+b+n-1,b+n-1-y]/(a+b+n)

HBottomByAltFormula[n_,a_,b_]:=Table[HBottomByAltFormula[n,a,b,y],{y,0,n-1}]

{HBottom[4,6,7],HBottomByWeight[4,6,7],HBottomByAltFormula[4,6,7],HBottom[4,-5,-2],HBottomByWeight[4,-5,-2]//Quiet,HBottomByAltFormula[4,-5,-2]//Quiet}

Entries in binomial transform and Lemma 2.2(iv)

Diff[λs_List,k_,y_]:=Sum[(-1)^e*Binomial[k,e]*λs[[y+e+1]],{e,0,k}]

Diff[k_,y_]:=Diff[λs[k+y+1],k,y]

Global anti-diagonal property and Ochiai et al proof

BottomRowMatrix[xs_]:=Module[{n=xs//Length},Table[Table[0,{y,0,n-1}],{x,0,n-2}]~Join~{xs}]

HE[n_]:=H[λs[n]]+BottomRowMatrix[Table[ϵ[y],{y,0,n-1}]]/.{ϵ[n-1]->0}/.{ϵ[0]->0}

PE[n_]:=HE[n].J[n]

PEX[n_]:=PE[n]-X*IdentityMatrix[n]

By choice that ϵ[0]=0, PE and P have the same trace (necessary condition), so we do not need a coefficient equation for n-1

By choice that ϵ[n-1]=0, PE and P have the same determinant, as obtained by expanding on the antidiagonal (necessary condition) so we do not need a coefficient equation for 0.

By multilinearity of the determinant, we can compute det PE - det P as the determinant of the single matrix obtained by subtracting the bottom row of P from PE, i.e. leaving just 0, ϵ[n-1], ..., ϵ[1], 0

NM[n_]:=PEX[n]-BottomRowMatrix[P[λs[n]][[(n-1)+1]]]

h[n_]:=Det[NM[n]]

ϵCoefficients[n_]:=Module[{es=Table[Coefficient[h[n],X^m],{m,1,n-2}]},Table[Table[Coefficient[e,ϵ[y]],{y,1,n-2}],{e,es}]]

This is the matrix C in proof

Check agrees with earlier approach, not using multilinearity of determinant

hOrig[n_]:=Det[PEX[n]]-Product[(-1)^x*λ[x]-X,{x,0,n-1}]

ϵCoefficientsOrig[n_]:=Module[{es=Table[Coefficient[hOrig[n],X^m],{m,1,n-2}]},Table[Table[Coefficient[e,ϵ[y]],{y,1,n-2}],{e,es}]]

{PEX[4]//MF,NM[4]//MF}

{ϵCoefficients[5]-ϵCoefficientsOrig[5]//MF,ϵCoefficients[5]//MF,ϵCoefficientsOrig[5]//MF}

Specialisation: to show that det ϵCoefficients is non-zero enough to show this with a specialisation of λ[0], ..., λ[n-1]. We take λ[x] -> x. (Might get smaller numbers but same pattern of non-zeros with λ[x]->n-1-x = x*.)

LSpec[n_]:=Table[λ[x]->x,{x,0,n-1}]

NMCut[n_]:=Table[NM[n][[k]][[2;;n-1]],{k,3,n}]

Example

{NM[7]/.LSpec[7]//MF,ϵCoefficients[7]/.LSpec[7]//MF}

{NMCut[7]/.LSpec[7]//MF,NMCut[8]/.LSpec[8]//MF}

The triangularity (after column rearrangement) in C is explained by rearranging the columns of NM so that it is tri-diagonal.

TridiagonalColumnOrder[n_,k_]:=If[EvenQ[k],k/2,(n-1)-(k-1)/2]

TridiagonalColumnOrder[n_]:=Table[TridiagonalColumnOrder[n,k],{k,0,n-1}]

TridiagonalRowOrder[n_,0]:=0

TridiagonalRowOrder[n_,k_]/;k<n-1:=If[OddQ[k],(k+1)/2,(n-1)-k/2]

TridiagonalRowOrder[n_,k_]/;k==n-1:=n-1

TridiagonalRowOrder[n_]:=Table[TridiagonalRowOrder[n,k],{k,0,n-1}]

NMR[n_]:=Module[{M=NM[n]},Table[Transpose[ Table[Transpose[M][[c+1]],{c,TridiagonalColumnOrder[n]}]][[r+1]],{r,TridiagonalRowOrder[n]}]]

NMRCut[n_]:=Table[NMR[n][[k]][[3;;n]],{k,3,n}]

Example in write-up

{TridiagonalRowOrder[7],TridiagonalColumnOrder[7]}

{NMR[7]/.LSpec[7]//MF,NMRCut[7]/.LSpec[7]//MF,ϵCoefficients[7]/.LSpec[7]//MF}

{TridiagonalRowOrder[8],TridiagonalColumnOrder[8]}

{NMR[8]/.LSpec[8]//MF,NMRCut[8]/.LSpec[8]//MF}

Explicit example of determinant: has degree (n-1)(n-2)/2

Det[ϵCoefficients[4]]//Factor

Example 2.5

H25:={{1,0,0,0},{1/3,2/3,0,0},{-1/12,5/6,1/4,0},{-9/20,11/10+4/5*τ,3/20+1/5*τ,1/5}}

L25:={{1,0,0,0,0},{1/2,1/2,0,0,0},{1/2-τ,2τ,1/2-τ,0,0},{1/2-2τ,3τ,0,1/2-τ,0},{1/2-4τ,5τ,0,τ,1/2-2τ}}

E25:={{1,-2,0,-4,-4τ(10τ-3)},{1,-5/4,0,-4,-28τ^2+11τ-1},{1,-1/2,1,-2τ-1,(4τ-1)(4τ^2-6τ+1)},{1,1/4,0,2-4τ,(4τ-1)(8τ^2-5τ+1)},{1,1,0,2-4τ,2τ(4τ-1)(10τ-3)}}

MF/@{H25,L25,E25}

{Eigenvalues[H25.J[4]],Eigenvalues[L25.J[5]]//Expand,Inverse[E25].L25.J[5].E25//Simplify//MF}

Section 4.1

InvariantDistribution[P[4,1]]

Check on Lemma 5.2

γ[a_,b_,y_,x_]:=Binomial[y+a,y]*Binomial[x-y+b,x-y]

Nγ[a_,b_,x_]:=Sum[γ[a,b,y,x],{y,0,x}]

CheckγNorm[n_,a_,b_]:=Table[{x,Nγ[a,b,x],Binomial[x+1+a+b,x]},{x,0,n-1}]

Table[CheckγNorm[5,a,b],{a,0,4},{b,0,4}]//TF

Discrete Jacobi: Section 5.4

v[n_,d_]:=Table[Binomial[x,d],{x,0,n-1}]

MatrixProduct[{M_}]:=M

MatrixProduct[{M_,Ms___}]:=M.MatrixProduct[Ms//List]

NormalizeVectorByMax[xs_]:=xs/Max@@Abs/@xs

Projector[n_,a_,b_,0]:=IdentityMatrix[n]

Projector[n_,a_,b_,d_]:=Module[{PM=P[n,a,b]},MatrixProduct[Table[PM-(-1)^e*λ[a,b,e]*IdentityMatrix[n],{e,0,d-1}]]]

w[n_,a_,b_,e_]:=v[n,e].Transpose[Projector[n,a,b,e]]

wNM[n_,a_,b_,e_]:=w[n,a,b,e]//NormalizeVectorByMax

wNMS[n_,a_,b_,e_]:=Table[{d/n,wNM[n,a,b,e][[d+1]]},{d,0,n-1}]

ws[n_,a_,b_]:=Table[w[n,a,b,e],{e,0,n-1}]//Transpose

wsNM[n_,a_,b_]:=Table[wNM[n,a,b,e],{e,0,n-1}]//Transpose

JacobiPN[d_,A_,B_,y_]:=Module[{f=JacobiP[d,A,B,z]},f/Abs[f/.{z->-1}]/.{z->y}]

Check on w(1)

Checkw1[n_,a_,b_]:={w[n,a,b,1]//Simplify,Table[(a+b+2)(n-1)-(2a+b+3)*x,{x,0,n-1}]*w[n,a,b,1][[1]]/(a+b+2)/(n-1)}

Table[Checkw1[5,a,b],{a,0,5},{b,0,5}]//TF

Check on Corollary 5.6

CheckwLast[n_,a_,b_]:={w[n,a,b,n-1]/w[n,a,b,n-1][[1]],Table[(-1)^x*Binomial[n-1,x]/Binomial[n-x+a-1,a]/Binomial[x+1+a+b,a+b+1],{x,0,n-1}]*Binomial[n+a-1,a]}

Table[CheckwLast[5,a,b],{a,0,5},{b,0,5}]//TF

CheckwLastZero[n_,b_]:={w[n,0,b,n-1]/w[n,0,b,n-1][[1]],Table[(-1)^x*Binomial[n+b,x+b+1]/Binomial[n+b,b+1],{x,0,n-1}]}

Table[CheckwLastZero[5,b],{b,0,5}]//TF

CheckwLastOne[n_,b_]:={w[n,1,b,n-1]/w[n,1,b,n-1][[1]],Table[(-1)^x*Binomial[n+b+2,x+b+2]/Binomial[n+b+2,b+2],{x,0,n-1}]}

Table[CheckwLastOne[5,b],{b,0,5}]//TF

Example

wNMS[20,1,1,3]

Show[ListPlot[{wNMS[10,0,0,3],wNMS[40,0,0,3]},PlotRange->{{0,1},{-0.4,0.4}},PlotStyle->{Gray,Black}],Plot[JacobiPN[3,0,1,2x-1],{x,0,1},PlotRange->{-0.4,0.4},PlotStyle->{Black,Dashed}],Ticks->{{0,{0.2,Null},{0.4,Null},{0.6,Null},{0.8,Null},{1,Null}},{{-0.4,Null},{-0.2,Null},{0.2,Null},{0.4,Null}}}]

Show[ListPlot[{wNMS[10,0,0,3],wNMS[40,0,0,3]},PlotRange->{{0,1},{-0.4,0.4}},PlotStyle->{Gray,Black}],Plot[JacobiPN[3,0,1,2x-1],{x,0,1},PlotRange->{-0.4,0.4},PlotStyle->{Black,Dashed}],Ticks->{{0,{0.25,Null},{0.5,Null},{0.75,Null},{1,Null}},{{-0.4,Null},{-0.2,Null},{0.2,Null},{0.4,Null}}}]

JacobiPN[3,0,1,2x-1]

Some binomial identities

See https : // wildonblog.wordpress.com/2020/12/28/signed-binomial-matrices-of-small-order/

It follows easily from Lemma 2.1 that B(n)^{-1}JB(n) = Diag(1,-1,...)JB(n)J and, harder, that Diag(1,-1,...)JB(n) has cube (-1)^(n-1)I(n) (elementary proof given in blog post)

DPM[n_]:=DM[Table[(-1)^x,{x,0,n-1}]]

BJRelations[n_]:=MF/@{Inverse[BM[n]].J[n].BM[n],DPM[n].J[n].BM[n].J[n],BM[n].J[n].Inverse[BM[n]],J[n].BM[n].J[n].DPM[n],MatrixPower[DPM[n].J[n].BM[n],3]};{BJRelations[5]}

We verify some identities stated without proof in the blog post on binomial matrices

BM[n_,a_]:=Table[Binomial[x+a,y+a],{x,0,n-1},{y,0,n-1}]

Rot[M_]:= Module[{n=Length[M]},J[n].M.J[n]]

TestPowerIdentity[n_,a_,b_,p_]:=MatrixPower[BM[n,a].DPM[n].Rot[BM[n,b]].DPM[n],p];{MF/@{TestPowerIdentity[5,0,0,3],TestPowerIdentity[5,0,1,12],TestPowerIdentity[5,0,1,6],TestPowerIdentity[5,0,5,3],TestPowerIdentity[5,1,1,3]},MF/@{TestPowerIdentity[6,0,0,3],TestPowerIdentity[6,0,1,6],TestPowerIdentity[6,0,1,12],TestPowerIdentity[6,0,6,3],TestPowerIdentity[6,1,1,6]}}

Final remark in blog post

{BM[6,0].DPM[6].Rot[BM[6,1]].DPM[6]//MF,MatrixPower[BM[6,0].DPM[6].Rot[BM[6,1]].DPM[6],12]//MF}

{BM[5,0].DPM[5].Rot[BM[5,4]].DPM[5].BM[5,0].DPM[5].Rot[BM[5,8]].DPM[5]//MF,MatrixPower[BM[5,0].DPM[5].Rot[BM[5,4]].DPM[5].BM[5,0].DPM[5].Rot[BM[5,8]].DPM[5],12]//MF}

Section 6: Stochastic P^λ

FMatrix: defined for weights and for general eigenvalues by inverse binomial transform

FMatrix[λs_]:=Module[{HM=H[λs],n=Length[λs]},Table[If[x≥z,HM[[x+1]][[z+1]]/Binomial[x,z],0],{x,0,n-1},{z,0,n-1}]]//Simplify

Example of F matrix and H with repeated eigenvalues

MF/@{FMatrix[{1,1/2,1/3,1/4}],H[5]/.{λ[0]->1,λ[3]->λ[2],λ[4]->λ[2]}}

Example with behaviour from Lemma 6.4

MF/@{H[{1,1/2,ν,ν,ν}]/.{ν->1/3},H[{1,1/2,ν,ν,ν,ν}]/.{ν->3/8}}

Example with several large zero triangles

Section 7: Reversibility of P^λ, under global reversibility hypothesis

μAB[a_,b_]:=λs[3,a,b][[2]]

νAB[a_,b_]:=λs[3,a,b][[3]]

Base case motivates definition of a and b as functions of μ and ν

H[3,a,b]//Simplify//MF

a[μ_,ν_]:=μ(μ-ν)/(ν-μ^2)-1

b[μ_,ν_]:=(1-μ)(μ-ν)/(ν-μ^2)-1

{(1+a[μ,ν])/(2+a[μ,ν]+b[μ,ν])//Simplify,(1+a[μ,ν])(2+a[μ,ν])/(2+a[μ,ν]+b[μ,ν])/(3+a[μ,ν]+b[μ,ν])//Simplify}

bP is the negation of b

bP[μ_,ν_]:=1+(1-μ)(μ-ν)/(μ^2-ν)

Critical values where bP[μ,ν] is a natural number

ν[m_,μ_]:=μ(m*μ-1)/(m-2+μ)

aPm[m_,μ_]:=((m-2)μ+1)/(1-μ)

Check on critical values

{-a[μ,ν[m,μ]],aPm[m,μ],-a[μ,ν[m,μ]]-aPm[m,μ],bP[μ,ν[m,μ]]}//Simplify

aPm[2,μ]

There is a constraint on the number of bands in terms of μ

mLeastWrong[n_,μ_]:=Ceiling[n(1-μ)/μ-2/μ]+4

mLeast[n_,μ_]:=Floor[(1-μ)/μ*(n-2)]+2

For instance if μ=1/2, this says m >= n, but can never have more than n bands, so the δ case is only relevant when μ>1/2. Or, another way, if μ<1/2 then μ =(a’-1)/((a’-1)+(m-1)) so we have m > a’, but always n <= Ceiling(a’), so doesn’t arise. (Included in revised build-up to statement of sharp result.)

Plot[{10,mLeast[10,μ]},{μ,0.4,1},PlotStyle->{Black,Blue}]

Example 7.4

Table[3.1-0.575*k,{k,0,5}]

{a[2/3,ν]//Simplify,b[2/3,ν]//Simplify}

Table[{m,ν[m,2/3]},{m,1,10}]//Transpose//TF

e:=1/100;Show[{Plot[{a[2/3,ν],b[2/3,ν],-1},{ν,0,2/3},PlotStyle->{Black,{Black,Dashed},{Black,Dotted}},PlotRange->{-6.5,3}],Plot[-2,{x,1/3-e,1/3+e},PlotStyle->Black],Plot[-3,{x,2/5-e,2/5+e},PlotStyle->Black],Plot[-4,{x,5/12-e,5/12+e},PlotStyle->Black],Plot[-5,{x,14/33-e,14/33+e},PlotStyle->Black],Plot[-6,{x,3/7-e,3/7+e},PlotStyle->Black]},Ticks->{{0,{1/3,Null},{2/5,Null},{5/12,Null},{14/33,Null},{3/7,Null}},{0,{-1,Null},{-2,Null},{-3,Null},{-4,Null},{-5,Null},{-6,Null}}}]

mLeast[n,2/3]

Derivative of bP showing that is strictly increasing on [0, μ^2)

D[bP[μ,ν],ν]//Simplify

Example

Plot[{bP[2/3,ν],bP[1/2,ν],bP[1/3,ν]},{ν,0,4/9},PlotStyle->{Red,Blue,Green}]

Discrete spectrum

DiscreteSpectrum[n_,μ_]:=Table[{m,ν[m,μ],-a[μ,ν[m,μ]],-b[μ,ν[m,μ]]},{m,Ceiling[mLeast[n,μ]],n-1}]

Example 7.4 and checks on lower bound on m

{DiscreteSpectrum[10,2/3]//TF,DiscreteSpectrum[11,2/3]//TF,17/39-10/23,DiscreteSpectrum[10,1/2]//TF,DiscreteSpectrum[10,2/3-1/100]//TF,DiscreteSpectrum[10,2/3+1/100]//TF}

DiscreteSpectrum[11,2/3]//TF

Thus can’t have 5 diagonal bands with μ=2/3-1/100, and correspondingly, Ceiling(a’) = 9, ruling out 10. Note negative entry in second matrix below

{-a[2/3-1/100,ν[5,2/3-1/100]],-a[2/3+1/100,ν[5,2/3+1/100]],P[10,-891/103,-5]//MF,P[10,-1083/103,-6]//MF,P[10,-909/97,-5]//MF}

Base case, including graph showing a’ is increasing

-b[μ,2μ-1+τ(1-μ)]-2//Simplify

-a[μ,2μ-1+τ(1-μ)]//Simplify

-a[μ,2μ-1+τ(1-μ)]-(1+(μ-μ*τ)/(1-μ-τ))//Simplify

Show[Plot[-a[μ,2μ-1+τ(1-μ)]/.{μ->2/3},{τ,0,1-2/3-0.1},PlotRange->{2,5}],Plot[-a[μ,2μ-1+τ(1-μ)]/.{μ->3/5},{τ,0,1-3/5-0.1}]]

Do we need global reversibility? Probably not

NormalizeToDistribution[xs_]:=Module[{T=Plus@@xs},#1/T&/@xs]

InvariantDistribution[P_]:=Eigenvectors[Transpose[P]][[1]]//NormalizeToDistribution//FullSimplify

ClearAll[IsReversibleFor]

IsReversibleFor[PM_,ps_]:=Module[{n=Length[PM]},And@@Table[And@@Table[ps[[x]]*PM[[x]][[y]]==ps[[y]]*PM[[y]][[x]],{y,1,n}],{x,1,n}]]

IsReversible[PM_]:=IsReversibleFor[PM,InvariantDistribution[PM]]

Kolmogorov condition on cycles: using 0 to n-1 labelling

KolmogorovProduct[P_,xs_]:=Module[{n=Length[P]},
Product[P[[xs[[i]]+1]][[xs[[i+1]]+1]],{i,1,Length[xs]-1}]*P[[xs[[-1]]+1]][[xs[[1]]+1]]]

KolmogorovProducts[P_,xs_]:={KolmogorovProduct[P,xs],KolmogorovProduct[P,Reverse[xs]]}

KolmogorovDifference[P_,xs_]:=KolmogorovProduct[P,xs]-KolmogorovProduct[P,Reverse[xs]]//Simplify

PWalks[n_,0]:=Table[{x},{x,0,n-1}]

PWalks[n_,l_]:=Join@@Table[Table[p~Join~{x},{x,n-1-p[[-1]],n-1}],{p,PWalks[n,l-1]}]

PPaths[n_,l_]:=Select[PWalks[n,l],DeleteDuplicates[#]==#&]

IsWalkUntilEnd[p_]:=DeleteDuplicates[p[[1;;-2]]]==p[[1;;-2]]

PCyclesWithStartAndEnd[n_,l_]:=Select[PWalks[n,l],IsWalkUntilEnd[#]~And~(#[[-1]]==#[[1]])&]

PCyclesAll[n_,l_]:=#[[1;;-2]]&/@PCyclesWithStartAndEnd[n,l]

PCycles[n_,l_]:=Select[PCyclesAll[n,l],Min[#]==#[[1]]&]

PCyclesUpToReverse[n_,l_]/;l≤2:=PCycles[n,l]

PCyclesUpToReverse[n_,l_]:=Select[PCycles[n,l],#[[-1]]>#[[2]]&]

AllKolmogorovProducts[P_,l_]:=Table[{xs,KolmogorovProducts[P,xs]},{xs,PCyclesUpToReverse[Length[P],l]}]

AllKolmogorovProducts[P_]:=Join@@Table[AllKolmogorovProducts[P,l],{l,3,P//Length}]

AllKolmogorovDifferences[P_,l_]:=Table[{xs,KolmogorovDifference[P,xs]},{xs,PCyclesUpToReverse[Length[P],l]}]

AllKolmogorovDifferences[P_]:=Join@@Table[AllKolmogorovDifferences[P,l],{l,3,P//Length}]

SolveKolmogorovDifferences[P_,l_,ys_]:=Module[{ds=AllKolmogorovDifferences[P,l]},Solve[Table[d[[2]]==0,{d,ds}],ys]]

SolveKolmogorovDifferences[P_,ys_]:=Module[{ds=AllKolmogorovDifferences[P]},Solve[Table[d[[2]]==0,{d,ds}],ys]]

SolveKolmogorovDifferencesAB[P_,l_,ys_]:=Evaluate[SolveKolmogorovDifferences[P,l,ys]]/.{μ->μAB[a,b],ν->νAB[a,b]}//Factor

SolveKolmogorovDifferencesAB[P_,ys_]:=Evaluate[SolveKolmogorovDifferences[P,ys]]/.{μ->μAB[a,b],ν->νAB[a,b]}//Factor

For inductive approach, suppose have n-1 known eigenvalues given by an (a,b)-weight, and get condition on last

KolmogorovProductAB[n_,a_,b_,xs_]:=Module[{PM=P[n,a,b]},Product[PM[[xs[[i]]+1]][[xs[[i+1]]+1]],{i,1,Length[xs]-1}]*PM[[xs[[-1]]+1]][[xs[[1]]+1]]]//Expand

KolmogorovProductsAB[n_,a_,b_,xs_]:={KolmogorovProductAB[n,a,b,xs],KolmogorovProductAB[n,a,b,Reverse[xs]]}

Example: two reversible, then one non-reversible

AllKolmogorovProducts[P[5,0,0],3]

AllKolmogorovProducts[P[6,0,0],4]

AllKolmogorovProducts[P[{1,1/2,1/3,1/5}],3]

Generic 4 by 4 binomial transform is reversible iff r4 vanishes: can pick any μ and ν and then only solution is the eigenvalue for the binomial weight

{AllKolmogorovDifferences[P[{1,μ,ν,θ}],3]/.{μ->1/2,ν->1/3},AllKolmogorovDifferences[P[{1,μ,ν,θ}],3]//Expand}

Collect[-3 θ μ+6 θ (μ^2)-3 θ ν-3 (μ^2) ν+6 (ν^2)-3 μ (ν^2),θ]

r4[μ_,ν_]:=θ/.SolveKolmogorovDifferences[P[{1,μ,ν,θ}],3,{θ}][[1]]//Simplify

{r4[μ,ν],r4[μAB[a,b],νAB[a,b]]//Factor}

SolveKolmogorovDifferencesAB[P[{1,μ,ν,θ}],3,{θ}]

Generic 5 by 5 binomial transform: no more stochastic solutions

AllKolmogorovDifferences[P[{1,μ,ν,θ,ϕ}]]//TF

AllKolmogorovDifferences[P[{1,μ,ν,θ,ϕ}]]/.{μ->1/2,ν->1/3}//TF

SolveKolmogorovDifferences[P[{1,1/2,1/3,θ,ϕ}],{θ,ϕ}]

SolveKolmogorovDifferencesAB[P[{1,μ,ν,θ,ϕ}],{θ,ϕ}]//TF

Non-stochastic example that is algebraically reversible: curiously with respect to an ‘invariant distribution’ that is a signed row of Pascal’s Triangle

{P5=P[{1,1/2,1/3,1/2,1}];P5//MF,Eigenvectors[Transpose[P5]]//MF,Eigenvalues[Transpose[P5]],{{1,-2,3,0,1}}.P5,IsReversibleFor[P5,{1,-2,3,0,1}+α{0,1/2,-3/4,1,0}]//Simplify,IsReversibleFor[P5,{1,-2,3,0,1}-4{0,1/2,-3/4,1,0}],{1,-2,3,0,1}-4{0,1/2,-3/4,1,0}}

Generic 6 by 6 binomial transform

AllKolmogorovDifferences[P[{1,μ,ν,θ,ϕ,ψ}]]//TF

SolveKolmogorovDifferencesAB[P[{1,μ,ν,θ,ϕ,ψ}],{θ,ϕ,ψ}]//TF

Generic 7 by 7 binomial transform

(*SolveKolmogorovDifferencesAB[P[{1,μ,ν,θ,ϕ,ψ,η}],{θ,ϕ,ψ,η}]//TF*)

(*P7=P[{1,μ,ν,θ,ϕ,ψ,η}]/.SolveKolmogorovDifferencesAB[P[{1,μ,ν,θ,ϕ,ψ,η}],{θ,ϕ,ψ,η}][[1]];Q7=P7/.{μ->μAB[a,b],ν->νAB[a,b]}/.{a->0,b->0}/.{θ->1/3};Q7//MF*)

Generic 8 by 8 binomial transform

(*SolveKolmogorovDifferencesAB[P[{1,μ,ν,θ,ϕ,ψ,η,ξ}],{θ,ϕ,ψ,η,ξ}]//TF*)

Directly from invariant distribution? Looks nasty, like the hunt for stoachastic GADEP matrices

(*InvariantDistribution[P[λs[5]]]//Simplify//TF*)

Continuous case

Lemma 9.1, showing that norm of κ[a,b] is a!b!/(a+b+1)! = Beta[a+1,b+1]

Integrate[w^a*(1-w)^b,{w,0,1}]

Integrate[w^a*(1-w)^b,{w,0,1}]/Beta[a+1,b+1]//FullSimplify

((Gamma[1+a] Gamma[1+b])/(Gamma[2+a+b]))-Beta[a+1,b+1]//FullSimplify

Integrate[w^a*(1-w)^b,{w,0,1}]/(a!*b!/(a+b+1)!)//FullSimplify

Integrate[y^a*(x-y)^b*(a+b+1)*Binomial[a+b,a]/x^(a+b+1),{y,0,x}]//FullSimplify

Normalization of invariant distribution

Integrate[(1-x)^a*x^(a+b+1),{x,0,1}]

Integrate[(1-x)^a*x^(a+b+1)/(a!*(a+b+1)!/(2a+b+2)!),{x,0,1}]//FullSimplify

Example: a = b = -1/2 (only allowed N0). Can see directly from (1/2)! = 1/2(-1/2)! and cancels with integral of 2.

Integrate[1/Sqrt[1-x],{x,0,1}]

Integrate[1/Sqrt[1-x],{x,0,1}]/((-1/2)!/(1/2)!)

Integrate[y*(x-y),{y,0,x}]

Trigonometric example: LH, LP acting on functions in z, returning functions in x

LHTrig[f_]:=Pi*Integrate[Sin[Pi*z]*f/(1-Cos[Pi*x]),{z,0,x}]

LPTrig[f_]:=Pi*Integrate[Sin[Pi*z]*f/(1-Cos[Pi*x]),{z,1-x,1}]

Table[LHTrig[Cos[d*Pi*x]],{d,0,5}]

Plot[Sin[π*x]*(1-Cos[π*x]),{x,0,1}]

Integrate[Pi*Sin[π*x]*(1-Cos[π*x])/2,{x,0,1}]

Lemma 7.6

L86[d_]:=Integrate[Sin[Pi*z]*Cos[Pi*z]^d,{z,1-x,1}]

Table[{d,L86[d]},{d,0,5}]//Transpose//TF

Section 9: left-multiplication operator in Theorem 9.6 and check on final remark in Section 9.3

LP[a_,b_,f_]:=(a+b+1)*Binomial[a+b,a]*Integrate[(1-z)^a*(x+z-1)^b/x^(a+b+1)*f,{z,1-x,1}]

LH[a_,b_,f_]:=(a+b+1)*Binomial[a+b,a]*Integrate[(1-z)^a*(x+z-1)^b/x^(a+b+1)*f,{z,0,x}]

Example

LH[0,0,z^d]

LP[0,0,z^d]

Final remark in Section 9.2

Integrate[1/Sqrt[y]/Sqrt[x-y],{y,0,x}]

Beta[1/2,1/2]