Подгонка кривой по интегралу, содержащему массив данных и функцию в Python

Question

У меня есть набор данных, описываемый интегралом с неизвестными константами, который я пытаюсь определить, используя python's curve_fit. Однако подынтегральная функция содержит функцию, умноженную на набор данных

def integrand(tm, Pm, args):
    dt, alpha1, alpha2 = args
    return Pm*(1-np.e**( -(alpha1 * (tm+dt))) )*np.e**(-(alpha2 * (tm+dt)))

Где Pm - это одномерный массив собранных данных импульсных характеристик, данных изображений импульсов и интегральной кривой

Оранжевая кривая представляет данные импульса, а синяя кривая - то, что интеграл должен оценить

tm - переменная интегрирования, а dt, alpha1, alpha2 - неизвестные константы, и границы интегрирования будут от 0 до tm.

Каков наилучший способ подбора кривой для интеграла такого типа или, возможно, других способов решения неизвестных констант?

Наборы данных связаны здесь

Answer 1

Из длины наборов данных кажется, что цель состоит в том, чтобы подогнать подынтегральное выражение (t) к выходному значению (t + dt). В модуле scipy optimize есть несколько функций, которые можно использовать для этого. Для простого примера мы покажем реализацию, использующую scipy.optimize.leastsqr (). Для получения дополнительной информации см. Учебник по scipy optimize

Основная схема заключается в создании функции, которая оценивает функцию модели по независимой координате и возвращает массив пустот, содержащий остатки, разницу между моделью и наблюдениями в каждой точке. Функция leastsq () находит значения набора параметров, который минимизирует сумму квадратов невязок.

В качестве предупреждения мы отмечаем, что подгонка может быть чувствительной к первоначальному предположению. Имитация отжига часто используется для нахождения вероятного глобального минимума и предоставления приблизительной оценки параметров подгонки перед уточнением подгонки. Значения, используемые здесь для первоначального предположения, приведены только для условных целей.

from scipy.optimize import leastsq
import numpy as np

# Read the data files    
Pm = np.array( [ float(v) for v in open( "impulse_data.txt", "r" ).readlines() ] )
print type(Pm), Pm.shape

tm = np.array( [ float(v) for v in open( "Impulse_time_axis.txt", "r" ).readlines() ] )
print type(tm), tm.shape

output = np.array( [ float(v) for v in open( "Output_data.txt", "r" ).readlines() ] )
print type(output), output.shape

tout = np.array( [ float(v) for v in open( "Output_time_axis.txt", "r" ).readlines() ] )
print type(tout), tout.shape

# Define the function that calculates the residuals
def residuals(  coeffs, output, tm ):
    dt, alpha1, alpha2 = coeffs
    res = np.zeros( tm.shape )
    for n,t in enumerate(tm):
        # integrate to "t"
        value = sum( Pm[:n]*(1-np.e**( -(alpha1 * (tm[:n]+dt))) )*np.e**(-(alpha2 * (tm[:n]+dt))) )
        # retrieve output at t+dt
        n1 = (np.abs(tout - (t+dt) )).argmin()
        # construct the residual
        res[n] =  output[n1] - value
    return res

# Initial guess for the parameters
x0 = (10.,1.,1.)

# Run the least squares routine
coeffs, flag = leastsq( residuals, x0, args=(output,tm) )

# Report the results
print( coeffs )
print( flag )