scikit-learn 回归基础
2016-07-16 21:54
501 查看
numpy ndarray属性:shape dtype
有切片运算
当对赋值后的容器部分元素进行修改,影响原ndarray。
ndarray Array creation routines:
可以用来构造ndarray的一些常用接口(Ones and zeros and so on)
一些ndarray的接口:
ndarray.T(转置)
ndarray.flat: 类似于指针,文档中将其类比于迭代器。(这当然在两种语言中可以)
flat具有赋值功能,index与向量初始化顺序有关。
ndarray.item: 行为与直接取值运算相同,结果不同在于返回python 内置类型而非array数据类型。
由于直接学习numpy比较没有针对性,现考虑利用numpy及scipy进行scikit-learn一些结果的搭建。
1.1.Generalized Linear Models
一般最小二乘:
模型的求解是求导
简单地讲,一般最小二乘要求设计阵非(渐近)线性相关。
岭回归:这部分diy代码与调用模块结果不太一样:
利用求导方法将一般最小二乘推广到岭回归,可以直接得到。
由于岭回归与约束条件最小二乘等价,模型选择是一个问题。
岭回归是一个具有超参数的方法,超参数的选择可以通过交叉验证完成,即所谓GCV( generalized Cross-Validation)
这种交叉验证有简单的形式 https://en.wikipedia.org/wiki/Cross-validation_(statistics)
Leave-p-out cross-validation 可以进行组合数的交叉验证 利用MSE
岭回归是在矩阵“可逆”与“精确性”间选择。这里默认地选择了p=1的形式。
下面要进行关于lasso回归的讨论,其单纯的是将岭回归的二范数约束变为一范数,由于尖点不可导,
应该讨论最优化问题,下面先对scipy中的优化方式进行简单介绍。
下面是一个无约束条件的最优问题的求解代码:
这里nelder-mead是指使用单纯形法进行求解。(求解线性规划的算法,线性目标函数,线性不等式约束)
虽然不知道为什么要对非线性问题使用。
options中xtol为停止迭代的自变量误差,disp是否显示拟合结果,x0初值点。
下面再见一个非线性规划的代码:
这里使用SLSQP(细节并不明确)jac是导函数。
下面尝试对lasso回归使用此种优化方法求解。
lasso回归代码:
这里对diy代码的约束进行了较大扩展,接近真值,但是否过拟合?
后面还对岭回归进行了非线性求解模拟,仅仅是norm不同。
对于lasso回归具有稀疏解的问题,其直观的理解为对二维等高线及约束的图像观察,大多数取约束尖点为
最优解的原因,为登高圆圆心以概率1与约束多边形中心不垂直。(http://blog.csdn.net/xidianzhimeng/article/details/20856047)
由于很可能通过lasso回归得到稀疏解,故其合适对高维数据进行特征选取。(高维因子模型降维)
具体细节以后展开。(L1-based feature selection)
lasso回归参数alpha称为稀疏度,其的选取类似于岭回归交叉验证,细节暂时不做讨论。(还有AIC BIC等选择方法)
将lasso与Ridge进行结合得到的Elastic Net方法。下面是实现:
这里predict值都很接近自变量,但是在参数方面不太统一。
但设定为不估计常数项后结果非常接近。
对于lasso回归应用于多元情况(因变量为矩阵),仅仅是将维数进行了
扩充,但在代码实现上要注意,最优化不等式约束仅仅能对向量进行,对于矩阵,要将其拉直为向量
形式(否则出错),下面的最优化diy代码对其进行了实践。
并将参数的差距与ols形式(无偏估计)进行了比较。
import numpy as np
from numpy.random import normal, poisson, multivariate_normal
from numpy.linalg import norm, inv
from functools import partial
from scipy.optimize import minimize
X = np.zeros([5, 5])
W = np.zeros([5, 5])
for i in range(5):
X[i:] = normal(10, 3, 5)
W[i:] = poisson(10, 5)
Y = np.dot(X, W) + multivariate_normal(np.ones(5), np.eye(5), 5)
def func(X, Y, w):
W = np.zeros([5, 5])
for i in range(len(w)):
first, second = divmod(i, 5)
W[first][second] = w[i]
return norm(np.dot(X, W) - Y) / (2 * X.shape[0])
penalty = 10
alpha = 0.01
def strain(penalty, alpha, w):
W = np.zeros([5, 5])
for i in range(len(w)):
first, second = divmod(i, 5)
W[first][second] = w[i]
return penalty - alpha * sum(np.diag(np.dot(W.T, W)) ** 0.5)
strain_0 = partial(strain, penalty)
strain_1 = partial(strain_0, alpha)
cons = ({'type':'ineq',
'fun': strain_1},)
func_0 = partial(func, X)
func_1 = partial(func_0, Y)
res = minimize(func_1, np.zeros(5 * 5), constraints = cons, method = 'SLSQP',\
options = {'disp': True})
W_bar = np.zeros([5, 5])
for i in range(len(res.x)):
first, second = divmod(i, 5)
W_bar[first][second] = res.x[i]
print norm(W)
print norm(W - W_bar)
print norm(W - np.dot(np.dot(inv(np.dot(X.T, X)), X.T), Y))
从结果来看alpha较小时diy multi-lasso 收敛到ols,与理论相符。
下面实现了其的sklearn版本及DIY版本,注意在自定义代码中常数项的加法,对数据阵的加法与一元情况相同,
对参数阵比数据阵第二维度大(可以推出)
利用np.resize()进行横向拉直拷贝。
有切片运算
当对赋值后的容器部分元素进行修改,影响原ndarray。
ndarray Array creation routines:
可以用来构造ndarray的一些常用接口(Ones and zeros and so on)
一些ndarray的接口:
ndarray.T(转置)
ndarray.flat: 类似于指针,文档中将其类比于迭代器。(这当然在两种语言中可以)
flat具有赋值功能,index与向量初始化顺序有关。
ndarray.item: 行为与直接取值运算相同,结果不同在于返回python 内置类型而非array数据类型。
由于直接学习numpy比较没有针对性,现考虑利用numpy及scipy进行scikit-learn一些结果的搭建。
1.1.Generalized Linear Models
一般最小二乘:
模型的求解是求导
from sklearn import linear_model clf = linear_model.LinearRegression() print clf.fit([[0, 2], [1, 1], [2, 0]], [0, 1, 2]) print clf.coef_ print clf.intercept_ print "\n" * 5 import numpy as np from numpy.linalg import inv from numpy import linalg x = np.array([[0, 2], [1, 1], [2, 0]]) y = np.array([0, 1, 2]) def LinearRegression(X, y): # X_bar is the matrix for [1 X] X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X if linalg.matrix_rank(X_bar) == min(X_bar.shape): X = X_bar return np.dot(np.dot(inv(np.dot(X.T, X)), X.T), y) print LinearRegression(x, y)
简单地讲,一般最小二乘要求设计阵非(渐近)线性相关。
岭回归:这部分diy代码与调用模块结果不太一样:
from sklearn import linear_model clf = linear_model.Ridge(alpha = 0.5) print clf.fit([[0, 0], [0, 0], [1, 1]], [0, 0.1, 1]) print clf.coef_ print clf.intercept_ import numpy as np from numpy.linalg import inv x = np.array([[0,0], [0,0], [1,1]]) y = np.array([0, 0.1, 1]) def Ridge(alpha ,X, y): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X X = X_bar return np.dot(np.dot(inv(alpha * np.eye(X.shape[0]) + np.dot(X.T, X)), X.T), y) print Ridge(0.5, x, y)
利用求导方法将一般最小二乘推广到岭回归,可以直接得到。
由于岭回归与约束条件最小二乘等价,模型选择是一个问题。
岭回归是一个具有超参数的方法,超参数的选择可以通过交叉验证完成,即所谓GCV( generalized Cross-Validation)
这种交叉验证有简单的形式 https://en.wikipedia.org/wiki/Cross-validation_(statistics)
Leave-p-out cross-validation 可以进行组合数的交叉验证 利用MSE
岭回归是在矩阵“可逆”与“精确性”间选择。这里默认地选择了p=1的形式。
from sklearn import linear_model clf = linear_model.RidgeCV(alphas = [0.01,0.1, 1.0, 10.0]) print clf.fit([[0,0], [0,0], [1,1]], [0, 0.1, 1]) print clf.alpha_ import numpy as np from numpy.linalg import inv, norm x = np.array([[0, 0], [0, 0], [1, 1]]) y = [0, 0.1, 1] alphas = [0.01,0.1, 1.0, 1] def RidgeCV(alphas, X, y): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X def MSE(alpha): para = np.dot(np.dot(inv(np.dot(X_bar.T, X_bar) + np.eye(X_bar.shape[1]) * alpha), X_bar.T), y) return norm(np.dot(X_bar, para) - y) min_alpha_MSE = [] for element in alphas: MSE_val = MSE(element) if not min_alpha_MSE: min_alpha_MSE = [element, MSE_val] elif MSE_val < min_alpha_MSE[1]: min_alpha_MSE = [element, MSE_val] return min_alpha_MSE print RidgeCV(alphas, x, y)
下面要进行关于lasso回归的讨论,其单纯的是将岭回归的二范数约束变为一范数,由于尖点不可导,
应该讨论最优化问题,下面先对scipy中的优化方式进行简单介绍。
下面是一个无约束条件的最优问题的求解代码:
import numpy as np from scipy.optimize import minimize def rosen(x): return sum(100*(x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2) x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2]) res = minimize(rosen, x0, method='nelder-mead', \ options = {'xtol':1e-8, 'disp': True}) print res.x
这里nelder-mead是指使用单纯形法进行求解。(求解线性规划的算法,线性目标函数,线性不等式约束)
虽然不知道为什么要对非线性问题使用。
options中xtol为停止迭代的自变量误差,disp是否显示拟合结果,x0初值点。
下面再见一个非线性规划的代码:
def func(x, sign = 1.0): return sign * (2 * x[0] * x[1] + 2 * x[0] - x[0]**2 - 2*x[1]**2) def func_deriv(x, sign = 1.0): dfdx0 = sign * (-2 * x[0] + 2 * x[1] + 2) dfdx1 = sign * (2 * x[0] - 4 * x[1]) return np.array([dfdx0, dfdx1]) cons = ({'type': 'eq', 'fun': lambda x:np.array([x[0]**3 - x[1]]), 'jac': lambda x:np.array([3.0 * (x[0]**2), -1.0]) }, {'type': 'ineq', 'fun': lambda x:np.array([x[1] - 1]), 'jac': lambda x:np.array([0.0, 1.0]) }) res = minimize(func, [-1.0, 1.0], args = (-1.0,), jac = func_deriv, \ method = 'SLSQP', options = {'disp': True, }) print res.x res = minimize(func, [-1.0, 1.0], args = (-1.0,), jac = func_deriv, \ constraints = cons ,method = 'SLSQP', options = {'disp': True, }) print res.x
这里使用SLSQP(细节并不明确)jac是导函数。
下面尝试对lasso回归使用此种优化方法求解。
lasso回归代码:
这里对diy代码的约束进行了较大扩展,接近真值,但是否过拟合?
后面还对岭回归进行了非线性求解模拟,仅仅是norm不同。
from sklearn import linear_model clf = linear_model.Lasso(alpha = 0.1) clf.fit([[0, 0], [1, 1]], [0, 1]) print clf.predict([[1, 1]]) print clf.predict([[0, 0]]) import numpy as np from numpy.linalg import norm from scipy.optimize import minimize import functools def func(X ,y, beta): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X return norm(np.dot(X_bar, beta) - y) / (2 * X.shape[0]) #set alpha to 0.1 cons = ({ 'type': 'ineq', 'fun': lambda beta:np.array([10 - norm(beta, 1)]) },) X = np.array([[0, 0], [1, 1]]) y = np.array([0, 1]) func_0 = functools.partial(func, X) func_1 = functools.partial(func_0, y) #func_1 has the only para beta res = minimize(func_1, [0, 0, 0], constraints = cons, method = 'SLSQP',\ options = {'disp': True}) beta = res.x print beta print np.dot(np.array([1, 1, 1]), beta) print np.dot(np.array([1, 0, 0]), beta) print #use this method to ridge regression cons = ({ 'type': 'ineq', 'fun': lambda beta:np.array([10 - norm(beta)]) },) res = minimize(func_1, [0, 0, 0], constraints = cons, method = 'SLSQP',\ options = {'disp': True}) beta = res.x print beta print np.dot(np.array([1, 1, 1]), beta) print np.dot(np.array([1, 0, 0]), beta) print clf = linear_model.Ridge(alpha = 0.1) clf.fit(X, y) print clf.intercept_ print clf.coef_ print clf.intercept_ + np.dot(clf.coef_, np.array([ 4000 1, 1])) print clf.intercept_ + np.dot(clf.coef_, np.array([0, 0]))
对于lasso回归具有稀疏解的问题,其直观的理解为对二维等高线及约束的图像观察,大多数取约束尖点为
最优解的原因,为登高圆圆心以概率1与约束多边形中心不垂直。(http://blog.csdn.net/xidianzhimeng/article/details/20856047)
由于很可能通过lasso回归得到稀疏解,故其合适对高维数据进行特征选取。(高维因子模型降维)
具体细节以后展开。(L1-based feature selection)
lasso回归参数alpha称为稀疏度,其的选取类似于岭回归交叉验证,细节暂时不做讨论。(还有AIC BIC等选择方法)
将lasso与Ridge进行结合得到的Elastic Net方法。下面是实现:
from sklearn import linear_model clf = linear_model.ElasticNet(l1_ratio = 0.1) clf.fit([[108, 33], [44, 10001]], [861, 444]) print clf.coef_ print clf.intercept_ print clf.predict([[108, 33]]) print clf.predict([[44, 10001]]) print import numpy as np from numpy.linalg import norm from scipy.optimize import minimize import functools def func(X, y, beta): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X return norm(np.dot(X_bar, beta) - y) / (2 * X.shape[0]) alpha = 1.0 l1_ratio = 0.1 penalty = 10 cons = ({ 'type':'ineq', 'fun': lambda beta: penalty - alpha * l1_ratio * norm(beta, 1) - alpha * (1 - l1_ratio) * norm(beta) / 2 },) X = np.array([[108, 33], [44, 10001]]) y = np.array([861, 444]) func_0 = functools.partial(func, X) func_1 = functools.partial(func_0, y) res = minimize(func_1, [0, 8, 1], constraints = cons, method = 'SLSQP',\ options = {'disp': True}) print res.x print np.dot([1 ,108, 33], res.x) print np.dot([1 ,44, 10001], res.x)
这里predict值都很接近自变量,但是在参数方面不太统一。
但设定为不估计常数项后结果非常接近。
对于lasso回归应用于多元情况(因变量为矩阵),仅仅是将维数进行了
扩充,但在代码实现上要注意,最优化不等式约束仅仅能对向量进行,对于矩阵,要将其拉直为向量
形式(否则出错),下面的最优化diy代码对其进行了实践。
并将参数的差距与ols形式(无偏估计)进行了比较。
import numpy as np
from numpy.random import normal, poisson, multivariate_normal
from numpy.linalg import norm, inv
from functools import partial
from scipy.optimize import minimize
X = np.zeros([5, 5])
W = np.zeros([5, 5])
for i in range(5):
X[i:] = normal(10, 3, 5)
W[i:] = poisson(10, 5)
Y = np.dot(X, W) + multivariate_normal(np.ones(5), np.eye(5), 5)
def func(X, Y, w):
W = np.zeros([5, 5])
for i in range(len(w)):
first, second = divmod(i, 5)
W[first][second] = w[i]
return norm(np.dot(X, W) - Y) / (2 * X.shape[0])
penalty = 10
alpha = 0.01
def strain(penalty, alpha, w):
W = np.zeros([5, 5])
for i in range(len(w)):
first, second = divmod(i, 5)
W[first][second] = w[i]
return penalty - alpha * sum(np.diag(np.dot(W.T, W)) ** 0.5)
strain_0 = partial(strain, penalty)
strain_1 = partial(strain_0, alpha)
cons = ({'type':'ineq',
'fun': strain_1},)
func_0 = partial(func, X)
func_1 = partial(func_0, Y)
res = minimize(func_1, np.zeros(5 * 5), constraints = cons, method = 'SLSQP',\
options = {'disp': True})
W_bar = np.zeros([5, 5])
for i in range(len(res.x)):
first, second = divmod(i, 5)
W_bar[first][second] = res.x[i]
print norm(W)
print norm(W - W_bar)
print norm(W - np.dot(np.dot(inv(np.dot(X.T, X)), X.T), Y))
从结果来看alpha较小时diy multi-lasso 收敛到ols,与理论相符。
下面实现了其的sklearn版本及DIY版本,注意在自定义代码中常数项的加法,对数据阵的加法与一元情况相同,
对参数阵比数据阵第二维度大(可以推出)
rom sklearn import linear_model clf = linear_model.MultiTaskLasso(alpha = 0.1) clf.fit([[0,0], [1,1], [2,2]], [[0,0], [1,1], [2,2]]) print clf.coef_ print clf.intercept_ print "\n" import numpy as np from numpy.linalg import norm from scipy.optimize import minimize from functools import partial import copy def func(X, Y, w): W = np.resize(w, [X.shape[1] + 1, Y.shape[1]]) X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X return norm(np.dot(X_bar, W) - Y) / (2 * X.shape[0]) penalty = 1 def strain(penalty , first_dim, second_dim ,w): W = np.resize(w ,[first_dim, second_dim]) return penalty - sum(np.diag(np.dot(W.T, W)) ** 0.5) X = np.array([[0,0], [1,1], [2,2]]) Y = copy.deepcopy(X) strain_0 = partial(strain, penalty) strain_1 = partial(strain_0, X.shape[1] + 1) strain_2 = partial(strain_1, Y.shape[1]) cons = ({'type': 'ineq', 'fun': strain_2},) func_0 = partial(func, X) func_1 = partial(func_0, Y) res = minimize(func_1, np.zeros([(X.shape[1] + 1) * Y.shape[1]]),\ method = 'SLSQP', constraints = cons, options = {'disp': True}) print res.x.reshape(X.shape[1] + 1, Y.shape[1])
利用np.resize()进行横向拉直拷贝。
相关文章推荐
- Python动态类型的学习---引用的理解
- Python3写爬虫(四)多线程实现数据爬取
- 垃圾邮件过滤器 python简单实现
- 下载并遍历 names.txt 文件,输出长度最长的回文人名。
- install and upgrade scrapy
- Scrapy的架构介绍
- Centos6 编译安装Python
- 使用Python生成Excel格式的图片
- 让Python文件也可以当bat文件运行
- [Python]推算数独
- Python中zip()函数用法举例
- Python中map()函数浅析
- Python将excel导入到mysql中
- Python在CAM软件Genesis2000中的应用
- 使用Shiboken为C++和Qt库创建Python绑定
- FREEBASIC 编译可被python调用的dll函数示例
- Python 七步捉虫法