Взвешенные наименьшие квадраты в Statsmodels против Numpy?

Я пытаюсь воспроизвести функциональность Statsmodels функция взвешивания наименьших квадратов (WLS) с обычный метод наименьших квадратов (OLS) (т. е. Numpy называет OLS просто методом наименьших квадратов).

Другими словами, я хочу вычислить WLS в Numpy. Я использовал эту запись Stackoverflow в качестве справки, но возникают совершенно другие значения R² при переходе от Statsmodel к Numpy.

Возьмите следующий пример кода, который повторяет это:

import numpy as np
import statsmodels.formula.api as smf
import pandas as pd

# Test Data
patsy_equation = "y ~ C(x) - 1" # Use minus one to get ride of hidden intercept of "+ 1"
weight = np.array([0.37, 0.37, 0.53, 0.754])
y = np.array([0.23, 0.55, 0.66, 0.88])
x = np.array([3, 3, 3, 3])
d = {"x": x.tolist(), "y": y.tolist()}
data_df = pd.DataFrame(data=d)

# Weighted Least Squares from Statsmodel API
statsmodel_model = smf.wls(formula=patsy_equation, weights=weight, data=data_df)
statsmodel_r2 = statsmodel_model.fit().rsquared      

# Weighted Least Squares from Numpy API
Aw = x.reshape((-1, 1)) * np.sqrt(weight[:, np.newaxis]) # Multiply two column vectors
Bw = y * np.sqrt(weight)
numpy_model, numpy_resid = np.linalg.lstsq(Aw, Bw, rcond=None)[:2]
numpy_r2 = 1 - numpy_resid / (Bw.size * Bw.var())

print("Statsmodels R²: " + str(statsmodel_r2))
print("Numpy R²: " + str(numpy_r2[0]))

После запуска такого кода я получаю следующие результаты:

Statsmodels R²: 2.220446049250313e-16
Numpy R²: 0.475486515775414

Здесь явно что-то не так! Кто-нибудь может указать здесь на мои недостатки? Я скучаю по формуле патси?


person Code Doggo    schedule 25.05.2018    source источник
comment
Вы убедились, что два метода дают одно и то же решение задачи взвешенных наименьших квадратов? Если нет, то нет смысла проверять R².   -  person Warren Weckesser    schedule 25.05.2018
comment
Все ваши значения x одинаковы! Как вы ожидаете подогнать линию к этим данным?   -  person Warren Weckesser    schedule 25.05.2018
comment
Если подумать, очевидно, вы хотите, чтобы линия проходила через начало координат. В этом случае несколько значений y при одном значении x дают четко определенную проблему. Вы можете вычислить наклон, взяв средневзвешенное значение этих значений y, а затем разделив его на x: np.average(y, weights=weight)/3 дает 0.21441370223978917, что согласуется с numpy_model.   -  person Warren Weckesser    schedule 25.05.2018
comment
@WarrenWeckesser Как вычисление наклона находит значение R²? Я вычислил неверное значение numpy_R²? Кроме того, вы сказали, что делите на x, однако вы только что разделили на моду вектора x (т.е. мода = наиболее повторяющееся число). Не могли бы вы немного пояснить свои расчеты?   -  person Code Doggo    schedule 25.05.2018
comment
@WarrenWeckesser Глядя на параметры модели Statsmodels WLS, я получаю один параметр 0.643241..., тогда как при использовании Numpy я получаю 0.2144137.... Модели не согласуются, поэтому я предполагаю, что неправильно вычисляю numpy_model?   -  person Code Doggo    schedule 25.05.2018
comment
Как при вычислении наклона определяется значение R²? Конечно, нет. Я просто исправлял свой предыдущий комментарий и вопрос о том, что все значения x одинаковы. На ваш вопрос о R² вам, вероятно, понадобится кто-то, более знакомый с реализацией WLS в statsmodels, чтобы дать ответ.   -  person Warren Weckesser    schedule 25.05.2018
comment
@WarrenWeckesser А, хорошо! Я буду продолжать изучать эту проблему! В любом случае спасибо!   -  person Code Doggo    schedule 25.05.2018
comment
@WarrenWeckesser Я вижу здесь проблему. В patsy_equation у меня есть C(x). Я рассматриваю свои значения x как категориальные, а не как числовые переменные. Если я просто поставлю x вместо C(x), то получу ту же модель. Вы хоть представляете, что происходит с переходом от числового к категориальному? Что еще более важно, как меняется математика?   -  person Code Doggo    schedule 25.05.2018
comment
(спустя годы) ваш Aw - это просто вектор, 4 x 1 ?? Также веса иногда ~ сигма, иногда ~ 1/сигма; попробовать веса все 1.   -  person denis    schedule 30.03.2020
comment
@denis Если я правильно помню, так оно и было. Мне пришлось использовать 1/weight, чтобы это начало работать. Кроме того, как указано в моем предыдущем комментарии, я использовал неправильный patsy_equation. Я рассматривал значения как категориальные, тогда как должен был рассматривать их как непрерывные числовые значения.   -  person Code Doggo    schedule 30.03.2020