KeyTo9_Fans 发表于 2024-4-2 18:26:56

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

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



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

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

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



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

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

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

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

Jack315 发表于 2024-4-2 23:02:00

本帖最后由 Jack315 于 2024-4-2 23:06 编辑

提供个思路。

下图为 Mathematica 所画的图:


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

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

说明:

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

KeyTo9_Fans 发表于 2024-4-3 00:45:05

还好他中断的时候,振幅就只剩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 发表于 2024-4-3 05:32:55

如图所示:

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

能晒晒模型残差的图吗?

Jack315 发表于 2024-4-4 18:50:57

本帖最后由 Jack315 于 2024-4-4 21:44 编辑

将 LZ 的模型改写成下列等效形式:
\(A(n)=0.9989539604497421125056471767972039295561+\)
\(10^{-28}exp[-0.01004512(n-11000)-1.05418]\cos\)
数据文件 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\)
即将 \(A(n)\) 减去常数项,再除去指数分量,得到余弦分量。两个余弦分量的差作为模型的残差 \(y_2\) 。
作出残差图:

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

Jack315 发表于 2024-4-4 20:45:42

将模型定义为:
\(\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^2\)

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

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

Jack315 发表于 2024-4-7 12:44:25

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

通过傅立叶变换后,画出频谱图:

这个频谱结果表明模型的形式应该是:
\(\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 发表于 2024-4-8 17:22:12

本帖最后由 Jack315 于 2024-4-8 20:56 编辑

拟合模型 \(\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 代码:


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

这个残差图有明显的模式,即与假设的正态随机分布不符,提示模型仍可进一步改进。
残差图的左边提示一个类似对数函数的形状:\(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 发表于 2024-4-8 20:29:00

【模型改进】
由 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)\) 的能力也有点可疑了:L

Jack315 发表于 2024-4-9 17:01:16

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

把幅度取对数。如果是简单的函数,这个数据经过若干次差分后应该能得到比较容易识别的函数形式。
不过现在还是看不出来是个啥。或许数据还是被振荡分量干扰了。
但用直线拟合后的残差有个固定的弯弓模式,不管数据起止位置取自何处。
页: [1] 2
查看完整版本: 振幅指数级衰减的数据拟合方法