読者です 読者をやめる 読者になる 読者になる

草薙の研究ログ

英語教育関係。でも最近は統計(特にR)ネタが中心。

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

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

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

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

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

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

f:id:kusanagik:20170111175015p:plain

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

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))

f:id:kusanagik:20170111174956p:plain

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

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

たとえば,

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

f:id:kusanagik:20170111180126p:plain

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

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

#ここまでの作図
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")

f:id:kusanagik:20170111180744p:plain


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

ちなみに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


ほい。