从flexmix对象（R）

Question

我一些数据拟合到两个高斯的混合物分布在flexmix预测：

data("NPreg", package = "flexmix") 
mod <- flexmix(yn ~ x, data = NPreg, k = 2, 
      model = list(FLXMRglm(yn ~ x, family= "gaussian"), 
         FLXMRglm(yn ~ x, family = "gaussian"))) 

模型拟合如下：

> mod 

Call: 
flexmix(formula = yn ~ x, data = NPreg, k = 2, model = list(FLXMRglm(yn ~ x, family = "gaussian"), 
    FLXMRglm(yn ~ x, family = "gaussian"))) 

Cluster sizes: 
    1 2 
74 126 

convergence after 31 iterations 

但是，如何我实际上从这个模型预测？

当我做

pred <- predict(mod, NPreg) 

我得到的预测名单从两个部件

为了得到一个预测，我必须在簇大小像这样加？

single <- (74/200)* pred$Comp.1[,1] + (126/200)*pred$Comp.2[,2] 

Answer 1

我以下列方式使用flexmix进行预测：

pred = predict(mod, NPreg) 
clust = clusters(mod,NPreg) 
result = cbind(NPreg,data.frame(pred),data.frame(clust)) 
plot(result$yn,col = c("red","blue")[result$clust],pch = 16,ylab = "yn") 

而混淆矩阵：

table(result$class,result$clust) 

为了获得预测值yn，我选择数据点所属的集群的组件值。

for(i in 1:nrow(result)){ 
    result$pred_model1[i] = result[,paste0("Comp.",result$clust[i],".1")][i] 
    result$pred_model2[i] = result[,paste0("Comp.",result$clust[i],".2")][i] 
} 

实际VS预测结果表明，该配合（这里只增加其中一个作为您的两种机型都是一样的，你可以使用pred_model2第二模型）。

qplot(result$yn, result$pred_model1,xlab="Actual",ylab="Predicted") + geom_abline() 

RMSE = sqrt(mean((result$yn-result$pred_model1)^2)) 

给出了一个根均的5.54方误差。

这个答案是基于许多SO答案，我通过flexmix工作时阅读。它适用于我的问题。

您可能也对可视化两个分布感兴趣。我的模型如下，它显示了一些重叠，因为组件的比例并不接近1。

Call: 
flexmix(formula = yn ~ x, data = NPreg, k = 2, 
model = list(FLXMRglm(yn ~ x, family = "gaussian"), 
      FLXMRglm(yn ~ x, family = "gaussian"))) 

     prior size post>0 ratio 
Comp.1 0.481 102 129 0.791 
Comp.2 0.519 98 171 0.573 

'log Lik.' -1312.127 (df=13) 
AIC: 2650.255 BIC: 2693.133 

我还生成一个密度分布与直方图visulaize这两个组件。这受到来自betareg维护者的SO answer的启发。

a = subset(result, clust == 1) 
b = subset(result, clust == 2) 
hist(a$yn, col = hcl(0, 50, 80), main = "",xlab = "", freq = FALSE, ylim = c(0,0.06)) 
hist(b$yn, col = hcl(240, 50, 80), add = TRUE,main = "", xlab = "", freq = FALSE, ylim = c(0,0.06)) 
ys = seq(0, 50, by = 0.1) 
lines(ys, dnorm(ys, mean = mean(a$yn), sd = sd(a$yn)), col = hcl(0, 80, 50), lwd = 2) 
lines(ys, dnorm(ys, mean = mean(b$yn), sd = sd(b$yn)), col = hcl(240, 80, 50), lwd = 2) 

# Joint Histogram 
p <- prior(mod) 
hist(result$yn, freq = FALSE,main = "", xlab = "",ylim = c(0,0.06)) 
lines(ys, p[1] * dnorm(ys, mean = mean(a$yn), sd = sd(a$yn)) + 
     p[2] * dnorm(ys, mean = mean(b$yn), sd = sd(b$yn))) 

Answer 2

你可以传递一个额外的参数，以您的通话预测。

pred <- predict(mod, NPreg, aggregate = TRUE)[[1]][,1]