振幅指数级衰减的数据拟合方法

KeyTo9_Fans

由于竞价实例被回收的缘故：



我正在运行的程序被腾讯云中断了，只拿到了中断前的输出数据：

n        A(n)
8751 0.9989539604497421124893224675482733328353
8752 0.9989539604497421134029256815780793109571
……
11975 0.9989539604497421125056471767972029556836

省略号里的数据在下面这个附件里：

data.zip (36.43 KB, 下载次数: 6)

2024-4-2 18:20 上传

点击文件名下载附件

如果我的程序没有被中断，那么当n=13750时，A(n)的值就能达到我需要的精度，不需要数据拟合

但现在我的程序被中断了，我想用数据拟合的方法，预测一下当n=13750时，A(n)的前40位小数是多少

（重新运行一次大约需要花费160块钱，我不想花这笔钱重新算一次，万一又中断，这钱还会打水漂）

特此向论坛里的大神们求助~

Jack315

提供个思路。

下图为 Mathematica 所画的图：


1. 可以建立振幅的模型，并预测第 n 个数据的振幅。

2. 通过前期数据，可以用均值估计小数点前若干位的“真值”。
而后期的数据随着精度的提高，提供了这个“真值”。因而可以得到估计值的误差。
如果振幅与上述误差也存在有意义的模型，则有可能预测出具有所需精度的“真值”。

说明：

• 0.99895396044974211250564717679720 是 32 位精度的“真值”。
• Matlab 是常规的单精度或双精度数据，不能满足要求。
• Mathematica 有任意精度浮点数；Python 有任意位数整数。应该都能用于这个问题的相关计算。
  

KeyTo9_Fans

还好他中断的时候，振幅就只剩10的-32次方了

小数点后的第33~40位我简单地估算了一下，大概是39295561，如下表所示：

n                 A(n)
8751        0.9989539604497421124893224675482733328353
8752        0.9989539604497421134029256815780793109571
……
11975        0.9989539604497421125056471767972029556836
……
$\infty$        0.9989539604497421125056471767972039295561

我是用这个模型估算出A(n)的极限值的：

(A(n)-0.998953960449742112505647176)*10^28 ≈ 7.972039295561 + exp(-0.01004512*(n-11000)-1.05418)*cos(0.25092733*(n-11000)-1.70536)

如下图所示：



我对第38~40位的估算结果是没有把握的，不知道还有没有更精准的模型，把第38~40位估算得准一些呢？

Jack315

如图所示：

注意到数据的振荡部分上尖下钝，应该不是简单的余弦函数。
或许改进下模型能得到更好的结果？

能晒晒模型残差的图吗？

Jack315

将 LZ 的模型改写成下列等效形式：
\(A(n)=0.9989539604497421125056471767972039295561+\)
\(10^{-28}exp[-0.01004512(n-11000)-1.05418]\cos[0.25092733(n-11000)-1.70536]\)
数据文件 Data.txt 中第一列数据对应 \(n\)，第二列数据对应 \(A(n)\)。
做下列运算：
\(y_1=\frac{A(n)-0.9989539604497421125056471767972039295561}{10^{-28}exp[-0.01004512(n-11000)-1.05418]}\)
\(y_2=y_1-\cos[0.25092733(n-11000)-1.70536]\)
即将 \(A(n)\) 减去常数项，再除去指数分量，得到余弦分量。两个余弦分量的差作为模型的残差 \(y_2\) 。
作出残差图：

上面一张图是全部数据的残差图，下面一张是 \(n\ge10750\) 的数据残差图。
残差图表明模型参数或有进一步优化的空间。

Jack315

将模型定义为：
\(\begin{matrix}y(n)=c+a\exp(kx)\cos(px+q),&x=n-805\end{matrix}\)

数据文件 Data.txt 中第一列数据对应变量 \(n\) ，第二列数据对应变量 \(A(n)\) 。
取中间 \(806\le n\le2418\) 一段数据建模：
自变量：\(x=1,2,\dots,1613\)
应变量：\(y_0=A(806),A(807),\dots,A(2418)\)

取数据文件 Data.txt 中最后一行的数据 \(A(11975)\) 作为模型参数 \(c\) 的初始值。
即取：\(c = 0.9989539604497421125056471767972029556836\) 。

减去常数项，并取绝对值和对数，即：\(y_1=\log|y_0-c|\) 。
拟合线性模型：\(y_1=kx+\beta\) 得 \(k=-0.0100399257488,\beta=51.708275431\) ，
取 \(k\) 的初始值 \(k=-0.0100399257488\) 。

再去除指数分量：\(y_2=10^{23}(y_0-c)\exp(-kx)\) 。
注：常数 \(10^{23}\) 将数据放大，以方便模型拟合计算。
拟合正弦模型：\(y_2=a\sin(px+q)\) 得 \(a=6.986,p=0.2509,q=1.701\)，
取 \(a,p,q\) 的初始值 \(a=6.986\times10^{-23},p=0.2509,q=1.701\) 。

至此得到模型参数的初始值：
\(c = 0.99895396044974211250564717679720295568360\)
\(a=6.986\times10^{-23}\)
\(k=-0.0100399257488\)
\(p=0.2509\)
\(q=1.701\)

从模型的初始参数出发，通过优化模型参数，使得模型残差最小。
模型残差的衡量指标为两个正弦分量的均方差 (MSE - Mean Squared Error) ：
\(MSE=\frac{1}{N}\sum_{i=1}^N[y_{observed}(i)-y_{predicted}(i)]^2\)

下图为初始模型参数下得到的残差图：

这个图说明模型参数有待优化。

Jack315

下图为去除常数项和指数分量后的数据：

通过傅立叶变换后，画出频谱图：

这个频谱结果表明模型的形式应该是：
\(\begin{matrix}y(n)=c+a\exp(kx)\cos(px+q)\sum_{t=1}^m\cos(tux+v),&x=n-805\end{matrix}\)
其中：
\(u<p\) ；
\(m=(1,30)\) ，视精度要求而定。
对于实际情况，这个模型不知道是否说得通。
可以先试试 \(m=1\) 的情况，一次加一个余弦分量。

Jack315

拟合模型 \(\begin{matrix}y(t)=\alpha+\exp(\beta_0+\beta_1t)\cos(\omega_ct+\varphi)&t=1,2,3,\cdots\end{matrix}\) 得到下列参数值（100 个数字精度）：
\(\alpha\) : 0.9989539604497421125056471767971999999999999999999999999999311015804376986097035328711763334810745355
\(\beta_0\) : -51.015619015721981240221401471857680531610703346762303492719086604684890126358315979396096699546571657
\(\beta_1\) : -0.01003982393965083431145124623483739033177033532910499176433880583718181839729084561842989035728673070
\(\omega_c\) : 0.2509267866097859372353970050172701274538115876708992946364332536654003902156089168959049451324005986
\(\varphi\) : 0.1301875971824696437222819241053059827064484189848839294069835642110835241946799272047714801647799732
其中：\(\alpha\) 取数据文件最后稳定的小数点后 32 位有效数字。其它参数通过拟合模型得到。

对此初始模型参数尝试进行优化，试图寻找使 RMSE 更小的参数，用了二种不同的方法均告失败。
详细过程参考下列 Mathematica 代码：
ModelParameter.rar (344.51 KB, 下载次数: 0)

2024-4-8 16:32 上传

点击文件名下载附件

【以下讨论可能有问题，暂留个记录。】
用上述模型参数作残差图如下：

这个残差图有明显的模式，即与假设的正态随机分布不符，提示模型仍可进一步改进。
残差图的左边提示一个类似对数函数的形状：\(f_1(t)\)。
如果在 \(t=1\) 处的残差能被 \(f_1(t)\) 调整到 0 附近，右边基本不动，且整个残差“平稳”，
则残差有可能提示一个 \(\begin{matrix}f_2(t)=\cos(\omega_ct+\varphi)\sum_{k=1}^n\cos(k\omega t+\varphi_k),&\omega<\omega_c\end{matrix}\) 的模型。

如果是 \(f_1(t)+f_2(t)\) 的形式，则改进的模型形式为：
\(y(t)=\alpha+\exp(\beta_0+\beta_1t)\cos(\omega_ct+\varphi)+f_1(t)+f_2(t)\)
\(=\alpha+[\exp(\beta_0+\beta_1t)+\sum_{k=1}^n\cos(k\omega t+\varphi_k)]\cos(\omega_ct+\varphi)+f_1(t)\)

如果是 \(f_1(t)f_2(t)\) 的形式，则改进的模型形式为：
\(y(t)=\alpha+\exp(\beta_0+\beta_1t)\cos(\omega_ct+\varphi)+f_1(t)f_2(t)\)
\(=\alpha+[\exp(\beta_0+\beta_1t)+f_1(t)\sum_{k=1}^n\cos(k\omega t+\varphi_k)]\cos(\omega_ct+\varphi)\)

还是老问题：对于实际情况，改进的模型不知道是否说得通。

Jack315

【模型改进】
由 8# 的残差图构造一个插值函数 \(f_1(t)\)，如下图红色曲线：

减去 \(f_1(t)\) 后的 \(f_2(t)\) 如下图所示：

\(f_2(t)\) 的频谱如下图所示：

这时的模型为：
\(y(t)=\alpha+\exp(\beta_0+\beta_1t)[\cos(\omega_ct+\varphi)-f_1(t)-f_2(t)]\)
感觉有点怪异，即使 \(f_1(t), f_2(t)\) 都很完美，模型预测 \(\alpha=y(\infty)\) 的能力也有点可疑了

Jack315

振荡的周期约为 25，以依次移动的 25 个数据点用离散傅立叶变换可得到振荡分量的幅度变化过程。
如果这个幅度变化过程可以被精确拟合，则剩下的部分同样用离散傅立叶变换可取得振荡分量的常数项。
微调幅度模型的参数，使得这个常数项在整个过程中小数点后 40 位保持不变，即可得到所需结果。

把幅度取对数。如果是简单的函数，这个数据经过若干次差分后应该能得到比较容易识别的函数形式。
不过现在还是看不出来是个啥。或许数据还是被振荡分量干扰了。
但用直线拟合后的残差有个固定的弯弓模式，不管数据起止位置取自何处。