Подгонка ступенчатой ​​функции

Question

Я пытаюсь установить пошаговую функцию, используя 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, а второе - значение для фактической оптимальной функции, которую он должен найти.

Answer 1

Оказывается, что подгонка будет намного лучше, если я добавлю аргумент epsfcn = в leastsq:

p1, success = leastsq(errfunc, p0, args=(x, y), epsfcn=10.)

и результат

[ 248.00000146   -8.8273455     0.40818216]

Мое основное понимание состоит в том, что первый свободный параметр должен быть перемещен больше, чем расстояние между соседними точками, чтобы повлиять на квадрат невязок, и epsfcn имеет какое-то отношение к тому, какие большие шаги нужно использовать, чтобы найти градиент, или к похоже.

Answer 2

Предлагаю аппроксимировать пошаговую функцию.Вместо бесконечного наклона в «точке изменения» сделайте его линейным на расстоянии одного 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.

Answer 3

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

Почему бы вам не использовать вместо этого приближение ряда Фурье? Вы всегда будете зависеть от явления Гиббса на разрыве, но остальная часть функции может быть аппроксимирована так, как вы и ваш процессор можете себе позволить.

Для чего именно вы собираетесь это использовать? Некоторый контекст может помочь.

Answer 4

Скорее всего, ваша оптимизация застряла в локальных минимумах . Я не знаю, как на самом деле работает leastsq, но если вы дадите ему начальную оценку (0, 0, 0), он тоже там застрянет.

Вы можете проверить градиент при первоначальной оценке численно (оценить в +/- epsilon для очень маленького эпсилона и разделить bei 2 * epsilon, взять разницу), и я держу пари, что он будет примерно равен 0.