例10.1
点击查看代码
import numpy as np
import statsmodels.api as sm
import pylab as plt
def check(d):
x0 = d[0]; y0 = d[1]; d = {'x':x0, 'y':y0}
re = sm.formula.ols('y~x', d).fit()
print(re.summary())
print(re.outlier_test())
print('残差的方差', re.mse_resid)
pre=re.get_prediction(d)
df = pre.summary_frame(alpha=0.05)
dfv = df.values; low, upp = dfv[:,4:].T
r = (upp-low)/2
num = np.arange(1, len(x0)+1)
plt.errorbar(num, re.resid, r, fmt='o')
plt.show()
a = np.loadtxt('data10_1.txt')
plt.rc('font, size=15');plt.plot(a[0],a[1],'o')
plt.figure();check(a)
a2 = a; a2 = np.delete(a2, 8, axis=1)
check(a2);a3 = a2
a3 = np.delete(a3, 4, axis=1);check(a3)
例10.2
点击查看代码
'''公式求解'''
import numpy as np
import statsmodels.api as sm
import pylab as plt
a = np.loadtxt('data10_2.txt')
plt.rc('text', usetex=False);plt.rc('font',size=16)
plt.plot(a[0], a[2], '*', label='$x_1$')
plt.plot(a[1],a[2],'o', label='$x_2$')
plt.legend(loc='upper left')
d = {'x1': a[0], 'x2': a[1], 'y': a[2]}
re = sm.formula.ols('y~x1+x2', d).fit()
print(re.summary())
yh = re.predict({'x1': [9,10], 'x2': [10,9]})
print('残差的方差：', re.mse_resid)
print('预测值：',yh);plt.show()
'''数组求解'''
import numpy as np
import statsmodels.api as sm
import pylab as plt
a = np.loadtxt('data10_2.txt')
plt.rc('text', usetex=False);plt.rc('font',size=16)
plt.plot(a[0], a[2], '*', label='$x_1$')
plt.plot(a[1],a[2],'o', label='$x_2$')
plt.legend(loc='upper left')
X = sm.add_constant(a[:2].T)
re = sm.OLS(a[2], X).fit()
print(re.summary())
yh = re.predict(np.array([[1,9,10],[1,10,9]]))
print('残差的方差：', re.mse_resid)
print('预测值：',yh);plt.show()
例10.3
点击查看代码
import numpy as np
import statsmodels.formula.api as smf
import pylab as plt
x = np.arange(17, 30 ,2)
a = np.loadtxt('data10_3.txt')
plt.rc('text', usetex=False);plt.rc('font',size=16)
plt.plot(x, a[0], '*', label='$y_1$')
plt.plot(x, a[1], 'o', label='$y_2$')
x = np.hstack([x,x]); d = {'y':a.flatten(),'x':x}
re = smf.ols('y~x+I(x ** 2)',d).fit()
print(re.summary())
print('残差的方差：', re.mse_resid)
plt.legend();plt.show()
例10.4
点击查看代码
import numpy as np
import statsmodels.formula.api as smf
import pylab as plt
a = np.loadtxt('data10_4.txt'); x1 = a[0]; x2 = a[1]; y = a[2]
plt.rc('text', usetex=False);plt.rc('font',size=16)
plt.plot(x1,y,'*',label='$x_1$');plt.plot(x2,y,'o',label='$x_2$')
d = {'y':y,'x1':x1,'x2':x2}
re1 = smf.ols('y~x1+x2', d).fit()
print('线性回归的残差方差：', re1.mse_resid)
re2 = smf.ols('y~x1+x2+I(x1 ** 2)+I(x2 ** 2)', d).fit()
print('纯二次的残差方差：', re2.mse_resid)
re3 = smf.ols('y~x1 * x2', d).fit()
print('交叉二次的残差方差：', re3.mse_resid)
re4 = smf.ols('y~x1 * x2+I(x1 ** 2)+I(x2 ** 2)', d).fit()
print('完全二次的残差方差：', re4.mse_resid)
print('预测值：', re2.predict({'x1':170,'x2':160}))
print(re2.summary());plt.legend;plt.show()
