作者：桂。
時間：2017-04-21 21:11:23
鏈接：http://www.cnblogs.com/xingshansi/p/6743780.html
前言
看到最近大家都在用Tensorflow，一查才發現火的不行。想着入門看一看，Tensorflow使用手冊第一篇是基於MNIST的手寫數字識別的，用到softmax regression，而這個恰好與我正在看的《統計信號處理》相關。本文借此梳理一下：
1）羅傑斯特回歸
2）Softmax Regression
3）基於Tensorflow的MNIST手寫數字識別框架
內容為自己的學習記錄，其中多有借鑒他人的地方，最后一並給出鏈接。
一、羅傑斯特回歸（Logistic Regression）
A-問題描述
一般說的Logistic回歸，是指二分類問題。
對於連續的變量X，如果說它服從Logistic回歸，則：
對應概率密度以及分布函數：
F(x)更像是一個threshold，f(x)的表達式也利於求導。
定義Logistic 回歸模型：
為輸入，
為對應輸出，二項Logistic回歸對應如下模型：
聯系到之前寫的感知機：
分類時，對應概率：
為正的概率大，為負的概率小，這個指數有點像把正/負→進行指數處理，這樣就容易理解了。
有時為了表示方便，記：
，則表達式簡化為：
為什了叫Logistic regression呢？
這樣明顯是一個線性回歸的表達式。
B-理論分析
對於參數的求解，同感知機一樣，就是求解w。利用最大似然求解：
定義sigmoid函數：
准則函數重新寫為：
可以借助梯度下降或者牛頓法，這幾個方法之前已經有過介紹。
利用求解得到的w，即可以進行概率判斷，哪個概率大就判給哪個類別，整個思路與感知機也是一樣的：
C-理論應用
用一個code說明其工作原理，還是借助之前感知器用的例子，因為Logistic回歸可以認為是Softmax回歸的特例，這里借用softmax框架給出對應代碼：
%Logistic regression
clear all; close all; clc;
%Every category samples
N = 100;
X = [randn(N,2)+2*ones(N,2);...
randn(N,2)-2*ones(N,2)];
figure;
subplot 121
label = [ones(N,1);2*ones(N,1)];
plot(X(label==1,1),X(label==1,2),'b.');hold on;
plot(X(label==2,1),X(label==2,2),'r.');hold on;
%%Train
numClass = 2;
parm.learningRate = 0.5;
pram.iter = 5000;
parm.precost = 0;
X = [ones(size(X,1),1), X]';
[nfeatures, nsamples] = size(X);
W = zeros(nfeatures, numClass);
% 1, 2, 3 -> 001, 010, 100.
groundTruth = full(sparse(label, 1:nsamples, 1));
for i = 1:pram.iter
rsp = W'*X;
rsp = bsxfun(@minus, rsp, max(rsp, [], 1));
rsp = exp(rsp);
prob = bsxfun(@rdivide, rsp, sum(rsp));
% compute gradient based on current probability
grad = - X * (groundTruth - prob)' / nsamples ;
% update W
W(:,2:end) = W(:,2:end) - parm.learningRate * grad(:,2:end);
% compute cost
logProb = log(prob);
idx = sub2ind(size(logProb), label', 1:size(logProb, 2));
parm.cost = - sum(logProb(idx)) / nsamples;
if i~=0 && abs(parm.cost - parm.precost) / parm.cost <= 1e-4
break;
end
parm.precost = parm.cost;
end
%Predict
rsp = W' * X;
[~, prediction] = max(rsp, [], 1);
subplot 122
X = X(2:end,:)';
plot(X(prediction==1,1),X(prediction==1,2),'b.');hold on;
plot(X(prediction==2,1),X(prediction==2,2),'r.');hold on;
x = -5:.1:5;
b = W(1,2);
y = -W(2,2)/W(3,2)*x-b/W(3,2);
plot(x,y,'k','linewidth',2);
對應結果圖：
二、Softmax 回歸
Logistic回歸分析的是兩類，對於多類的情形，通常稱為：Softmax 回歸。基於概率所以soft，概率中最大值所以Max。
其實這個式子可以理解為Y=K時，$w_k*x = 0$，上面式子可以表述為，為了與下面表述一致，這里w改為$\theta$表示：
假設樣本個數為m,同logistic回歸類似，Softmax准則函數可以寫為：
進一步寫為：
思路同Logistic回歸，同樣利用梯度下降/牛頓-Rapson可以求解。求解出w：
其中，
分類的思路與Logistic回歸也是完全一致。另外在求解的時候，可以發現：
求解的時候，都減去一個數，是不影響結果的。前面提到Y=K時，$w_k*x = 0$，其實有一個為零就行，哪一類是無所謂的，所以優化只優化K-1，給出代碼：
%Test Softmax Regression
clear all; close all; clc;
N = 100;
X = [randn(N,2)+5*ones(N,2);...
randn(N,2)-1*ones(N,2);...
randn(N,2)-4*ones(N,2)];
figure;
subplot 121
label = [ones(N,1);2*ones(N,1);3*ones(N,1)];
plot(X(label==1,1),X(label==1,2),'b.');hold on;
plot(X(label==2,1),X(label==2,2),'r.');hold on;
plot(X(label==3,1),X(label==3,2),'k.');hold on;
title('Orignal Data')
%%Train
numClass = 3;
learningRate = 0.5;
X = [ones(size(X,1),1), X]';
[nfeatures, nsamples] = size(X);
W = zeros(nfeatures, numClass);
% 1, 2, 3 -> 001, 010, 100.
groundTruth = full(sparse(label, 1:nsamples, 1));
iter = 5000;
precost = 0;
for i = 1:iter
rsp = W'*X;
rsp = bsxfun(@minus, rsp, max(rsp, [], 1));
rsp = exp(rsp);
prob = bsxfun(@rdivide, rsp, sum(rsp));
% compute gradient based on current probability
grad = - X * (groundTruth - prob)' / nsamples ;
% update W
W(:,2:end) = W(:,2:end) - learningRate * grad(:,2:end);
% compute cost
logProb = log(prob);
idx = sub2ind(size(logProb), label', 1:size(logProb, 2));
cost = - sum(logProb(idx)) / nsamples;
if i~=0 && abs(cost - precost) / cost <= 1e-4
break;
end
precost = cost;
end
%Predict
rsp = W' * X;
[~, prediction] = max(rsp, [], 1);
subplot 122
X = X(2:end,:)';
plot(X(prediction==1,1),X(prediction==1,2),'b.');hold on;
plot(X(prediction==2,1),X(prediction==2,2),'r.');hold on;
plot(X(prediction==3,1),X(prediction==3,2),'k.');hold on;
title('Prediction')
對應結果（線性可分時，結果正確，線性不可分時，邊界點容易出現錯誤，如下圖）：
而實際應用中，我們不願將任何一個$\theta$直接置零，而是保留所有
，但此時我們需要對代價函數做一個改動：加入權重衰減。准則函數變為：
這個時候只要將上面對應的code按如下修改即可：
grad = - X * (groundTruth - prob)' / nsamples +.02 * W;%lamdba = 0.02
% grad = - X * (groundTruth - prob)' / nsamples ;
% update W
% W(:,2:end) = W(:,2:end) - learningRate * grad(:,2:end);
W = W - learningRate * grad;
三、基於Tensorflow的MNIST手寫數字識別
MNIST數據庫是一個手寫數字識別數據庫，
MNIST(Mixed National Institude of Standards and Technology database)，由幾萬張圖片組成，其中每個照片像素為28*28（在讀取時已經被拉成28*28=784的向量了），訓練數據集train：55000個樣本，驗證數據集validation：5000個樣本，測試數據集：10000個樣本。
首先分別讀取訓練和測試數據：
import numpy as np
import matplotlib.pyplot as plt
from tensorflow.examples.tutorials.mnist import input_data
mnist = input_data.read_data_sets("MNIST_data/", one_hot=True)
先看看數據長什么樣，任意讀取一張圖片：
#imshow data imgTol = mnist.train.images img = np.reshape(imgTol[1,:],[28,28]) plt.imshow(img)
大概就是這個樣子（其實是個gray圖）.
下面開始對數據進行基於Softmax Regression的處理。
首先創建一個會話（session），為什么需要創建session可以參考Tensorflow基本用法：
import tensorflow as tf sess = tf.InteractiveSession()
通過為輸入圖像和目標輸出類別創建節點，來開始構建計算圖：
sess = tf.InteractiveSession()
x = tf.placeholder("float", shape=[None, 784])
y_ = tf.placeholder("float", shape=[None, 10])
其中'None'表示該維度大小不指定。
我們現在為模型定義權重W和偏置b。可以將它們當作額外的輸入量，但是TensorFlow有一個更好的處理方式：變量。一個變量代表着TensorFlow計算圖中的一個值，能夠在計算過程中使用，甚至進行修改。在機器學習的應用過程中，模型參數一般用Variable來表示：
W = tf.Variable(tf.zeros([784,10])) b = tf.Variable(tf.zeros([10]))
變量需要通過seesion初始化后，才能在session中使用。這一初始化步驟為，為初始值指定具體值（本例當中是全為零），並將其分配給每個變量,可以一次性為所有變量完成此操作：
sess.run(tf.initialize_all_variables())
初始化完成后，接下來可以實現Softmax Regression了：
y = tf.nn.softmax(tf.matmul(x,W) + b)
softmax是tf.nn下的一個函數，tf.nn包含了很多神經網絡的組件，tf.matmul在tensorflow里的意義是矩陣相乘。給出損失函數（交叉熵）：
交叉熵沒聽過其實也不要緊，Maximum Likelihood（ML, 最大似然）估計應該都知道吧？
其中N為樣本總數，且
我們頻數以及頻率是觀測得到的，如何估計概率呢？我們通常用最大似然：
如果變形一下呢？同除以N：
最大似然估計不就是最小交叉熵估計？另一方面，交叉熵對誤差的衡量類似聯合概率密度的作用。給出對應code:
cross_entropy = -tf.reduce_sum(y_*tf.log(y))
reduce.sum為求和，reduce.mean為求取均值。有了准則函數以及理論框架，利用梯度下降進行訓練：
train_step = tf.train.GradientDescentOptimizer(0.01).minimize(cross_entropy)
這里利用小批量梯度mini-batch下降,batch = 50個訓練樣本:
for i in range(1100):
batch = mnist.train.next_batch(50)
train_step.run(feed_dict={x: batch[0], y_: batch[1]})
訓練后進行預測並打印：
correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_,1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, "float"))
print(sess.run(accuracy, feed_dict={x: mnist.test.images, y_: mnist.test.labels}))
給出總的code:
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow.examples.tutorials.mnist import input_data
mnist = input_data.read_data_sets("MNIST_data/", one_hot=True)
#imshow data
#imgTol = mnist.train.images
#img = np.reshape(imgTol[1,:],[28,28])
#plt.imshow(img)
sess = tf.InteractiveSession()
x = tf.placeholder("float", shape=[None, 784])
y_ = tf.placeholder("float", shape=[None, 10])
W = tf.Variable(tf.zeros([784,10]))
b = tf.Variable(tf.zeros([10]))
sess.run(tf.global_variables_initializer())
y = tf.nn.softmax(tf.matmul(x,W) + b)
cross_entropy = -tf.reduce_sum(y_*tf.log(y))
train_step = tf.train.GradientDescentOptimizer(0.01).minimize(cross_entropy)
for i in range(1100):
batch = mnist.train.next_batch(50)
train_step.run(feed_dict={x: batch[0], y_: batch[1]})
correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_,1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, "float"))
print(sess.run(accuracy, feed_dict={x: mnist.test.images, y_: mnist.test.labels}))
預測結果為：91.56%.
參考：
- Softmax regression:http://deeplearning.stanford.edu/wiki/index.php/Softmax_Regression
- 李航：《統計學習方法》
- Tensorflow官網文檔：https://www.tensorflow.org/get_started/mnist/beginners
- http://tensorfly.cn/tfdoc/tutorials/mnist_pros.html