Как аппроксимировать cos^2(x) гауссианами в SciPy и настроить curve_fit

Как аппроксимировать cos^2(x) суммой гауссиан в SciPy: фиксируем число пиков через find_peaks и задаем стартовые оценки, чтобы curve_fit сходился устойчиво.

Подгонка сигнала cos^2(x) суммой гауссовых пиков кажется простой задачей, однако наивное применение scipy.optimize.curve_fit может привести к болезненно плохому результату. Проблема не в данных и не в выборе функции, а в том, как заданы и как инициализированы параметры модели.

Постановка задачи

Цель — аппроксимировать cos^2(kx) набором гауссиан, равномерно смещённых на постоянный шаг. Приведённый ниже код формирует целевую кривую, базис из гауссиан и пытается одновременно подобрать все параметры, включая количество пиков.

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

freq = 1
mult = 3
x_lo = -mult * np.pi / (2 * freq)
x_hi =  mult * np.pi / (2 * freq)
grid = np.linspace(x_lo, x_hi, 1000)

target = np.cos(freq * grid) ** 2

def bell(xv, center, sigma):
    return np.exp(- (xv - center) ** 2 / (2 * sigma ** 2))

def bell_stack(xv, center, sigma, spacing, count):
    total = 0
    for idx in range(int(count)):
        pos = center + idx * spacing
        total += bell(xv, pos, sigma)
    return total

est = curve_fit(bell_stack, grid, target)[0]

c0 = est[0]
sg = est[1]
gap = est[2]
num = int(est[3])

approx = 0
for idx in range(num):
    pos = c0 + idx * gap
    approx += bell(grid, pos, sg)

fig, ax = plt.subplots()
ax.plot(grid, target, label="cosine function")
ax.plot(grid, approx, label="gaussian approximation to cosine")
ax.set_xlabel("x")
ax.set_ylabel("Amplitude")
plt.grid()
plt.legend()
plt.show()

Почему подгонка срывается

Два фактора работают против решателя. Во‑первых, модель включает параметр, отвечающий за число гауссиан в сумме. Этот параметр по своей природе дискретный и недифференцируемый. Если оптимизатор изменит его на сколь угодно малую величину, выход не изменится. Локальные методы, такие как curve_fit, опираются на локальную градиентную информацию, поэтому не могут получить осмысленное обновление для такого параметра. Во‑вторых, у функции ошибки много локальных минимумов. Локальный оптимизатор легко застревает, и старт с типовой инициализации единицами делает это особенно болезненным для количества пиков и их шага.

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

Практическое решение: убрать недифференцируемый параметр и задать осмысленные начальные значения

Идея — исключить дискретный параметр из подгонки и задать информированные начальные оценки. Надёжный способ здесь — обнаружить количество и позиции пиков в целевом сигнале. С помощью scipy.signal.find_peaks можно определить, сколько пиков использовать, оценить расстояние между ними и выбрать положение первого пика. Затем передайте в curve_fit модель с фиксированным числом пиков и позвольте ей оптимизировать только непрерывные параметры.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
from functools import partial

freq = 1
mult = 3
x_lo = -mult * np.pi / (2 * freq)
x_hi =  mult * np.pi / (2 * freq)
grid = np.linspace(x_lo, x_hi, 1000)

target = np.cos(freq * grid) ** 2

def bell(xv, center, sigma):
    return np.exp(- (xv - center) ** 2 / (2 * sigma ** 2))

def bell_stack(xv, center, sigma, spacing, count):
    total = 0
    for idx in range(int(count)):
        pos = center + idx * spacing
        total += bell(xv, pos, sigma)
    return total

peak_idx, _ = find_peaks(target)
num = len(peak_idx)
init_gap = np.mean(np.diff(grid[peak_idx]))
init_c0 = grid[peak_idx[0]]

model = partial(bell_stack, count=num)
fit_params = curve_fit(model, grid, target, p0=(init_c0, 1, init_gap))[0]

c0, sg, gap = fit_params
approx = model(grid, *fit_params)

fig, ax = plt.subplots()
ax.plot(grid, target, label="cosine function")
ax.plot(grid, approx, label="gaussian approximation to cosine")
ax.set_xlabel("x")
ax.set_ylabel("Amplitude")
plt.grid()
plt.legend()
plt.show()

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

Почему это важно

Даже при разумной модели оптимизация срывается, если параметризация нарушает предпосылки метода. curve_fit — локальный алгоритм; ему нужны градиенты и приличные начальные оценки. Введение дискретных выборов в градиентную процедуру или надежда на значения по умолчанию ведут к плохой сходимости. Напротив, предварительная оценка структурных свойств сигнала и ограничение подгонки непрерывными переменными дают алгоритму то, что он действительно умеет оптимизировать.

Выводы

При аппроксимации структурированных сигналов вроде cos^2(x) суммами гауссиан важно учитывать возможности оптимизатора. Сохраняйте параметрическое пространство дифференцируемым — фиксируйте дискретные количества заранее — и инициализируйте подгонку сведениями из данных: числом пиков, их шагом и координатой первого пика. С этими поправками знакомые инструменты начинают работать предсказуемо и дают нужную аппроксимацию.