Решение проблемы Пончика с помощью логистической регрессии
В данной теме мы поговорим об одной практической проблеме, касающейся логистической регрессии. Она называется «проблемой пончика». Вы поймёте смысл этого названия, после того как мы получим графическое её отображение.
В этом примере мы используем куда больше данных, чтобы видимый эффект был более существенным.
import numpy as np
import matplotlib.pyplot as plt
N = 1000
D = 2
Суть в том, что у нас есть два радиуса – внутренний и внешний, равные 5 и 10 соответственно. Половина данных равномерно распределяется по внутреннему радиусу, при этом полярные координаты мы сразу преобразуем в прямоугольные.
R_inner = 5
R_outer = 10
R1 = np.random.randn(N/2) + R_inner
theta = 2*np.pi*np.random.random(N/2)
X_inner = np.concatenate([[R1 * np.cos(theta)], [R1 * np.sin(theta)]]).T
То же самое делаем с внешним радиусом, также с равномерно распределёнными данными.
R2 = np.random.randn(N/2) + R_outer
theta = 2*np.pi*np.random.random(N/2)
X_outer = np.concatenate([[R2 * np.cos(theta)], [R2 * np.sin(theta)]]).T
Теперь определим общее X и, конечно, целевую переменную.
X = np.concatenate([ X_inner, X_outer ])
T = np.array([0]*(N/2) + [1]*(N/2))
И отобразим это графически, чтобы вы могли увидеть, как всё это выглядит.
plt.scatter(X[:,0], X[:,1], c=T)
plt.show()
На диаграмме вы можете наглядно увидеть «проблему пончика». Линейный разделитель логистической регрессии здесь мало применим, поскольку не существует прямой, способной разделить эти классы. Но я покажу вам, что на самом деле эта проблема вполне разрешима.
Как обычно, создадим столбец из единиц для члена смещения.
ones = np.array([[1]*N]).T
Фокус в обходе «проблемы пончика» в том, что мы создаём ещё один столбец, характеризующий радиус точки. Благодаря этому вы сможете увидеть, как сделать наши данные линейно разделяемыми.
r = np.zeros((N,1))
for i в xrange(N):
r[i] = np.sqrt(X[i,:].dot(X[i,:]))
Xb = np.concatenate((ones, r, X), axis=1)
Установим опять случайным образом весовые коэффициенты.
w = np.random.rand(D+2)
Остальной же код из предыдущих лекций, вроде функции сигмоиды, остаётся без изменений, поэтому не будем его переписывать.
z = Xb.dot(w)
def sigmoid(z):
return 1/(1 + np.exp(-z))
Y = sigmoid(z)
def cross_entropy(T, Y):
E = 0
for i in xrange(N):
if T[i] == 1:
E -= np.log(Y[i])
else:
E -= np.log(1 – Y[i])
return E
Установим специально подобранное мною значение коэффициента обучения, а также количество итераций. В общем же случае нужно немного поэкспериментировать, чтобы найти правильные значения, или же использовать что-то вроде перекрёстной проверки. Кроме того, давайте покажем, как весовые коэффициенты изменяются с течением времени при проходе каждой итерации. Зададим выдавать результат каждые 100 проходов и используем градиентный спуск с регуляризацией.
learning_rate = 0.0001
error = []
for i in xrange(5000):
e = cross_entropy(T, Y)
error.append(e)
if i % 100 == 0:
print e
w += learning_rate * ( np.dot((T – Y).T, Xb) – 0.01*w)
Y = sigmoid(Xb.dot(w))
И выведем на экран изменение ошибки с течением времени, окончательные весовые коэффициенты и рейтинг классификации.
plt.plot(error)
plt.title(‘’Cross-entropy’’)
print ‘’Final w’’, w
print ‘’Final classification rate:’’, 1 – np.abs(T – np.round(Y)). sum / N
Запустим программу. Итак, мы вновь видим «пончик» и прохождение итераций.
В конце концов мы получаем очень неплохой коэффициент классификации. Значения весовых коэффициентов также весьма интересны – коэффициенты при x и y почти равны нулю. Таким образом, классификация на самом деле не зависит от прямоугольной системы координат, что и показала наша модель.
Но мы обнаруживаем, что она зависит от члена смещения, поэтому если мы поставим слишком малый радиус, то получим отрицательное значение члена смещения, что подталкивает классификацию к нулю. Если же радиус чересчур велик, то это подталкивает значение коэффициента классификации к единице.
Вот таким образом можно решить проблему, используя логистическую регрессию.