必ず起きることは、いつまでに起きる?3

  • n個の変数が0-1の値をとるとき、その和の最小値は0、最大値はn。その和の確率密度分布を引いてみよう
  • P^{(n)}(T)=\int_{T-1}^{T} P^{(n-1)}(t) \times (P^{(1)}(T-t) dt
  • これは、前の記事に記載した通り、0-1,1-2,2-3,...ごとに区切る必要がある
  • K \le T \le K+1なる整数Kがあるとき
  • P^{(n)}=\int_{T-1}^K P^{(n-1)}(t)dt + \int_{K}^T P^{(n-1)}(t)dtとすることで、すべてのP^{(n)}(T)について、式が漸化式で立てられる
  • これをRで実行するには、パッケージpolynomを用いると、その中には、多項式の係数で多項式を作ったり係数を取り出したりする関数、微分積分をする関数があるので、定積分を伴ったり、変数変換に対応したりする処理が可能となる
  • n=1では水平、n=2では三角屋根、n=3で丸みを帯びて、どんどんベル型になっていく

  • n=3において、構成3関数のすべてを描けば次のようになる。T=1,2での移行具合もわかる

  • 関数が実際どんな風になっているかといえば
> Ps
[[1]]
[[1]][[1]]
1 


[[2]]
[[2]][[1]]
x 

[[2]][[2]]
2 - x 


[[3]]
[[3]][[1]]
0.5*x^2 

[[3]][[2]]
-1.5 + 3*x - x^2 

[[3]][[3]]
4.5 - 3*x + 0.5*x^2 


[[4]]
[[4]][[1]]
0.1666667*x^3 

[[4]][[2]]
0.6666667 - 2*x + 2*x^2 - 0.5*x^3 

[[4]][[3]]
-7.333333 + 10*x - 4*x^2 + 0.5*x^3 

[[4]][[4]]
10.66667 - 8*x + 2*x^2 - 0.1666667*x^3 


[[5]]
[[5]][[1]]
0.04166667*x^4 

[[5]][[2]]
-0.2083333 + 0.8333333*x - 1.25*x^2 + 0.8333333*x^3 - 0.1666667*x^4 

[[5]][[3]]
6.458333 - 12.5*x + 8.75*x^2 - 2.5*x^3 + 0.25*x^4 

[[5]][[4]]
-27.29167 + 32.5*x - 13.75*x^2 + 2.5*x^3 - 0.1666667*x^4 

[[5]][[5]]
26.04167 - 20.83333*x + 6.25*x^2 - 0.8333333*x^3 + 0.04166667*x^4 


[[6]]
[[6]][[1]]
0.008333333*x^5 

[[6]][[2]]
0.05 - 0.25*x + 0.5*x^2 - 0.5*x^3 + 0.25*x^4 - 0.04166667*x^5 

[[6]][[3]]
-3.95 + 9.75*x - 9.5*x^2 + 4.5*x^3 - x^4 + 0.08333333*x^5 

[[6]][[4]]
36.55 - 57.75*x + 35.5*x^2 - 10.5*x^3 + 1.5*x^4 - 0.08333333*x^5 

[[6]][[5]]
-91.45 + 102.25*x - 44.5*x^2 + 9.5*x^3 - x^4 + 0.04166667*x^5 

[[6]][[6]]
64.8 - 54*x + 18*x^2 - 3*x^3 + 0.25*x^4 - 0.008333333*x^5 
library(polynom)
# n個の変数があると、n個の関数が必要
# それを格納するリスト
Ps<-list()
n<-10
Ps[[1]]<-list()
# n=1の場合は指定可能
Ps[[1]][[1]]<-as.polynomial(c(1))
# あとは、漸化式
for(i in 2:n){
	Ps[[i]]<-list()
	Ps[[i]][[1]]<-integral(Ps[[i-1]][[1]])
	tmpcoef<-coef(Ps[[i]][[1]])
	# 0<=T<=1と n-1<=T<=nの部分は、2つの関数にまたがらないので特別
	# この両端の対称性を使う
	Ps[[i]][[i]]<-as.polynomial(tmpcoef*((-1)^(0:(length(tmpcoef)-1))))
	# 原点移動関数を使う
	Ps[[i]][[i]]<-change.origin(Ps[[i]][[i]],-i)
	# 両端以外の関数を作る
	if(i > 2){
		for(j in 2:(i-1)){
			tmp1<-integral(Ps[[i-1]][[j-1]])
			tmp1.2<-change.origin(tmp1,-1)
			tmp2<-integral(Ps[[i-1]][[j]])
			# 2つの定積分のうちTに依存する部分を出して足し合わせる
			coef1<-coef(tmp1)
			coef1.2<--coef(tmp1.2)
			print(coef1.2)
			coef2<-coef(tmp2)
			print(coef2)
			coef<-coef1.2+coef2
			# 2つの定積分の定数部分
			C1<-sum((j-1)^(0:(length(coef1)-1))*coef1)
			C2<-sum((j-1)^(0:(length(coef2)-1))*coef2)
			coef[1]<-coef[1]-C2+C1
			Ps[[i]][[j]]<-as.polynomial(coef)
		}
	}
}
Ps


for(i in 1:length(Ps)){
	t<-seq(from=0,to=i,length=100)
	X<-matrix(0,length(Ps[[i]]),length(t))
	X2<-rep(0,length(t))
	for(j in 1:length(Ps[[i]])){
		coef<-coef(Ps[[i]][[j]])
		for(k in 1:length(t)){
			X[j,k]<-sum(t[k]^(0:(length(Ps[[i]][[j]])-1))*coef)
		}
		tmpt<-intersect(which(t>=(j-1)) , which(t<=j))
		X2[tmpt]<-X[j,tmpt]
	}
	matplot(t,t(X),type="l")
	plot(t,X2,type="l")
	abline(h=0)
	Sys.sleep(1)
}
library(animation)
saveGIF({
for(i in 1:length(Ps)){
	t<-seq(from=0,to=i,length=100)
	X<-matrix(0,length(Ps[[i]]),length(t))
	X2<-rep(0,length(t))
	for(j in 1:length(Ps[[i]])){
		coef<-coef(Ps[[i]][[j]])
		for(k in 1:length(t)){
			X[j,k]<-sum(t[k]^(0:(length(Ps[[i]][[j]])-1))*coef)
		}
		tmpt<-intersect(which(t>=(j-1)) , which(t<=j))
		X2[tmpt]<-X[j,tmpt]
	}
	#matplot(t,t(X),type="l")
	plot(t,X2,type="l")
	abline(h=0)
	Sys.sleep(1)
}
},interval=0.5)