理论部分
理解参考博客:逻辑回归的本质
假设函数hθ:
g(z)称为S型函数或逻辑函数(“Sigmoid Function,” also called the “Logistic Function”),图形如下:
这个假设函数hθ的含义是:输入x和参数θ后,y=1的概率。概率论上的数学表达式是:
损失函数:
这个Cost我们刚开始肯定想的是和线性回归一样,用误差平方和来当代价函数:
但是他的图像是这样的(左图):
这样的话梯度下降就找不到全局最低点,我们希望的图像应该是右边这样的。
于是我们就找其他的Cost的表达式,根据极大似然估计(还没学到,这里先占个坑),得出下面的Cost表达式:
将两个cost写在一起得到:
把这个cost表达式代入上面的Jθ,得到最终的损失函数的表达式:
这个式子的向量形式表示为:(写出这个形式是为了好编代码)
梯度下降:
还记得线性回归时候用的梯度下降是:
这里使用同样的形式,把Jθ代入上面的这个式子,求偏导之后得到:
向量形式为:
代码部分
初始参数设置
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
#显示中文(根据自己电脑设置)
plt.rcParams['font.sans-serif'] = ['SimHei']
#显示负数
plt.rcParams['axes.unicode_minus'] = False
导入数据
data = np.genfromtxt('data1.txt',delimiter=',')
print(data[:5])
x,y = data[:,:-1],data[:,-1]
数据样式(前五行):
数据预处理
#数据预处理
def preProcess(x,y):
#标准化特征缩放
x-=np.mean(x,axis = 0)
x/=np.std(x,axis=0,ddof=1)
#ddof eans Delta Degrees of Freedom. The divisor used in calculations is N - ddof, where N represents the number of elements. By default ddof is zero.
#在x前加上一列1 方便代码中的矩阵操作
x = np.c_[np.ones(len(x)),x]
y = np.c_[y]
#np.c_是按行连接两个矩阵,就是把两矩阵左右相加,要求行数相等。
return x,y
x,y = preProcess(x,y)
p.s. 这里的数据处理会直接使x和y发生改变,在稍后绘制决策边界的时候,不能直接使用这里的xy,而应该使用‘导入数据’模块里的原始数据。
定义函数
#定义模型
def model(x, theta):
z = np.dot(x,theta)
return 1/(1+np.exp(-z))
#损失函数
def costFunc(h,y):
m = len(y)
J = 1.0/m*np.sum(-np.dot(y.T,np.log(h))-np.dot((1-y).T,np.log(1-h)))
# J = (-1.0/m)*np.sum(y*np.log(h)+(1-y)*np.log(1-h))
return J
# 梯度下降
alpha = 0.01
iter_num = 20000
def graDesc(x,y,alpha,iter_num):
m,n = x.shape
#n行1列theta参数
theta = np.zeros((n,1))
#初始化代价值
J_history = np.zeros(iter_num)
#迭代更新
for i in range(iter_num):
h = model(x,theta)
J_history[i] = costFunc(h,y)
#梯度下降算法,更新θ
deletaTheta = 1.0/m*np.dot(x.T,h-y)
theta-=alpha*deletaTheta
return J_history,theta
跑起来
#调用函数
J_history,theta = graDesc(x,y,alpha,iter_num)
#得出预测值
h = model(x,theta)
画图
损失函数图
plt.title("cost function chart")
plt.plot(range(iter_num),J_history)
plt.xlabel('迭代次数')
plt.ylabel('cost')
plt.show()
决策边界曲线图
def plotDescisionBoundary(X,y,theta):
# 补充画当前样本散点图的代码
# 分类决策面 theta0+theta1*x1+theta2*x2 = 0
pos = np.where(y == 1)
#pos是y中数据等于1的下标索引
# print(pos)
neg = np.where(y==0)
#neg是y中数据等于0的下标索引
# print(neg)
x1 = np.arange(min(X[:, 1]), max(X[:, 1]), 0.1)
x2 = -(theta[1]*x1+theta[0])/theta[2]
plt.plot(x1,x2)
#python中数据可视化函数scatter(数据的横坐标向量,数据的纵坐标向量,marker='0'数据以点的形式显示,c='b'数据点是blue颜色)
plt.scatter(X[pos,0],X[pos, 1],marker='o', c='b')
plt.scatter(X[neg,0],X[neg, 1],marker='x', c='r')
plt.show()
plotDescisionBoundary(data[:,:-1],data[:,-1],theta)
检查迭代效果
这个图是我自己调参时,为了看迭代效果明显一点画的图,y轴的值为y(真实值)与h(预测值)的差,这个差离0越近,说明预测得效果越好。
plt.scatter(range(len(y)),y-h,label='预测值')
plt.plot(range(len(y)),[0 for _ in range(len(y))],label='参考线',color='red')
plt.xlabel('样本')
plt.ylabel('真实值与预测值的差值')
plt.legend()
plt.show()
总结
- 数学的欠缺让我推导公式有些困难,尤其是矩阵变换上,公式代码是照着吴恩达给的向量公式敲的,没有理解充分;
- 边界曲线绘制容易搞混,红点和蓝点是标签(y)==1和0时候的图,不是预测的结果,预测结果体现在决策边界(那条蓝色的直线)上;
- 边界曲线代码中
X[pos,0],X[pos, 1]
意思是标签(y)==1,然后x0和x1两个特征值。