如何将高斯曲线拟合到这些数据？ [复制]

Question

我是 R 的新手，我正试图得到一条曲线来拟合这个散布数据，这给了我一个高斯曲线。我真的很感激任何帮助。 数据：

library(tidyverse)
MK20 <- tribble(~X.Intensity,    ~Average,
             0.400,  0.0000000,
             0.463,  0.0000000,
             0.536,  0.000000,
             0.621,  0.0000000,
             0.719,  0.0000000,
             0.833,  0.0000000,
             0.965,  0.0000000,
             1.120,  0.0000000,
             1.290,  0.0000000,
             1.500,  0.0000000,
             1.740,  0.0000000,
             2.010,  0.0000000,
             2.330,  0.0000000,
             2.700,  0.0000000,
             3.120,  0.0000000,
             3.620,  0.0000000,
             4.190,  0.0000000,
             4.850,  0.0000000,
             5.610,  0.0000000,
             6.500,  0.0000000,
             7.530,  0.0000000,
             8.720,  0.0000000,
             10.100,  0.0000000,
             11.700,  0.0000000,
             13.500,  0.0000000,
             15.700,  0.0000000,
             18.200,  0.0000000,
             21.000,  0.0000000,
             24.400,  0.0000000,
             28.200,  0.0000000,
             32.700,  0.0000000,
             37.800,  0.0000000,
             43.800,  0.7023333,
             50.700,  3.3700000,
             58.800,  7.3933333,
             68.100, 11.4666667,
             78.800, 14.3666667,
             91.300, 15.4000000,
             106.000, 14.5000000,
             122.000, 12.0000000,
             142.000,  8.6433333,
             164.000,  5.2200000,
             190.000,  2.4500000,
             220.000,  0.7580000,
             255.000,  0.1306667,
             295.000,  0.0000000,
             342.000,  0.0000000,
             396.000,  0.0000000,
             459.000,  0.0000000,
             531.000,  0.0000000,
             615.000,  0.0000000,
             712.000,  0.0000000,
             825.000,  0.0000000,
             955.000,  0.0000000,
             1110.000,  0.0000000,
             1280.000,  0.0000000,
             1480.000,  0.0000000,
             1720.000,  0.0000000,
             1990.000,  0.0000000,
             2300.000,  0.0000000,
             2670.000,  0.0000000,
             3090.000,  0.0000000,
             3580.000,  0.0000000,
             4150.000,  0.0000000,
             4800.000,  0.0000000,
             5560.000,  0.0000000,
             6440.000,  0.0000000,
             7460.000,  0.0000000,
             8630.000,  0.0000000)

我用来绘图的代码是：

plot(log10(MK20$X.Intensity), MK20$Average, col=1, xlim=c(-0.5,4), ylim=c(0,20), xlab="Log(Average diameter)", ylab="Intensity", xaxt='n')

我正在使用 minor.tick.axis 函数在对数 x 轴上添加次要刻度。我想向该数据添加一条高斯曲线（最适合）。我尝试在绘图上添加type='l'，但曲线并不平滑，我不希望曲线必须触及每个数据点，而是最适合的曲线。

如果解决方案很简单，但我无法弄清楚，我很抱歉。

Answer 1

接触每个点的曲线肯定会最适合您的数据。 :)

除此之外，您可以尝试包含平滑曲线，例如

plot(log10(MK20$X.Intensity), MK20$Average, col=1, xlim=c(-0.5,4), ylim=c(0,20), 
     xlab="Log(Average diameter)", ylab="Intensity", xaxt='n', type='n')
lines(lowess(MK20$Average ~ log10(MK20$X.Intensity), f=0.3))

您可以在（0 和 1）之间改变 f= 参数以更改平滑级别。

这是 f=0.3 时的输出。

Answer 2

在这种情况下，我们不能使用通常的fitdistr 方法来拟合正态分布，因为我们没有原始数据。看起来“平均”列是某种类型的密度估计。如果它是 pdf，那么它应该集成到 1 但它没有。

f <- approxfun(x = log10(MK20$X.Intensity), y= MK20$Average)
integrate(f, lower = log10(0.4), upper = log10(8630))

#6.142134 with absolute error < 0.00043

因此我们可以通过将其缩小约 6.14 来将其转换为 pdf，然后尝试找到与该 pdf 匹配的均值和标准差。

这是对简单高斯拟合的首次尝试。首先，我选择了平均值 2（通过查看密度最大的位置）、k = 6.14（积分值）的比例因子，然后使用 sd 直到有一个合理的拟合。

m=2
s=0.15
k=6.14
x_seq = seq(1,3,length.out = 100)
df <- tibble(x_seq = x_seq, dens = dnorm(x_seq, m, s))


MK20 %>% 
  mutate(log_intensity = log10(X.Intensity)) %>% 
  ggplot(aes(log_intensity, Average/k)) +
  geom_point() +
  geom_line(data = df, aes(x_seq, dens)) 

接下来我使用 optimx 通过最小化拟合和数据之间的平方和来拟合 3 个参数（k = 比例因子、m = 平均值、s = 标准偏差）。

目标函数（拟合和数据之间差异的平方和）

f <- function(x) {
  k = x[1]
  m = x[2]
  s = x[3]
  MK20 %>% 
  mutate(log_intensity = log10(X.Intensity)) %>%
  mutate(fit = dnorm(log_intensity, m, s)) %>% 
  summarise(sum((fit - Average/k)^2)) %>% pull
}

使用 optimx 查找参数（最小平方和） 参数的初始值取自眼睛拟合。

library(optimx)    
optimx(par = c(6.14, 2, 0.15), fn = f )

#k = 6.294696 m = 1.971488 s= 0.1583936 

让我们用拟合的参数重新绘制

# points for a gaussian
x_seq = seq(1,3,length.out = 100) 
df <- tibble(x_seq = x_seq, dens = dnorm(x_seq, m, s))


MK20 %>% 
  mutate(log_intensity = log10(X.Intensity)) %>% 
  ggplot(aes(log_intensity, Average/k)) +
  geom_point() +
  geom_line(data = df, aes(x_seq, dens))