R 中高斯混合模型的 Matlab 'fit' 等价物？
Equivalent of Matlab's 'fit' for Gaussian mixture models in R?
我有一些时间序列数据，如下所示：
x <- c(0.5833, 0.95041, 1.722, 3.1928, 3.941, 5.1202, 6.2125, 5.8828,
4.3406, 5.1353, 3.8468, 4.233, 5.8468, 6.1872, 6.1245, 7.6262,
8.6887, 7.7549, 6.9805, 4.3217, 3.0347, 2.4026, 1.9317, 1.7305,
1.665, 1.5655, 1.3758, 1.5472, 1.7839, 1.951, 1.864, 1.6638,
1.5624, 1.4922, 0.9406, 0.84512, 0.48423, 0.3919, 0.30773, 0.29264,
0.19015, 0.13312, 0.25226, 0.29403, 0.23901, 0.000213074755156413,
5.96565965097398e-05, 0.086874, 0.000926808687858284, 0.000904641782399267,
0.000513042259030044, 0.40736, 4.53928073402494e-05, 0.000765719624469057,
0.000717419263673946)
我想使用 1 到 5 个高斯分布的混合来拟合此数据的曲线。在 Matlab 中，我可以执行以下操作：
fits{1} = fit(1:length(x),x,fittype('gauss1'));
fits{2} = fit(1:length(x),x,fittype('gauss2'));
fits{3} = fit(1:length(x),x,fittype('gauss3'));
...等等。
在 R 中，我很难找到类似的方法。
dat <- data.frame(time = 1:length(x), x = x)
fits[[1]] <- Mclust(dat, G = 1)
fits[[2]] <- Mclust(dat, G = 2)
fits[[3]] <- Mclust(dat, G = 3)
...但这似乎并没有真正做同样的事情。例如，我不确定如何使用 Mclust 解决方案计算拟合曲线与原始数据之间的 R^2。
在 base R 中是否有更简单的替代方法来使用混合高斯拟合曲线？
函数
使用下面给出的代码，如果幸运地找到了好的初始参数，您应该能够对数据进行高斯曲线拟合。
在函数fit_gauss中，目标是y ~ fit_gauss(x)，要使用的高斯数由参数初始值的长度决定: a, b, d 所有这些都应该等长
我已经演示了将 OP 数据曲线拟合到三个高斯分布。
指定初始值
这几乎是我用 nls 完成的大部分工作（感谢 OP）。所以，我不太确定select初始值的最佳方法是什么。当然，它们取决于峰的高度 (a)、峰周围 x 的平均值和标准差（b 和 d）。
一个选项是对于给定数量的高斯分布，尝试使用多个起始值，并根据残差标准误差 fit$sigma.
找到最适合的那个
为了找到初始参数我费了一番功夫，但我敢说参数和
三个高斯模型的情节看起来很稳固。
将一、二和高斯拟合到示例数据
ind <- 1 : length(x)
# plot original data
plot(ind, x, pch = 21, bg = "blue")
# Gaussian fit
fit_gauss <- function(y, x, a, b, d) {
p_model <- function(x, a, b, d) {
rowSums(sapply(1:length(a),
function(i) a[i] * exp(-((x - b[i])/d[i])^2)))
}
fit <- nls(y ~ p_model(x, a, b, d),
start = list(a=a, b = b, d = d),
trace = FALSE,
control = list(warnOnly = TRUE, minFactor = 1/2048))
fit
}
单高斯
g1 <- fit_gauss(y = x, x = ind, a=1, b = mean(ind), d = sd(ind))
lines(ind, predict(g1), lwd = 2, col = "green")
两个高斯
g2 <- fit_gauss(y = x, x = ind, a = c(coef(g1)[1], 1),
b = c(coef(g1)[2], 30),
d = c(coef(g1)[1], 2))
lines(ind, predict(g2), lwd = 2, col = "red")
三个高斯
g3 <- fit_gauss(y = x, x = ind, a=c(5, 4, 4),
b = c(12, 17, 11), d = c(13, 2, 2))
lines(ind, predict(g3), lwd = 2, col = "black")
三个高斯拟合的总结
summary(g3)
# Formula: x ~ p_model(ind, a, b, d)
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# a1 5.9307 0.5588 10.613 5.93e-14 ***
# a2 3.5689 0.7098 5.028 8.00e-06 ***
# a3 -2.2066 0.8901 -2.479 0.016894 *
# b1 12.9545 0.5289 24.495 < 2e-16 ***
# b2 17.4709 0.2708 64.516 < 2e-16 ***
# b3 11.3839 0.3116 36.538 < 2e-16 ***
# d1 11.4351 0.8568 13.347 < 2e-16 ***
# d2 1.8893 0.4897 3.858 0.000355 ***
# d3 1.0848 0.6309 1.719 0.092285 .
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.7476 on 46 degrees of freedom
#
# Number of iterations to convergence: 34
# Achieved convergence tolerance: 8.116e-06
我有一些时间序列数据，如下所示：
x <- c(0.5833, 0.95041, 1.722, 3.1928, 3.941, 5.1202, 6.2125, 5.8828,
4.3406, 5.1353, 3.8468, 4.233, 5.8468, 6.1872, 6.1245, 7.6262,
8.6887, 7.7549, 6.9805, 4.3217, 3.0347, 2.4026, 1.9317, 1.7305,
1.665, 1.5655, 1.3758, 1.5472, 1.7839, 1.951, 1.864, 1.6638,
1.5624, 1.4922, 0.9406, 0.84512, 0.48423, 0.3919, 0.30773, 0.29264,
0.19015, 0.13312, 0.25226, 0.29403, 0.23901, 0.000213074755156413,
5.96565965097398e-05, 0.086874, 0.000926808687858284, 0.000904641782399267,
0.000513042259030044, 0.40736, 4.53928073402494e-05, 0.000765719624469057,
0.000717419263673946)
我想使用 1 到 5 个高斯分布的混合来拟合此数据的曲线。在 Matlab 中，我可以执行以下操作：
fits{1} = fit(1:length(x),x,fittype('gauss1'));
fits{2} = fit(1:length(x),x,fittype('gauss2'));
fits{3} = fit(1:length(x),x,fittype('gauss3'));
...等等。
在 R 中，我很难找到类似的方法。
dat <- data.frame(time = 1:length(x), x = x)
fits[[1]] <- Mclust(dat, G = 1)
fits[[2]] <- Mclust(dat, G = 2)
fits[[3]] <- Mclust(dat, G = 3)
...但这似乎并没有真正做同样的事情。例如，我不确定如何使用 Mclust 解决方案计算拟合曲线与原始数据之间的 R^2。
在 base R 中是否有更简单的替代方法来使用混合高斯拟合曲线？
函数
使用下面给出的代码，如果幸运地找到了好的初始参数，您应该能够对数据进行高斯曲线拟合。
在函数fit_gauss中，目标是y ~ fit_gauss(x)，要使用的高斯数由参数初始值的长度决定: a, b, d 所有这些都应该等长
我已经演示了将 OP 数据曲线拟合到三个高斯分布。
指定初始值
这几乎是我用 nls 完成的大部分工作（感谢 OP）。所以，我不太确定select初始值的最佳方法是什么。当然，它们取决于峰的高度 (a)、峰周围 x 的平均值和标准差（b 和 d）。
一个选项是对于给定数量的高斯分布，尝试使用多个起始值，并根据残差标准误差 fit$sigma.
为了找到初始参数我费了一番功夫，但我敢说参数和 三个高斯模型的情节看起来很稳固。
将一、二和高斯拟合到示例数据
ind <- 1 : length(x)
# plot original data
plot(ind, x, pch = 21, bg = "blue")
# Gaussian fit
fit_gauss <- function(y, x, a, b, d) {
p_model <- function(x, a, b, d) {
rowSums(sapply(1:length(a),
function(i) a[i] * exp(-((x - b[i])/d[i])^2)))
}
fit <- nls(y ~ p_model(x, a, b, d),
start = list(a=a, b = b, d = d),
trace = FALSE,
control = list(warnOnly = TRUE, minFactor = 1/2048))
fit
}
单高斯
g1 <- fit_gauss(y = x, x = ind, a=1, b = mean(ind), d = sd(ind))
lines(ind, predict(g1), lwd = 2, col = "green")
两个高斯
g2 <- fit_gauss(y = x, x = ind, a = c(coef(g1)[1], 1),
b = c(coef(g1)[2], 30),
d = c(coef(g1)[1], 2))
lines(ind, predict(g2), lwd = 2, col = "red")
三个高斯
g3 <- fit_gauss(y = x, x = ind, a=c(5, 4, 4),
b = c(12, 17, 11), d = c(13, 2, 2))
lines(ind, predict(g3), lwd = 2, col = "black")
三个高斯拟合的总结
summary(g3)
# Formula: x ~ p_model(ind, a, b, d)
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# a1 5.9307 0.5588 10.613 5.93e-14 ***
# a2 3.5689 0.7098 5.028 8.00e-06 ***
# a3 -2.2066 0.8901 -2.479 0.016894 *
# b1 12.9545 0.5289 24.495 < 2e-16 ***
# b2 17.4709 0.2708 64.516 < 2e-16 ***
# b3 11.3839 0.3116 36.538 < 2e-16 ***
# d1 11.4351 0.8568 13.347 < 2e-16 ***
# d2 1.8893 0.4897 3.858 0.000355 ***
# d3 1.0848 0.6309 1.719 0.092285 .
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.7476 on 46 degrees of freedom
#
# Number of iterations to convergence: 34
# Achieved convergence tolerance: 8.116e-06