可变输入函数scipy.curve_fit

问题描述

我正在处理我的实验数据的峰值反卷积，我想生成一个Python脚本，在其中我可以轻松地改变非线性曲线拟合/峰值反卷积的方程。使用高斯曲线和线性偏移量，scipy.Optimize.curve_fit可以很好地与以下代码配合使用：

def Combined(x,*params):
    off = Linear(x,params[0],params[1])
    peak1 = Gaussian(x,params[2],params[3],params[4])
    peak2 = Gaussian(x,params[5],params[6],params[7])
    peak3 = Gaussian(x,params[8],params[9],params[10])
    return off + peak1 + peak2 + peak3

popt, pcov = opt.curve_fit(Combined, data[10][0], data[10][1], method='lm', check_finite=True, p0=[0.1, 0.1, 115, 508.33, 7.1,130, 508.33, 7.1, 165.84, 508.33, 7.1])

所有方程式都预先在一个函数中定义：

def ZeroOrder(x,a):
    return a

def Linear(x,a,b):
    return a+b*x

def SecondOrder(x,a,b,c):
    return a+b*x+c*x**2

我想将函数创建为Combine(x,baseline='ZeroOrder',peak1='Gaussian',peak2='Gaussian',peak3='Gaussian')，其中我可以轻松地分配不同的峰值函数，而不必为每个组合创建特定的函数。然而，在我的理解中，curve_fit函数非常严格，需要Combined(x,*params)这样的输入函数。如何更改代码以使其在所需函数中工作？

---

解决方案

对于这种情况，我会使用更通用的名称，比如来自scipy的f_min。 以下是如何为实验数据创建通用回归变量的示例。

我创建抽象类只是为了定义Shape函数的接口，这有助于创建更通用的方法

在下面的示例中，求解器获得的系数与原始函数不同，但对于确定的域，结果函数满足收敛条件，这一定是因为我创建的函数有很多自由度。

from abc import ABC, abstractmethod
from typing import List

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fmin

from pprint import pprint

class ShapeFunction(ABC):
    def __init__(self, params_count: int):
        self.params_count = params_count
        self.params = [0] * params_count

    @property
    def params(self):
        return self.__params

    @params.setter
    def params(self, params: List[float]):
        if len(params) != self.params_count:
            raise ValueError(f"params count must be {self.params_count}")
        self.__params = params

    @abstractmethod
    def evaluate(self, t: List[float]) -> List[float]:
        pass


class Gaussian(ShapeFunction):
    def __init__(self):
        params_count = 3
        super().__init__(params_count)

    def evaluate(self, t: List[float]) -> List[float]:
        return self.Gaussian(t, *self.params)

    @staticmethod
    def Gaussian(x, mu, sigma, scale):
        return scale/(sigma*np.sqrt(2*np.pi))*np.exp(-0.5*((x-mu)/sigma)**2)


class ZeroOrder(ShapeFunction):
    def __init__(self):
        params_count = 1
        super().__init__(params_count)

    def evaluate(self, t: List[float]) -> List[float]:
        return self.ZeroOrder(t, *self.params)

    @staticmethod
    def ZeroOrder(x, a):
        return a


class Linear(ShapeFunction):
    def __init__(self):
        params_count = 2
        super().__init__(params_count)

    def evaluate(self, t: List[float]) -> List[float]:
        return self.Linear(t, *self.params)

    @staticmethod
    def Linear(x, a, b):
        return a+b*x


class Quadratic(ShapeFunction):
    def __init__(self):
        params_count = 3
        super().__init__(params_count)

    def evaluate(self, t: List[float]) -> List[float]:
        return self.Quadratic(t, *self.params)

    @staticmethod
    def Quadratic(x, a, b, c):
        return a+b*x+c*x**2


f1 = Gaussian()
f1.params = [0, 1, 10]

f2 = ZeroOrder()
f2.params = [1]

f3 = Linear()
f3.params = [1, 1]

f4 = Quadratic()
f4.params = [1, 1, 1]


shape_functions = [f1, f2, f3, f4]
coefs_count = sum([func.params_count for func in shape_functions])

original_coefs = []
for func in shape_functions:
    original_coefs.extend(func.params)

x = np.linspace(-5, 5, 100)
sample = f1.evaluate(x) + f2.evaluate(x) + f3.evaluate(x) + f4.evaluate(x)


def cost_function(coefs: list[float], *params):
    funcs: List[ShapeFunction] = params[0]
    parametric_values: List[float] = params[1]
    obj_y: List[float] = params[2]

    y = [0] * len(parametric_values)

    for func in funcs:
        func.params = coefs[:func.params_count]
        y += func.evaluate(parametric_values)

        coefs = coefs[func.params_count:]

    return np.sum((y - obj_y)**2)


x0 = [1] * coefs_count
solution = fmin(cost_function, x0, args=(shape_functions, x, sample))


for sol, coef in zip(solution, original_coefs):
    print(f"Original: {coef:.2f} -> solution: {sol:.2f}")


y_sol = np.sum([func.evaluate(x) for func in shape_functions], axis=0)

plt.plot(x, sample, label="Original")
plt.plot(x, y_sol, label="Solution", alpha=0.5, marker=".")

plt.legend()

plt.show()