■

• こちらで、漸化式の話題がある
• Wikiの漸化式の記事のとおりに線形代数処理をRで行う
• １イニングあたりの得点がx=0,1,2,...である確率
• tイニングでの得点がx=0,1,2,...である確率
• 
• この関係を行列式でを列ベクトルとしてと表すとき、は正方行列で

• を用いて
• 行列のべき(パワー)は、固有値分解して…と行うわけだが(こちら)、もちろん、パッケージを使えば、その手間もいらないわけで、パッケージはexpmパッケージ、で次のようにする(検索したサイト)

M%^% t

• 地道に計算するのと、とで同じ結果になることを確かめよう

maxX<-100 # 考慮する最大得点
maxT<-5 # イニング数
maxXininng<-5 # １イニングあたりの最大得点

Ps<-matrix(0,maxT,maxX+1)
library(MCMCpack)
Ps[1,]<-rdirichlet(1,c(rep(1,maxXininng+1),rep(0,maxX-maxXininng)))

for(i in 2:maxT){
	for(j in 1:(maxX+1)){
		for(k in 1:j){
			Ps[i,j]<-Ps[i,j]+Ps[i-1,k]*Ps[1,(j-k+1)]
		}

	}
}

apply(Ps,1,sum) # t イニングでの得点確率の和(１になっている(maxXが十二分に大きい値に設定できているから)

ylim<-c(0,max(Ps))
plot(Ps[1,],type="l",ylim=ylim) # １イニングあたりの得点確率
par(new=TRUE)
plot(Ps[maxT,],type="l",col=2,ylim=ylim) # maxTイニングの得点確率

######### 線形代数で
P1<-Ps[1,]
M<-matrix(0,maxX+1,maxX+1)
for(i in 1:(maxX+1)){
	for(j in 1:i){
		M[i,j]<-P1[i-j+1]
	}
}
t<-9
library(expm)
Pt<-M%^%t %*% P1
plot(Ps[t+1,],Pt)