详细原理推导先看这篇https://blog.csdn.net/To_be_to_thought/article/details/81511529
与线性回归的梯度下降法类似。代码如下:
import numpy as np
import math
from matplotlib import pyplot as plt
from sklearn.linear_model import LogisticRegression#用现成的库进行对比试验
from sklearn.preprocessing import StandardScaler
from sklearn import datasets
def cost_function(data,labels,positive_probability):
n,p=data.shape
a=labels*positive_probability
m=1-positive_probability
b=(1-labels)*np.log(m)
c=a+b
return -1/n*np.sum(c,axis=0)[0]
#sigmoid函数计算事件为正的概率
def sigmoid(data,beta):
inner_product=np.dot(data,beta)
probability=1/(1+np.exp(-inner_product))
return probability
#随机梯度下降求解logistic回归参数
#迭代终止条件: 1:达到最大迭代次数 2:前后两次梯度变化小于一个极小值 3:迭代前后损失函数值变化极小
def LogisticRegByBGD(data,labels,eps1,eps2,max_iter,alpha):
n,p = data.shape
labels = labels.reshape((n,1))
beta = np.random.rand(p+1,1) #参数初始化,(p+1)*1 #np.random.rand(p+1,1)
data = np.column_stack((data,np.ones(n))) #数据集最后一列添加一列1,n*(p+1)
postive_probability = sigmoid(data,beta) #每个样本为正样本的概率向量 (p+1)*1
iter=0
cost=0
while True:
new_beta=beta.copy()
for i in range(p+1):
y1=data[:,i].reshape((n,1))*(postive_probability-labels)
b = alpha * 1/n * np.sum(y1,axis = 0)[0]
new_beta[i,0] = new_beta[i,0] - b
new_cost=cost_function(data,labels,postive_probability)
gradient_change=np.sum(np.abs(new_beta-beta),axis = 0)[0]
iter += 1
print("迭代次数: "+str(iter))
print(new_beta)
print("梯度变化量: "+str(gradient_change))
print("损失函数变化量:"+str(abs(new_cost-cost)))
if iter > max_iter or abs(new_cost-cost)<eps1 or gradient_change < eps2:
break
beta=new_beta
postive_probability = sigmoid(data,beta)
return beta
#根据预测的正事件概率计算负事件概率,比较两者大小确定类别
def classification(positive_probability):
negative_probability=1-positive_probability
probability=np.column_stack((negative_probability,positive_probability))
return np.argmax(probability,axis=1)
测试代码使用MINIST数据集进行,代码如下:
digits=datasets.load_digits()
X, y = digits.data, digits.target
X = StandardScaler().fit_transform(X)
y = (y > 4).astype(np.int)#将大于4和小于4的数字分两组并二值化
max_iter=20000
eps=1e-1
eps1=1e-6
eps2=1e-6
alpha=1e-2
beta=LogisticRegByBGD(X , y , eps1 , eps2 , max_iter , alpha)
X=np.column_stack((X,np.ones(X.shape[0])))
predict=classification(sigmoid(X,beta))
print("训练集分类精度"+str((predict==y).astype(np.int).mean()))
结果如下:
再来看成熟的sklearn库的代码:
L1_LR=LogisticRegression(C=0.1, penalty='l1', tol=0.01)
L1_LR.fit(X,y)
L1_LR_predict=L1_LR.predict(X)
print((L1_LR_predict==y).astype(np.int).mean())
精度为0.903171953255,看上去差不多,但我写的收敛速度比较慢,因为它用了算法复杂度更加低的优化算法