凉凉 如果发现拟合效果一般, 或者还想要稍微微调一些参数去重新拟合怎么办?
请修改参数列表: {{mu1, 0.1}, ...} 以 {dummyVariable, initValue} 的形式来表示数据.
换用其他的函数: LinearModelFit, NoneLinearModelFit 等等
关于 MMA 拟合的更多内容
NonlinearModelFit
和 FindFit
在 MMA 内部是一样的,区别只在于参数的形式不同,但最终拟合的结果不会更好或更坏。
LinearModelFit
是真的线性拟合,你需要指定 {f1[x], f2[x], .., fk[x]}
,MMA 会帮你线性组合这几个函数,达到平方误差最小的拟合。
如何得到更好的拟合结果
一方面可以手动调初值,拟合基本上是从初值开始(默认似乎是 0 ),使用的是局部最小值的算法。因此在很多情况下,局部最小就需要调整到一个合适的初值,使得从该初值出发,诸如梯度下降等算法可以达到最小值点。
但如果你不想手动调初值,或者你想写一个稍微通用些的算法,那就需要尝试从局部最优跳出来,到达全局最优。
可以参考MMA官方的关于全局最优的教程:数值非线性全局最优化 其中也提到了一些不同算法。下面主要以退火算法 Simulated Annealing 为例。
我们先生成一个较为复杂的函数(这是在官方文档中介绍退火算法的函数),并生成 200 个样本点。
f[x_, y_] :=
20 Sin[\[Pi]/2 (x - 2 \[Pi])] +
20 Sin[\[Pi]/2 (y - 2 \[Pi])] + (x - 2 \[Pi])^2 + (y - 2 \[Pi])^2;
Plot3D[f[x, y], {x, 0, 10}, {y, 0, 10}]
data = Table[
With[{x = RandomReal[]*10, y = RandomReal[]*10}, {x, y,
f[x, y]}], {i, 0, 200}];
我们如果想获取一个函数的最小值,需要使用 NMinimize
函数,如下:
In[] = NMinimize[f[x, y], {x, y},
Method -> {"SimulatedAnnealing", "PerturbationScale" -> 7}]
Out[] = {-38.0779, {x -> 5.32216, y -> 5.32216}}
可以看到,NMinimize 找到了一个局部最优解后,就误以为自己找到了全局最优,在很远处就结束掉了。
我们需要使用退火算法,让 NMinimize 更好的寻找结果。(退火算法的详细介绍自己看wiki或者上计算物理课)
In[] = NMinimize[f[x, y], {x, y}]
Out[] = {8.0375, {x -> 1.48098, y -> 1.48098}}
那么问题来了,这只是寻找函数的最小值,和 NonlinearModelFit
有什么关系呢?有的,因为 NonlinearModelFit
内部就需要调用 Minimize 算法,找到一组参数,使均方误差最小。我们只需要告诉 MMA ,寻找最小值时,使用退火算法来寻找就可以了。
先看一下如果没有任何算法,会怎样:
nlm = NonlinearModelFit[data,
a1 Sin[k1 x + b1] +
a2 Sin[k2 y + b2] + (x - c1)^2 + (y - c2)^2, {a1, a2, k1, k2, b1,
b2, c1, c2}, {x, y}];
nlm["BestFitParameters"]
Plot3D[nlm[x, y], {x, 0, 10}, {y, 0, 10}]
{a1 -> -1.52271, a2 -> 3.39773, k1 -> 10.3748, k2 -> 8.19673,
b1 -> -52.2501, b2 -> -44.6849, c1 -> 6.26415, c2 -> 6.28012}
可以看出来,拟合似乎是把二次函数的中心位置找到了,但是正弦函数的拟合就很糟糕。
下面简单的引入退火算法,需要让 NonlinearModelFit
在内部使用 NMinimize
算法并且极小值算法为 SimulatedAnnealing
。这可以使用 Method 设置,Method 的进阶的参数形式如下
{name1, Method->{name2, ...}} 用方法和子方法
因此使用MMA代码,使用 NMinimize
方法和它的子方法 "SimulatedAnnealing":
nlm = NonlinearModelFit[data,
a1 Sin[k1 x + b1] +
a2 Sin[k2 y + b2] + (x - c1)^2 + (y - c2)^2, {a1, a2, k1, k2, b1,
b2, c1, c2}, {x, y},
Method -> {"NMinimize", Method -> "SimulatedAnnealing"}];
nlm["BestFitParameters"]
Plot3D[nlm[x, y], {x, 0, 10}, {y, 0, 10}]
{a1 -> -20.717, a2 -> -3.5992, k1 -> 1.57375, k2 -> -4.21244*10^-8,
b1 -> -0.389284, b2 -> 2.37695, c1 -> 6.2613, c2 -> 6.37348}
可以发现,我们成功地找到了一组的正弦,但是没有找到全部的最优解。这个时候可以调整更多的参数来完成:
nlm = NonlinearModelFit[data,
a1 Sin[k1 (x - b1)] +
a2 Sin[k2 (y - b2)] + (x - c1)^2 + (y - c2)^2, {a1, a2, k1, k2,
b1, b2, c1, c2}, {x, y},
Method -> {"NMinimize",
Method -> {"SimulatedAnnealing", "PerturbationScale" -> 5}}];
nlm["BestFitParameters"]
Plot3D[nlm[x, y], {x, 0, 10}, {y, 0, 10}]
{a1 -> 20., a2 -> -20., k1 -> 1.5708, k2 -> -1.5708, b1 -> -5.71682,
b2 -> 6.28319, c1 -> 6.28319, c2 -> 6.28319}
非常棒!现在所有的参数都达到预期了。
更多更复杂的函数的退火可能需要更多的退火的参数,还是看上面给出的参考文档的链接就可以了。它里面也给出了退火算法可以设置的方法
选项名 | 默认值 | 说明 |
"BoltzmannExponent" | Automatic | 概率函数的指数 |
"InitialPoints" | Automatic | 初始点集 |
"LevelIterations" | 50 | 停留在某一给定点的最大迭代次数 |
"PenaltyFunction" | Automatic | 用于约束条件以惩罚无效点的函数 |
"PerturbationScale" | 1.0 | 随机跳转的尺度 |
"PostProcess" | Automatic | 是否利用局部搜索方法进行后处理 |
"RandomSeed" | 0 | 随机数字生成器的初始值 |
"SearchPoints" | Automatic | 初始点的数目 |
"Tolerance" | 0.001 | 接受对约束违反程度的宽限 |
如果简单的退火不能达到要求,调整一下退火的具体方法就差不多可以了。