7

元のデータから取得した 100 万レコードのサンプルがあります。(参考までに、ほぼ同様の分布を生成する可能性があるこのダミー データを使用することができます。

b <- data.frame(matrix(rnorm(2000000, mean=c(8,17), sd=2)))
c <- b[sample(nrow(b), 1000000), ]

) 私はヒストグラムが 2 つの対数正規分布の混合であると信じており、次のコードを使用して EM アルゴリズムを使用して合計分布を適合させようとしました。

install.packages("mixtools")
lib(mixtools)
#line below returns EM output of type mixEM[] for mixture of normal distributions
c1 <- normalmixEM(c, lambda=NULL, mu=NULL, sigma=NULL) 
plot(c1, density=TRUE)

最初のプロットは対数尤度プロットで、2 番目のプロット (もう一度 Return キーを押した場合) は、次の密度曲線のようになります。

混合モデル密度曲線

前述したように、c1 は mixEM[] 型であり、plot() 関数はそれに対応できます。密度曲線を色で塗りつぶしたい。これは ggplot2() を使用して簡単に実行できますが、ggplot2() は mixEM[] 型のデータをサポートしておらず、次のメッセージがスローされます。

ggplot は、クラス mixEM のデータを処理する方法を知りません

この問題に対して私が取ることができる他のアプローチはありますか?

4

2 に答える 2

5

geom_ploygon(...)への複数回の呼び出しの代わりに使用する、わずかに異なるアプローチを次に示しますstat_function(...)。問題の 1 つstat_function(...)は、パラメーターを使用して渡される 2 次引数 (この例では mu、sigma、および lambda) をargs=list(...)美的マッピングに含めることができないため、stat_function(...)@Spacedman のソリューションのように複数の呼び出しを行う必要があることです。 .

このアプローチは、ggplot の外部で PDF を構築し、geom_polygon(...). その結果、混合物内の任意の数の分布に対して変更なしで機能します。

# ggplot mixture plot
gg.mixEM <- function(EM) {
  require(ggplot2)
  x       <- with(EM,seq(min(x),max(x),len=1000))
  pars    <- with(EM,data.frame(comp=colnames(posterior), mu, sigma,lambda))
  em.df   <- data.frame(x=rep(x,each=nrow(pars)),pars)
  em.df$y <- with(em.df,lambda*dnorm(x,mean=mu,sd=sigma))
  ggplot(data.frame(x=EM$x),aes(x,y=..density..)) + 
    geom_histogram(fill=NA,color="black")+
    geom_polygon(data=em.df,aes(x,y,fill=comp),color="grey50", alpha=0.5)+
    scale_fill_discrete("Component\nMeans",labels=format(em.df$mu,digits=3))+
    theme_bw()
}

library(mixtools)
# two components
set.seed(1)    # for reproducible example
b <- rnorm(2000000, mean=c(8,17), sd=2)
c <- b[sample(length(b), 1000000) ]
c2 <- normalmixEM(c, lambda=NULL, mu=NULL, sigma=NULL) 
gg.mixEM(c2)


# three components
set.seed(1)
b <- rnorm(2000000, mean=c(8,17,30), sd=c(2,3,5))
c <- b[sample(length(b), 1000000) ]
library(mixtools)
c3 <- normalmixEM(c, k=3, lambda=NULL, mu=NULL, sigma=NULL) 
gg.mixEM(c3)

于 2014-08-14T18:43:27.207 に答える