numpy ndarray属性:shape dtype
有切片运算
当对赋值后的容器部分元素进行修改,影响原ndarray。
有切片运算
当对赋值后的容器部分元素进行修改,影响原ndarray。
ndarray Array creation routines:
可以用来构造ndarray的一些常用接口(Ones and zeros and so on)
一些ndarray的接口:
ndarray.T(转置)
ndarray.flat: 类似于指针,文档中将其类比于迭代器。(这当然在两种语言中可以)
flat具有赋值功能,index与向量初始化顺序有关。
ndarray.item: 行为与直接取值运算相同,结果不同在于返回python 内置类型而非array数据类型。
可以用来构造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
岭回归是一个具有超参数的方法,超参数的选择可以通过交叉验证完成,即所谓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初值点。
虽然不知道为什么要对非线性问题使用。
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回归使用此种优化方法求解。
lasso回归代码:
这里对diy代码的约束进行了较大扩展,接近真值,但是否过拟合?
后面还对岭回归进行了非线性求解模拟,仅仅是norm不同。
这里对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([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)
最优解的原因,为登高圆圆心以概率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()进行横向拉直拷贝。