我对大量2D数据点进行高斯拟合，以测试它们是否(至少大致)遵循高斯曲线。然而，我意识到，当我将y维度标准化为其最大值时，我得到的结果与不标准化时不同。这里有一个例子，两个向量分别包含y和x数据：(我知道这些数据实际上没有遵循高斯分布，这是我想通过这个测试得到的答案)
y = c(0, 4551.2783203125, 17204.81640625, 18554.16015625, 36956.65625,
37979.30859375, 41470.2265625, 61240.8359375, 106161.9609375,
87546.9375, 100634.7109375, 186276.765625, 221888.4375, 248687.84375,
252164.890625, 343520.21875, 280060.625, 442446.15625, 424090.8125,
442001.40625, 375247.46875, 458099.25, 607202.25, 452370.15625,
361559.75, 481408.28125, 323888.53125, 500188.5625, 409107.1875,
380734.96875, 312657.8125, 507054.5625, 423043.59375, 465331.0625,
567257.375, 365184.125, 685140.5625, 459672.34375, 845416.5,
455453.3125, 521206.03125, 613202.625, 477479.625, 579983.3125,
760272.375, 514784.03125, 652986.5625, 413048.75, 533935.9375,
686196.9375, 767804.9375, 665105.9375, 525717.6875, 631246.3125,
545236.5, 621491.8125, 891600.8125, 870612.625, 570333.8125,
562794.125, 680060.5625, 492349.90625, 492722.09375, 439363.0625,
793801.1875, 629333.25, 697334.625, 441465.53125, 389707.71875,
415848.25, 418587.0625, 527000.0625, 472959.34375, 508528.96875,
679527.25, 655779.1875, 498504.0625, 376315.875, 429004.71875,
328198.3125, 441643.53125, 406398.90625, 453388.03125, 349939.59375,
348471.34375, 629093.5, 325516.03125, 285678.53125, 298847.8125,
375509.875, 392465.15625, 367188.46875, 419100.9375, 311449.59375,
260993.9375, 290824.6875, 236561.0625, 265879.5625, 379404.25,
311229.6875, 307151.96875, 379062.46875, 280802.5, 457000.78125,
260519.28125, 349163.875, 291823.40625, 260145.359375, 255300.3125,
334640.84375, 306605.34375, 425454.46875, 257264.5, 220655.140625,
223242.28125, 327066.5625, 300993.46875, 234662.15625, 262443.5,
230376.796875, 227217.25, 236800.0625, 214616.421875, 260351.578125,
226784.921875, 282392.28125, 200394.671875, 256462.5, 182574.59375,
161353.78125, 165745.484375, 199003.859375, 175615.828125, 167471.859375,
204727.078125, 207417.140625, 202296.46875, 183818.984375, 247653.640625,
163297.9375, 171750.921875, 161632.78125, 201405.53125, 149500,
123130.8125, 144252.359375, 173929.453125, 164804.953125, 144984.1875,
140006.96875, 126611.0859375, 131078.140625, 222015.546875, 124387.859375,
112429.4453125, 185341.9375, 83172.6640625, 142822.765625, 131457.234375,
122272.4921875, 99884.0546875, 128589.4765625, 110691.6328125
)
x = c(7.99422121047974, 7.99860048294067, 8.00297927856445, 8.00735855102539,
8.0117359161377, 8.01611709594727, 8.02049446105957, 8.02487277984619,
8.02925491333008, 8.03363132476807, 8.038010597229, 8.04239177703857,
8.04676914215088, 8.05114841461182, 8.05552768707275, 8.05990695953369,
8.06428337097167, 8.06866455078125, 8.07304191589355, 8.07742118835449,
8.08180236816406, 8.08617973327637, 8.0905590057373, 8.09494018554688,
8.09931659698486, 8.1036958694458, 8.10807704925537, 8.11245441436768,
8.11683368682861, 8.12121295928955, 8.12559032440186, 8.12996959686279,
8.13435077667236, 8.13872718811035, 8.14310646057129, 8.14748764038085,
8.15186500549316, 8.1562442779541, 8.16062545776367, 8.16500282287598,
8.16937923431396, 8.17376136779785, 8.17813777923583, 8.18251705169678,
8.18689823150635, 8.19127559661865, 8.19565486907958, 8.20003604888916,
8.20441341400146, 8.20879173278809, 8.21317386627197, 8.21755027770996,
8.2219295501709, 8.22630882263184, 8.23068618774414, 8.23506546020508,
8.23944664001465, 8.24382305145264, 8.24820232391357, 8.25258350372314,
8.25696086883544, 8.26134014129639, 8.26572132110596, 8.27009868621826,
8.27447795867919, 8.27885723114014, 8.28323554992676, 8.28761291503906,
8.29199409484863, 8.29637145996094, 8.30075073242188, 8.30513191223145,
8.30950927734375, 8.31388854980469, 8.31826972961426, 8.32264614105225,
8.32702541351318, 8.33140659332275, 8.33578395843506, 8.340163230896,
8.34454250335693, 8.34891986846924, 8.35329818725585, 8.35768032073975,
8.36205673217773, 8.36643600463867, 8.37081718444824, 8.37519454956055,
8.37957382202148, 8.38395500183105, 8.38833236694335, 8.39271068572998,
8.39708995819092, 8.40146923065186, 8.40584659576416, 8.41022777557373,
8.41460514068604, 8.41898441314697, 8.42336559295654, 8.42774295806885,
8.43212127685547, 8.43650245666504, 8.44087982177734, 8.44525909423828,
8.44963836669922, 8.45401763916016, 8.45839500427246, 8.46277618408203,
8.46715259552001, 8.47153186798096, 8.47591304779053, 8.48029041290283,
8.48466968536376, 8.48905086517333, 8.49342823028564, 8.49780750274658,
8.50218677520751, 8.50656509399414, 8.51094245910645, 8.51532363891602,
8.51970100402832, 8.52408027648926, 8.52846145629883, 8.53283882141113,
8.53721714019775, 8.54159927368164, 8.54597568511963, 8.55035495758057,
8.5547342300415, 8.55911350250244, 8.56349086761475, 8.56787204742432,
8.5722484588623, 8.57662773132324, 8.58100891113281, 8.58538627624512,
8.58976554870605, 8.59414672851563, 8.59852409362792, 8.60290336608887,
8.60728454589844, 8.61166095733643, 8.61604022979736, 8.6204195022583,
8.62479686737061, 8.62917613983154, 8.63355731964111, 8.63793468475342,
8.64231395721436, 8.64669513702393, 8.65107154846191, 8.65545082092285,
8.65983200073242, 8.66420936584473, 8.66858863830566, 8.6729679107666,
8.67734718322754, 8.68172359466553, 8.68610572814941, 8.6904821395874,
8.69486141204833, 8.69924259185791, 8.70361995697021)
这些函数用于创建高斯轮廓、计算均方根误差和优化高斯轮廓的三个参数：
GaussCurve <- function(rt.Vector,par) #generate Gauss profile
{
m <- par[1]
sd <- par[2]
k <- par[3]
Fct.V <- k * exp(-0.5 * ((rt.Vector - m)/sd)^2)
Fct.V
}
RMSE <- function(par) #calculate root mean square error
{
Fct.V <- GaussCurve(rt,par)
sqrt(sum((signal - Fct.V)^2)/length(signal))
}
signal <- y
rt <- x
#optimization
result <- optim(c(rt[which.max(signal)], unname(quantile(rt)[4]-quantile(rt)[2]), max(signal)),
lower = c(min(rt), -Inf, 0.1*max(signal)),
upper = c(max(rt), Inf, max(signal)),
RMSE, method="L-BFGS-B", control=list(factr=1e7))
result
#plot of result
plot(rt,signal,xlab="RT/min",ylab="I")
lines(seq(min(rt),max(rt),length=1000),GaussCurve(seq(min(rt),max(rt),length=1000),result$par),col=2)
当我使用这段代码时，我现在得到了一些结果，当然这与我对这些数据的期望不太吻合。然而，当我用运行相同的代码时
signal <- y/max(y)
合身看起来完全不同。为什么？我知道我不能很好地拟合这些数据，但从我对优化过程的理解来看，我希望在这两种情况下都能得到同样的拟合。在其中一个案例中，我被困在当地的最低限度吗？在这种情况下，这不会是一个问题，因为无论如何，拟合都很糟糕，但我想确保在拟合有效的情况下不会发生这种情况。那么，我可以调整代码中的任何内容来避免这种情况吗？
根据我的经验提出的一些建议：
- 找到更好的目标函数
- 试着让你的函数尽可能"平滑"(关于参数)。添加错误术语可以在此处提供帮助
- 比起RMSE，我更喜欢基于可能性的方法，我很少发现后者有用
- 试着让你的函数尽可能"平滑"(关于参数)。我会添加一个"error"参数
- 根据需要转换参数，而不是将限制传递给
optim
结合这些，我会让你的目标函数做：
optfn <- function(par) {
Fct.V <- GaussCurve(rt,par)
se2 <- exp(par[4])
-sum(dnorm(signal, Fct.V, se2, log=TRUE))
}
result <- optim(c(
rt[which.max(signal)],
quantile(rt)[[4]] - quantile(rt)[[2]],
max(signal) * 0.7,
log(max(signal) * 0.7)
), optfn)
这有点帮助，但没有多大帮助。
蒙特卡罗方法给了我更好的结果，我得到了8.248 <= par[0] <= 8.276和0.170 <= par[1] <= 0.205的95%置信区间。当我进行y/max(y)变换时，我得到了8.249 <= par[0] <= 8.275和0.171 <= par[1] <= 0.204，这似乎是适当一致的。它要复杂得多，所以取决于你想花多少钱学习