Ошибка нелинейного метода наименьших квадратов синусоидальной функции с экспоненциальным затуханием

Мои нелинейные данные аппроксимируются методом наименьших квадратов по формуле Asin(wt+phase)exp(-decay*t), при этом omega(w) остается постоянным. Я пробовал несколько подходов без успеха.

Ниже мой код

import numpy as np
from numpy import loadtxt
import lmfit
np.random.seed(2)
x = np.linspace(0, 10, 101)
decay = 4.5
shift = 0
amp = 0.0015
y = amp * np.sin(x*5+shift) * np.exp(-x*decay)
yn = y + np.random.normal(size=y.size, scale=0.450)


def resid(params, x, ydata):
    decay = params['decay'].value
    shift = params['shift'].value
    amp = params['amp'].value

    y_model = amp * np.sin(x*5+shift) * np.exp(-x*decay)
    return y_model - ydata
params = lmfit.Parameters()
params.add('shift', 0.0, min=-np.pi, max=np.pi)
params.add('amp', 0.0015, min=0, max=0.02)
params.add('decay', 4.0, min=0, max=10.0)
fit = lmfit.minimize(resid, params, args=(x, yn), method='differential_evolution')
print("\n\n# Fit using differential_evolution:")
lmfit.report_fit(fit)
plt.plot(x, y, 'ko', lw=2)
plt.plot(x, yn+fit.residual, 'b--', lw=2)
plt.legend(['data', 'leastsq', 'diffev'], loc='upper left')
plt.show()

python lmfit
person Prakash Gyawali    schedule 22.10.2020    source источник
comment
Прочтите, как создать MCVE: stackoverflow.com/help/minimal-reproducible-example. помогите здесь, в Stackoverflow, но у большинства из нас не так много времени, чтобы помочь. С MCVE вы упрощаете для нас воспроизведение ваших проблем, чтобы мы могли уделять больше времени тому, что у вас есть.   -  person Joooeey    schedule 22.10.2020
comment
Пожалуйста, предоставьте точный ввод и точный ожидаемый результат. Также ссылка на теорию.   -  person Gulzar    schedule 22.10.2020
comment
Глядя на ваши данные: ваш распад огромен по сравнению с амплитудой, что приводит к проблемам с небольшими значениями. У вас есть около 10 точек данных, которые ведут себя как ваша ожидаемая формула, но 90 точек, которые в основном равны нулю, и вдобавок к этому вы маскируете эти значения случайным шумом. Ваш подход работает так же хорошо для меньших значений затухания и больших значений усиления.   -  person Mr. T    schedule 22.10.2020

Ответы (1)


arrow_upward
0
arrow_downward

В комментариях указано, что распад у вас очень тяжелый. Но кроме того, шум, который вы добавляете к yn, огромен по сравнению с вашей затухающей синусоидой. Вы и другие, возможно, пропустили это, потому что вы не строите массив yn, который вам подходит, а строите y без добавленного шума data.

Если вы нарисуете данные, которые вы фактически используете в подгонке:

plt.plot(x, yn, 'ko', lw=2)
plt.show()

вы увидите это:

введите здесь описание изображения

Вы также помечаете самый плохой вариант, который вы получаете, как leastsq, хотя на самом деле вы использовали differential_evolution.

Если вы уменьшите шум и затухание в своих тестовых данных и действительно соответствуете leastsq, вы можете получить что-то вроде этого:

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

np.random.seed(2)
x = np.linspace(0, 10, 101)

decay = 0.6
shift = 0
amp = 0.025
y = amp * np.sin(x*5+shift) * np.exp(-x*decay)
yn = y + np.random.normal(size=y.size, scale=0.001)

def resid(params, x, ydata):
    decay = params['decay'].value
    shift = params['shift'].value
    amp = params['amp'].value

    y_model = amp * np.sin(x*5+shift) * np.exp(-x*decay)
    return y_model - ydata

params = lmfit.Parameters()
params.add('shift', 0.0, min=-np.pi, max=np.pi)
params.add('amp', 0.01)
params.add('decay', 0.2, min=0, max=50)
fit = lmfit.minimize(resid, params, args=(x, yn))


print("\n\n# Fit using differential_evolution:")
lmfit.report_fit(fit)
plt.plot(x, yn, 'ko', lw=2)
plt.plot(x, yn+fit.residual, 'b--', lw=2)
plt.legend(['data', 'leastsq'], loc='upper right')
plt.show()

который даст отчет о

# Fit using differential_evolution:
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 21
    # data points      = 101
    # variables        = 3
    chi-square         = 1.0894e-04
    reduced chi-square = 1.1117e-06
    Akaike info crit   = -1381.71974
    Bayesian info crit = -1373.87437
[[Variables]]
    shift:  0.00796328 +/- 0.01999031 (251.03%) (init = 0)
    amp:    0.02448871 +/- 7.5756e-04 (3.09%) (init = 0.01)
    decay:  0.59826922 +/- 0.02587107 (4.32%) (init = 0.2)
[[Correlations]] (unreported correlations are < 0.100)
    C(amp, decay)   =  0.725
    C(shift, amp)   = -0.147
    C(shift, decay) = -0.105

и график введите здесь описание изображения

Я думаю, что вывод состоит в том, что ваша функция невязки и настройка проблемы в порядке, но ваш очень шумный набор тестовых данных не очень хорошо представлен этой затухающей синусоидой.

person M Newville    schedule 23.10.2020