最尤推定した母数のもとでの確率密度曲線をヒストグラムに描き足す

えっと，リクエストがあったのでここに書く。（SPSSならそれっぽい曲線もつけてくれるのにRはそんなこともできないのか？といわれた）

まずはこんなデータがあるとしよう。形状母数が3，尺度母数1のガンマ分布にしたがう300個。

set.seed(0)
dat<-rgamma(300,3,1)

これをヒストグラムにする。

hist(dat,col="lightblue",main="")

まずはガンマ分布を仮定して，母数を最尤推定する。

library(MASS)
fitdistr(dat,densfun="gamma")

ま，3.01と1.03とかが推定されたと。
この母数のもとでの関数をヒストグラムに描き足したい。

まずはhist関数に，freq=Fを入れる。上に少し余裕を持たせるためにylimを指定しよう。

hist(dat,col="lightblue",main="",freq=F,ylim=c(0,.4))

すると縦軸が密度になる。

この上に関数を描けばいい。

たとえば，

x<-seq(0,9,.001)
lines(x,dgamma(x,3.01,1.03),col="blue")

はい。上手く当てはまっている感ある。

これで，別の関数も試してみたらいいかもだ。ワイブルと正規分布いこう。

#ここまでの作図
hist(dat,col="lightblue",main="",freq=F,ylim=c(0,.4))
x<-seq(0,9,.001)
lines(x,dgamma(x,3.01,1.03),col="blue")

#ワイブルからパラメーターを取る
w.shape<-fitdistr(dat,densfun="weibull")$estimate[[1]]
w.scale<-fitdistr(dat,densfun="weibull")$estimate[[2]]

#ワイブルを描き足す
lines(x,dweibull(x,w.shape,w.scale),col="red")

#正規分布からパラメーターを取る
mu<-fitdistr(dat,densfun="normal")$estimate[[1]]
sigma<-fitdistr(dat,densfun="normal")$estimate[[2]]

#正規分布を描き足す
lines(x,dnorm(x,mu,sigma),col="green")

緑がダメ，みたいな感じがわかる。

ちなみにAICを比較してみたり…

aics<-numeric(3)
aics[1]<-AIC(fitdistr(dat,densfun="gamma"))
aics[2]<-AIC(fitdistr(dat,densfun="weibull"))
aics[3]<-AIC(fitdistr(dat,densfun="normal"))
aics

ほい。