使用 python(或其他有用的语言)在积分内拟合参数

Fitting parameter inside an integral using python (or another useful language)

我有一组数据,基本上是f(x)作为x的函数的信息,以及x本身。我从问题的理论中知道我正在研究 f(x) 的格式,它由下面的表达式给出:

本质上我是想用这组数据求出参数a和b。我的问题是:我该怎么做?我应该使用什么图书馆?我想要使​​用 Python 的答案。但 R 或 Julia 也可以。

从我到目前为止所做的一切来看,我已经从 SciPy 库中读到了一个名为 curve fit 的功能,但是我遇到了一些问题,我将以哪种形式编写代码long 我的 x 变量位于积分限制之一。

为了更好地解决问题,我还有以下资源:

A sample set,我知道我正在寻找的参数。对于这个集合,我知道 a = 2 和 b = 1(并且 c = 3)。在它提出一些关于我如何知道这些参数的问题之前:我知道它们是因为我使用这些参数从上面的方程的积分创建了这个样本集,只是为了使用样本来研究我如何找到它们并获得参考。

I also have this set,对此我唯一的信息是 c = 4 并且想找到 a 和 b。

我还想指出:

i) 现在我这里没有 post 的代码,因为我不知道如何编写一些东西来解决我的问题。但我很乐意在阅读你们可以提供给我的任何答案或帮助后编辑和更新问题。

ii) 我首先要寻找一个我不知道 a 和 b 的解决方案。但如果它太难了,我会很高兴看到一些解决方案,我认为 a 或 b 是已知的。

编辑 1: 我想参考 this question 任何对这个问题感兴趣的人,因为它是对这里面临的问题的并行但重要的讨论

它们是三个不独立的变量a,b,c。如果我们想通过回归计算另外两个,则必须给出其中一个。对于给定的 c,求解 a、b 很简单:

下面的数值微积分的例子是用一个小数据(n=10)做的,以便于检查。

请注意,当数据分散时,函数 t(y) 的回归与 y(x) 的回归并不完全相同(如果没有分散,结果相同)。

如果绝对有必要对 y(x) 进行回归,则需要非线性回归。这涉及从对 a,b 的足够好的初始猜测开始的迭代过程。上面的微积分给出了很好的初始值。

另外:

与此同时,Andrea 发布了一个中肯的答案。当然用他的方法拟合更好,因为这是非线性回归,而不是上面注释中已经指出的线性回归。

尽管如此,尽管 (a=1.881 ; b=1.617) 与 (a=2.346, b=-0.361) 的值不同,但下面绘制的各个曲线相距不远:

蓝色曲线:来自线性回归(上述方法)

绿色曲线:来自非线性回归 (Andrea's)

第二组数据的大小写

https://mega.nz/#!echEjQyK!tUEx0gpFND7gucvsTONiB_wn-ewBq-5k-pZlfLxmfvw

回归失败,因为假设 c=3 为假。

在c=0的情况下,积分的解析与上面不同:

我会使用纯数值方法,即使不能直接求解积分也可以使用它。这是一个仅适合 a 参数的剪接器:

import numpy as np
from scipy.optimize import curve_fit
import pandas as pd
import matplotlib.pyplot as plt

def integrand(x, a):
    b = 1
    c = 3
    return 1/(a*np.sqrt(b*(1+x)**3 + c*(1+x)**4))

def integral(x, a):
    dx = 0.001
    xx = np.arange(0, x, dx)
    arr = integrand(xx, a)
    return np.trapz(arr, dx=dx, axis=-1)

vec_integral = np.vectorize(integral)

df = pd.read_csv('data-with-known-coef-a2-b1-c3.csv')
x = df.domin.values
y = df.resultados2.values
out_mean, out_var = curve_fit(vec_integral, x, y, p0=[2])

plt.plot(x, y)
plt.plot(x, vec_integral(x, out_mean[0]))
plt.title(f'a = {out_mean[0]:.3f} +- {np.sqrt(out_var[0][0]):.3f}')
plt.show()

vec_integral = np.vectorize(integral)

当然,您可以降低dx的值以获得所需的精度。虽然只拟合 a,但当您也尝试拟合 b 时,拟合不会正确收敛(在我看来,因为 ab 是强相关的) .这是你得到的:

def integrand(x, a, b):
    c = 3
    return 1/(a*np.sqrt(np.abs(b*(1+x)**3 + c*(1+x)**4)))

def integral(x, a, b):
    dx = 0.001
    xx = np.arange(0, x, dx)
    arr = integrand(xx, a, b)
    return np.trapz(arr, dx=dx, axis=-1)

vec_integral = np.vectorize(integral)

out_mean, out_var = sp.optimize.curve_fit(vec_integral, x, y, p0=[2,3])
plt.title(f'a = {out_mean[0]:.3f} +- {np.sqrt(out_var[0][0]):.3f}\nb = {out_mean[1]:.3f} +- {np.sqrt(out_var[1][1]):.3f}')

plt.plot(x, y, alpha=0.4)
plt.plot(x, vec_integral(x, out_mean[0], out_mean[1]), color='green', label='fitted solution')
plt.plot(x, vec_integral(x, 2, 1),'--', color='red', label='theoretical solution')
plt.legend()
plt.show()

如您所见,即使形成拟合的 ab 参数是 "not good",情节也非常相似。