第十章

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("D:/学习/Python数学建模算法与应用/第十章/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)

'''公式求解'''
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()

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()

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()