#!/usr/bin/python # -*- coding:utf-8 -*- import numpy as np from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV from sklearn.preprocessing import PolynomialFeatures import matplotlib.pyplot as plt from sklearn.pipeline import Pipeline import matplotlib as mpl import warnings def xss(y, y_hat): y = y.ravel() y_hat = y_hat.ravel() # Version 1 tss = ((y - np.average(y)) ** 2).sum() rss = ((y_hat - y) ** 2).sum() ess = ((y_hat - np.average(y)) ** 2).sum() r2 = 1 - rss / tss # print 'RSS:', rss, '\t ESS:', ess # print 'TSS:', tss, 'RSS + ESS = ', rss + ess tss_list.append(tss) rss_list.append(rss) ess_list.append(ess) ess_rss_list.append(rss + ess) # Version 2 # tss = np.var(y) # rss = np.average((y_hat - y) ** 2) # r2 = 1 - rss / tss corr_coef = np.corrcoef(y, y_hat)[0, 1] return r2, corr_coef if __name__ == "__main__": warnings.filterwarnings("ignore") # ConvergenceWarning np.random.seed(0) np.set_printoptions(linewidth=1000) N = 9 x = np.linspace(0, 6, N) + np.random.randn(N) x = np.sort(x) y = x ** 2 - 4 * x - 3 + np.random.randn(N) x.shape = -1, 1 y.shape = -1, 1 print '--' # 形如 1.00000000e+02 表示1.000*10^2 即1.00乘以10的2次幂 print np.logspace(-3, 2, 5) print np.logspace(-2, 9, 5) print np.logspace(0, 0, 5) print np.logspace(0, 1, 5) print '--' # 线性回归的目的是要得到输出向量Y和输入特征X之间的线性关系,求出线性回归系数θ,也就是Y=Xθ, # 其中Y的维度为mx1,X的维度为mxn,而θ的维度为nx1,m代表样本个数,n代表样本特征的维度 # 损失函数:损失函数是用来评价模型的预测值f(x)与真实值Y的不一致程度,它是一个非负实值函数 通常用L(Y,f(x))表示,损失函数越小,模型的性能就越好 # 正则化项:为了防止损失函数过拟合的问题,一般会在损失函数中加上正则化项,增加模型的泛化能力 models = [Pipeline([ ('poly', PolynomialFeatures()), # 损失函数:J(θ)=1/2(Xθ−Y)T(Xθ−Y) 优化方法:梯度下降和最小二乘法,scikit中采用最小二乘 # 使用场景:只要数据线性相关,LinearRegression是我们的首选,如果发现拟合或者预测的不够好,再考虑其他的线性回归库 ('linear', LinearRegression(fit_intercept=False))]), # Ridge回归(岭回归)损失函数的表达形式:J(θ)=1/2(Xθ−Y)T(Xθ−Y)+1/2α||θ||22(线性回归LineaRegression的损失函数+L2(2范式的正则化项)) # a为超参数 alphas=np.logspace(-3, 2, 50) 从给定的超参数a中选择一个最优的,logspace用于创建等比数列 本例中 开始点为10的-3次幂,结束点10的2次幂,元素个数为 # 50,并且从这50个数中选择一个最优的超参数 # linspace创建等差数列 # Ridge回归中超参数a和回归系数θ的关系,a越大,正则项惩罚的就越厉害,得到的回归系数θ就越小,最终趋近与0 # 如果a越小,即正则化项越小,那么回归系数θ就越来越接近于普通的线性回归系数 # 使用场景:只要数据线性相关,用LinearRegression拟合的不是很好,需要正则化,可以考虑使用RidgeCV回归, # 如何输入特征的维度很高,而且是稀疏线性关系的话, RidgeCV就不太合适,考虑使用Lasso回归类家族 Pipeline([ ('poly', PolynomialFeatures()), ('linear', RidgeCV(alphas=np.logspace(-3, 2, 50), fit_intercept=False))]), Pipeline([ ('poly', PolynomialFeatures()), # 损失函数:J(θ)=1/2m(Xθ−Y)T(Xθ−Y)+α||θ||1 线性回归LineaRegression的损失函数+L1(1范式的正则化项)) # Lasso回归可以使得一些特征的系数变小,甚至还使一些绝对值较小的系数直接变为0,从而增强模型的泛化能力 # 使用场景:对于高纬的特征数据,尤其是线性关系是稀疏的,就采用Lasso回归,或者是要在一堆特征里面找出主要的特征,那么 # Lasso回归更是首选了 ('linear', LassoCV(alphas=np.logspace(-3, 2, 50), fit_intercept=False))]), Pipeline([ ('poly', PolynomialFeatures()), # 损失函数:J(θ)=1/2m(Xθ−Y)T(Xθ−Y)+αρ||θ||1+α(1−ρ)/2||θ||22 其中α为正则化超参数,ρ为范数权重超参数 # alphas=np.logspace(-3, 2, 50), l1_ratio=[.1, .5, .7, .9, .95, .99, 1] ElasticNetCV会从中选出最优的 a和p # ElasticNetCV类对超参数a和p使用交叉验证,帮助我们选择合适的a和p # 使用场景:ElasticNetCV类在我们发现用Lasso回归太过(太多特征被稀疏为0),而Ridge回归也正则化的不够(回归系数衰减太慢)的时候 ('linear', ElasticNetCV(alphas=np.logspace(-3, 2, 50), l1_ratio=[.1, .5, .7, .9, .95, .99, 1], fit_intercept=False))]) ] mpl.rcParams['font.sans-serif'] = [u'simHei'] mpl.rcParams['axes.unicode_minus'] = False np.set_printoptions(suppress=True) plt.figure(figsize=(18, 12), facecolor='w') d_pool = np.arange(1, N, 1) # 阶 m = d_pool.size clrs = [] # 颜色 for c in np.linspace(16711680, 255, m): clrs.append('#%06x' % int(c)) line_width = np.linspace(5, 2, m) titles = u'线性回归', u'Ridge回归', u'LASSO', u'ElasticNet' tss_list = [] rss_list = [] ess_list = [] ess_rss_list = [] for t in range(4): model = models[t] plt.subplot(2, 2, t + 1) plt.plot(x, y, 'ro', ms=10, zorder=N) for i, d in enumerate(d_pool): model.set_params(poly__degree=d) model.fit(x, y.ravel()) lin = model.get_params('linear')['linear'] output = u'%s:%d阶,系数为:' % (titles[t], d) if hasattr(lin, 'alpha_'): idx = output.find(u'系数') output = output[:idx] + (u'alpha=%.6f,' % lin.alpha_) + output[idx:] if hasattr(lin, 'l1_ratio_'): # 根据交叉验证结果,从输入l1_ratio(list)中选择的最优l1_ratio_(float) idx = output.find(u'系数') output = output[:idx] + (u'l1_ratio=%.6f,' % lin.l1_ratio_) + output[idx:] print output, lin.coef_.ravel() x_hat = np.linspace(x.min(), x.max(), num=100) x_hat.shape = -1, 1 y_hat = model.predict(x_hat) s = model.score(x, y) r2, corr_coef = xss(y, model.predict(x)) # print 'R2和相关系数:', r2, corr_coef # print 'R2:', s, '\n' z = N - 1 if (d == 2) else 0 label = u'%d阶,$R^2$=%.3f' % (d, s) if hasattr(lin, 'l1_ratio_'): label += u',L1 ratio=%.2f' % lin.l1_ratio_ plt.plot(x_hat, y_hat, color=clrs[i], lw=line_width[i], alpha=0.75, label=label, zorder=z) plt.legend(loc='upper left') plt.grid(True) plt.title(titles[t], fontsize=18) plt.xlabel('X', fontsize=16) plt.ylabel('Y', fontsize=16) plt.tight_layout(1, rect=(0, 0, 1, 0.95)) plt.suptitle(u'多项式曲线拟合比较', fontsize=22) plt.show() y_max = max(max(tss_list), max(ess_rss_list)) * 1.05 plt.figure(figsize=(9, 7), facecolor='w') t = np.arange(len(tss_list)) plt.plot(t, tss_list, 'ro-', lw=2, label=u'TSS(Total Sum of Squares)') plt.plot(t, ess_list, 'mo-', lw=1, label=u'ESS(Explained Sum of Squares)') plt.plot(t, rss_list, 'bo-', lw=1, label=u'RSS(Residual Sum of Squares)') plt.plot(t, ess_rss_list, 'go-', lw=2, label=u'ESS+RSS') plt.ylim((0, y_max)) plt.legend(loc='center right') plt.xlabel(u'实验:线性回归/Ridge/LASSO/Elastic Net', fontsize=15) plt.ylabel(u'XSS值', fontsize=15) plt.title(u'总平方和TSS=?', fontsize=18) plt.grid(True) plt.show()
转载于:https://my.oschina.net/marjeylee/blog/1517895