代码实现
牛顿法的代码实现:
# 定义h函数
def h(x1, x2):
h = theta[0]+theta[1]*x1+theta[2]*x2
return 1/(1+math.exp(-h))
# 定义迭代次数
cnt = 0
# 定义f函数
def f(x1,x2):
h = theta[0]+theta[1]*x1+theta[2]*x2
return 1/((1+math.exp(h))**2)
def newton_method(x1, x2 , y):
global cnt
global error0
while True:
cnt += 1
diff = [0]*3
# 初始化hessian矩阵
Hessian = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
Hessian_mat = np.mat(Hessian,dtype=float)
# 计算l(theta)导数
for i in range(m):
diff[0] += y[i]-h(x1[i], x2[i])
diff[1] += (y[i]-h(x1[i], x2[i]))*x1[i]
diff[2] += (y[i]-h(x1[i], x2[i]))*x2[i]
# 计算hessian矩阵的每一个元素
for k in range(m):
Hessian_mat[(0, 0)] += (-f(x1[k], x2[k]))
Hessian_mat[(0, 1)] += (-f(x1[k], x2[k])) * x1[k]
Hessian_mat[(0, 2)] += (-f(x1[k], x2[k])) * x2[k]
Hessian_mat[(1, 0)] += (-f(x1[k], x2[k])) * x1[k]
Hessian_mat[(1, 1)] += (-f(x1[k], x2[k])) * x1[k] * x1[k]
Hessian_mat[(1, 2)] += (-f(x1[k], x2[k])) * x1[k] * x2[k]
Hessian_mat[(2, 0)] += (-f(x1[k], x2[k])) * x2[k]
Hessian_mat[(2, 1)] += (-f(x1[k], x2[k])) * x2[k] * x1[k]
Hessian_mat[(2, 2)] += (-f(x1[k], x2[k])) * x2[k] * x2[k]
Hessian_mat_I = Hessian_mat.I
theta[0] -= Hessian_mat_I[0]*np.mat(diff).T
theta[1] -= Hessian_mat_I[1]*np.mat(diff).T
theta[2] -= Hessian_mat_I[2]*np.mat(diff).T
error1 = 0
for i in range(m):
error1 += (h(x1[i], x2[i]) - y[i]) ** 2 / 2
print("cnt=%d,Hessian_mat_I=%f,error1=%f" % (cnt, np.linalg.det(Hessian_mat_I), error1))
if (abs(error1 - error0) < epsilon):
return (theta,cnt)
else:
error0 = error1
softmax回归的代码实现:
import numpy as np
import matplotlib.pyplot as plt
# 迭代次数
iter_num = 3500
# 设置学习率
lr = 0.0001
# 未使用正则化
lamba = 0
# 加载数据集,返回的特征集和标签集都是矩阵
def load_data(filename):
dataset = []
labelset = []
fr = open(filename)
for row in fr.readlines():
row = row.strip().split()
dataset.append([float(row[0]), float(row[1]), float(row[2])])
labelset.append([int(row[3])])
return np.mat(dataset), np.mat(labelset)
# 归一化特征集
def normalization(X):
m = X.shape
print(m)
Xmin = np.min(X, axis=0)
Xmax = np.max(X, axis=0)
Xmean = np.mean(X, axis=0)
X_norm = (X-Xmean)/(Xmax-Xmin)
return X_norm
# sigmoid函数
def sigmoid(x):
h = 1.0 / (1 + np.exp(-x))
return h
# 定义一个propab函数
def propab(i, j, x, theta):
divider_up = np.exp(theta[j] * x[i].T)[0][0, 0]
divider_down = np.sum(np.exp(theta * x[i].T))
return divider_up/divider_down
# 定义决策函数
def predict(w,x):
vector = w*x.T
y_ = np.argmax(vector)+1
return y_
def test(w, x, y):
m1, n1 = x.shape
cnt = 0
for i in range(m1):
y_ = predict(w, x[i])
if(y_ == y[i, 0]):
cnt += 1
accurancy = cnt/m1
return accurancy
# 定义一个loss函数
def calculate_loss(w, y, x, k):
loss = 0
error = 0
for i in range(m):
for j in range(k):
value = propab(i, j, x, w)
error = int(y[i, 0] == j+1) * np.log(value)
w_ravel = w.ravel()
loss += -error/m + lamba * (w_ravel * w_ravel.T)[0, 0] / m
return loss
# 定义一个softmax回归函数
def softmax_reg(datamat, labelmat, k, Iternum, alpha):
m, n = datamat.shape
print("m=", m, " n=", n)
# 记录正确率和损失值的变化
accur_list = []
loss_list = []
iter_list = []
# 初始化权重
w = np.mat(np.ones((k, n)))
diff = np.mat(np.ones((k, n)))
iter = 0
while(iter < Iternum):
for j in range(k):
for i in range(m):
value = propab(i, j, datamat, w)
diff[j] += -np.mat((int(labelmat[i, 0] == j+1) - value)*datamat[i])/m + lamba * w[j]
w[j] -= alpha*diff[j]
iter += 1
if(iter % 100 == 0):
print(iter)
accur = test(w, datamat, labelmat)
accur_list.append(accur)
loss = calculate_loss(w, labelmat, datamat, k)
loss_list.append(loss)
print("accurancy = ", accur)
print("loss = ", loss)
iter_list.append(iter)
return w, accur_list, loss_list, iter_list
if __name__ == "__main__":
data_mat, label_mat = load_data("web_yuehui.txt")
print(data_mat.shape)
print(label_mat.shape)
data_norm = normalization(data_mat)
m, n = data_norm.shape
offset = np.ones((m, 1))
trainmat = np.c_[offset, data_norm]
w_value, accurancy, loss, iter = softmax_reg(trainmat, label_mat, 3, iter_num, lr)
# 画图
plt.plot(iter, loss)
plt.xlabel("iter")
plt.ylabel("loss")
plt.savefig("loss_changing.png")
plt.show()