本帖最后由 Jack315 于 2023-8-28 13:47 编辑
实现搜索平均意义上拟合效果最佳区域的方法：
1. 任意给定一个搜索起点。
2. 在这个点为中心的一个区域内、大致均匀的各个点上。
分别计算代表拟合优劣的一组数值。
3. 根据这组数值建立一个被平均了的、光滑的曲面模型。
4. 用这个模型计算出梯度，这个梯度就指向
极值（最大或最小）点变化最快的方向。
5. 如果搜索区域内已经包含了最小值，迭代结束；
否则根据梯度方向，朝向最小值方向移动，继续搜索的过程。
暂且称这个搜索为全局搜索。
这个全局搜索提高了找到全局最优解的可能。
可以尝试几个不同的搜索起始点，然后再作比较和判断。
在全局搜索完成后，再用软件命令找出局部最优解。
本帖最后由 笨笨 于 2023-8-30 08:55 编辑
Ickiverar 发表于 2023-8-18 06:23
对于你 https://bbs.emath.ac.cn/thread-19045-1-1.html 帖子中加入一个新参数之后的情况，我用 kabcd 这五 ...
先生还在吗，楼主没看明白这一步，为啥引进红线部分函数：Sin，它有什么特殊之处，而不是其它函数？原理是啥？
xlist = N@Table, {x, 2/5, 1, 1/100}];
另外请问下面这些代码，前辈是现场直接手动一个一个打出输入到Mathematica中的？
还是说在其它软件里提前编辑好，然后把转化成的代码拷贝到Mathematica中后作小浮动修改？
比如我们在Mathtype里编辑好公式，然后直接拷贝到论坛，那么公式自然就转化成latex代码了，请问编辑下面代码有啥技巧？或者有其它软件作编程辅助，求告之。
grad :=
Module[{n = Length, x = x0, dlist, l = l0, y, oldy, oldd,
step = 0, nkernel = Max]], ymin = Infinity,
xmin},
oldy = f @@ x;
oldd = None;
While[True,
++step;
dlist =
Table[f @@ (x + #) - f @@ (x - #) &@
Table, {j, n}], {i, n}]/(2 d);
dlist /= Norm;
(*x-=l dlist;
y=f@@x;*)
y = Table[{i, #, f @@ #} &, {i, nkernel}];
y = SortBy[];
l *= 1.5^(y[] - 1 - Round[(nkernel - 1)/2]);
x = y[];
y = y[];
If;
If == 1,
Pause; Print[{{ymin, xmin}, {y, x}, l, step}]];
If[! (oldd === None), l *= 0.85 + 0.35 dlist.oldd];
oldd = dlist;
oldy = y;
]
];
本帖最后由 Jack315 于 2023-8-30 07:12 编辑
先说个直觉，这个补偿函数可能有个天花板。
直觉来自图形，不是证明。先用下列代码看下图……
首先是函数的定义：
target := Hypergeometric2F1[-(1/2), -(1/2), 1, x^2]
ramanujan := 1 + (3 x^2)/(10 + Sqrt)
compensation := ((3 x^10)/2^17 (1 + ((79/48 + (\ - 1 - 79/48) x) x^2)/(1 + x^a (1 - x^b)^c)^d)) /. \ -> 2^17/3 (4/\ - 14/11)
err := target - ramanujan - compensation
接下来是个函数拟合优劣的评估函数：
sseMean :=
Module[{a = aVal, b = bVal, c = cVal, d = dVal, \, data, sse},
(*计算数据点时采用的步长*)
\ = 0.01;
(*产生数据点平方*)
data = Table[(err)^2, {x, \, 1, \}];
(*计算 SSE*)
sse = Total;
(*返回 SSE 的平均值*)
sse \
]
略微说明下评估拟合的原理：
如果能100%拟合，则两曲线将会重合；
不然的话，两个曲线之间所夹的面积就代表了失拟（曲线形状失配）的程度。
这个面积应该用积分来计算。但计算积分比较耗时，
所以改成离散点上误差平方和的平均值。
与计算两曲线间的最大误差功能类似。但计算速度上可能会更快一点。
下面是用这个函数画图的代码。函数前面加了负号，
使得最小值成为最大值，就是把曲面翻了个个，便于观察。
第一幅图：
Plot3D[{-sseMean,
-sseMean,
-sseMean,
-sseMean},
{a, 1, 100}, {b, 1, 100},
MeshFunctions -> {#3 &},
BoundaryStyle -> Thick,
AxesLabel -> Automatic,
PlotLegends -> {"c,d=10,20", "c,d=20,10", "c,d=20,40", "c,d=40,20"}]
第二幅图：
Plot3D[{-sseMean,
-sseMean,
-sseMean,
-sseMean},
{b, 1, 100}, {c, 1, 100},
MeshFunctions -> {#3 &},
BoundaryStyle -> Thick,
AxesLabel -> Automatic,
PlotLegends -> {"a,d=10,20", "a,d=20,10", "a,d=20,40", "a,d=40,20"}]
第三幅图：
Plot3D[{-sseMean,
-sseMean,
-sseMean,
-sseMean},
{c, 1, 100}, {d, 1, 100},
MeshFunctions -> {#3 &},
BoundaryStyle -> Thick,
AxesLabel -> Automatic,
PlotLegends -> {"a,b=10,20", "a,b=20,10", "a,b=20,40", "a,b=40,20"}]
第四幅图：
Plot3D[{-sseMean,
-sseMean,
-sseMean,
-sseMean},
{a, 1, 100}, {d, 1, 100},
MeshFunctions -> {#3 &},
BoundaryStyle -> Thick,
AxesLabel -> Automatic,
PlotLegends -> {"b,c=10,20", "b,c=20,10", "b,c=20,40", "a,b=40,20"}]
从图上可以看出参数的不同只影响最佳参数点的位置，并不影响最佳参数点的值。
这是直觉的由来。若要证这一点，或可从误差函数的级数展开式入手。
注：全局搜索的结果就是到了这个参数最佳点，但这个点上拟合效果并不好。
再放两个演示参数变化如何影响拟合误差的动画。
选的演示点就在全局最佳参数的区域。
第一动画 err ~ a, b ：
Animate,
{x, 0, 1},
GridLines -> Automatic,
PlotRange -> All,
Frame -> True],
{a, 1, 100},
{b, 1, 100},
AnimationRunning -> False]
第二个动画 err ~ c, d ：
Animate,
{x, 0, 1},
GridLines -> Automatic,
PlotRange -> All,
Frame -> True],
{c, 1, 100},
{d, 1, 100},
AnimationRunning -> False]
在动画演示的过程中还可以看到计算发生困难的情况。
这是因为在补偿函数里采用了指数函数的原因：
\[\begin{matrix}f(x)=a^x,&a>0且a\ne1\end{matrix}\]
虽然有这个直觉，但并不能否定有可以满足拟合精度要求的解。
不过即使有，找到这个点的概率可能也不太大，就是要有点运气的意思。
鉴于成功前景不佳，所以只能另寻其它方法了。
Jack315 发表于 2023-8-30 06:51
先说个直觉，这个补偿函数可能有个天花板。
直觉来自图形，不是证明。先用下列代码看下图……
可以参靠一下楼上前辈上传的这个文件包：
Jack315 发表于 2023-8-30 06:51
先说个直觉，这个补偿函数可能有个天花板。
直觉来自图形，不是证明。先用下列代码看下图……
可以看看52楼吗？代码来源54楼压缩包
本帖最后由 Ickiverar 于 2023-8-30 10:55 编辑
笨笨 发表于 2023-8-30 00:21
先生还在吗，楼主没看明白这一步，为啥引进红线部分函数：Sin，它有什么特殊之处，而不是其它函数 ...
因为显然Er在x接近1的时候变化非常剧烈，需要采样更密集才能捕捉上升、下降和极值点。sinx起到让等距的x在1附近采样更密的条件。
函数是直接在Mathematica里敲的。这个是以前用过的函数所以直接粘上来了。
Mathematica是一种编程语言，不是纯粹的计算器：
它里面那些内置的优化、求解函数，FindMaximum，Solve，FindRoot什么的，都有内在的缺陷。简单函数还可以凑合用用，但不保证对复杂问题可以求解，或者解的正确性。因为它是编程语言，你需要在对问题有清晰认识的情况下，自己思考一个逻辑出来，用更基本的控制语句实现一个算法，才能充分利用Mathematica求解所有问题。
这是前辈区域的图：
这是我找到的一个区域的图：
（不能上图了:L ）
上画图代码……(参考53#)
前辈的：
ContourPlot,
{c, 0.5, 2}, {d, 0.5, 2},
Contours -> Function[{min, max}, Range],
PlotLegends -> Automatic
]
Plot3D,
{c, 0.5, 2}, {d, 0.5, 2},
MeshFunctions -> {#3 &},
BoundaryStyle -> Thick,
AxesLabel -> Automatic,
PlotLegends -> {"a,b,c,d=-2.63,11.47,0.94,1.36"}]
这是俺的：
ContourPlot,
{c, 0.1, 10}, {d, 1, 20},
Contours -> Function[{min, max}, Range],
PlotLegends -> Automatic
]
Plot3D,
{c, 0.5, 7.5}, {d, 1, 10},
MeshFunctions -> {#3 &},
BoundaryStyle -> Thick,
AxesLabel -> Automatic,
PlotLegends -> {"a,b,c,d=9,16,3.6,10"}]
注意两者数量级是一样的，只是周围地势不一样。
Jack315 发表于 2023-8-30 06:51
先说个直觉，这个补偿函数可能有个天花板。
直觉来自图形，不是证明。先用下列代码看下图……
53#
注：全局搜索的结果就是到了这个参数最佳点，但这个点上拟合效果并不好。
这就是说，需要在这个区域用前辈的方法仔细找最佳点。
55#（点评）
还有最后一个希望：在区域打概率，因为路面不平有望中奖的。
打概率的意思就是在【最佳点区域】内，随机选取 (a,b,c,d)，
然后看看拟合效果最佳的里面，有没有符合要求的。
俺找到的区域地势比较险峻，中奖的可能性或许会更大一点。
本帖最后由 笨笨 于 2023-8-30 12:41 编辑
Ickiverar 发表于 2023-8-30 10:51
因为显然Er在x接近1的时候变化非常剧烈，需要采样更密集才能捕捉上升、下降和极值点。sinx起到让等距的x ...
发贴失误……
本帖最后由 Ickiverar 于 2023-8-30 12:30 编辑
笨笨 发表于 2023-8-30 11:44
下面的作差效果要好点：
首先Hypergeometric的定义就是那个sum式，你这个行为和 x - x 没有区别，结果当然是恒0。
而EllipticE定义的F是定义不同的代数等价形式，它们的数值是完全精确相等的，只是因为计算式不同，因为浮点精度限制有所区别，所以画出来是毛毛糙糙的10^-16级别的误差。你可以用WorkingPrecision指定更高精度再试。
你进行任何数值计算，都会受到这个至少10^-16精度的影响。我画出10^-16级别的误差图就已经足够说明它们实际上完全一致了。
Plot[Evaluate[
FullSimplify@
Sum/n!)^2, {n, 0,
Infinity}] - F], {x, 0, 1}, WorkingPrecision -> 30]