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

Jack315

任取一个周期（25 个点）的数据，取模型中常数项初值为 \(c=1.0\)，初始步长 \(h=0.1\)，进行下列迭代：

• 在 \([c-h,c+h]\) 范围内均匀取 11 个点 \(c_i\)。
• 在每个点 \(c_i\) 上用 DFT 计算 \(y-c_i\) 的常数项。
• 以误差绝对值最小的 \(c_i\) 作为新的迭代参数 \(c\)，同时步长 \(h\) 缩小到原来的 1/10 。
  

重复上述迭代过程，直到误差绝对值最小的 \(c_i\) 满足指定的数度要求。

【Mathematica 代码】
加载数据

1. ClearAll["Global`*"]
   
2. SetDirectory["C:\"];
   
3. tbl = Import["Data.txt", "Table"];
   
4. n = Table[tbl[[i]][[1]], {i, Length[tbl]}];
   
5. y = SetPrecision[Table[tbl[[i]][[2]], {i, Length[tbl]}], 100];

复制代码

常数分量计算函数

1. calConstant[y_] :=
   
2. (*=======================================================*)
   
3. (*计算 y 中的常数分量。*)
   
4. (*=======================================================*)
   
5. (*【输入参数】*)
   
6. (*y: {Subscript[y, 1],Subscript[y, 2],...,Subscript[y, 25]} 原始数据。*)
   
7. (*【输出数据】*)
   
8. (*常数分量。*)
   
9. Module[{ferr, optimal, step, offset, err},
   
10.   (*-------------------------------------------------------*)
    
11.   (*初始化。*)
    
12.   (*-------------------------------------------------------*)
    
13.   ferr[data_, bias_] := 1/5 Fourier[data - bias][[1]];
    
14.   optimal = SetPrecision[1, 100];
    
15.   step = SetPrecision[0.1, 100];
    
16.   (*-------------------------------------------------------*)
    
17.   (*迭代。*)
    
18.   (*-------------------------------------------------------*)
    
19.   While[step > 10^-48,
    
20.    offset = Subdivide[optimal - step, optimal + step, 10];
    
21.    err = Abs[Table[ferr[y, offset[[i]]], {i, Length[offset]}]];
    
22.    optimal = Part[offset, Position[err, Min[err]][[1]][[1]]];
    
23.    step /= 10;
    
24.    ]; (*end while*)
    
25.   (*-------------------------------------------------------*)
    
26.   (*返回结果。*)
    
27.   (*-------------------------------------------------------*)
    
28.   optimal
    
29.   ] (*end Module *)

复制代码

计算并保存结果：

1. result = Table[{n[[i]], calConstant[y[[i - 24 ;; i]]]}, {i, 25, Length[y], 25}];
   
2. Export["Constant.txt", result, "TSV"]

复制代码

【结果】
小数点后稳定的数字比原数据多了 1~2 位。
这是小数点后稳定的 33 位数字：0.998953960449742112505647176797204

如果数据中的常数项信息都已基本被提取了出来，
则意味着这个数据并没有提供预测未来某个时间步常数项小数点后 40 位数字所需的信息。

如果假设前提不成立，还有什么方法来提取更多的信息？

Jack315

原数据减去 0.99895396044974211 再乘以 1e+18 作图看到了与 8# 形状相同的趋势。
因此模型应该修改为：
\(\begin{matrix}y(t)=\frac{p_1t^2+p_2t+p_3}{t^2+q_1t+q_2}+\exp(\beta_0+\beta_1t)\cos(\omega_ct+\varphi)&t=1,2,3,\cdots\end{matrix}\)
据此有：\(y(\infty)=p_1\)

建模数据取自中心左右各四分之一的部分，模型参数：
\(p_1\): 0.998953960449742112505647162231022030432
\(p_2\): 61.92624846613714685707263143501205487642
\(p_3\): 7240.718073573050381480974427123042069732
\(q_1\): 61.99109360180837709884047211271530449650
\(q_2\): 7248.300082131097003177689864784403599322
\(\beta_0\): -51.01143750080362404210707554419013802146
\(\beta_1\): -0.01004441927896885604364847479566826014080
\(\omega_c\): 0.2509257588303796389235689874214851411490
\(\varphi_c\): 0.1300513307914659036485461634070844812088

残差图：


离真相应该不太远了。
估计模型再加上 \(f(t)\cos[\omega(t)+\varphi]\) 这一项可能会有更理想的结果。

【讨论】
如果 \(f(x)\) 是一个衰减指数或幂函数，则 \(y(\infty)=p_1=\) 0.998953960449742112505647162231022030432...
需要经过非常长的时间才能得到这个极限，否则真实的有效数字只有小数点后 26~27 位的样子。

如果 \(f(x)\) 也是一个有理分式函数，则 \(y(\infty)\) 很有可能根本就不存在。

【更新】
非常可能：
\(f_(x)\) 是一个衰减的幂函数。
\(\omega(t)\) 是一个自然对数函数。
\(y(\infty)=p_1\)

KeyTo9_Fans

由于我要的是振幅衰减到无限小后，整个数列的极限值，所以我是优先考虑尾部的数据的，因为尾部的数据振幅最小，它们最接近极限值。

我把尾部的100个数据喂给文心大模型4.0后，他返回给我的是一段python代码：

https://yiyan.baidu.com/share/FUanbmPf8z



我以为他是随便写的伪代码，但神奇的是，我把y值填进去之后，这段python代码竟然可以运行出结果：



现在的大语言模型都这么厉害了吗？他竟然能够根据我的问题描述，编写出一段能解决这个问题的程序代码，

而且代码的运行结果的第1个数值和我在3楼里给出的预测值39295561几乎一样！

Jack315

这是残差的数据：
Residuals.rar (12.57 KB, 下载次数: 1)

2024-4-12 07:47 上传

点击文件名下载附件

喂点给无所不能的 AI 不知道会出个啥
会不会是类似这样的函数：\(x^{-0.5}\cos(5\ln x)\)

Jack315

接 12#...

模型简记为：\(f_1(t)+f_2(t)\cdot f_3(t)\)
继续寻找 \(f_1(t)\) 的模型感到有点难度了。

理论上只有完全确定准确的模型形式、
或者三个函数分量用傅立叶变换后在频率上能完全分离，
才可能得到 \(y(\infty)\) 比较准确的估计。

以我自己目前的结果，或许只能得到小数点后约 34 位有效数字的估计：
\(y(\infty)\approx\) 0.9989 5396 0449 7421 1250 5647 1767 9720 39

Jack315

指数衰减：

从图形上看这是一条直线，但查看细节，会有所发现。
首选用傅立叶变换取出幅度值：

1. amp = Log[Abs[Table[1/5 Fourier[y[[i - 24 ;; i]]][[2]], {i, 25, Length[y]}]]];

复制代码

然后用下列函数建立多项式模型：

1. buildAmpPolyModel[data_, order_] :=
   
2. (*=======================================================*)
   
3. (*生成并返回多项式模型参数。*)
   
4. (*=======================================================*)
   
5. (*【输入参数】*)
   
6. (*data: {Subscript[y, 1],Subscript[y, 2],...} 数据。*)
   
7. (*order: 多项式阶数。*)
   
8. (*【输出数据】*)
   
9. (*建立的模型。*)
   
10. Module[{fun, var},
    
11.   (*-------------------------------------------------------*)
    
12.   (*初始化。*)
    
13.   (*-------------------------------------------------------*)
    
14.   fun = Plus @@ Table[ToExpression["\[Beta]" <> ToString[i]] t^i, {i, 0, order}];
    
15.   var = Table[ToExpression["\[Beta]" <> ToString[i]], {i, 0, order}];
    
16.   (*-------------------------------------------------------*)
    
17.   (*拟合并返回模型。*)
    
18.   (*-------------------------------------------------------*)
    
19.   NonlinearModelFit[data, fun, var, t]
    
20.   ] (*end Module *)
    
21. 
    
22. mdlAmp = buildAmpPolyModel[amp, 20]

复制代码

模型残差如下图所示：

在不知道准确模型的情况下，无法知道这个残差
究竟是数据生成过程本身的特性，还是因为耦合了振荡所造成。

取最后 25 个数据点拟合下列模型：
\(y(t)=\alpha+\exp(\sum_{i=0}^n\beta_it^i)\cos(\omega t+\varphi)\)
当 \(n=1\) 时，得到 \(\alpha\) 95% 置信区间为：
置信下限：0.9989 5396 0449 7421 1250 5647 1767 9720 3929 5351
置信上限：0.9989 5396 0449 7421 1250 5647 1767 9720 3929 5712
小数点后 37 位有效数字。

增加 \(n\) 收效甚微。这是当 \(n=4\) 时的残差图：

不知道还有什么好的办法能获得更理想的结果。

Ickiverar

所以你具体在算什么？
也许有更好的算法。而且我有空闲算力，如果问题有趣的话也可以重算一遍。

Jack315

用时间序列的 AI 来做预测，这个或许还可以再玩一下

KeyTo9_Fans

Ickiverar 发表于 2024-4-15 04:43
所以你具体在算什么？
也许有更好的算法。而且我有空闲算力，如果问题有趣的话也可以重算一遍。 ...

我算的是下面这个问题的答案在n趋于无穷大时的极限值：
https://bbs.emath.ac.cn/thread-1819-1-1.html

这个问题早在1963年就已经有人给出0.593这3位数的近似结果了，我（This paper）在2011年发表了9位小数的计算结果，但在2015年被13位小数的计算结果赶超了，如下表所示：



我正在计算的是极限值的第14位小数及以后的数据，企图在2024年赶超2015年的计算结果。