Подгонка ступенчатой функции
Я пытаюсь соответствовать шаговой функции, используя scipy.optimize.leastsq. Рассмотрим следующий пример:
import numpy as np
from scipy.optimize import leastsq
def fitfunc(p, x):
y = np.zeros(x.shape)
y[x < p[0]] = p[1]
y[p[0] < x] = p[2]
return y
errfunc = lambda p, x, y: fitfunc(p, x) - y # Distance to the target function
x = np.arange(1000)
y = np.random.random(1000)
y[x < 250.] -= 10
p0 = [500.,0.,0.]
p1, success = leastsq(errfunc, p0, args=(x, y))
print p1
параметры - местоположение шага и уровень с обеих сторон. Странно то, что первый свободный параметр никогда не меняется, если вы запустите, то scipy даст
[ 5.00000000e+02 -4.49410173e+00 4.88624449e-01]
когда первый параметр будет оптимальным при установке на 250, а второй на -10.
У кого-нибудь есть понимание того, почему это может не сработать и как заставить его работать?
Если я бегу
print np.sum(errfunc(p1, x, y)**2.)
print np.sum(errfunc([250.,-10.,0.], x, y)**2.)
Я нахожу:
12547.1054663
320.679545235
где первое число - это то, что находит leastsq, а второе - значение для фактической оптимальной функции, которую он должен найти.
5 ответов
Оказывается, что подгонка будет намного лучше, если я добавлю аргумент epsfcn= в leastsq:
p1, success = leastsq(errfunc, p0, args=(x, y), epsfcn=10.)
и результат
[ 248.00000146 -8.8273455 0.40818216]
Мое основное понимание состоит в том, что первый свободный параметр должен быть перемещен больше, чем расстояние между соседними точками, чтобы повлиять на квадрат невязок, и epsfcn имеет какое-то отношение к тому, какие большие шаги нужно использовать, чтобы найти градиент, или что-то подобное.
Предлагаю аппроксимировать шаговую функцию. Вместо бесконечного наклона в "точке изменения" сделайте его линейным на расстоянии одного x (1,0 в примере). Например, если параметр x, xp, для функции определен как средняя точка на этой линии, тогда значение в xp-0.5 является более низким значением y, а значение в xp+0.5 является более высоким значением y, а промежуточные значения функции в интервал [хр-0,5; xp+0.5] - линейная интерполяция между этими двумя точками.
Если можно предположить, что ступенчатая функция (или ее приближение) переходит от более низкого значения к более высокому значению, то я думаю, что первоначальное предположение для последних двух параметров должно быть наименьшим значением y и наибольшим значением y соответственно вместо 0,0 и 0.0.
У меня есть 2 исправления:
1) np.random.random () возвращает случайные числа в диапазоне от 0,0 до 1,0. Таким образом, среднее значение составляет +0,5, а также является значением третьего параметра (вместо 0,0). И тогда второй параметр равен -9,5 (+0,5 - 10,0) вместо -10,0.
таким образом
print np.sum(errfunc([250.,-10.,0.], x, y)**2.)
должно быть
print np.sum(errfunc([250.,-9.5,0.5], x, y)**2.)
2) В исходной функции fitfunc() одно значение y становится 0.0, если x точно равно p[0]. Таким образом, в этом случае это не ступенчатая функция (больше похоже на сумму двух ступенчатых функций). Например, это происходит, когда начальное значение первого параметра равно 500.
Я не думаю, что метод наименьших квадратов - это способ приблизиться к шагу. Я не верю, что это даст вам удовлетворительное описание разрыва. Наименьшие квадраты не были бы моей первой мыслью при атаке на эту проблему.
Почему бы вам не использовать вместо этого приближение ряда Фурье? Вы всегда будете застревать с феноменом Гиббса на разрыве, но остальная часть функции может быть аппроксимирована так, как вы и ваш процессор можете себе позволить.
Для чего именно вы собираетесь это использовать? Некоторый контекст может помочь.
Скорее всего, ваша оптимизация застряла в локальных минимумах. Я не знаю, как на самом деле работает leastsq, но если вы дадите ему начальную оценку (0, 0, 0), он тоже там застрянет.
Вы можете проверить градиент при первоначальной оценке численно (оценить в +/- эпсилон для очень маленького эпсилона и разделить на 2* эпсилон, взять разницу), и я держу пари, что это будет что-то около 0.
используйте statsmodel ols. ols использует обычный метод наименьших квадратов для подбора кривой