Метод наименьших квадратов для суммы функций - PullRequest
Новая
       39

Метод наименьших квадратов для суммы функций

0 голосов
/ 26 ноября 2018

Я бы хотел использовать функцию curve_fit из модуля scipy.optimize, чтобы определить амплитуды, частоты, фазы суммы синусоидальных функций (и одну y0).Это легко сделать, когда я знаю несколько синусов для использования.Например, когда я знаю две частоты из ДПФ (дискретного преобразования Фурье): 1.152 и 0.432, я могу определить функцию: 

def func(x, amp1, amp2, freq1 , freq2, phase1, phase2, y0):
    return amp1*np.sin(freq1*x + phase1) + amp2*np.sin(freq2*x + phase2) + y0

Затем, используя curve_fit и ограничив интервалы частотЯ могу найти подходящий пример: 

param, _ = curve_fit(func, t, data, bounds=([-np.inf, -np.inf, 1.14, 0.43, -np.inf, -np.inf, -np.inf], [np.inf, np.inf, 1.16, 0.44, np.inf, np.inf, np.inf]))

Выглядит отлично: enter image description here

Но в этом случае я подготовил данные и знал,ряд частот.Знаете ли вы, как определить func только один раз и обрабатывать все случаи (например, пять функций синуса)?Я попытался поместить параметры в списки, например, amp = [amp1, amp2, ... ], и я перебрал их длину.Но существует проблема определения bounds для списков параметров.bounds очень важно для обеспечения реальности модели.

Решение не должно основываться на curve_fit.

1 Ответ

0 голосов
/ 27 ноября 2018

Предполагая, что вы знаете частоты заранее, проблема проста.Вы можете установить нижнюю границу на 0 и установить верхнюю границу на 2 * pi * freq для частоты.Для усилителей задайте любое число (или np.inf, если вы не хотите границы).

Вы можете сформулировать функцию в виде lambda x, amp1, phase1, amp2, phase2... : y, curve_fit может принимать функцию неопределенного числа аргументов до тех пор, покапо мере того как вы вводите правильное начальное предположение.

Пример кода для пяти частот: 

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

x = np.linspace(0,10,60)
w = [1,2,3,4,5]

a = [1,4,2,3,0.1]
x0 = [0,1,0,1,0.5]

y = np.sum(a_i * np.sin(w_i * x - x0_i) for w_i, a_i, x0_i in zip(w,a, x0))  #base_data
yr = y + np.random.normal(0,0.5, size=x.size)   #noisy data

def func(x, *args):
    """ function of the form lambda x, amp1, phase1, amp2, phase2...."""
    return np.sum(a_i * np.sin(w_i * (x-x0)) for w_i, a_i, x0
                  in zip(w,args[::2], args[1::2]))

ubounds = np.zeros(len(w) * 2)
ubounds[::2] = 10   #setting amp max value to 10 (arbitrary)
ubounds[1::2] = np.asarray(w) * 2 * np.pi
p0 = [0] * 10   # note p0 size
popt, pcov = curve_fit(func, x, yr, p0, bounds=(0, ubounds))
amps, phases = popt[::2], popt[1::2]

plt.plot(x,func(x, *popt))
plt.plot(x,yr, 'go')
...