线性回归和 logistic回归算法原理推导

线性回归是机器学习的入门基础。在机器学习中,主要分类两类:回归和分类。而线性回归属于回归,虽然logistics回归名字带有回归,其实这个模型完成的分类任务。简单的理解回归和分类,其实就是回归的输出是一个具体的数值,而回归的输出是一个特定的类别。

线性回归:

  • y_{i} = \theta ^{T}\mathbf{x_{i}}+ \varepsilon    (1)

这里的theta 是个矩阵,代表的是样本特征的系数,X代表的是样本的特征。\varepsilon代表的是误差项。这个误差项独立同分布的。服从标准正态分布。

  • 独立:这里的样本之间是没有联系的,也就是说一个样本的数据是不会影响另一个样本的数据的
  • 同分布:因为这是按照同一套数据进行划分的,也就是这套数据的使用场景是同一个场景。
  • 高斯分布:照现实生活的情况来估计的,因为在一个现实生活中,误差总是有正有负的,而且正负的比例是差不多的。所以服从误差为0.方差为\Theta ^{2}的标准正态分布。

误差为\varepsilon的概率密度为 p(\epsilon ) = \frac{1}{\sqrt{2\pi }\sigma} exp(-\frac{\varepsilon ^{2}}{2\sigma ^{2}})               (2)

现在将(1)式带入(2)式。得到p(y^{i}|x^{i}:\theta ) = \frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y^{i}-\theta ^{T}x)^{2}}{2\sigma ^{2}})    (3)

然后根据(3)式可以得到要求的似然函数:L(\theta ) = \prod_{1}^{m}p(y^{i}|x^{i}:\theta ) = \prod_{1}^{m}\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y^{i}-\theta ^{T}x)^{2}}{2\sigma ^{2}})   (4)

这个式子的来源就是因为误差项是独立同分布的,所以概率密度相乘的话越大也就是代表参数选择使得我们的预测值和真实值越接近。这个我是这么理解的,因为对于标准正态分布,在均值处取得最大值,离均值越远,值越小。而现实生活中那些特别离谱的值是很少存在的,所以只要调整参数,使误差项在均值附近就可以使L(\Theta )取得最大值。

因为乘法的求导是很麻烦的,所以利用对数似然函数。

Log L(\theta ) = Log\prod_{1}^{m}p(y^{i}|x^{i}:\theta ) =Log \prod_{1}^{m}\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y^{i}-\theta ^{T}x)^{2}}{2\sigma ^{2}})     (5)

展开后得到:mLog\frac{1}{\sqrt{2\pi }\sigma } - \frac{1}{\sigma ^{2}}\cdot \frac{1}{2}\sum_{i=1}^{m}(y^{i}-\theta ^{T}x)^{2}             (6)

所以现在我们只需要把    \sum_{i=1}^{m}(y^{i}-\theta ^{T}x)^{2}   (最小二乘法)      (7)   这个值最小就可以了。

也就是令其偏导数为0就可以了  得到: \theta = (X^{T}X)^{-1}X^{T}y

这个求解算法虽然被证明也是对的,但是这完全不是机器学习的套路,机器学习的套路是用一堆现有的数据,然后告诉我怎么样,想让我达到什么样的水平,然后我一点点的去学习。这就是用到了下面的梯度下降算法。(最小二乘法只是个巧合,正好能用公式算出来,但是并不是所有的数据都能适用这种方法)

梯度,这里其实就是数学上的对各个方向上的偏导数,当函数沿着梯度的方向改变变量的时候,函数的变化最快。但是按我的理解,是先有了变化最快的方向,然后人们给其起了个名字,把它叫做梯度。所以,当我们按着梯度的反方向更新参数的时候,就可以最快的达到最小值。

目标函数:J(\theta ) = \frac{1}{2m}\sum_{i=1}^{m}(y^{i}-h_{\theta }(x^{i}))^{2}    (8)

  • 批量梯度下降:\frac{J(\partial \theta )}{\partial \Theta_{j}} = -\frac{1}{m}\sum_{i=1}^{m}(y^{i}-h_{\theta }(x^{i}))x_{j}^{i}   (9)     \Theta _{j} = \theta _{j} +\frac{1}{m}\sum_{i=1}^{m}(y^{i}-h_{\theta }(x^{i}))x_{j}^{i}     (10)

这里是用了所有的样本一起进行了计算,然后为了达到最后的效果,只需要进行迭代次数的修改。此方法计算的时间很慢,因为一次要计算所有的样本值

  • 随机梯度下降:\Theta _{j} = \theta _{j} +(y^{i}-h_{\theta }(x^{i}))x_{j}^{i}

这里是用了一个样本进行参数的更新,但是这样有个不好的效果,那就是会产生震荡,因为数据里总会有些噪声。

  • 小批量梯度下降:\Theta _{j} = \theta _{j} +\frac{1}{m}\sum_{i=1}^{m}(y^{i}-h_{\theta }(x^{i}))x_{j}^{i}

这种办法是每次采用一批量的数据进行数据的更新。也是使用最多的方法。

  • logistic回归

logistic是经典的二分类算法,在机器学习中通常被用来先用逻辑回归进行分类,然后实在不行,再用复杂的分类算法进行分类。这样比较符合人们正常的思维,舍难而取易者也。这里的逻辑回归虽然是二分类算法,但是也可以进行多分类,softmax函数就可以实现这样的效果。softmax在神经网络里很常用。逻辑回归的边界是可以是非线性的,因为只要数据是高阶的,出来的边界就可以是非线性的。

这里要引入sigmoid函数:g(z) = \frac{1}{1+e^{- z}}        (11)

代码实现:

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

这个函数的定义域是(-\infty ,+\infty),值域是 [0,1]。

在这里我们做个映射  h_{\theta }(x) = g(\theta ^{T}x) = \frac{1}{1+e^{-\theta ^{T}x}}     (12)  现在线性函数的输出值变成了在区间0,1上的值了。

p(y=1|x:\theta ) = h_{\theta }(x)        (12)

p(y=0|x:\theta ) = 1-h_{\theta }(x)    (13)

将两者合在一起比较好些似然函数  得到:p(y|x:\theta ) = (h_{\theta }(x))^{y}(1-h_{\theta }(x))^{1-y}     (14)

似然函数:L(\theta ) = \prod_{i=1}^{m} (h_{\theta }(x))^{y}(1-h_{\theta }(x))^{1-y}     (15)

对数似然:log L(\theta ) =\sum_{i=1}^{m} (y^{i}log(h_{\theta }(x))+(1-y^{i})log(1-h_{\theta }(x)))     (16)

按照机器学习的套路,是求目标函数的最小值。所以这里改为求  -log(L(\theta )).

这个函数就可以理解为损失函数,我们的目的就是更新这个函数,使其达到最大值。但是还是按照机器学习的套路来,加个负号求其最小值。代码实现为:

#model :返回预测值
def gfunc(x,t):
    return sigmoid(np.dot(x,t))
print((gfunc(data_x,theta)).shape)
#cost : 根据参数计算损失
def cost(x,y):
    s = np.multiply(y,np.log(x)) + np.multiply((1-y),np.log(1-x))
    sum =-np.sum(s)/len(x)
    return sum

下面是对其求偏导的过程:(公式实在是太难打了,就直接给个结果吧。)

\frac{1}{m}\sum_{i=1}^{m}(h_{\theta }(x) - y^{i})x_{i}^{j}

所以最后的参数更新为:\theta _{j} = \theta _{j} - \alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\theta }(x) - y^{i})x_{i}^{j}   (17)

代码实现为:

def gradient(data,x,y,t):
    c = (x - y).ravel()
    for i in range(len(t)):
        s = np.sum(np.multiply(c,data[:,i]))/len(x)
        t[i] = t[i] - 0.001*s
for j in range(5000):
    gradient(data_x,gfunc(data_x,theta),data_y,theta)

对于怎么做分类,只需要指定阈值。因为最后输出的值是个概率,默认的值是0.5。但是这个值是可以改变的,这时候要考虑到实际的需求,比如要考虑召回率高的时候,就可以适当的降低这个阈值,但是此时带来的也是将反例错分类的代价。所以这个是要根据实际需要进行设定的。下面附录整个Python实现逻辑回归和梯度下降策略的选择。(训练数据没有上传)(环境:Python3.6+Spyder+Numpy+pandas+matplotlib)

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


#首先进行数据的读取
data  = pd.read_csv('LogiReg_data.txt',header = None,names = ['num1','num2','class'])
#print(data)         #打印数据


#首先增加偏置项
data.insert(0,'num0',1)
#拿到数据后,首先观看一下数据的分布特征

class_1 = data[data['class'] == 1]
class_0 = data[data['class'] == 0]

print('class_1的样本数量:     %s' %(len(class_1)))  #打印类别1的值
print('class_0的样本数量:     %s' %(len(class_0)))    #打印类别0的值


#用图显示数据的分布情况
plt.figure()
plt.scatter(class_1['num1'],class_1['num2'],c = 'r',marker = 'o',label = 'class_1')
plt.scatter(class_0['num1'],class_0['num2'],c = 'g',marker = 'x',label = 'class_0')
plt.legend(loc = 'best')
plt.show()

#为了取数据方便,将其改变为ndarray格式的数据
data = np.array(data)

data_x = data[:,0:3]
data_y = data[:,3:4]

print(data_x[:5])
print(data_y[:5])

#decent : 更新参数
#accuracy : 计算精度


#定义sigmoid函数

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

theta = np.zeros([3,1])
print(theta)


#model :返回预测值
def gfunc(x,t):
    return sigmoid(np.dot(x,t))
print((gfunc(data_x,theta)).shape)
#cost : 根据参数计算损失
def cost(x,y):
    s = np.multiply(y,np.log(x)) + np.multiply((1-y),np.log(1-x))
    sum =-np.sum(s)/len(x)
    return sum

#查看一下没有更新参数之前的函数损失值
print(cost(gfunc(data_x,theta),data_y))

#gradient:  计算每个参数的梯度方向
def gradient(data,x,y,t):
    c = (x - y).ravel()
    for i in range(len(t)):
        s = np.sum(np.multiply(c,data[:,i]))/len(x)
        t[i] = t[i] - 0.001*s
#for j in range(5000):
#    gradient(data_x,gfunc(data_x,theta),data_y,theta)
co = np.zeros([2,1])
print(co)
index = 0
while(True):
    gradient(data_x,gfunc(data_x,theta),data_y,theta)
    index += 1
    co[index%2,0] = cost(gfunc(data_x,theta),data_y)
    if index % 2 == 0:
        b = np.abs(co[0,0] - co[1,0])
        if b<= 1e-6:
            break
print(index)
print(co)

print(cost(gfunc(data_x,theta),data_y))

实验结果就是当迭代执行109902时,损失值由之前的0.69降到了0.37。

 

真正实践的时候,是不需要自己写这些函数的,写这些函数的目的主要是为了帮助理解内部的实现的原理。在真正使用的时候,我们所关注的是数据的处理。因为现实中的数据不是拿来就可以用的,如果不对数据进行人工处理,得到的结果肯定会不如意的,下面就为大家展示一下如果样本类别数据不一致的时候怎么处理,以及处理后的对比效果。

#对于这个程序就是为了试验对于样本不统一的时候选择不同的采样方法对logistic回归结果的影响。

#1:下采样

#2:过采样

本程序我选了下采样,也就是使类别多的数据向类别少的数据方向进行靠拢

# -*- coding: utf-8 -*-
"""
Created on Fri Aug 10 15:39:32 2018

@author: Fighting_Hou
"""

import numpy as np 
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.model_selection import train_test_split,KFold
from sklearn.preprocessing import StandardScaler
import warnings
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import recall_score,accuracy_score

warnings.filterwarnings(action = 'ignore')
data = pd.read_csv('creditcard.csv',header = 0)
#过程为:
#首先查看不同类别的数据个数
#然后选择不同的策略对数据进行处理
#最后顺便试验了一下惩罚因子的影响


#将数据分为特征和标记数据
X = data.ix[:,data.columns != 'Class']       #这个取值的方法还是不错的,想取哪列取哪列。
y = data.ix[:,data.columns == 'Class']
print(y.shape)
#对Amount数据进行均一化处理
scaler = StandardScaler().fit(X.loc[:,'Amount'].reshape(-1,1))      #对于归一化  首先是fit一个模型
X['Amount_new'] = scaler.transform(X.loc[:,'Amount'].reshape(-1,1)) #然后调用模型的 fit_transform
X = X.drop(['Time','Amount'],axis = 1)    #去除dataframe里不想要的列  传入list数组
#print(X.head())

#统计下不同类别的数量
positive_data = data[data.Class == 1]        #这里就是用到了  .列名。  很高级,竟然还可以用.属性
negative_data = data[data.Class == 0]
print('正例的数量为:  %s' % (len(positive_data)))
print('负例的数量为:  %s' % (len(negative_data)))
#首先采用下采样进行模型的训练
positive_index =np.array( positive_data.index)     #有必要 array一下  因为之后还要取值用
#print(positive_index)
negative_index = np.array(np.random.choice(negative_data.index,size = len(positive_data),replace = False))
print(len(negative_index))
connect = np.concatenate((positive_index,negative_index))
data_new = data.iloc[connect,:]   #以后知道index后,就是用iloc这个函数取值了
print(data_new.shape)
X1 = data_new.ix[:,data_new.columns != 'Class']
y1 = data_new.ix[:,data_new.columns == 'Class']
print(y1.shape)
#交叉验证
#首先将数据分为训练集和测试集
x_train,x_test,y_train,y_test = train_test_split(X1,y1,test_size = 0.3,random_state = 1)
print(len(x_train))
print(len(x_test))
param = [0.01,0.1,1,10,100]
flod = KFold(n_splits=5,shuffle = False)    #k折交叉验证
for para in range(len(param)):
    print('--------')
    print('惩罚因子为:  %s ' % (param[para]))
    cost = []  
    for index_x,index_y in flod.split(x_train):        
        model = LogisticRegression(C = param[para],penalty = 'l1')       
        model.fit(x_train.iloc[index_x,:],y_train.iloc[index_x,:].values.ravel())       
        y_hat = model.predict(x_train.iloc[index_y,:])       
        cost.append(recall_score(y_hat,y_train.iloc[index_y,:].values))       
        print('inter  %s的召回率为:%s ' % (para,recall_score(y_hat,y_train.iloc[index_y,:])))   
    print('平均  : %s' % (np.mean(cost)))       
#### 展示一下不处理数据时的召回率   证明这个召回率是明显低于 处理后的情况的    
x_train,x_test,y_train,y_test = train_test_split(X,y,test_size = 0.3,random_state = 1)    
model = LogisticRegression(C = 0.01,penalty = 'l1')
model.fit(x_train, y_train)
y_hat = model .predict(x_test)
print(recall_score(y_hat,y_test))
        
        
        
  

最后的输出结果是:   使用处理后的数据的召回率(衡量算法好坏的一种方法)为97%

                                      没有使用数据处理的召回率为:85%  而且识别的错误的数据很多。

  • 4
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值