Python机器学习笔记 Logistic Regression

可能我们最常听到的是梯度下降算法，它与这里的梯度上升算法是一样的，只是公式中的
加法需要变成减法，梯度上升算法用来求函数的最大值，而梯度下降算法是用来求函数的最小值

每个回归系数初始化为1重复R次： 计算整个数据集的梯度 使用alpha * gradient 更新回归系数的向量 返回回归系数

#_*_coding:utf-8_*_
from numpy import *# 读取数据
def loadDataSet(filename):
'''
对于testSet.txt，每行前两个值分别是X1和X2，第三个值数据对应的类别标签
而且为了设置方便，该函数还将X0的值设置为1.0
:return:
'''
dataMat = []
labelMat = []
fr = open(filename)
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat,labelMatdef sigmoid(inX):
return 1.0/(1+exp(-inX))def gradAscent(dataMatIn,classLabels):
'''
:param dataMatIn: 是一个2维Numpy数组，每列分别代表每个不同的特征
每行则代表每个训练样本。
:param classLabels: 是类别标签，是一个1*100的行向量，为了便于矩阵运算，需要将行向量
转换为列向量，就是矩阵的转置，再将其赋值与labelMat。
:return:
'''
dataMatrix = mat(dataMatIn)
labelMat = mat(classLabels).transpose()
# labelMat = mat(classLabels).T
m,n = shape(dataMatrix)
# alpha是向目标移动的步长
alpha = 0.001
# 迭代次数
maxCycles = 500
weights = ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix*weights)
error = (labelMat-h)
weights = weights + alpha*dataMatrix.transpose()*error
return weights

if __name__ == '__main__':
filename = 'testSet.txt'
dataArr,labelMat = loadDataSet(filename)
weights_res = gradAscent(dataArr,labelMat)
print(weights_res)'''
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
'''

def plotBestFit(wei):
import matplotlib.pyplot as plt
weights = wei.getA()
dataMat,labelMat = loadDataSet(filename)
dataArr = array(dataMat)
n = shape(dataArr)[0]
xcord1 = []
ycord1 = []
xcord2 = []
ycord2 = []
for i in range(n):
if int(labelMat[i]) ==1:
xcord1.append(dataArr[i,1])
ycord1.append(dataArr[i,2])
else:
xcord2.append(dataArr[i, 1])
ycord2.append(dataArr[i, 2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1,ycord1,s=30,c='red',marker='s')
ax.scatter(xcord2,ycord2,s=30,c='green')
x = arange(-3.0,3.0,0.1)
y = (-weights[0]-weights[1] * x) / weights[2]
ax.plot(x,y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()

if __name__ == '__main__':
filename = 'testSet.txt'
dataArr,labelMat = loadDataSet(filename)
weights_res = gradAscent(dataArr,labelMat)
print(weights_res)
plotBestFit(weights_res)

所有回归系数初始化为1对数据集中每个样本 计算该样本的梯度 使用alpha*gradient 更新回归系数值返回回归系数值

# 随机梯度上升算法
def stocGradAscent0(dataMatrix,classLabels):
m,n = shape(dataMatrix)
alpha = 0.01
weights = ones(n)
for i in range(m):
h = sigmoid(sum(dataMatrix[i]*weights))
error = classLabels[i] - h
weights = weights + alpha*error*dataMatrix[i]
return weights

if __name__ == '__main__':
filename = 'testSet.txt'
dataArr,labelMat = loadDataSet(filename)
weights_res = stocGradAscent0(array(dataArr),labelMat)
print(weights_res)
plotBestFit(weights_res)

# 改进的随机梯度上升算法
def stocGradAscent1(dataMatrix,classLabels,numIter=150):
m,n = shape(dataMatrix)
weights = ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.01
randIndex = int(random.uniform(0,len(dataIndex)))
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights + alpha *error*dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights

if __name__ == '__main__':
filename = 'testSet.txt'
dataArr,labelMat = loadDataSet(filename)
weights_res = stocGradAscent1(array(dataArr),labelMat)
print(weights_res)
plotBestFit(weights_res)

weights_res = stocGradAscent1(array(dataArr),labelMat,500)

weights = weights + alpha *error*dataMatrix[randIndex]

weights = weights

#-*_coding:utf-8_*_
import math
from numpy import *# Sigmoid函数的计算
def sigmoid(inX):
return 1.0/(1+exp(-inX))# 改进的随机梯度上升算法
def stocGradAscent1(dataMatrix,classLabels,numIter=150):
m,n = shape(dataMatrix)
weights = ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.01
randIndex = int(random.uniform(0,len(dataIndex)))
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights + alpha *error*dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights# Logistic回归分类函数
def classifyVetor(inX,weights):
prob = sigmoid(sum(inX*weights))
if prob > 0.5:
return 1.0
else:
return 0.0def colicTest(filetrain,filetest):
frTrain = open(filetrain)
frTest = open(filetest)
trainingSet = []
trainingLabeles = []
for line in frTrain.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabeles.append(float(currLine[21]))
trainWeights = stocGradAscent1(array(trainingSet),trainingLabeles,500)
errorCount = 0
numTestVec = 0
for line in frTest.readlines():
numTestVec += 1.0
currLine = line.strip().split('\t')
lineArr =[]
for i in range(21):
lineArr.append(float(currLine[i]))
if int(classifyVetor(array(lineArr),trainWeights)) != int(currLine[21]):
errorCount += 1
errorRate = (float(errorCount)/numTestVec)
print('the error rate of this test is %f'%errorRate)
return errorRatedef multTest(filetrain,filetest):
numTests = 6
errorSum = 1.0
for k in range(numTests):
errorSum += colicTest(filetrain,filetest)
print('after %d iterations the average error rate is %f'%(numTests,errorSum/float(numTests)))if __name__ == '__main__':
filetrain = 'horseColicTraining.txt'
filetest = 'horseColicTest.txt'
multTest(filetrain,filetest)

the error rate of this test is 0.373134
the error rate of this test is 0.358209
the error rate of this test is 0.402985
the error rate of this test is 0.432836
the error rate of this test is 0.462687
the error rate of this test is 0.343284
after 6 iterations the average error rate is 0.562189

import numpy as np
import pandas as pd
import matplotlib.pyplot as pltpath = 'LogiReg_data.txt'
pdData = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])
print(pdData.head(10))

positive = pdData[pdData['Admitted'] == 1]
negative = pdData[pdData['Admitted'] == 0]
fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(positive['Exam 1'], positive['Exam 2'], s=30, c='b', marker='o', label='Admitted')
ax.scatter(negative['Exam 1'], negative['Exam 2'], s=30, c='r', marker='x', label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')

def sigmoid(z):
return 1 / (1 + np.exp(-z))# 如果想查看Sigmoid函数，可以运行下面代码
nums = np.arange(-10, 10, step=1)
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(nums, sigmoid(nums), 'r')

# 增加一行1 在第0列
pdData.insert(0, 'Ones', 1)
print(pdData.head())
# 设置X和Y 即训练数据和训练变量
orig_data = pdData.values
cols = orig_data.shape[1] # (100, 4)
X = orig_data[:, 0:cols-1] # 等价于 np.matrix(X.values)
Y = orig_data[:, cols-1:cols] # 等价于 np.matrox(data.iloc[:, 3:4].value# 传递numpu矩阵 并初始化参数矩阵 theta
theta = np.zeros([1, 3]) # [[0. 0. 0.]]

def cost(X, y, theta):
left = np.multiply(-y, np.log(model(X, theta)))
right = np.multiply(1-y, np.log(1 - model(X, theta)))
return np.sum(left - right) / (len(X))

def gradient(X, y, theta):
grad = np.zeros(theta.shape)
error = (model(X, theta)-y).ravel()
for j in range(len(theta.ravel())): # for each parmeter
term = np.multiply(error, X[:, j])
grad[0, j] = np.sum(term) / len(X)
return grad

# 比较三种不同梯度下降方法
STOP_ITER = 0
STOP_COST = 1
STOP_GRAD = 2def stopCriterion(type, value, threshold):
# 设定三种不同的停止策略
if type == STOP_ITER:
return value > threshold
elif type == STOP_COST:
return abs(value[-1]-value[-2]) < threshold
elif type == STOP_GRAD:
return np.linalg.norm(value) < threshold# 洗牌
def shuffleData(data):
np.random.shuffle(data)
cols = data.shape[1]
X = data[:, 0:cols-1]
y = data[:, cols-1:]
return X, ydef descent(data, theta, batchSize, stopType, thresh, alpha):
# 梯度下降求解
init_time = time.time()
i = 0 # 迭代次数
k = 0 # batch
X, y = shuffleData(data)
grad = np.zeros(theta.shape) # 计算的梯度
costs = [cost(X, y, theta)] # 损失值 while True:
grad = gradient(X[k:k+batchSize], y[k:k+batchSize], theta)
k += batchSize # 取batch数量个数据
if k>=n:
k = 0
X, y = shuffleData(data) # 重新洗牌
theta = theta = alpha*grad # 参数更新
costs.append(cost(X, y, theta)) # 计算新的损失
i += 1 if stopType == STOP_ITER:
value = i
elif stopType == STOP_COST:
value = costs
elif stopType == STOP_GRAD:
value = grad
if stopCriterion(stopType, value, thresh):
break
return theta, i-1, costs, grad, time.time()-init_time

def runExpe(data, theta, batchSize, stopType, thresh, alpha):
theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh, alpha)
name = 'Original' if (data[:, 1] > 2).sum() > 1 else 'Scaled'
name += 'data - learning rate: {} -'.format(alpha)
if batchSize == n:
strDescType = 'Gradient'
elif batchSize == 1:
strDescType = 'Stochastic'
else:
strDescType = 'Min-batch({})'.format(batchSize)
name += strDescType + 'descent - Stop: '
if stopType == STOP_ITER:
strStop = "{} iterations".format(thresh)
elif stopType == STOP_COST:
strStop = "costs change < {}".format(thresh)
else:
strStop = 'gradient norm < {}'.format(thresh)
name += strStop
print("***{}\nTheta: {} - Iter:{}-Last cost: {:03.2f}--Duration: {:03.2f}s".format(
name, theta, iter, costs[-1], dur
))
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(np.arange(len(costs)), costs, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title(name.upper() + '- Error vs. Iteration')
return theta

if __name__ == '__main__':
# 设定迭代次数 选择的梯度下降方法是基于所有样本的
n = 100
# runExpe(orig_data, theta, n, STOP_ITER, thresh=5000, alpha=0.000001)
# 根据损失值停止 设定阈值 1e-6 差不多需要110 000 次迭代
# runExpe(orig_data, theta, n, STOP_COST, thresh=0.000001, alpha=0.001)
# 根据梯度变换停止 设定阈值 0.05，差不多需要 40000次迭代
# runExpe(orig_data, theta, n, STOP_GRAD, thresh=0.05, alpha=0.001)
# 对比不同的梯度下降方法 stochastic descent 这个不收敛
# runExpe(orig_data, theta, 1, STOP_ITER, thresh=5000, alpha=0.001)
# 我们将学习率调小，提高thresh再试试，发现收敛了
# runExpe(orig_data, theta, 1, STOP_ITER, thresh=15000, alpha=0.000001)
# Mine-batch descent 这个也不收敛 ，我们将alpha调小 0.000001 发现收敛了
runExpe(orig_data, theta, 16, STOP_ITER, thresh=15000, alpha=0.001)

scaled_data = orig_data.copy()
scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3])
runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.01)

scaled_data = orig_data.copy()
scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3])
# runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.01)
runExpe(scaled_data, theta, n, STOP_GRAD, thresh=0.02, alpha=0.01)

runExpe(scaled_data, theta, 1, STOP_GRAD, thresh=0.002/5, alpha=0.01)

# 设置阈值
def predict(x, theta):
return [1 if x >= 0.5 else 0 for x in model(X, theta)]if __name__ == '__main__':
# 设定迭代次数 选择的梯度下降方法是基于所有样本的
n = 100
scaled_data = orig_data.copy()
scaled_X = scaled_data[:, :3]
y = scaled_data[:, 3]
predictions = predict(scaled_X, theta)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for a, b in zip(predictions, y)]
accuracy = (sum(map(int, correct)) % len(correct))
print('accuracy = {0}%'.format(accuracy))
print(len(correct))

#_*_ codingLutf-8_*_from sklearn.linear_model import LogisticRegressiondef colicSklearn(filetrain,filetest):
frTrain = open(filetrain)
frTest = open(filetest)
trainingSet = []
trainingLabels = []
testSet = []
testLabels = []
for line in frTrain.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[-1]))
for line in frTest.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
testSet.append(lineArr)
testLabels.append(float(currLine[-1]))
classifier = LogisticRegression(solver='liblinear',max_iter=20).fit(trainingSet,trainingLabels)
test_accurcy = classifier.score(testSet,testLabels)*100
print("正确率为%s%%"%test_accurcy)if __name__ == '__main__':
filetrain = 'horseColicTraining.txt'
filetest = 'horseColicTest.txt'
colicSklearn(filetrain,filetest)

正确率为73.13432835820896%

#_*_coding:utf-8_*_
import numpy as np
import matplotlib.pyplot as pltdef sigmoid(z):
return 1.0/(1.0 + np.exp(-z))z = np.arange(-7, 7, 0.1)
print(z)
phi_z = sigmoid(z)def cost_1(z):
return -np.log(sigmoid(z))def cost_0(z):
return -np.log(1-sigmoid(z))z = np.arange(-10, 10, 0.1)
phi_z = sigmoid(z)c1 = [cost_1(x) for x in z]
plt.plot(phi_z, c1, label='J(w) if y=1')c0 = [cost_0(x) for x in z]
plt.plot(phi_z, c0, linestyle='--', label='J(w) if y=0')plt.ylim(0.0, 5.1)
plt.xlim([0, 1])
plt.xlabel('$\phi$(z)')
plt.ylabel('J(w)')
plt.legend(loc='best')
plt.tight_layout()
# plt.savefig('./figures/log_cost.png', dpi=300)
plt.show()

from IPython.display import Image# Added version check for recent scikit-learn 0.18 checks
from distutils.version import LooseVersion as Version
from sklearn import __version__ as sklearn_versionfrom sklearn import datasets
import numpy as npiris = datasets.load_iris()
# http://scikit-learn.org/stable/auto_examples/datasets/plot_iris_dataset.html
X = iris.data[:, [2, 3]]
print(X.shape)
y = iris.target # 取species列，类别if Version(sklearn_version) < '0.18':
from sklearn.cross_validation import train_test_split
else:
from sklearn.model_selection import train_test_split# train_test_split方法分割数据集
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.3, random_state=0)from sklearn.preprocessing import StandardScalersc = StandardScaler() # 初始化一个对象sc去对数据集作变换
sc.fit(X_train) # 用对象去拟合数据集X_train，并且存下来拟合参数
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)from sklearn.linear_model import LogisticRegressiondef sigmoid(z):
return 1.0 / (1.0 + np.exp(-z))lr = LogisticRegression(C=1000.0, random_state=0)
# http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression
lr.fit(X_train_std, y_train)# 计算该预测实例点属于各类的概率
lr.predict_proba(X_test_std[0, :].reshape(1, -1))
# Output:array([[ 2.05743774e-11, 6.31620264e-02, 9.36837974e-01]])# 验证predict_proba的作用
c = lr.predict_proba(X_test_std[0, :].reshape(1, -1))
# c[0, 0] + c[0, 1] + c[0, 2]
# Output:0.99999999999999989# 查看lr模型的特征系数
lr = LogisticRegression(C=1000.0, random_state=0)
# http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression
lr.fit(X_train_std, y_train)
print(lr.coef_)
# Output:[[-7.34015187 -6.64685581]
# [ 2.54373335 -2.3421979 ]
# [ 9.46617627 6.44380858]]# 验证predict_proba工作原理
Zz = np.dot(lr.coef_, X_test_std[0, :].T) + lr.intercept_
np.array(sigmoid(Zz)) / sum(np.array(sigmoid(Zz)))
# Output:array([ 2.05743774e-11, 6.31620264e-02, 9.36837974e-01])
# 此结果就是预测实例点各类的概率

from sklearn import datasets
from numpy import *
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_splitdef colicSklearn():
iris = datasets.load_iris()
X = iris.data
Y = iris.target
trainingSet,testSet,trainingLabels,testLabels = train_test_split(X,Y,test_size=0.25,random_state=40)
classifier = LogisticRegression(solver='sag', max_iter=5000).fit(trainingSet, trainingLabels)
test_accurcy = classifier.score(testSet, testLabels) * 100
print("正确率为%s%%" % test_accurcy)if __name__ == '__main__':
colicSklearn()

正确率为100.0%

#_*_coding:utf-8_*_
'''
下面这个例子，从数据产生，到数据提取，数据标准化
模型训练和评估来说明各个API的调用过程
'''
print(__doc__)
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# import matplotlib as mpl
import pylab as mpl# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False# 定义目标函数通过改函数产生对应的y
# y = 1 *x[0] + 2 *x[1] + ... (n+1) *x[n]
def l_model(x):
params = np.arange(1, x.shape[-1] +1)
y = np.sum(params *x) + np.random.randn(1) * 0.1
return y# 定义数据集
x = pd.DataFrame(np.random.rand(500,6))
# print(x)
y = x.apply(lambda x_rows:pd.Series(l_model(x_rows)),axis=1)
# print(y)# 划分数据集
x_train, x_test, y_train, y_test = train_test_split(x,y,test_size=0.3,random_state=2)# 数据标准化
ss = StandardScaler()
x_train_s = ss.fit_transform(x_train)
x_test_s = ss.fit_transform(x_test)# 输出下元数据的标准差和平均数
print(ss.scale_)
print(ss.mean_)# 训练模型
lr = LinearRegression()
lr.fit(x_train_s , y_train)
# 训练后的输入端模型系数，如果label有两个，即y值有两列，那么是一个2D的array
print(lr.coef_)
# 截距
print(lr.intercept_)# 用模型预测,并计算得分
y_predict = lr.predict(x_test_s)
test_accuracy = lr.score(x_test_s, y_test)print("正确率为%s%%" % test_accuracy)# 预测值和实际值画图比较
t = np.arange(len(x_test_s))
# 建一个画布，facecolor是背景色
plt.figure(facecolor='W')
plt.plot(t, y_test, 'r-', linewidth= 2, label = '真实值')
plt.plot(t, y_predict, 'b-', linewidth= 1, label = '预测值')
# 显示图例，设置图例的位置
plt.legend(loc= 'upper left')
plt.title("线性回归预测真实值之间的关系", fontsize = 20)
# 加网格
plt.grid(b = True)
plt.show()

下面这个例子，从数据产生，到数据提取，数据标准化
模型训练和评估来说明各个API的调用过程[0.28299035 0.3138177 0.29522932 0.30125841 0.28402977 0.28989605]
[0.53282546 0.51045318 0.48925819 0.5202345 0.4989613 0.54335117]
[[0.29202621 0.5787407 0.87280326 1.12786483 1.50043119 1.73271735]]
[10.38441274]正确率为0.9706613427158242%