3

我将一些数据拟合到两个高斯的混合分布中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]
4

2 回答 2

7

我通过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")

NPreg 中的集群

和混淆矩阵:

table(result$class,result$clust)

NPreg 的混淆矩阵

为了获得 的预测值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]
}

实际结果与预测结果显示了拟合(此处仅添加其中一个,因为您的两个模型相同,您将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.

该答案基于我在使用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 

我还使用直方图生成密度分布,以可视化这两个组件。是受到betareg.

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

在此处输入图像描述

于 2016-02-11T20:55:26.833 回答
0

您可以将附加参数传递给您的预测调用。

pred <- predict(mod, NPreg, aggregate = TRUE)[[1]][,1]
于 2017-06-08T08:05:55.163 回答