机器学习笔记（三）：线性回归大解剖（代码部分）
这里，让我手把手教你如何用逻辑回归分析数据
根据学生分数预测是否录取：
#必备3个库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
让我们读入数据：
import os
path = "data" + os.sep + "LogiReg_data.txt"; # os.sep 表示的是文件分割符 “ / ” 因为win和linux分隔符不同所以 为了跨平台就不直接写
pdData = pd.read_csv(path,header=None,names=['Exam1','Exam2','Admitted']); #header默认使用第一行做标题行，这里改为自己设置
pdData.shape; #数据的纬度 (100 , 3)
pdData.head() #不加参数表示默认显示5行
前两列数据表示分数，最后一列表示是否录取：
让我们简单看一下，录取和没录取在不同分数上的显示情况
positive = pdData[pdData['Admitted' ] == 1] # 录取的集合
nagative = pdData[pdData['Admitted' ] == 0] # 没录取集合
fig, ax=plt.subplots(figsize=(10,5)) # ax前面必须有fig figsize设置子图的宽和高
#分1 做x 分2做y
ax.scatter(positive['Exam1'],positive['Exam2'],s=30,c='blue',marker='o',label='Admitted') #s=点大小 c='颜色' marker=标记
ax.scatter(nagative['Exam1'],nagative['Exam2'],s=30,c='r',marker='x',label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam1 score')
ax.set_ylabel('Exam2 score'
plt.show()
可以看到有明显的分界线：
接下来就是正式的分析了
我们的目标是：建立分类器 ，求出3个参数（分数1的参数，分数2的参数，误差）
设定阈值，根据阈值判断是否录取（这里>=0.5就当作录取）
要完成的模块：
#1.sigmoid：映射到概率（0-1）的函数 g(z) = 1/(1+e^-z)
#2.model: 返回预测结果值
#3.cost: 根据参数计算损失
#4.gradient: 计算每个参数的梯度方向
#5.descent: 进行参数更新
#6.accuracy: 计算精度
1.sigmoid 函数
值域（0-1） g(-∞) = 0 g(+∞) = 1 g(0)=0.5
def sigmoid(z):
return 1/(1+np.exp(-z))
2.model
就是用数据作为sigmoid的输入，从而构造出预测函
数据 (1,x1,x2) * (θ0 , θ1 , θ2)的转置 = θ0+θ1X1+θ2X2
然后传到 sigmoid中得到结果
def model(X,theta):
return sigmoid( np.dot(X,theta.T) )
#结果中 θ0需要和1相乘，数据中缺少1这1列，所以要补进
pdData.insert(0,'Ones',1); #在第0列（第一列之前，索引之后 插入列名为‘ones’，所有数据为1的一列
data = pdData.as_matrix() #将pdData从DataFrame 转为矩阵，方便进行矩阵乘法
cols = data.shape[1] #数据有几列 4 ones exam1 exam2 admitte
X = data[:,0:cols-1] #数据矩阵X 1，分数1，分数2
y = data[:,cols-1:cols]; # 结果（是否选上
theta=np.zeros([1,3]) # 3个参数 参数个数=数据个数+1
3.损失函数
这里就是对数似然函数 * -1/m [加-是因为要把梯度上什转为梯度下降，1/m是为了求平均] 似然函数 = 每个 x 对应的 P(y|x;θ) 的连乘 ， 然后取对数得到对数似然函数
损失函数（loss function）是用来衡量模型的预测值f(x)与真实值Y的不一致程度，它是一个非负实值函数，损失函数越小，模型越优（还需考虑过拟合等问题）
这里选择用 似然函数*-1/m 作为损失函数，是不是因为 似然函数越大越好 ， 加一个 负号就越小越好 ， 所以只要平均似然函数的值越来越小了，就说明越来越好了?
损失函数用： J(θ) = -1/m * sum( (y*log(hθ(X)) + (1-y)log(1-hθ(x)) ) )
def cost(X,y,theta):
tt = y-model(X,theta)
left = np.multiply(tt,tt)
return np.sum(left)/(len(X))
用 1/m * sum( (y - hθ(X) )^2 ) 做损失函数 效果相同
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)
4.求梯度
多维图像的梯度，也就是多维图像在某点各个方向的斜率中斜率最大的那个
算法就是对J(θ) 求偏导 dJ(θ)/dθ = -1/m * sum( (yi - hθ(xi))xij )
def gradient(X,y,theta):
grad = np.zeros(theta.shape) #有几个θ，就找几个梯度
# dJ(θ)/dθ = -1/m * sum( (yi - hθ(xi))xij )
error = (model(X,theta) - y).ravel(); # yi - hθ(xi) 将-号提取进来 ravel将多维降为一维矩阵 ,因为原error是列向量，所以作用和专T类似
#求每个θi的梯度
for j in range (theta.shape[1]):
term = np.multiply(error,X[:,j]) #Xij 表示的是 X的每一个样本，针对每一个参数θj 都只取和θj有关的那一列即X（ ， j）行不限
grad[0,j] = np.sum(term) / len(X)
return grad
4 （附属） 停止策略
梯度更新需要有个头，有个结束点
#设计了3种不同的停止策略
STOP_ITER = 0 #按照迭代次数停止 (一次迭代即更新一次梯度值)
STOP_COST = 1 #目标函数两次迭代之间的损失（差异）很小
STOP_GRAD = 2 #两次迭代之间梯度没什么变化
def stopCriterion(type,value,threshold):
#设定3种不同的停止策略
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; #求2范数 sqrt(x1^2+x2^2+.....)
4 （附属） 对数据进行洗牌
打乱（防止自己收集数据的时候按照某种规律收集）
def shuffleData(data):
np.random.shuffle(data) #洗牌数据
#重新获取数据
cols = data.shape[1] #数据有几列 4 ones exam1 exam2 admitted
X = data[:,0:cols-1] #数据矩阵X 1，分数1，分数2
y = data[:,cols-1:cols]; # 结果（是否选上）
return X,y
5，6重点！
查看不同梯度下降策略对时间的影响
import time
def descent(data, theta, batchSize, stopType ,thresh, alpha): #batchSize = 1 随机梯度下降 =样本个数 梯度下降 =1-样本个数 小批量样本夏婧
init_time = time.time();#开始时候的时间 #thresh 阈值 alpha学习率
i = 0;#迭代次数
k = 0;#batch ，即当前已经消耗的样本数量
X,y = shuffleData(data)
grad = np.zeros(theta.shape) #梯度
costs = [cost(X,y,theta)] #损失值
#***5.descent: 进行参数更新***
while True:
grad = gradient(X[k:k+batchSize],y[k:k+batchSize],theta)
#更新参数 ， 让θ 不停的靠近梯度为0的点
theta = theta - alpha*grad;
costs.append(cost(X,y,theta)) #用新的 θ 计算出新的损失值
i += 1; #迭代次数 +1
k+=batchSize
if k >= len(X):
k=0;
X,y = shuffleData(data); #重新洗牌
value = 0;
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):
#import pdb; pdb.set_trace();
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 = "Mini-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,8))
ax.plot(np.arange(len(costs)), costs, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title(name.upper() + ' - Error vs. Iteration')
plt.show()
return theta
数据标准化
（（a-均值）/方差）
from sklearn import preprocessing as pp
scaled_data = data.copy()
scaled_data[:, 1:3] = pp.scale(data[:, 1:3])
n=100
#损失值 < 0.000001时停
theta = runExpe(scaled_data, theta, n, STOP_COST, thresh=0.000001, alpha=0.001)
6.accuracy: 计算精度
def predict(X,theda): #设定阈值 预测 大于0.5就是可以
return [1 if x>=0.5 else 0 for x in model(X,theta)]
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))