p<-c(0.6,0.4)
q<-c(0.5,0.5)
z<-c(0.7,0.3)
p*log(p/(p+q))
sum(p*log(2*p/(p+z))+z*log(2*z/(p+z)))

#表3の計算
x<-c(120,80,60,55,20)
y<-c(24,28,75,80,42)
xy<-c(20,20,20,20,20)
round(log2((xy/10000)/((x/10000)*(y/10000))),3)


#データの準備
bs<-matrix(0,10,7)
bs[,1]<-1:10
bs[,2]<-c(53,1988,3538,1927,999,409,189,53,34,10)
bs[,3]<-bs[,2]/sum(bs[,2])
Su<-sum(bs[,2])

#ポアソンモデル
M1mean<-sum(bs[,1]*bs[,3])
ramuda<-M1mean-1
fact <- function(n) {if(n<=1) 1 else prod(1:n)}; #階乗を求める
for(i in 1:10){bs[i,5]<-exp(-ramuda)*ramuda^(i-1)/fact(i-1)}
#ポアソン関数dpoisを用いてもよい
bs[,5]<-dpois(bs[,1]-1, ramuda)
bs[,4]<-round(Su*bs[,5])

#対数正規モデル
x<-bs[,1]
xln<-log(x)
xlmean<-sum(xln*bs[,3])
xlvar<-sum(xln^2*bs[,3])-xlmean^2
bs[,7]<-(2*pi*xlvar*x^2)^(-1/2)*exp(-(xln-xlmean)^2/(2*xlvar))
#対数正規関数dlnormを用いてもよい
bs[,7]<-dlnorm(bs[,1],xlmean,sqrt(xlvar))
bs[,6]<-round(Su*bs[,7])

round(bs,3)

#ポアソンモデルのAIC
lm1<-sum(bs[,2]*log(bs[,5]))
-2*lm1+2

#対数正規モデルのAIC
lm2<-sum(bs[,2]*log(bs[,7]))
-2*lm2+2*2

par(mar=c(2,4,1,1))
matplot(bs[,c(3,5,7)],type="b",ylab="相対度数", pch = 1:3, col =c(1,2,4))
legend(6, max(bs[,3]), c("実測値","ポアソンモデル","対数正規モデル"),col=c(1,2,4),lty=1:3,pch = 1:3)