最大エントロピー法（メディア情報処理の基礎数理 尾関和彦先生）

packageを使用しないで最大エントロピー法に挑んでみました。
この方法は、エネルギー量とそのエネルギー量を持つ粒子の数に対して、それらの総エネルギー量が決められている場合に、各エネルギー量をもつ粒子の頻度を最適化するパラメータ(ボルツマン-ギブス分布)を決定することでエントロピーを最大化する方法です。ここでは最後にエントロピーの区間推定を行っています。

#最大エントロピー法

library(dplyr)

classes=100

N_vec=sample(classes)
E_vec=sample(classes)
N=sum(N_vec);E=sum(E_vec*N_vec)


if(round((1/N)*(sum(log(c(1:N)))-sum(log((factorial(N_vec))))))==round(-sum((N_vec/N)*log(N_vec/N)))){
  print("Result is right!")
}

ave_energy=sum(E_vec*(N_vec/N))

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=0
for(j in 1:nrow(newton_raphson)){
f<-sum((E_vec-ave_energy)*exp(-value*E_vec))
df<-sum((E_vec-ave_energy)*exp(-value*E_vec)*(-E_vec))

  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value

  newton_raphson$value[j]=value
}

solution=newton_raphson$value[nrow(newton_raphson)]

prob=array(0,dim=c(1,length(E_vec)))

for(j in 1:length(prob)){
  prob[1,j]=exp(-solution*E_vec[j])/sum(exp(-solution*E_vec))
}

if(signif(sum(prob),10)==1){
  print("Result is right!")
}


print("probabilities are")
print(prob)

entropy=-signif(sum(prob*log(prob)),10)

#interval estimation for the entropy values
n=length(E_vec);m=1
diff_H=(1/N)*qchisq(1-0.05,n-m-1)
print(paste0("Interval is [",entropy-diff_H," , ",entropy,"]"))