- n個の変数が0-1の値をとるとき、その和の最小値は0、最大値はn。その和の確率密度分布を引いてみよう
- これは、前の記事に記載した通り、ごとに区切る必要がある
- なる整数Kがあるとき
- とすることで、すべてのについて、式が漸化式で立てられる
- これをRで実行するには、パッケージpolynomを用いると、その中には、多項式の係数で多項式を作ったり係数を取り出したりする関数、微分・積分をする関数があるので、定積分を伴ったり、変数変換に対応したりする処理が可能となる
- では水平、では三角屋根、で丸みを帯びて、どんどんベル型になっていく
- において、構成3関数のすべてを描けば次のようになる。での移行具合もわかる
> 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)
Ps<-list()
n<-10
Ps[[1]]<-list()
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]])
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]])
coef1<-coef(tmp1)
coef1.2<--coef(tmp1.2)
print(coef1.2)
coef2<-coef(tmp2)
print(coef2)
coef<-coef1.2+coef2
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]
}
plot(t,X2,type="l")
abline(h=0)
Sys.sleep(1)
}
},interval=0.5)