Сигмоидальная регрессия с scipy, numpy, python и т.д.

Ответ 1

Использование scipy.optimize.leastsq:

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

def sigmoid(p,x):
    x0,y0,c,k=p
    y = c / (1 + np.exp(-k*(x-x0))) + y0
    return y

def residuals(p,x,y):
    return y - sigmoid(p,x)

def resize(arr,lower=0.0,upper=1.0):
    arr=arr.copy()
    if lower>upper: lower,upper=upper,lower
    arr -= arr.min()
    arr *= (upper-lower)/arr.max()
    arr += lower
    return arr

# raw data
x = np.array([821,576,473,377,326],dtype='float')
y = np.array([255,235,208,166,157],dtype='float')

x=resize(-x,lower=0.3)
y=resize(y,lower=0.3)
print(x)
print(y)
p_guess=(np.median(x),np.median(y),1.0,1.0)
p, cov, infodict, mesg, ier = scipy.optimize.leastsq(
    residuals,p_guess,args=(x,y),full_output=1,warning=True)  

x0,y0,c,k=p
print('''\
x0 = {x0}
y0 = {y0}
c = {c}
k = {k}
'''.format(x0=x0,y0=y0,c=c,k=k))

xp = np.linspace(0, 1.1, 1500)
pxp=sigmoid(p,xp)

# Plot the results
plt.plot(x, y, '.', xp, pxp, '-')
plt.xlabel('x')
plt.ylabel('y',rotation='horizontal') 
plt.grid(True)
plt.show()

дает

с сигмовидными параметрами

x0 = 0.826964424481
y0 = 0.151506745435
c = 0.848564826467
k = -9.54442292022

Обратите внимание, что для более новых версий scipy (например, 0.9) есть функция scipy.optimize.curve_fit, которая проще в использовании, чем leastsq, Соответствующее обсуждение сигмоидов, использующих curve_fit, можно найти здесь.

Изменить: добавлена ​​функция resize, чтобы исходные данные можно было масштабировать и сдвигать, чтобы они соответствовали любому требуемому ограничивающему прямоугольнику.

"ваше имя появляется как писатель от скудной документации"

ОТКАЗ ОТ ОТВЕТСТВЕННОСТИ: Я не писатель скупой документации. Я просто пользователь, и новичок в этом. Большая часть того, что я знаю о leastsq, происходит от чтения этого учебника, написанного Трэвисом Олифантом.

1.) Вызывает ли lesssq() остатки(), которые затем возвращают разницу между входным вектором y и y-вектор, возвращаемый сигмоидом() функция?

Да! точно.

Если да, то как это объясняет разница в длинах ввода y-вектора и y-вектора, возвращаемого функция sigmoid()?

Длины одинаковы:

In [138]: x
Out[138]: array([821, 576, 473, 377, 326])

In [139]: y
Out[139]: array([255, 235, 208, 166, 157])

In [140]: p=(600,200,100,0.01)

In [141]: sigmoid(p,x)
Out[141]: 
array([ 290.11439268,  244.02863507,  221.92572521,  209.7088641 ,
        206.06539033])

Одна из замечательных вещей о Numpy заключается в том, что она позволяет писать "векторные" уравнения, которые работают на всех массивах.

y = c / (1 + np.exp(-k*(x-x0))) + y0

может выглядеть так, как будто он работает с float (на самом деле это было бы), но если вы сделаете x массив numpy и c, k, x0, y0 float, то уравнение определяет y, чтобы быть массивом numpy той же формы, что и x. Таким образом, sigmoid(p,x) возвращает массив numpy. Более полное объяснение того, как это работает в numpybook (требуется чтение для серьезных пользователей numpy).

2.) Похоже, я могу назвать lesssq() для любого математического уравнения, если я доступ к этому математическому уравнению через функция остатков, которая, в свою очередь, вызывает математическую функцию. Это правда?

True. leastsq пытается минимизировать сумму квадратов остатков (разностей). Он ищет пространство параметров (все возможные значения p), ища p, который минимизирует эту сумму квадратов. x и y, отправленные в residuals, являются вашими необработанными значениями данных. Они исправлены. Они не меняются. Это p (параметры в сигмоидной функции), которые leastsq пытается свести к минимуму.

3.) Кроме того, я замечаю, что p_guess имеет такое же количество элементов, что и p. Есть ли это означает, что четыре элемента p_guess соответствуют по порядку, соответственно, с возвращенными значениями через x0, y0, c и k?

Точно так! Как и метод Ньютона, leastsq требует начального предположения для p. Вы указываете его как p_guess. Когда вы видите

scipy.optimize.leastsq(residuals,p_guess,args=(x,y))

вы можете думать, что как часть алгоритма наименьшего квадрата (действительно алгоритм Левенбурга-Марквардта) в качестве первого прохода, наименьший вызов вызывает residuals(p_guess,x,y). Обратите внимание на визуальное сходство между

(residuals,p_guess,args=(x,y))

и

residuals(p_guess,x,y)

Это может помочь вам запомнить порядок и смысл аргументов leastsq.

residuals, например sigmoid возвращает массив numpy. Значения в массиве квадратичны, а затем суммируются. Это номер, который нужно бить. p_guess изменяется тогда как leastsq ищет набор значений, который минимизирует residuals(p_guess,x,y).

4.) Является ли p, который отправляется как аргумент остаткам() и sigmoid() выполняет ту же функцию p, что будет выводиться с помощью наименьших квадратов(), а Функция minsq() использует это значение p внутренне, прежде чем возвращать его?

Ну, не совсем. Как вы знаете, p_guess изменяется, поскольку leastsq выполняет поиск значения p, которое минимизирует residuals(p,x,y). p (er, p_guess), который отправляется на leastsq, имеет ту же форму, что и p, которая возвращается leastsq. Очевидно, что значения должны быть разными, если вы не догадываетесь:)

5.) Может p и p_guess иметь любое количество элементов, в зависимости от сложность используемого уравнения в качестве модели, если число элементы в p равны числу элементов в p_guess?

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

Ответ 2

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

Я не программист на Python, поэтому я не знаю, имеет ли numpy более общую кривую рутина. Если вам нужно катиться самостоятельно, возможно, эта статья в Логистическая регрессия даст вам несколько идей.

Ответ 3

Для логистической регрессии в Python scikits-learn предоставляет высокопроизводительный код подстановки:

http://scikit-learn.sourceforge.net/modules/linear_model.html#logistic-regression