如何选择最佳的拟合方法？

Question

我想将参数 （代码中的

xs

）拟合到函数

fp

，由

给出

其中

v

和

L

是给定参数。通过向量化并使用一些统一过滤，可以相对快速地计算

fp

。以下代码实现了这个优化问题，其中函数

cfit

（通过

fitf

）通过在每一步计算

xs

并将其与真实数据进行比较来缓慢调整

fp

的值 (

list

)在代码中），根据

问题：有没有办法改进这个特定函数的拟合方法以及我计算它的方式？我尝试使用

scipy.optimize

，但一切似乎都比这个简单的乘法梯度方法慢。

这是代码实现，其中跟踪均方误差

import numpy as np
import math

def fitfunction(list, maxiter, error_threshold):
    
    v = 1.4
    st = 5
    exp_v = np.exp(-1/v)
    x00 = np.array([(math.pi/(4*v))*i**(-2) for i in list])
    lm = 1000
    
    def mse(y_true, y_pred):
        mse_value = sum((yt - yp) ** 2 for yt, yp in zip(y_true, y_pred)) / len(y_true)
        return mse_value
    
    def fast_roll_add(dst, src, shift):
        dst[shift:] += src[:-shift]
        dst[:shift] += src[-shift:]
    
    # Main function
    def fp(x, L, v):
        n = len(x)
        y = np.zeros(n)
        last_exp_2_raw = np.zeros(n)
        last_exp_2 = np.ones(n)
        unitary = x.copy()
        for k in range(L+1):
            if k != 0:
                fast_roll_add(unitary, x, k)
                fast_roll_add(unitary, x, -k)
            exp_1_raw = last_exp_2_raw
            exp_1 = last_exp_2
            exp_2_raw = exp_1_raw + unitary / v
            exp_2 = np.exp(-exp_2_raw)
            y += (exp_1 - exp_2) / unitary
            last_exp_2_raw = exp_2_raw
            last_exp_2 = exp_2
        return y

    # Fitting functions
    def fitf(time, lst, x0, j):
        return x0[j] * (lst[j] / time[j])**2
    def cfit(time, lst, x0):
        result = np.empty_like(x0)
        for j in range(len(x0)):
            if fitf(time, lst, x0, j) < 10**(-20):
                result[j] = 10**(-20)
            elif abs(time[j] - lst[j]) < .5:
                result[j] = x0[j]
            else:
                result[j] = fitf(time, lst, x0, j)
        return result
    
    xs = x00
    ys = fp(xs, len(xs)//st, v)
    err = 10**10
    
    for j in range(maxiter):
        if err > error_threshold:
            xs = cfit(list, ys, xs) # Fitting
            ys = fp(xs, len(xs)//st, v) # Update function values
            err = mse(list[lm:-lm], ys[lm:-lm])
            print(str(j+1) + '/' + str(maxiter) + ' err: ' + str('{:.20f}'.format(err)), end="\r")            
        else:
            break

    fire_rates = ['{:.20f}'.format(i) for i in xs]
    time_sim = ys
    
    return [fire_rates, time_sim]

以及一些带有生成数据的工作示例

from scipy.ndimage import gaussian_filter1d

# Data generation
np.random.seed(18)
data = np.random.normal(loc=0, scale=10, size=10000).cumsum()
data = (data - data.min()) / (data.max() - data.min()) * 500

# Gaussian smoothing
data = gaussian_filter1d(data, sigma=100)

data_fit = fitfunction(data, 100, 2)

并可视化它

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.plot(data, label='Data')
plt.plot(data_fit[1], label='Data Fit')
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Data and Data Fit')
plt.legend()
plt.show()

简要评论：

• 自然，这个拟合“看起来”相当好并且相对较快，但我计划在更大的范围内（10^6点）进行，并且我也想保证尽可能最佳的拟合。我的乘法方法可能有限制。

• fp

  的实施本质上是周期性的，这会导致最终出现大量错误。因此，计算误差

  err

  时不包括这些区域。
• 贴合功能提供更好的贴合度

def fitf(time, lst, x0, j):
    return x0[j]**(np.log(time[j]) / np.log(lst[j]))

但是，当这些值太接近真实值时，它似乎非常不稳定。

Answer 1

此答案仅实现中等速度增益 (10%)。

但是，我删除了不必要的循环和数组副本。我还重命名了一些变量，如果它们在代码的其他部分具有不同的名称（但在其他方面相同），并删除了其他多余甚至未使用的变量。因此，代码更短并且更具可读性。 我希望它可以帮助其他人找到更实质性的改进。

import numpy as np

def fitfunction(data, maxiter, error_threshold):

    v = 1.4
    maximum_shift = len(data) // 5
    margin = 1000

    def mse(y_true, y_pred):
        return np.sum((y_true - y_pred)**2)

    def fast_roll_add(dst, src, shift):
        dst[shift:] += src[:-shift]
        dst[:shift] += src[-shift:]

    # Main function
    def fp(x, maximum_shift, v):
        unitary = x.copy()
        exponent = unitary / v
        new_power = np.exp(-exponent)
        y = (np.ones_like(x) - new_power) / unitary
        for k in range(1, maximum_shift+1):
            fast_roll_add(unitary, x, k)
            fast_roll_add(unitary, x, -k)
            exponent += unitary / v
            old_power = new_power
            new_power = np.exp(-exponent)
            y += (old_power - new_power) / unitary
        return y

    # Fitting function
    def cfit(y_true, y_pred, x):
        result = np.clip(x * (y_pred / y_true)**2, 1e-20, None)
        residuals = np.abs(y_true - y_pred)
        return np.where(residuals < .5, x, result)

    xs = np.pi / (4 * v * data**2)
    ys = fp(xs, maximum_shift, v)
    err = mse(data[margin:-margin], ys[margin:-margin])
    print("1" + '/' + str(maxiter) + ' err: ' + str('{:.20f}'.format(err)), end="\r")

    for j in range(1, maxiter):
        if err > error_threshold:
            xs = cfit(data, ys, xs) # Fitting
            ys = fp(xs, maximum_shift, v) # Update function values
            err = mse(data[margin:-margin], ys[margin:-margin])
            print(str(j+1) + '/' + str(maxiter) + ' err: ' + str('{:.20f}'.format(err)), end="\r")
        else:
            break
    print()

    return xs, ys


if __name__ == "__main__":

    import time
    import matplotlib.pyplot as plt

    from scipy.ndimage import gaussian_filter1d

    # Data generation
    np.random.seed(18)
    data = np.random.normal(loc=0, scale=10, size=10_000).cumsum()
    data = (data - data.min()) / (data.max() - data.min()) * 500

    # Gaussian smoothing
    data = gaussian_filter1d(data, sigma=100)

    tic = time.time()
    data_fit = fitfunction(data, 100, 2)
    toc = time.time()
    print(f"Time elapsed: {toc-tic:.2f} seconds")

    plt.figure(figsize=(10, 6))
    plt.plot(data, label='Data')
    plt.plot(data_fit[1], label='Data Fit')
    plt.xlabel('Index')
    plt.ylabel('Value')
    plt.title('Data and Data Fit')
    plt.legend()
    plt.show()