python加权最小二乘法_Python Pandas中加权最小二乘的意外标准错误

这是我的示例数据集,其中一些导入使用Pandas并直接使用scikits.statsmodels WLS：import pandas

import numpy as np

from statsmodels.regression.linear_model import WLS

# Make some random data.

np.random.seed(42)

df = pd.DataFrame(np.random.randn(10, 3), columns=['a', 'b', 'weights'])

# Add an intercept term for direct use in WLS

df['intercept'] = 1

# Add a number (I picked 10) to stabilize the weight proportions a little.

df['weights'] = df.weights 10

# Fit the regression models.

pd_wls = pandas.ols(y=df.a, x=df.b, weights=df.weights)

sm_wls = WLS(df.a, df[['intercept','b']], weights=df.weights).fit()

我使用％cpaste在IPython中执行此操作,然后打印两个回归的摘要：In [226]: %cpaste

Pasting code; enter '--' alone on the line to stop or use Ctrl-D.

:import pandas

:import numpy as np

:from statsmodels.regression.linear_model import WLS

:

:# Make some random data.

np:np.random.seed(42)

:df = pd.DataFrame(np.random.randn(10, 3), columns=['a', 'b', 'weights'])

:

:# Add an intercept term for direct use in WLS

:df['intercept'] = 1

:

:# Add a number (I picked 10) to stabilize the weight proportions a little.

:df['weights'] = df.weights 10

:

:# Fit the regression models.

:pd_wls = pandas.ols(y=df.a, x=df.b, weights=df.weights)

:sm_wls = WLS(df.a, df[['intercept','b']], weights=df.weights).fit()

:--

In [227]: pd_wls

Out[227]:

-------------------------Summary of Regression Analysis-------------------------

Formula: Y ~

Number of Observations: 10

Number of Degrees of Freedom: 2

R-squared: 0.2685

Adj R-squared: 0.1770

Rmse: 2.4125

F-stat (1, 8): 2.9361, p-value: 0.1250

Degrees of Freedom: model 1, resid 8

-----------------------Summary of Estimated Coefficients------------------------

Variable Coef Std Err t-stat p-value CI 2.5% CI 97.5%

--------------------------------------------------------------------------------

x 0.5768 1.0191 0.57 0.5869 -1.4206 2.5742

intercept 0.5227 0.9079 0.58 0.5806 -1.2567 2.3021

---------------------------------End of Summary---------------------------------

In [228]: sm_wls.summ

sm_wls.summary sm_wls.summary_old

In [228]: sm_wls.summary()

Out[228]:

"""

WLS Regression Results

==============================================================================

Dep. Variable: a R-squared: 0.268

Model: WLS Adj. R-squared: 0.177

Method: Least Squares F-statistic: 2.936

Date: Wed, 17 Jul 2013 Prob (F-statistic): 0.125

Time: 15:14:04 Log-Likelihood: -10.560

No. Observations: 10 AIC: 25.12

Df Residuals: 8 BIC: 25.72

Df Model: 1

==============================================================================

coef std err t P>|t| [95.0% Conf. Int.]

------------------------------------------------------------------------------

intercept 0.5227 0.295 1.770 0.115 -0.158 1.204

b 0.5768 0.333 1.730 0.122 -0.192 1.346

==============================================================================

Omnibus: 0.967 Durbin-Watson: 1.082

Prob(Omnibus): 0.617 Jarque-Bera (JB): 0.622

Skew: 0.003 Prob(JB): 0.733

Kurtosis: 1.778 Cond. No. 1.90

==============================================================================

"""

注意不匹配的标准误差：Pandas声称标准误差是[0.9079,1.0191],而statsmodels说[0.295,0.333].

首先,您可以看到标准错误是函数的报告：def _std_err_raw(self):

"""Returns the raw standard err values."""

return np.sqrt(np.diag(self._var_beta_raw))

所以看着self._var_beta_raw我发现：def _var_beta_raw(self):

"""

Returns the raw covariance of beta.

"""

x = self._x.values

y = self._y.values

xx = np.dot(x.T, x)

if self._nw_lags is None:

return math.inv(xx) * (self._rmse_raw ** 2)

else:

resid = y - np.dot(x, self._beta_raw)

m = (x.T * resid).T

xeps = math.newey_west(m, self._nw_lags, self._nobs, self._df_raw,

self._nw_overlap)

xx_inv = math.inv(xx)

return np.dot(xx_inv, np.dot(xeps, xx_inv))

在我的用例中,self._nw_lags将永远是None,所以这是令人费解的第一部分.由于xx只是回归矩阵的标准乘积：x.T.dot(x),我想知道权重是如何影响这一点的.术语self._rmse_raw直接来自于在OLS的构造函数中拟合的statsmodels回归,因此最明确地包含权重.

这提示了以下问题：

>为什么报告的标准错误是在RMSE部分中应用了权重,而不是在回归量变量中应用.

>如果您想要“非变换”变量(那么您是否也想要非变换后的RMSE？)这是标准做法吗？有没有办法让Pandas回馈标准错误的完全加权版本？

>为什么所有的误导？在构造函数中,计算完整的statsmodels拟合回归.为什么绝对不是每个汇总统计数据都直接来自那里？为什么它混合和匹配,以便一些来自statsmodels输出,一些来自熊猫家庭烹饪计算？

看起来我可以通过执行以下操作来协调Pandas输出：In [238]: xs = df[['intercept', 'b']]

In [239]: trans_xs = xs.values * np.sqrt(df.weights.values[:,None])

In [240]: trans_xs

Out[240]:

array([[ 3.26307961, -0.45116742],

[ 3.12503809, -0.73173821],

[ 3.08715494, 2.36918991],

[ 3.08776136, -1.43092325],

[ 2.87664425, -5.50382662],

[ 3.21158019, -3.25278836],

[ 3.38609639, -4.78219647],

[ 2.92835309, 0.19774643],

[ 2.97472796, 0.32996453],

[ 3.1158155 , -1.87147934]])

In [241]: np.sqrt(np.diag(np.linalg.inv(trans_xs.T.dot(trans_xs)) * (pd_wls._rmse_raw ** 2)))

Out[241]: array([ 0.29525952, 0.33344823])

我对这种关系感到非常困惑.这是统计学家常见的问题：将权重与RMSE部分联系起来,然后在计算系数的标准误差时选择是否对变量进行加权？如果是这种情况,为什么系数本身也不会在Pandas和statsmodel之间有所不同,因为那些类似地是从statsmodels首先转换的变量中得到的？

作为参考,这里是我的玩具示例中使用的完整数据集(如果np.random.seed不足以使其可重现)：In [242]: df

Out[242]:

a b weights intercept

0 0.496714 -0.138264 10.647689 1

1 1.523030 -0.234153 9.765863 1

2 1.579213 0.767435 9.530526 1

3 0.542560 -0.463418 9.534270 1

4 0.241962 -1.913280 8.275082 1

5 -0.562288 -1.012831 10.314247 1

6 -0.908024 -1.412304 11.465649 1

7 -0.225776 0.067528 8.575252 1

8 -0.544383 0.110923 8.849006 1

9 0.375698 -0.600639 9.708306 1

解决方法:

这里没有直接回答你的问题,但是,一般来说,你应该更喜欢将statsmodels代码用于建模的pandas.最近在statsmodels中发现了一些WLS问题,现在已经修复了. AFAIK,它们也固定在熊猫中,但大多数情况下,熊猫建模代码没有得到维护,中期目标是确保大熊猫中的所有可用内容都已弃用并且已移至statsmodels(下一版本为0.6.0,用于statsmodels)应该这样做).

为了更清楚一点,熊猫现在是statsmodels的依赖.您可以将DataFrame传递给statsmodel或在statsmodels中使用公式.这是未来的预期关系.来源：https://www.icode9.com/content-1-429651.html