curve_fit 函数求和

Question

我的学士论文需要拟合广义麦克斯韦函数。函数如下：

我从 .csv 文件中获取数据 (x,y) 并将其用于 curve_fit。

目前我正在研究一阶拟合（所以我用 N = 1 填充了公式，这样对我自己来说更容易）。我不知道如何将求和的额外部分添加到我的函数中并适合所有额外参数。我知道 N 的最大值为 10.

def first_order(G_0, G_1, t_1, omega):
    return G_0 * ((G_1*t_1*omega)/(1+(t_1**2)*(omega**2)))


def calc_gmm(dframe):
    array_omega = np.array(dframe['Angular Frequency']).flatten()
    array_G = np.array(dframe['Loss Modulus']).flatten()

    print(array_omega)
    print(array_G)
    variables, _ = curve_fit(first_order, array_omega, array_G, p0=[1,5,0.5])
    print(variables)
    print(_)

    plt.figure(figsize=(12,8))
    plt.plot(array_omega, array_G)
    plt.plot(array_omega, first_order(variables[0], variables[1], variables[2], array_omega))
    plt.show()

这是一些示例数据。

这是固定的结果（在SO的帮助下）。

Answer 1

您可以将函数参数设为 g 和 tau 数组，然后使用 sum.

def gmm(G_0, g, tau, omega):
    return G_0 * ((g*tau*omega)/(1+(tau**2)*(omega**2))).sum()

示例：

import numpy as np
np.random.seed(1)
g = np.random.rand(10)
tau = np.random.rand(10)
gmm(1, g, tau, 0.5)  # returns 0.6530812207319884

Answer 2

流变学！我就此完成了我的博士论文，并且我拟合了许多麦克斯韦模型。这是我给你的建议。

首先，您对 G' 和 G'' 都感兴趣，还是只对 G'' 感兴趣？我通常必须同时满足两者，并且为了获得更好的结果，G' 和 G'' 上的弛豫时间和模量必须相同，所以我认为您必须改变您的方法来考虑这一点。

其次，我认为像 lmfit 这样的包更适合这样做，因为您可以更好地控制最小化功能。

第三，由于 n 是一个整数，我认为你必须在 n=1、n=2、...、n=10 评估你的模型并检查参数的标准误差。太多是过拟合，太少是欠拟合。我认为不能真正实现自动化。

让我们先构建一些玩具数据。

import matplotlib.pyplot as plt
import numpy as np
import lmfit

def G2Prime(g_i, t_i, w):  # G''
    return g_i * (t_i * w) / (1 + t_i ** 2 * w **2)

def GPrime(g_i, t_i, w):   # G'
    return g_i * (t_i * w)**2 / (1 + t_i ** 2 * w **2)

# Generate a sample model with 3 components
omegas = np.logspace(-2, 1)
# G0 = 1
test_data_GPrime = 1 + GPrime(1, 1, omegas) + GPrime(1, 10, omegas) + GPrime(1, 30, omegas)
test_data_G2Prime = 1 + G2Prime(1, 1, omegas) + G2Prime(1, 10, omegas) + G2Prime(1, 30, omegas)

这是图表。

接下来，让我们创建要使用的参数 lmfit。

params = lmfit.Parameters()  # Creates a parameter object
params.add('n', value=2, vary=False, min=1, max=10) # start with n=2, so it's not exact
params.add('G0', value=1, min=0)
for i in range(params['n'].value):   # Adds the relaxation times and moduli separately
    params.add(f't_{i}', value=1, min=0)
    params.add(f'g_{i}', value=1, min=0)

然后，让我们定义同时考虑 G' 和 G'' 的最小化函数。

def min_function(params, x, data_GPrime, data_G2Prime):
    n = int(params['n'].value)
    G0 = params['G0']
    # Calculate the first component
    model_GPrime = G0 + GPrime(params['g_0'], params['t_0'], x)
    model_G2Prime = G0 + G2Prime(params['g_0'], params['t_0'], x)
    for i in range(1, n): # Go through the other components
        model_GPrime += GPrime(params[f'g_{i}'], params[f't_{i}'], x)
        model_G2Prime += G2Prime(params[f'g_{i}'], params[f't_{i}'], x)
    # return the total residual of both G' and G''.
    return (model_GPrime - data_GPrime) + (model_G2Prime - data_G2Prime)

最后，让我们调用最小化函数。使用这种方法，您不能使用变化 n，因此您必须自己进行变化。

res = lmfit.minimize(min_function, params, args=(omegas, test_data_GPrime, test_data_G2Prime))

让我们看看 n=2 的结果。

plt.plot(omegas, test_data_GPrime)
plt.plot(omegas, test_data_GPrime + res.residual, c='r', ls='--')
plt.plot(omegas, test_data_G2Prime)
plt.plot(omegas, test_data_G2Prime + res.residual, c='r', ls='--')
plt.xscale('log')
plt.yscale('log')

n=3是绝配，就不展示了。这是拟合的输出报告，其中 lmfit.report_fit(res).

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 72
    # data points      = 50
    # variables        = 5
    chi-square         = 0.04825415
    reduced chi-square = 0.00107231
    Akaike info crit   = -337.164818
    Bayesian info crit = -327.604702
[[Variables]]
    n:    2 (fixed)
    G0:   1.10713874 +/- 0.01190976 (1.08%) (init = 1)
    t_0:  1.11030322 +/- 0.02837998 (2.56%) (init = 1)
    g_0:  1.07272282 +/- 0.01532421 (1.43%) (init = 1)
    t_1:  16.6536979 +/- 0.34791430 (2.09%) (init = 1)
    g_1:  1.71017461 +/- 0.02099472 (1.23%) (init = 1)
[[Correlations]] (unreported correlations are < 0.100)
    C(G0, g_1)  = -0.769
    C(G0, t_1)  = -0.731
    C(g_0, t_1) =  0.699
    C(t_0, g_0) =  0.497
    C(t_0, t_1) =  0.493
    C(G0, g_0)  = -0.442
    C(t_1, g_1) =  0.263
    C(t_0, g_1) = -0.255
    C(G0, t_0)  = -0.231
    C(g_0, g_1) = -0.157

现在，您必须遍历其他可能 n，检查拟合参数并确定哪个是理想的。

curve_fit 函数求和

curve_fit of a summation of functions

python

model-fitting