一文打盡：線性回歸和邏輯斯蒂線性回歸

01-24

最近忙著弄畢業論文實驗和一個競賽，好久沒發文章了，想著該發一篇了吧。。。正好最近也在看Hands-On Machine Learning with Scikit-Learn and TensorFlow這本書，於是將看的內容和之前的知識做了一個總結，打算陸續寫一個tensorflow+機器學習的系列。本文是第一篇。

線性模型是機器學習中最基礎的模型，它形式簡單、易於建模，許多非線性模型可以在線性模型的基礎上通過引入層級結構或高維映射而得。如下圖所示，紅色的線條就是一條最簡單的線性模型。給定一些點集D={(x1,y1),(x2,y2),(x3,y3),...,(xN,YN)}，我們希望通過這些點去學習一個線性預測模型，當給出 $x_{N+1}$ 時，我們可以通過這個模型預測出 $y_{N+1}$ 的值。

上圖的紅色線條可以用 $y=ax+b$ 來表示，要學習這個線性模型，就得學習a和b這兩個參數。這是一個二元一次方程，而我們的線性模型可以是多元一次方程或者多元高次方程。

1、線性模型的基本形式

在上述例子中，我們的一個點就是一個坐標，表示的是訓練集中的實例有且只有1個屬性。不失一般性，如果實例有d個屬性，將實例描述為 $x=(x_1;x_2;...;x_d)$ ，其中xi是x在第i個屬性上的取值，線性模型試圖學得一個通過屬性的線性組合來進行預測的函數，即：

$f(x)=w_1x_1+w_2x_2+...+w_dx_d+b$

上式就是一個多元一次線性方程，因此也叫多元線性回歸。一般使用向量形式寫成：

$f(x)=w^Tx+b$

其中 $w=(w_1;w_2;...;w_d)$ ，w和b就是模型需要學習的參數。

2、線性回歸模型 linear regression

假定有數據集 $D={(x_1,y_1),(x_2,y_2),(x_3,y_3),...,(x_N,y_N)}$ ,其中 $x_i=(x_{i1};x_{i2};...;x_{id})$ , $y_i in R$ 。線性回歸模型就要試圖學得一個線性模型來儘可能準確的預測實值輸出標記。對於屬性 $x_i$ ，它的取值可能是連續的或者離散的變數，如果是離散的變數，若屬性間存在『序』的關係，可通過連續化將其轉化為連續值，比如{高，矮}轉化為{1.0，0.0}；對於不存在序關係的離散變數可進行獨熱編碼轉化，如{黃瓜，冬瓜，西瓜}可以轉化為（0，0，1），（0，1，0），（1，0，0）。對於數據集D，線性模型試圖學得 $f(x_i)=wx_i+b$ ，使得 $f(x_i)simeq y_i$

那麼怎麼學呢？我們知道，線性模型學習的就是參數w和b，將所有實例x的值代入上式，我們可以得到N個 $f(x_i)$ 值，為了讓得到的 $f(x_i)$ 和 $y_i$ 儘可能的接近，我們需要最小化 $|f(x_i)-y_i|$ 的累積和，用數學公式表達如下：

$(w^*,b^*)=arg minsum_{i=1}^{N}{(f(x_i)-y_i)^2}=arg minsum_{i=1}^{N}{(y_i-wx_i-b)^2}$

上式中 $(f(x_i)-y_i)^2$ 表示均方誤差，它對應著歐幾里德距離，我們的目的就是讓所有點到我們學習到的線條的歐式距離之和最小，也就是誤差最小化。因此上式也是線性模型的損失函數。求解 $w^*,b^*$ 使得損失函數最小化的過程稱為線性回歸模型的最小二乘「參數估計」，將上式分別對w和b求導，得到

$frac {partial E}{partial w} = 2(wsum_{i=1}^{N}{x_i ^2}-sum_{i=1}^{N}{(y_i-b)x_i })$

$frac {partial E}{partial b} = 2(Nb-sum_{i=1}^{N}{(y_i-wx_i) })$

然後令以上兩式為0，既可以得到w和b的最優閉式解。

對於多元線性回歸模型，變數比較多，數據集D可以表示成一個N×(d+1)的矩陣X。當 $X^TX$ 為滿秩矩陣或者正定矩陣失，通過最小二乘的估算，我們可以得到

$w^*=(X^TX)^{-1}X^Ty$

此時w*有且僅有一個解。那麼學得的線性回歸模型為

$f(x_i)=x_i^T(X^TX)^{-1}X^Ty$

然而現實中， $X^TX$ 往往不是滿秩矩陣，許多任務中的實例可能有大量變數，甚至超過實例的總數目，導致X的列數多於行數，此時顯然不是滿秩矩陣，此時可以解出多個w*，他們搜能使均方誤差最小，那麼選擇哪一個解作為輸出呢？這由演算法的學習偏好決定，長用的方法是加入正則化項。

3、對數線性回歸模型

線性模型的目的就是讓預測值逼近真實標記y，那麼可否令預測值逼近y的衍生物呢？假設我們認為訓練集中的示例的輸出標記y是在指數尺度上變化的，那就可以將輸出標記的對數作為線性模型的目標，即：

$ln y=w^T+b$

這就是「對數線性回歸」，更一般的，考慮單調可微函數g(.)，令

$y=g^{-1}(w^T+b)$

這樣得到的模型稱為「廣義線性模型」，顯然對數線性回歸是廣義線性模型在g(.)=ln(.)時的特例。

4、邏輯斯蒂線性回歸

線性回歸模型進行的是回歸學習，也就是輸出值的預測問題，那麼如果要做的任務是分類呢？只需要找到一個單調可微函數將分類任務的真實標記y與線性回歸模型的預測值聯繫起來就可以了。因此邏輯斯蒂回歸雖然名字叫回歸，但實際上它是一個分類演算法。

考慮二分類任務，假設輸出標記 $yin {0,1}$ ,對於線性回歸模型產生的預測值 $z=w^Tx+b$ ,z是一個實值，我們需要將它轉化為0或者1。在邏輯斯蒂回歸模型中，使用的轉化函數是sigmod函數，如下圖所示：

上圖橫軸是z，縱軸是y，可以看到sigmod函數非常漂亮，它是連續可微的，而且當z大於0時y趨近於1,小於0時趨近於0，因此我們可以設定當線性回歸模型的預測值大於0時我們將實例判定為正例，反之為反例。因此邏輯斯蒂回歸模型可以表達如下：

$y=frac {1}{1+e^{-(w^Tx+b)}}$

上式可以解讀為，用線性回歸模型的預測結果去逼近真實標記的對數幾率，因此，也叫對數幾率回歸模型（logistic regression）。

將上式做變化如下：

$ln frac{y}{1-y}=w^Tx+b$

若將y視為樣本x作為正例的可能性，則1-y是其反例可能性，兩者的比值稱為幾率，反映了x作為整例的相對可能性。對幾率取對數澤的到對數幾率（log odds）。

邏輯斯蒂回歸方法有許多優點，它是直接對分類可能性進行建模，無需事先假設數據分布，這樣就避免了假設分布不準確所帶來的問題（和樸素貝葉斯相比）；它不僅預測出類別，而是可以得到近似概率預測，這對於許多需要利用概率輔助決策的任務很有用；此外，對率函數是任意階可導函數，具有很好的數學性質，有許多數值優化演算法可以用於求最優解。

邏輯斯蒂回歸模型參數的估計

在上文中，我們將y視為樣本x作為正例的可能性，將1-y視為其反例的可能性，那麼可以有如下表達式：

$y=P(Y=1|x)$

$1-y=P(Y=0|x)$

結合等式 $y=frac {1}{1+e^{-(w^Tx+b)}}$ ，我們有

$P(Y=1|x)=frac{e^{(wx+b)}}{1+e^{wx+b}}$

$P(Y=0|x)=frac{1}{1+e^{wx+b}}$

有時為了方便，將權值向量和輸出向量加以擴充，仍記為w,x，即 $w=(w^1,w^2,..,w^n,b)^T$ ， $x=(x^1,x^2,...,x^n,1)^T$ ,這時邏輯斯蒂模型改寫如下：

$P(Y=1|x)=frac{e^{(wx)}}{1+e^{wx}}$

$P(Y=0|x)=frac{1}{1+e^{wx}}$

我們設 $P(Y=1|x)=pi (x),P(Y=0|x)=1-pi (x)$ ,那麼似然函數為

$prod_{i=1}^{N}[pi(x_i)]^{y_{i}}[1-pi (x_i)]^{1-y_i}$

對數似然函數為

$L(w)=sum_{i=1}^{N}{[y_ilog(x_i)-(1-y_i)log(1-pi(x_i))]}=sum_{i=1}^{N}{[y_i(wx_i)-log(1+exp(wx_i))]}$

對L(w)求極小值，得到w的估計值。

5、使用tensorflow實現邏輯斯蒂回歸演算法

我們使用sklearn裡面的moons數據集來做示例：

from sklearn.datasets import make_moonsm = 1000X_moons, y_moons = make_moons(m, noise=0.1, random_state=42)

這個數據集是個二分類數據集，只有兩個屬性x1,x2，標記為1和0，我們取了1000個樣例，在二維平面上這1000個樣例分布如下圖所示：

將輸入向量X_moons加一列全是1的列進行擴充，

X_moons_with_bias = np.c_[np.ones((m, 1)), X_moons]

X_moons_with_bias數據形式如下

array([[ 1. , -0.05146968, 0.44419863], [ 1. , 1.03201691, -0.41974116], [ 1. , 0.86789186, -0.25482711], [ 1. , 0.288851 , -0.44866862], [ 1. , -0.83343911, 0.53505665]])

將標籤向量重塑，使之每一個標籤與數據實例意義對應

y_moons_column_vector = y_moons.reshape(-1, 1)

劃分訓練集和測試集

test_ratio = 0.2test_size = int(m * test_ratio)X_train = X_moons_with_bias[:-test_size]X_test = X_moons_with_bias[-test_size:]y_train = y_moons_column_vector[:-test_size]y_test = y_moons_column_vector[-test_size:]

在邏輯斯提回歸模型的參數估計中，我們使用梯度下降法來對所要學習的參數進行優化，為了節約資源，一般是使用小批量的數據進行迭代優化，而不是每次優化都將所有訓練數據放進去計算優化，因此我們需要實現一個隨機生成批量訓練數據的函數。

def random_batch(X_train, y_train, batch_size): rnd_indices = np.random.randint(0, len(X_train), batch_size) X_batch = X_train[rnd_indices] y_batch = y_train[rnd_indices] return X_batch, y_batch

上述函數有放回的每次隨機的選擇batch_size個訓練樣本，這樣做可能有些樣本始終都不會被選到，而有些樣本可能經常被選到，也就是訓練集中可能只有三分之二的數據被用於訓練，但是在實踐中這是無關緊要的。

使用tensorflow就必須要建立計算圖

n_inputs = 2X = tf.placeholder(tf.float32, shape=(None, n_inputs + 1), name="X")y = tf.placeholder(tf.float32, shape=(None, 1), name="y")theta = tf.Variable(tf.random_uniform([n_inputs + 1, 1], -1.0, 1.0, seed=42), name="theta")logits = tf.matmul(X, theta, name="logits")y_proba = 1 / (1 + tf.exp(-logits))

上述代碼建立了一個X節點和y節點，theta是隨機初始化的範圍在-1和1之間的三行一列的矩陣，logits表示z， $z=w^Tx+b$ ，而1 / (1 + tf.exp(-logits))就是我們的邏輯斯蒂回歸模型了，當然我們也可以直接使用tensorflow中的專門sigmod函數來改寫這一句代碼：

y_proba = tf.sigmoid(logits)

接下來實現邏輯斯蒂模型的對數似然L(w)，也就是損失函數

epsilon = 1e-7 loss = -tf.reduce_mean(y * tf.log(y_proba + epsilon) + (1 - y) * tf.log(1 - y_proba + epsilon))

epsilon為一個常數因子，為避免計算對數時，所計算的值為0。另外損失函數也可以使用tensorflow中的原生函數

loss = tf.losses.log_loss(y, y_proba)

使用梯度下降優化函數來最小化損失函數

learning_rate = 0.01optimizer = tf.train.GradientDescentOptimizer(learning_rate=learning_rate)training_op = optimizer.minimize(loss)

learning_rate即為學習率，也就是步長，GradientDescentOptimizer為梯度下降優化器。

下面進行訓練

init = tf.global_variables_initializer()n_epochs = 1000batch_size = 50n_batches = int(np.ceil(m / batch_size))with tf.Session() as sess: sess.run(init) for epoch in range(n_epochs): for batch_index in range(n_batches): X_batch, y_batch = random_batch(X_train, y_train, batch_size) sess.run(training_op, feed_dict={X: X_batch, y: y_batch}) loss_val = loss.eval({X: X_test, y: y_test}) if epoch % 100 == 0: print("Epoch:", epoch, " Loss:", loss_val) y_proba_val = y_proba.eval(feed_dict={X: X_test, y: y_test})precision_score(y_test, y_pred)recall_score(y_test, y_pred)

進行1000次迭代，輸出結果

Epoch: 0 Loss: 0.792602Epoch: 100 Loss: 0.343463Epoch: 200 Loss: 0.30754Epoch: 300 Loss: 0.292889Epoch: 400 Loss: 0.285336Epoch: 500 Loss: 0.280478Epoch: 600 Loss: 0.278083Epoch: 700 Loss: 0.276154Epoch: 800 Loss: 0.27552Epoch: 900 Loss: 0.2749120.862745098039215730.88888888888888884

從結果來看效果並不好，準確率和召回率都只有88%左右，將分類結果可視化可以看到，邏輯斯蒂回歸模型僅僅學到了一條二元一次方程直線。

下面嘗試學習一條多元三次方程曲線來建立回歸模型。

X_train_enhanced = np.c_[X_train, np.square(X_train[:, 1]), np.square(X_train[:, 2]), X_train[:, 1] ** 3, X_train[:, 2] ** 3]X_test_enhanced = np.c_[X_test, np.square(X_test[:, 1]), np.square(X_test[:, 2]), X_test[:, 1] ** 3, X_test[:, 2] ** 3]

這裡為每個樣例添加了4個屬性 $x_1^2,x_2^2,x_1^3,x_2^3$ ，這樣就可以建立一個高次方的方程了。

訓練數據形式如下：

array([[ 1.00000000e+00, -5.14696757e-02, 4.44198631e-01, 2.64912752e-03, 1.97312424e-01, -1.36349734e-04, 8.76459084e-02], [ 1.00000000e+00, 1.03201691e+00, -4.19741157e-01, 1.06505890e+00, 1.76182639e-01, 1.09915879e+00, -7.39511049e-02], [ 1.00000000e+00, 8.67891864e-01, -2.54827114e-01, 7.53236288e-01, 6.49368582e-02, 6.53727646e-01, -1.65476722e-02], [ 1.00000000e+00, 2.88850997e-01, -4.48668621e-01, 8.34348982e-02, 2.01303531e-01, 2.41002535e-02, -9.03185778e-02], [ 1.00000000e+00, -8.33439108e-01, 5.35056649e-01, 6.94620746e-01, 2.86285618e-01, -5.78924095e-01, 1.53179024e-01]])

根據上文步驟，我們建立一個函數集成一些操作

def logistic_regression(X, y, initializer=None, seed=42, learning_rate=0.01): n_inputs_including_bias = int(X.get_shape()[1]) with tf.name_scope("logistic_regression"): with tf.name_scope("model"): if initializer is None: initializer = tf.random_uniform([n_inputs_including_bias, 1], -1.0, 1.0, seed=seed) theta = tf.Variable(initializer, name="theta") logits = tf.matmul(X, theta, name="logits") y_proba = tf.sigmoid(logits) with tf.name_scope("train"): loss = tf.losses.log_loss(y, y_proba, scope="loss") optimizer = tf.train.GradientDescentOptimizer(learning_rate=learning_rate) training_op = optimizer.minimize(loss) loss_summary = tf.summary.scalar("log_loss", loss) with tf.name_scope("init"): init = tf.global_variables_initializer() with tf.name_scope("save"): saver = tf.train.Saver() return y_proba, loss, training_op, loss_summary, init, saver

進行訓練

n_inputs = 2 + 4logdir = log_dir("logreg")X = tf.placeholder(tf.float32, shape=(None, n_inputs + 1), name="X")y = tf.placeholder(tf.float32, shape=(None, 1), name="y")y_proba, loss, training_op, loss_summary, init, saver = logistic_regression(X, y)n_epochs = 10001batch_size = 50n_batches = int(np.ceil(m / batch_size))checkpoint_path = "/tmp/my_logreg_model.ckpt"checkpoint_epoch_path = checkpoint_path + ".epoch"final_model_path = "./my_logreg_model"with tf.Session() as sess: if os.path.isfile(checkpoint_epoch_path): # if the checkpoint file exists, restore the model and load the epoch number with open(checkpoint_epoch_path, "rb") as f: start_epoch = int(f.read()) print("Training was interrupted. Continuing at epoch", start_epoch) saver.restore(sess, checkpoint_path) else: start_epoch = 0 sess.run(init) for epoch in range(start_epoch, n_epochs): for batch_index in range(n_batches): X_batch, y_batch = random_batch(X_train_enhanced, y_train, batch_size) sess.run(training_op, feed_dict={X: X_batch, y: y_batch}) loss_val, summary_str = sess.run([loss, loss_summary], feed_dict={X: X_test_enhanced, y: y_test}) file_writer.add_summary(summary_str, epoch) if epoch % 500 == 0: print("Epoch:", epoch, " Loss:", loss_val) saver.save(sess, checkpoint_path) with open(checkpoint_epoch_path, "wb") as f: f.write(b"%d" % (epoch + 1)) saver.save(sess, final_model_path) y_proba_val = y_proba.eval(feed_dict={X: X_test_enhanced, y: y_test}) os.remove(checkpoint_epoch_path)y_pred = (y_proba_val >= 0.5)print("precision_score",precision_score(y_test, y_pred))print("recall_score",recall_score(y_test, y_pred))

這裡有一些代碼是對模型進行保存，不詳細說了。最終準確率和找回率輸出結果如下

0.979797979797979780.97979797979797978

可以看到效果提升了非常多。最後將分類結果進行可視化輸出

參考資料：

1、周志華《機器學習》
2、李航《統計學習方法》
3、Aurélien Géron《Hands-On Machine Learning with Scikit-Learn and TensorFlow》

未經允許，禁止轉載！