1# -*- coding: utf-8 -*- 2"""Examples for Regression Plots 3 4Author: Josef Perktold 5 6""" 7 8import numpy as np 9import statsmodels.api as sm 10import matplotlib.pyplot as plt 11 12from statsmodels.sandbox.regression.predstd import wls_prediction_std 13import statsmodels.graphics.regressionplots as smrp 14from statsmodels.graphics.tests.test_regressionplots import TestPlot 15 16#example from tut.ols with changes 17#fix a seed for these examples 18np.random.seed(9876789) 19 20# OLS non-linear curve but linear in parameters 21# --------------------------------------------- 22 23nsample = 100 24sig = 0.5 25x1 = np.linspace(0, 20, nsample) 26x2 = 5 + 3* np.random.randn(nsample) 27X = np.c_[x1, x2, np.sin(0.5*x1), (x2-5)**2, np.ones(nsample)] 28beta = [0.5, 0.5, 1, -0.04, 5.] 29y_true = np.dot(X, beta) 30y = y_true + sig * np.random.normal(size=nsample) 31 32#estimate only linear function, misspecified because of non-linear terms 33exog0 = sm.add_constant(np.c_[x1, x2], prepend=False) 34 35# plt.figure() 36# plt.plot(x1, y, 'o', x1, y_true, 'b-') 37 38res = sm.OLS(y, exog0).fit() 39#print res.params 40#print res.bse 41 42 43plot_old = 0 #True 44if plot_old: 45 46 #current bug predict requires call to model.results 47 #print res.model.predict 48 prstd, iv_l, iv_u = wls_prediction_std(res) 49 plt.plot(x1, res.fittedvalues, 'r-o') 50 plt.plot(x1, iv_u, 'r--') 51 plt.plot(x1, iv_l, 'r--') 52 plt.title('blue: true, red: OLS') 53 54 plt.figure() 55 plt.plot(res.resid, 'o') 56 plt.title('Residuals') 57 58 fig2 = plt.figure() 59 ax = fig2.add_subplot(2,1,1) 60 #namestr = ' for %s' % self.name if self.name else '' 61 plt.plot(x1, res.resid, 'o') 62 ax.set_title('residuals versus exog')# + namestr) 63 ax = fig2.add_subplot(2,1,2) 64 plt.plot(x2, res.resid, 'o') 65 66 fig3 = plt.figure() 67 ax = fig3.add_subplot(2,1,1) 68 #namestr = ' for %s' % self.name if self.name else '' 69 plt.plot(x1, res.fittedvalues, 'o') 70 ax.set_title('Fitted values versus exog')# + namestr) 71 ax = fig3.add_subplot(2,1,2) 72 plt.plot(x2, res.fittedvalues, 'o') 73 74 fig4 = plt.figure() 75 ax = fig4.add_subplot(2,1,1) 76 #namestr = ' for %s' % self.name if self.name else '' 77 plt.plot(x1, res.fittedvalues + res.resid, 'o') 78 ax.set_title('Fitted values plus residuals versus exog')# + namestr) 79 ax = fig4.add_subplot(2,1,2) 80 plt.plot(x2, res.fittedvalues + res.resid, 'o') 81 82 # see http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/partregr.htm 83 fig5 = plt.figure() 84 ax = fig5.add_subplot(2,1,1) 85 #namestr = ' for %s' % self.name if self.name else '' 86 res1a = sm.OLS(y, exog0[:,[0,2]]).fit() 87 res1b = sm.OLS(x1, exog0[:,[0,2]]).fit() 88 plt.plot(res1b.resid, res1a.resid, 'o') 89 res1c = sm.OLS(res1a.resid, res1b.resid).fit() 90 plt.plot(res1b.resid, res1c.fittedvalues, '-') 91 ax.set_title('Partial Regression plot')# + namestr) 92 ax = fig5.add_subplot(2,1,2) 93 #plt.plot(x2, res.fittedvalues + res.resid, 'o') 94 res2a = sm.OLS(y, exog0[:,[0,1]]).fit() 95 res2b = sm.OLS(x2, exog0[:,[0,1]]).fit() 96 plt.plot(res2b.resid, res2a.resid, 'o') 97 res2c = sm.OLS(res2a.resid, res2b.resid).fit() 98 plt.plot(res2b.resid, res2c.fittedvalues, '-') 99 100 # see http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/ccpr.htm 101 fig6 = plt.figure() 102 ax = fig6.add_subplot(2,1,1) 103 #namestr = ' for %s' % self.name if self.name else '' 104 x1beta = x1*res.params[1] 105 x2beta = x2*res.params[2] 106 plt.plot(x1, x1beta + res.resid, 'o') 107 plt.plot(x1, x1beta, '-') 108 ax.set_title('X_i beta_i plus residuals versus exog (CCPR)')# + namestr) 109 ax = fig6.add_subplot(2,1,2) 110 plt.plot(x2, x2beta + res.resid, 'o') 111 plt.plot(x2, x2beta, '-') 112 113 114 #print res.summary() 115 116doplots = 1 117if doplots: 118 fig1 = smrp.plot_fit(res, 0, y_true=None) 119 smrp.plot_fit(res, 1, y_true=None) 120 smrp.plot_partregress_grid(res, exog_idx=[0,1]) 121 smrp.plot_regress_exog(res, exog_idx=0) 122 smrp.plot_ccpr(res, exog_idx=0) 123 smrp.plot_ccpr_grid(res, exog_idx=[0,1]) 124 125tp = TestPlot() 126tp.test_plot_fit() 127 128fig1 = smrp.plot_partregress_grid(res, exog_idx=[0,1]) 129#add lowess 130ax = fig1.axes[0] 131y0 = ax.get_lines()[0]._y 132x0 = ax.get_lines()[0]._x 133lres = sm.nonparametric.lowess(y0, x0, frac=0.2) 134ax.plot(lres[:,0], lres[:,1], 'r', lw=1.5) 135ax = fig1.axes[1] 136y0 = ax.get_lines()[0]._y 137x0 = ax.get_lines()[0]._x 138lres = sm.nonparametric.lowess(y0, x0, frac=0.2) 139ax.plot(lres[:,0], lres[:,1], 'r', lw=1.5) 140 141#plt.show() 142