ggplotでF分布と棄却域を描く自作関数
RでF分布と棄却域を描画したいとき,片所強先生が公開されている関数がたいへん便利である。
引数に"norm", "t", "chisq", "F"のいずれかを指定するだけで簡単に分布曲線が描かれます。描こうとする分布によっては、x軸の値が適当でないこともあるので、その場合は引数myself = TRUEとすることによって分布を描くための関数自体が返されます。*1
分布のプロット(正規分布、t分布、F分布、χ2分布)
片所先生の関数を使えば,例えば,df(2,100)のF分布と有意水準を5%とした場合の棄却域は,次のように書けば良い。
#※予め関数の定義を読ませておく d <- my.dist("F", myself=TRUE) #自由度を聞かれるので2と100を入力する plot(d, 0, 10) #第2,3引数は,x軸の描画範囲 my.paint(d, qf(1-0.05, 2, 100), 10, col="grey") #第2,3引数は,塗りつぶしに対応するx軸の範囲
これで十分であるが,場合によっては,もう少し丁寧に描画したいことがある。
今回は,ggplotを使った関数を書いた。
以下に関数の定義コードを示す。
F.dist <- function(df1,df2,alpha=0.05,fmax=20){ my.F <- data.frame(X=x<-seq(0,fmax,len=10000), Y=df(x,df1,df2)) #F分布の確率密度関数データフレーム p <- ggplot() p <- p + geom_path(data=my.F, aes(x=X, y=Y), size=1.0) p <- p + labs(x = "F value", y = "Probability density") #以下,棄却域 upper <- qf(1-alpha, df1, df2) #臨界値 ymin <- par("usr")[3] ymax <- par("usr")[4] pos1 <- ((ymin+ymax)/100)*20*0.70 #矢印↓のy座標上限 pos2 <- ((ymin+ymax)/100)*20 #テキストのy座標 p <- p + geom_ribbon(data=subset(my.F, X>=upper), aes(x=X, ymin=0, ymax=Y), fill="gray") p <- p + annotate("text",size=4.5,label=format(upper,digits=2,nsmall=2,scientific=F),x=upper,y=pos2, fill="white") p <- p + annotate("segment",size=1, x=upper, xend=upper, y=pos1, yend=0.0,arrow=arrow()) p <- p + scale_x_continuous(expand = c(0, 0)) p <- p + scale_y_continuous(expand = c(0, 0)) p <- p + annotate("text", x = Inf, y = Inf, label = paste("F(",df1,",",df2,")\nα=",alpha,sep=""), hjust=1.1, vjust=1.1, col="black", cex=6, fontface = "bold", alpha = 0.8) #以下,書式の調整 p <- p + theme(axis.title.x = element_text(vjust = 0), #vjustでラベル位置調整 axis.title.y = element_text(vjust = 2, angle = 90), #vjustでラベル位置調整 panel.background = element_rect(fill = "white", colour = NA), plot.background = element_rect(colour = "white"), axis.line = element_line(colour = "#000000", size = 1, linetype = "solid", lineend = "round"), axis.text = element_text(size = rel(1)), axis.text.x = element_text(vjust = 1,colour="#000000"), axis.text.y = element_text(hjust = 1,colour="#000000"), axis.ticks = element_line(size = 1,colour="#000000"), axis.ticks.length = unit(0.2, "cm"), axis.ticks.margin = unit(0.2, "cm"), panel.grid = element_blank(), text = element_text(face = "plain", colour = "black", size = 12, hjust = 0.5, vjust = 0.5, angle = 0, lineheight = 0.9)) plot(p) return(paste("F分布(",df1,",",df2,")における,α=",alpha,"の棄却域はF>=",format(upper,digits=2,nsmall=2,scientific=F),"です。",sep="")) }
この関数を使って,df(2,100)のF分布と有意水準を5%とした場合の棄却域を描くには,次のようにすれば良い。
#※予め関数の定義を読ませておく library(ggplot2) library(grid) F.dist(2,100,fmax=10) #fmaxはx軸の描画上限値(省略時20)
有意水準を変えたい場合は,引数にalphaを追加すれば良い。
F.dist(2, 100, fmax=10, alpha=0.1)
一行で描画できるため,次のようなことを理解する際に良い。
次回は,F分布に形状が近い,χ2乗分布の関数を作成する予定である。
*1:2015年12月19日閲覧