【Python数据挖掘】回归模型与应用

 

线性回归 ( Linear Regression )

　　线性回归中，只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示，这种回归称为一元线性回归。

　　如果回归分析中包括两个或两个以上的自变量，且因变量和自变量之间是线性关系，则称为多元线性回归。

　　在监督学习中，学习样本为 D = { (x⁽ⁱ⁾, y⁽ⁱ⁾)；i =1, . . . , m } ，预测的结果y⁽ⁱ⁾为连续值变量，需要学习映射 f：X → Y ，并且假定输入X和输出Y之间有线性相关关系。

给出一组数据：

　　

其中x是实数域中的二维向量。比如，xⁱ₁是第i个房子的居住面积，xⁱ₂是这个房子的房间数。 

为了执行监督学习，我们需要决定怎样在计算机中表示我们的函数/假设。我们可以近似地使用线性函数来表示。 

（矩阵形式）

现在，有了训练数据，我们怎么挑选，或者说得知θ的值呢？一个可信的方法是使得h(x)与y更加接近，至少对于我们训练的例子来说是这样的。

于是，我们定义一个损失函数 / 成本函数( loss function / cost function )： 

我们把 x 到 y 的映射函数 f 记作 θ 的函数 h_θ(x)

损失函数有很多种类型，根据需求进行选择。

然后进行最小化损失函数，将函数优化成凸函数 (往往只会有一个全局最优解，不用过多担心算法收敛到局部最优解) 。

梯度下降 ( Gradient Descent algorithm )

　　最快的速度最小化损失函数，比作如何最快地下山，也就是每一步都应该往坡度最陡的方向往下走，而坡度最陡的方向就是损失函数相应的偏导数。

因此算法迭代的规则是：

假设现在有n个特征、或者变量x_j (j=1…n)

其中α是算法的参数learning rate，α越大每一步下降的幅度越大，速度也会越快，但过大有可能反复震荡，导致算法不准确。

欠拟合与过拟合（Underfitting and Overfitting）

　　欠拟合问题：特征值少，模型过于简单不足与支撑。

　　过拟合问题：有非常多特征，模型很复杂, 我们的假设函数曲线可以对原始数据拟合得非常好， 但丧失了一般性， 从而导致对新给的待预测样本，预测效果差。

正则项、正则化

　　通过正则项控制参数幅度。

正则项有多种方式选择，常采用的有：

　　L1正则：|θ_j|

　　L2正则：θ_j²

 

Logistic 回归（Logistic Regression）

采用线性回归解决分类问题时，往往会遇到模型健壮性低，遇到噪声时，受干扰严重。

我们可以对旧的线性回归算法来进行适当的修改来得到我们想要的函数。

引入sigmoid 函数：

对原函数h_θ(x)进行改写得到：

观察函数图像发现：当x大于0时，y的值大于0.5，根据这特性可以将线性回归得到的预测值压缩在0~1范围内。

1.线性判定边界：

假设线性函数为：，

当 h_θ(x) > 0 时，g(h_θ(x)) 的值为大于 0.5；

当 h_θ(x) < 0 时，g(h_θ(x)) 的值为小于 0.5；

2.非线性判定边界：

假设函数为：

当θ₀=0，θ₁=0，θ₂=0，θ₃=1，θ₄=1，得到函数g(x₁²+x₂²-1)，边界为一个圆，圆内点的值小于0

定义损失函数：

该函数为非凸函数，有局部最小值，应选择其他函数。

定义损失函数为：

该函数的图像如下：

 

我们可以发现该函数在：

y=1的正样本中，h_θ(x)趋向于0.99~9 ，此时我们希望得到的代价越小，而当得到的预测值是0.00~1时，我们希望它的代价越大；

y=0的负样本中，h_θ(x)趋向于0.00~1 ，此时我们希望得到的代价越小，而当得到的预测值是0.99~9时，我们希望它的代价越大；

损失函数可以改写成：

加入正则项：

二分类与多分类

one vs one

one vs rest

方法一：

1.先对三角形与叉叉进行分类，得到分类器C1，以及概率值P_c1(x) 和 1-P_c1(x)

2.然后对三角形与正方形进行分类，得到分类器C2，以及概率值P_c2(x) 和 1-P_c2(x)

3.最后对正方形与叉叉进行分类，得到分类器C3，以及概率值P_c3(x) 和 1-P_c3(x)

得到通过3个分类器,6个概率值，概率值最大的判断为相应的类型！

方法二：

1.先对三角形进行分类，判断是否为三角形，得到分类器C1，以及概率值P_c1(x)

2.然后对正方形进行分类，判断是否为正方形，得到分类器C2，以及概率值P_c2(x)

3.最后对叉叉叉进行分类，判断是否为叉叉叉，得到分类器C3，以及概率值P_c3(x)

得到3个分类器，3个概率值，概率值最大的判断为相应的类型！

 

应用一：( Linear Regression )

1、导入相应的包，设置画图格式：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
 
import seaborn as sns
sns.set_context('notebook')
sns.set_style('white')
plt.figure(figsize=(8,6))

2、数据准备：                                                                           ?下载地址

df = pd.read_csv(r'..\..\input\data1.txt',header=None,names=['x','y'])
 
xdata = np.c_[np.ones(df.shape[0]),df['x']]           # 转换成numpy 形式
ydata = np.c_[df['y']]

3、画出数据的散点图：

plt.scatter(xdata[:,1],ydata,s=50,marker='x',color='r')
plt.xlim(4,24)
plt.xlabel('Population of City in 10,000s')
plt.ylabel('Profit in $10,000s');
plt.show()

4、定义Loss function：

def computeCost(x,y,theta=[[0],[0]]):
    m = y.size
    h = x.dot(theta)
    j = 0
    j = 1.0/(2*m)*(np.sum(np.square(h-y)))
    return j

5、定义梯度下降算法：

def gradientDescent(x,y,theta=[[0],[0]],alpha=0.01, num_iters=1500):
    m = y.size
    J_history = np.zeros(num_iters)
    
    for i in range(num_iters):
        h = x.dot(theta)
        theta = theta - alpha * (1.0/m) * (x.T.dot(h-y))
        J_history[i]  = computeCost(x,y,theta)
    return (theta,J_history)

6、计算&画图：

theta,Cost_J = gradientDescent(xdata,ydata)
print('theta: ',theta.ravel())                       # theta:  [-3.63029144  1.16636235]
 
plt.figure(figsize=(8,6))
plt.plot(Cost_J)
plt.ylabel('Cost J')
plt.xlabel('Iterations');
plt.show()

得到θ的变化曲线：

通过计算我们能够得到 θ 的两个值，由此我们能画出一条与预测直线。

xx = np.arange(5,23)              # 生成直线上的点
yy = theta[0]+theta[1]*xx         # 计算出直线的方程表达式

7、散点图中画出直线 与 sklearn得到的结果进行对比：

# 画出我们自己写的线性回归梯度下降收敛的情况
plt.figure(figsize=(8,6))
plt.scatter(xdata[:,1],ydata,s=50,marker='x',color='r')
plt.plot(xx,yy, label='Linear regression (Gradient descent)',c='b')
 
# 和Scikit-learn中的线性回归对比一下 
from sklearn.linear_model import LinearRegression
regr = LinearRegression()
regr.fit(xdata[:,1].reshape(-1,1), ydata.ravel())
plt.plot(xx, regr.intercept_+regr.coef_*xx, label='Linear regression (Scikit-learn GLM)',c='g')
 
plt.xlim(4,24)
plt.xlabel('Population of City in 10,000s')
plt.ylabel('Profit in $10,000s');
plt.legend(loc=4);
plt.show()

可以看出两条直线是非常接近的！

8、进行数据的预测：

# 预测一下人口为35000和70000的城市的结果
print(theta.T.dot([1, 3.5])*10000)
print(theta.T.dot([1, 7])*10000)
 
[ 4519.7678677]
[ 45342.45012945]

补充：

numpy.c_              # 将切片对象转换为沿第二轴的连接。
>>> np.c_[np.array([[1,2,3]]), 0, 0, np.array([[4,5,6]])]   
array([[1, 2, 3, 0, 0, 4, 5, 6]])        # 1行8列!
 
#######################################
 
numpy.ravel           # 返回一个连续的扁平数组。
>>> x = np.array([[1, 2, 3], [4, 5, 6]])
>>> print(np.ravel(x))
[1 2 3 4 5 6]
 
#######################################
 
numpy.reshape(-1,1)   # 为数组提供新形状，而不更改其数据。
此为: n行1列

　　

应用二：( Logistic Regression )

1、导入相应的包，设置画图格式：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

2、数据准备

data = np.loadtxt(file, delimiter=delimeter)
 
X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]
y = np.c_[data[:,2]]

3、画出样本的散点图

def plotData(data, label_x, label_y, label_pos, label_neg, axes=None):
    # 获得正负样本的下标(即哪些是正样本，哪些是负样本)
    neg = data[:,2] == 0
    pos = data[:,2] == 1
    
    if axes == None:
        axes = plt.gca()
    axes.scatter(data[pos][:,0], data[pos][:,1], marker='+', c='k', s=60, linewidth=2, label=label_pos)
    axes.scatter(data[neg][:,0], data[neg][:,1], c='y', s=60, label=label_neg)
    axes.set_xlabel(label_x)
    axes.set_ylabel(label_y)
    axes.legend(frameon= True, fancybox = True);
 
plt.figure(figsize=(8,6))
plotData(data, 'Exam 1 score', 'Exam 2 score', 'Pass', 'Fail')

4、定义sigmoid函数

def sigmoid(z):
    return 1 / ( 1 + np.exp(-z) )

5、定义损失函数

def costFunction(theta, X, y):
    m = y.size
    h = sigmoid(X.dot(theta))
    
    J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y))
    
    if np.isnan(J[0]):
        return(np.inf)
    return J[0]

6、求解剃度函数：

def gradient(theta, X, y):
    m = y.size
    h = sigmoid(X.dot(theta.reshape(-1,1)))
    grad =(1.0/m)*X.T.dot(h-y)
    return(grad.flatten())

7、最小化损失函数：

initial_theta = np.zeros(X.shape[1])
res = minimize(costFunction, initial_theta, args=(X,y), jac=gradient, options={'maxiter':400})

8、预测：

sigmoid(np.array([1, 45, 85]).dot(res.x.T))

9、画出决策边界

plt.scatter(45, 85, s=60, c='r', marker='v', label='(45, 85)')
plotData(data, 'Exam 1 score', 'Exam 2 score', 'Admitted', 'Not admitted')
x1_min, x1_max = X[:,1].min(), X[:,1].max(),
x2_min, x2_max = X[:,2].min(), X[:,2].max(),
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))
h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(res.x))
h = h.reshape(xx1.shape)
plt.contour(xx1, xx2, h, [0.5], linewidths=1, colors='b');

10、添加正则项：

def costFunctionReg(theta, reg, *args):
    m = y.size
    h = sigmoid(XX.dot(theta))
    
    J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y)) + (reg/(2.0*m))*np.sum(np.square(theta[1:]))
    
    if np.isnan(J[0]):
        return(np.inf)
    return(J[0])
 
def gradientReg(theta, reg, *args):
    m = y.size
    h = sigmoid(XX.dot(theta.reshape(-1,1)))
      
    grad = (1.0/m)*XX.T.dot(h-y) + (reg/m)*np.r_[[[0]],theta[1:].reshape(-1,1)]
        
    return(grad.flatten())

11、最优化：

# C = [0.0, 1.0, 100.0] 
res2 = minimize(costFunctionReg, initial_theta, args=(C, XX, y), jac=gradientReg, options={'maxiter':3000})

　　

Kaggle比赛应用：

San Francisco Crime Classification

Predict the category of crimes that occurred in the city by the bay

import pandas as pd
import numpy as np

train = pd.read_csv('\\train.csv')
test = pd.read_csv('\\test.csv')

all_address = list(np.array(train['Address'].tolist() + test['Address'].tolist()))

from sklearn.feature_extraction.text import CountVectorizer
stop_words = ['dr', 'wy', 'bl', 'av', 'st', 'ct', 'ln', 'block', 'of']
vectorizer = CountVectorizer(max_features=300, stop_words=stop_words)
features = vectorizer.fit_transform(all_address).toarray()

X = features[:train.shape[0]]
y = train.Category

#分成80%的训练集和20%的验证集
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=44)

from sklearn.linear_model import LogisticRegression
log_model = LogisticRegression()
log_model.fit(X_train,y_train)
results = log_model.predict_proba(X_test)

from sklearn.metrics import log_loss
log_loss_score = log_loss(y_test, results)
print('log loss score: {0}'.format(round(log_loss_score, 3)))
# log loss score: 2.511

# 使用整个train
log_model = LogisticRegression().fit(X=features[:train.shape[0]], y=train.Category)
results = log_model.predict_proba(features[train.shape[0]:])

submission = pd.DataFrame(results)
submission.columns = sorted(train.Category.unique())
submission.set_index(test.Id)
submission.index.name="Id"
submission.to_csv('py_submission_logreg_addr_300.csv')