（四）回归模型

1 回归模型的引入

由于客观事物内部规律的复杂性及人们认识程度的限制，无法分析实际对象内在的因果关系，建立合乎机理规律的数学模型。所以在遇到有些无法用机理分析建立数学模型的时候，通常采取搜集大量数据的办法，基于对数据的统计分析去建立模型，其中用途最为广泛的一类随即模型就是统计回归模型。

回归模型确定的变量之间是相关关系，在大量的观察下，会表现出一定的规律性，可以借助函数关系式来表达，这种函数就称为回归函数或回归方程。
 

2 用回归模型解题的步骤

回归模型解题步骤主要包括两部分，一：确定回归模型属于那种基本类型，然后通过计算得到回归方程的表达式；二：是对回归模型进行显著性检验。

一：①根据试验数据画出散点图；

    ②确定经验公式的函数类型；

③通过最小二乘法得到正规方程组；

④求解方程组，得到回归方程的表达式。

二：①相关系数检验，检验线性相关程度的大小；

②F检验法（这两种检验方法可以任意选）；

③残差分析；

④对于多元回归分析还要进行因素的主次排序；

    如果检验结果表示此模型的显著性很差，那么应当另选回归模型了。
 

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
d=datasets.load_boston()
print(d.data)
print(d.DESCR)
print(d.feature_names)
print(d.data[:,5])
x=d.data[d.target<50]
y=d.target[d.target<50]

#1-2使用多元线性回归法对其进行训练和预测
from sklearn.linear_model import LinearRegression   #引入多元线性回归算法模块进行相应的训练
simple2=LinearRegression()
from sklearn.model_selection import train_test_split
x_train,x_test,y_train,y_test=train_test_split(x,y,random_state=666)
simple2.fit(x_train,y_train)
print(simple2.coef_)               #输出多元线性回归的各项系数
print(simple2.intercept_)          #输出多元线性回归的常数项的值
y_predict=simple2.predict(x_test)

#1-3利用sklearn里面的merics模块导出三大评价指标进行评价，直接进行调用计算
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score                #直接调用库函数进行输出R2
print(mean_squared_error(y_test,y_predict))
print(mean_absolute_error(y_test,y_predict))
print(r2_score(y_test,y_predict))
print(simple2.score(x_test,y_test))
print(simple2.coef_)               #输出多元回归算法的各个特征的系数矩阵
print(np.argsort(simple2.coef_))  #输出多元线性回归算法各个特征的系数排序，可以知道各个特征的影响度
print(d.feature_names[np.argsort(simple2.coef_)])  #输出各个特征按照影响系数从小到大的顺序

二、逐步回归模型

(一)逐步回归

逐步回归的基本思想是将变量逐个引入模型，每引入一个解释变量后都要进行F检验，并对已经选入的解释变量逐个进行t检验，当原来引入的解释变量由于后面解释变量的引入变得不再显著时，则将其删除。以确保每次引入新的变量之前回归方程中只包含显著性变量。这是一个反复的过程，直到既没有显著的解释变量选入回归方程，也没有不显著的解释变量从回归方程中剔除为止。以保证最后所得到的解释变量集是最优的。

实例

前提:当自变量x1,x2,x3…xn过多时,希望进行简化,找到对因变量贡献相对较大的自变量

需要计算出相关系数矩阵

然后根据自变量的贡献系数

找出贡献最大的自变量,与贡献最小的自变量,再进一步与F检验的Fin,Fout两个临界值比较,

判断是否能被选入,是否被剔除

再选出一个自变量之后,该自变量变为主元,对相关系数矩阵进行变化

在进行下一轮判断,直到所有的自变量都被确定

import pandas as pd
import numpy as np
import statsmodels.api as sm
def stepwise(X,y,alpha_in=0.1,alpha_out=0.15):
    '''X为所有可能的自变量构成的DataFrame,
      y为响应变量，
      alpha_in为引入变量的显著性水平，
      alpha_out为剔除变量的显著性水平'''
    included=[]
    while True:
        changed=False
        excluded=list(set(X.columns)-set(included))
        p_val=pd.Series(index=excluded)
        for new_column in excluded:
            model=sm.OLS(y,sm.add_constant(X[included+[new_column]])).fit()
            p_val[new_column]=model.pvalues[new_column]
        min_pval=p_val.min()
        #forward step
        if min_pval < alpha_in:
            changed=True
            add_feature=p_val.idxmin()
            included.append(add_feature)
            print("Add {:20} with p_value   {:.6}".format(add_feature,min_pval))
        #backward step
        model=sm.OLS(y,sm.add_constant(X[included])).fit()
        p_val=model.pvalues.iloc[1:]
        max_pval=p_val.max()
        if max_pval > alpha_out:
            changed=True
            drop_feature=p_val.idxmax()
            included.remove(drop_feature)
            print("Drop {:20} with p_value   {:.6}".format(drop_feature,max_pval))
        if not changed:
            break
    return included


data=pd.read_csv("/Users/mac/Desktop/test.csv",sep=",")
columns=['x1', 'x2', 'x3', 'x4']
X=data[columns]
y=data["y"]
result = stepwise(X,y)
print('resulting features:')
print(result)
####在找到了最优自变量时，下面就使用该变量建立回归方程
new_colmuns=result
model=sm.OLS(y,sm.add_constant(X[new_colmuns])).fit()
model.summary2()


#根据最大调整R方前向逐步回选模型
def  forward_selection(X,y):
    included=[]
    current_R,best_R=0.0,0.0
    while True:
        changed=False
        excluded=list(set(X.columns)-set(included))
        adj_R=pd.Series(index=excluded)
        for new_column in excluded:
            model=sm.OLS(y,sm.add_constant(X[included+[new_column]])).fit()
            adj_R[new_column]=model.rsquared_adj
        best_R=adj_R.max()
        add_feature=adj_R.idxmax()
        if best_R > current_R:
            changed=True
            included.append(add_feature)
            current_R=best_R
            print("Add {}".format(add_feature))
        if not changed:
            break
    return included      

三、岭回归模型

岭回归一般可以用来解决线性回归模型系数无解的两种情况，一方面是自变量间存在高度多重共线性，另一方面则是自变量个数大于等于观测个数

岭回归就是在矩阵xTx上增加一项使得矩阵非奇异，从而能够对其求逆。从上面两端代码我们可以看到，在之前对xTx求逆时都需要先判断xTx是否可以求逆，而岭回归就是解决这个问题的。岭回归的回归系数计算公式为：

 

import numpy as np
import matplotlib.pyplot as plt
 
def ridgeRegres(xMat,yMat,lam=0.2):
    xTx = xMat.T*xMat
    denom = xTx + np.eye(np.shape(xMat)[1])*lam
    if np.linalg.det(denom) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = denom.I * (xMat.T*yMat)
    return ws
 
def ridgeTest(xArr,yArr):
    xMat = np.mat(xArr); yMat=np.mat(yArr).T
    yMean = np.mean(yMat) # 数据标准化
    # print(yMean)
    yMat = yMat - yMean
    # print(xMat)
    #regularize X's
    xMeans = np.mean(xMat,0)
    xVar = np.var(xMat,0)
    xMat = (xMat - xMeans) / xVar #（特征-均值）/方差
    numTestPts = 30
    wMat = np.zeros((numTestPts,np.shape(xMat)[1]))
    for i in range(numTestPts): # 测试不同的lambda取值，获得系数
        ws = ridgeRegres(xMat,yMat,np.exp(i-10))
        wMat[i,:]=ws.T
    return wMat
 
 
# import data
ex0 = np.loadtxt('abalone.txt',delimiter='\t')
xArr = ex0[:,0:-1]
yArr = ex0[:,-1]
# print(xArr,yArr)
ridgeWeights = ridgeTest(xArr,yArr)
# print(ridgeWeights)
plt.plot(ridgeWeights)
plt.show()

• Ridge 回归：
  • 执行 L2 正则化，即增加相当于系数幅度的平方的惩罚
  • 最小化目标函数：LS Obj + a*（系数平方和）
• Lasso 回归：
  • 执行 L1 正则化，即增加相当于系数幅度的绝对值的惩罚
  • 最小化目标函数：LS Obj + a*（系数绝对值之和）

        惩罚λ的引入能够减少不重要的参数，从而更好滴理解数据。参数λ的选择，需要将原训练数据分成训练数据和测试数据，在训练数据上训练回归系数，然后在测试数据上测试性能，通过选取不同的λ来重复上述过程，最终选取使得预测误差最小的λ。
岭回归缺陷：
1.主要靠目测选择岭参数
2.计算岭参数时，各种方法结果差异较大
所以一般认为，岭迹图只能看多重共线性，却很难做变量筛选

 

四、Lasso回归模型

lasso回归通过构造一个一阶惩罚函数获得一个精炼的模型；通过最终确定一些指标（变量）癿系数为零（岭回归估计系数等于0癿机会微乎其微，造成筛选变量困难），解释力很强。
擅长处理具有多重共线性的数据，筛选变量，与岭回归一样是有偏估计。lasso回归则是与岭回归略有不同，使用了L1正则

import numpy as np
# 计算残差
def rss(X,y,w):
	return (y-X*w).T*(y-X*w)
	
# lasso 回归	
def lassoRegression(X,y,a):
	threshold = 0.1
	# init para
	m,n = X.shape

	w = np.matrix(np.zeros((n,1)))
	r = rss(X,y,w)
	z = [0 for i in range(n)]

	tempw = [0 for i in range(n)]
	# 记录迭代次数
	iteration = 0
	# 迭代以计算参数
	while True:
		iteration += 1
		for k in range(n):
			temp_X = X[:, k].reshape(-1,1)
			z[k] = (temp_X.T*temp_X)[0, 0]
			p = [0 for i in range(n)]
			p_part = [0 for i in range(n)]
			for i in range(m):
				for j in range(n):
					if k != j:
						# print w[j,0]
						p_part[j] = X[i, j]*w[j, 0]
				p[k] += X[i, k] * (y[i, 0] - sum(p_part))
			if p[k] < -a/2:
				tempw[k] = (p[k] + a/2)/z[k]
			elif p[k] > a/2:
				tempw[k] = (p[k] - a/2)/z[k]
			else:
				tempw[k] = 0
			w[k, 0] = tempw[k]
		tempr = rss(X, y, w)
		delta = abs(tempr - r)[0, 0]
		r = tempr

		if delta < threshold:
			print 'iteration', iteration
			print 'delta', delta
			break
	return w
	
def load_data(filename):
	X, Y = [], []
	with open(filename, 'r') as f:
		for line in f:
			splited_line = [float(i) for i in line.split()]
			x, y = splited_line[: -1], splited_line[-1]
			X.append(x)
			Y.append(y)
	X, Y = np.matrix(X), np.matrix(Y).T
	return X, Y
def std(X):
	std_deviation = np.std(X, 0)
	mean = np.mean(X, 0)
	return (X - mean)/std_deviation
if '__main__'==__name__:
	# 使用鲍鱼数据集
	X, y = load_data('abalone.txt')
	X, y = std(X), std(y)
	enum = 30
    w_all = lassoTrace(X, y, enum)

    lambdas = [i-10 for i in range(enum)]

    data = pd.DataFrame({
        'x1': w_all[:, 0],
        'x2': w_all[:, 1],
        'x3': w_all[:, 2],
        'x4': w_all[:, 3],
        'x5': w_all[:, 4],
        'x6': w_all[:, 5],
        'x7': w_all[:, 6],
        'x8': w_all[:, 7],
    })
    plt.figure(12)

    plt.subplot(221)
    data.boxplot()

    plt.subplot(222)
    plt.plot(lambdas, w_all)
    plt.show()

 Lasso和岭回归