用`statsmodels`绘制掩码值的残差答案

【问题标题】：Plotting residuals of masked values with `statsmodels`用`statsmodels`绘制掩码值的残差
【发布时间】：2019-01-17 13:07:39
【问题描述】：

我正在使用statsmodels.api 来计算两个变量之间 OLS 拟合的统计参数：

def computeStats(x, y, yName):
    '''
    Takes as an argument an array, and a string for the array name.
    Uses Ordinary Least Squares to compute the statistical parameters for the
    array against log(z), and determines the equation for the line of best fit.
    Returns the results summary, residuals, statistical parameters in a list, and the 
    best fit equation.
    '''

    #   Mask NaN values in both axes
    mask = ~np.isnan(y) & ~np.isnan(x)
    #   Compute model parameters
    model = sm.OLS(y, sm.add_constant(x), missing= 'drop')
    results = model.fit()
    residuals = results.resid

    #   Compute fit parameters
    params = stats.linregress(x[mask], y[mask])
    fit = params[0]*x + params[1]
    fitEquation = '$(%s)=(%.4g \pm %.4g) \\times redshift+%.4g$'%(yName,
                    params[0],  #   slope
                    params[4],  #   stderr in slope
                    params[1])  #   y-intercept
    return results, residuals, params, fit, fitEquation

函数的第二部分（使用stats.linregress）可以很好地处理掩码值，但statsmodels 不能。当我尝试使用plt.scatter(x, resids) 绘制残差与 x 值时，尺寸不匹配：

ValueError: x and y must be the same size

因为有 29007 个 x 值和 11763 个残差（这是通过掩蔽过程的 y 值的数量）。我尝试将model 变量更改为

model = sm.OLS(y[mask], sm.add_constant(x[mask]), missing= 'drop')

但这没有任何效果。

如何将残差与它们匹配的 x 值进行散点图？

【问题讨论】：

标签： python-3.x statistics regression statsmodels

【解决方案1】：

嗨 @jim421616 由于 statsmodels 丢弃了一些缺失值，您应该使用模型的 exog 变量来绘制散点图，如图所示。

plt.scatter(model.model.exog[:,1], model.resid)

完整的虚拟示例供参考

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

#generate data
x = np.random.rand(1000)
y =np.sin( x*25)+0.1*np.random.rand(1000)

# Make some as NAN
y[np.random.choice(np.arange(1000), size=100)]= np.nan
x[np.random.choice(np.arange(1000), size=80)]= np.nan


# fit model
model = sm.OLS(y, sm.add_constant(x) ,missing='drop').fit()
print model.summary()

# plot 
plt.scatter(model.model.exog[:,1], model.resid)
plt.show()

【讨论】：