Использование максимальной вероятности для получения коэффициентов регрессии в Python

Question

В качестве обучающего упражнения я пытаюсь оценить параметры регрессии, используя метод MLE в Python.

Из того, что я собрал, функция логарифмического правдоподобия сводится к максимизации следующего :

введите описание изображения здесь

Итак, мне нужно взять частные производные относительно точки пересечения и наклона, установив для каждой ноль, и это должно дать мне коэффициенты.

Я пытался приблизиться это с использованием sympy следующим образом:

from sympy import *
b = Symbol('b')# beta1
a = Symbol('a')# intercept
x = Symbol('x', integer=True)
y = Symbol('y', integer=True)
i = symbols('i', cls=Idx)

x_values = [2,3,2]
y_values = [1,2,3]
n = len(x_values)-1

function = summation((Indexed('y',i) - a+b*Indexed('x',i))**2, (i, 0, n))

partial_intercept = function.diff(a)

print(partial_intercept)
# 6*a - 2*b*x[0] - 2*b*x[1] - 2*b*x[2] - 2*y[0] - 2*y[1] - 2*y[2]

intercept_f = lambdify([x, y], partial_intercept)

inter = solve(intercept_f(x_values, y_values), a)

print(inter)
# [7*b/3 + 2]

Я ожидал бы одно значение для наклона, такое, что переменная 'b' исчезнет. Однако я вижу, что это было бы невозможно, учитывая, что переменная b все еще присутствует в моем производном уравнении.

Есть ли у кого-нибудь совет, где я ошибаюсь?

Спасибо!

Редактировать: Исправлена ​​опечатка в кодовом блоке

Answer 1

Выражение 7*b/2 + 2 в конце говорит вам, что мы должны удовлетворить a = 7*b/2 + 2, это зависит от количества b.

Вы должны решить для a и b в качестве системы одновременно.

В следующем коде я нахожу соотношение, которое a и b должны удовлетворять, и решайте их одновременно.

from sympy import *

b = Symbol('b')# beta1
a = Symbol('a')# intercept
x = Symbol('x', integer=True)
y = Symbol('y', integer=True)
i = symbols('i', cls=Idx)

x_values = [2,3,2]
y_values = [1,2,3]
n = len(x_values)-1

function = summation((Indexed('y',i) - a+b*Indexed('x',i))**2, (i, 0, n))

partial_intercept = function.diff(a)

print(partial_intercept)
# 6*a - 2*b*x[0] - 2*b*x[1] - 2*b*x[2] - 2*y[0] - 2*y[1] - 2*y[2]

intercept_f = lambdify([x, y], partial_intercept)

inter = solve(intercept_f(x_values, y_values), a)

print(inter)
#[7*b/3 + 2]

partial_gradient = function.diff(b)

print(partial_gradient)
# 6*a - 2*b*x[0] - 2*b*x[1] - 2*b*x[2] - 2*y[0] - 2*y[1] - 2*y[2]

intercept_f = lambdify([x, y], partial_gradient)

inter2 = solve(intercept_f(x_values, y_values), b)
print(inter2)

ans = solve([a-inter[0], b-inter2[0]])
print(ans)

Вот результаты:

6*a - 2*b*x[0] - 2*b*x[1] - 2*b*x[2] - 2*y[0] - 2*y[1] - 2*y[2]
[7*b/3 + 2]
2*(-a + b*x[0] + y[0])*x[0] + 2*(-a + b*x[1] + y[1])*x[1] + 2*(-a + b*x[2] + y[2])*x[2]
[7*a/17 - 14/17]
{a: 2, b: 0}

a должно быть установлено как 2, а b должно быть установлено как 0.