分类

为什么不能直接用线性回归？-多分类不适用，penalize samples that are too correct 目标是接近1，但可能远大于1，影响了最优分类面的确定

概率生成式模型

分类问题核心是对于输入$x$，找到使得$p(C|x)$最大的类别$C$

按照全概率公式，上面的计算可以变成$p(C|x)=\frac{p(x|C)p(C)}{p(x|C_1)p(C_1)+...+p(x|C_n)p(C_n)}$

​为什么被称为生成式模型，因为为了计算$p(C|x)$，已经得到了全部的先验和似然度，模型描述了每一个$x$ 的概率。如何计算$p(C_i)$ 和$p(x|C_i)$? 前者很简单，只需要根据数据的种类统计一下占比即可。后者的计算一般是假设每个类别的数据服从某个分布（比如Guassian），这样如果有$|C|$类就有$|C|$个类别函数。再使用极大似然估计利用该类别的数据计算分布的参数（如果是Guassian，他的极大似然刚好就是$\mu=\frac{1}{n}\sum x$，$\Sigma=\frac{1}{n}\sum(x-\mu)(x-\mu)^T$）

如果假设所有的feature都是独立的，那么分布就变成连乘，Naive Bayes Method

判别式模型：逻辑回归

对于二分类问题，$P(C_1|x)=\frac{P(x|C_1)P(C_1)}{P(x|C_1)P(C_1) + P(x|C_2)P(C_2)}=\frac{1}{1 + \frac{P(x|C_2)P(C_2)}{P(x|C_1)P(C_1)}}=\frac{1}{1+e^{-z}}$。 定义sigmoid函数$\sigma(z)=\frac{1}{1+e^{-x}}$，那么分类的概率就是$\sigma(z)$，其中$z = \ln \frac{P(x|C_1)P(C_1)}{P(x|C_2)P(C_2)}$。

接下来分析一下这个$z$

所以重点就是这个z，如果直接把两个概率看作一个整体训练出来，那相当于直接学习了$P(C|x)$，这就是判别式模型。

如果假设$z$的各个属性服从联合高斯分布，并且不同的类别协方差一样（这是一个为了简化计算的假设），那么有

经过化简，二次项没有了，只剩下一次和0次项，因此可以看作是一个线性模型，所以直接用$w^Tx+b$拟合。 所以$P(C_1|x)=\sigma(w^Tx+b)$，这就是逻辑回归。

逻辑回归只能求解二分类问题，这是因为他用了线性回归决定的。逻辑二字也说明他只能0-1。如果多分类就只能1对1或者1对1。

loss function $L = \frac{1}{n} \sum_i (y_i \log(\sigma(w'x_i + b)) + (1-y_i)\log(1 - \sigma(w'x_i + b)))$。很直观，sigma表示输出分类1的概率，所以交叉熵就是这样的。

求梯度，还是表示成增广的形式，[b, w]，数据每条前面增加1。 这样loss变成了 $L = -\frac{1}{n}\sum_i(y_i\log(\sigma(w^Tx_i)) + (1 - y_i)\log(1 - \sigma(w^Tx_i)))$

求导（这里有个小技巧，$[\sigma(z)]'=\sigma(z)(1-\sigma(z))$），得到 $\frac{\partial L}{\partial w_k}=-\frac{1}{n}\sum_i (y_i (1 - \sigma(w^Tx_i))x_{ik}- (1 - y_i)\sigma(w^Tx_i)x_{ik}) = -\frac{1}{n}\sum_ic_i x_{ik}=\frac{1}{n} [c_1, c_2, ..., c_n]\begin{bmatrix}x_{1k}\\x_{2k}\\ ...\\x_{nk}\end{bmatrix}$ ，

其中 $c_i = \sigma(w^Tx_i) - y_i$

现在回到全体参数 $\frac{\partial L}{\partial w^T}=\frac{1}{n}[c_1, c_2, ..., c_n]\begin{bmatrix} 1 & x_{11} & x_{12} & ... & x_{1k} \\ 1 & x_{21} & x_{22} & ... & x_{2k} \\ ... & ... & ... &... & ... \\ 1 & x_{n1} & ... &... & x_{nk}\end{bmatrix}$

虽然两个模型set一样，但即使用同样的数据也会训练出不同的模型，因为生成式模型假设了先验的分布，一般来说，在数据量少的时候生成式模型会比较有优势。

对于线性不可分问题，可以使用级联多个分类器的方法，实现feature transformation，这就是神经网络了

softmax的推导过程

思考一个多分类问题，他的输出应该是什么样子，显然直接用0，1，2，3标号是不合适的，因为这会让模型误以为1和3的差距要比1和2大，所以并不能只输出一个数字，而应该输出一个向量，有多少类，向量就有多长，模型输出对每个类的分数，分数最高的类就是模型的分类。

在明确输出以后，在思考输出和标签之间的关系，如何构造loss function。一个直观的方法是$L=\max\{\sum_{i\ne y} (z_i - z_y), 0\}$ 其中y是正确的分类在向量中的位置，这个loss的意思是让所有不正确分类的分数都尽可能小。而如果分类正确了，也不能让其他位置的值无限制的降低，所以只要$z_i-z_y$小于0了loss就变成0不再鼓励优化了。这样的一个小问题是梯度可能很大，因为求导的时候所有的$z_i$都会有一个1的梯度，为了避免梯度太大的问题，重新写成$L=\max\{\max_{i\ne y}z_i - z_y, 0\}$ 这样含义变成了$z_y$比其他最大的更大，意思没变，好处是可以避免梯度过大，因为只会有两项的偏导是1或-1。坏处由两个：1 只产生两个梯度学习速度太慢，2 max函数不能求导

解决上述两个问题的思路是，利用平滑的函数近似max函数，从而使得他可以求导，并且由于平滑的过程需要使用所有$z_i$表示，但是梯度的大小仍然可控。接下来的问题是如何找到max函数的smooth近似

首先考虑两个任意大小变量$x, y$。那么$\max\{x,y\}=\frac{1}{2}(x+y+|x-y|)$ 。问题变成了如何找绝对值函数的近似。核心思想是迂回一下，先找到绝对值函数的导数的近似，再求不定积分就得到了绝对值函数的近似。绝对值函数的导数除了0位置不可导，其他地方是$f(x)=\begin{cases} 1, x > 0 \\-1, x < 0 \end{cases}$，这可以从另一个函数$\delta(x)=\begin{cases}1, x > 0 \\ 0, x < 0 \end{cases}$ 转化过来，$f(x)=2\delta(x)-1$

而$\delta(x)$物理学早已给出一个相当好的近似$a(x)=\lim_{k\rightarrow\infty} \frac{1}{1+e^{-kx}}$ 。$k$这个参数其实控制了平滑的程度，k越大越陡峭。所以$f(x)$就可以写成$f(x)=\lim_{k\rightarrow\infty}\frac{2}{1+e^{-kx}}-1$ 。通分凑分母分子齐次可以转化成$f(x)=\lim_{k\rightarrow\infty}\frac{-2e^{-kx}}{1+e^{-kx}} + 1$ ，配齐次是为了方便后面积分。得到$f(x)$对其积分就得到了绝对值函数的近似$F(x)=\lim_{k\rightarrow\infty} \frac{2}{k}\ln(1+e^{-kx}) + x$。统一形式后面的$x$变成$\frac{1}{k}\ln e^{kx}$，和前面的ln合并一起乘，注意前一项系数的2会变成平方，最后变成$F(x)=\lim_{k\rightarrow\infty} \frac{1}{k}\ln(e^{-kx}+2+e^{kx})$ 。在$k\rightarrow \infty$ 过程中，ln里面的+2和都可以忽略掉，所以$F(x)=\lim_{k\rightarrow\infty}\frac{1}{k}\ln(e^{kx} + e^{-kx})$。

在得到了绝对值函数的smooth形式以后，回到max函数，可以得到

$\max(x,y)=\frac{1}{2}(x+y+|x-y|)=\frac{1}{2}(\frac{1}{k}\ln e^{kx} + \frac{1}{k}\ln e^{ky} +\frac{1}{k}\ln(e^{k(x-y)} + e^{-k(x-y)})) = \frac{1}{2k}\ln(e^{2kx}+e^{2ky})$

再考虑前面的极限，其实一般情况下是取k是一个常数就足够了，所以这里取k=1/2，就得到了max函数的一个平滑近似 $\max(x,y)=\ln(e^x + e^y)$ 。

更进一步，可以推广到多元的最大值，因为$\max(x,y,z)=\max(\max(x,y), z)=\ln(e^{\ln(e^x+e^y)}+e^z)=\ln(e^x+e^y+e^z)$ 。这个函数又被称为logSumExp，是max函数

现在回到前面的loss function，$L=\max\{\max_{i\ne y}z_i - z_y, 0\}=\ln(e^{\ln(\sum_{i\ne y}e^{z_i}) - z_y}+1)=-\ln(\frac{e^{z_y}}{\sum_i e^{z_i}})$。为了表示y的取值的多样性，所以需要变成$L=-\sum_{c}\Pr[c=y]\ln(\frac{e^{z_y}}{\sum_i e^{z_i}})$

现在结合交叉熵来看一下$H(p,q)=-\sum_x p(x)\log q(x)$。它衡量了两个分布p和q的相似程度。现在我们求出来的这个loss，Pr其实就是正确标签的分布，他是一个独热码，后面的ln就可以看作是我们输出的概率分布。这就是为什么softmax函数可以看作输出的是一种概率，而softmax的输入$z$ 则可以称之为logit，这其实是log-it的意思，对概率求log，就得到了z。这也是为什么多分类的最后一层要在z的基础上sofmax一下。

因此我们发现，这个loss function的本质其实是从优化参数之间的大小关系出发，通过把max函数进行平滑近似，得到的结果恰好就是交叉熵。即希望把输出的这个分布和预期的one-hot尽可能相同，妙啊。

再重新结合二分类问题看，对比softmax和sigmoid能看出一些问题，如果二分类的输出不再是逻辑回归的单一值，同样也输出一个长度维2的向量，那么softmax的输出就是$[\frac{e^{z_1}}{e^{z_1}+e^{z_2}}, \frac{e^{z_2}}{e^{z_1}+e^{z_2}}] = [\frac{1}{1+e^{-(z_1-z_2)}}, \frac{1}{1+e^{-(z_2-z_1)}}]=[\frac{1}{1+e^{-(z_1-z_2)}}, 1 - \frac{1}{1+e^{-(z_1-z_2)}}]$。所以可以看出，二分类如果使用逻辑回归方法，其实他是把两个类别的feature的差值作为logit。而回顾一下最大似然估计，进一步证明了这个位置是对概率求了log。

生成模型与判别模型代码实现

from types import new_class
import numpy as np


def gen_data_linear(n_entries=100):
    w = np.array([1, 2, 3, 4, 5, 6]).reshape(6, 1)
    b = 1
    X = (np.random.randn(w.shape[0] * n_entries) * 5).reshape(n_entries, w.shape[0])
    y = np.matmul(X, w) + b + np.random.randn(n_entries).reshape(n_entries, 1)
    return X, y

def gen_class_data(n_entries=100, n_classes=2):
    X, y = gen_data_linear(n_entries)
    y_logits = np.zeros(shape=(n_entries, n_classes))
    max_y, min_y = np.max(y), np.min(y)
    seg = (max_y - min_y) / n_classes
    for i in range(n_entries):
        c = max(0, min(int((y[i][0] - min_y) / seg), n_classes - 1))
        y_logits[i][c] = 1
    return X, y_logits


def sigmoid(z):
    # sigma(z) = 1 / (1 + e^(-z)), z is a vector [dim * 1] or scalar
    n = z.shape[0]
    z.reshape(n)
    return 1. / (1. + np.exp(-z))


def softmax(z):
    # softmax(z) = e^z{z_j} / sum{e^z_i}, z: [dim * 1]
    return (np.exp(z) / np.sum(np.exp(z))).reshape(n, 1)

class Probalistic_Generative_Model:

    def __init__(self):
        self.prior = None
        self.shared_var = None
        self.mus = None
        self.n_features = 0
        self.n_classes = 0
        self.guassian_coefficient = 1 / np.sqrt(2 * np.pi)

    def fit(self, data, labels):
        n_samples, n_features = data.shape[0], data.shape[1]
        n_classes = labels.shape[1]
        self.n_features, self.n_classes = n_features, n_classes

        class_cnt = np.sum(labels, axis=0)
        self.prior = class_cnt / n_samples
        # assuming each property is Guassian distribution and share the same covariance, and independent
        mus = np.zeros(shape=(n_classes, n_features))
        for i in range(n_samples):
            c = np.argmax(labels[i])
            mus[c] += data[i]
        for c in range(n_classes):
            mus[c] /= class_cnt[c]
        self.mus = mus

        var = np.zeros(shape=(n_classes, n_features))
        for i in range(n_samples):
            c = np.argmax(labels[i])
            var[c] += (data[i] - mus[c]) * (data[i] - mus[c])
        for c in range(n_classes):
            var[c] /= class_cnt[c]
        shared_var = np.zeros(n_features)
        for c in range(n_classes):
            shared_var += class_cnt[c] / n_samples * var[c]
        self.shared_var = shared_var

    def predict(self, samples):
        prediction = np.zeros(shape=(samples.shape[0], self.n_classes))
        for i in range(samples.shape[0]):
            posterior_numerator = np.ones(self.n_classes)
            for c in range(self.n_classes):
                for k in range(self.n_features):
                    posterior_numerator[c] *= (1. / np.sqrt(self.shared_var[k]) * np.exp(-np.square(samples[i][k] - self.mus[c][k]) / (2 * self.shared_var[k])))
                posterior_numerator[c] *= self.prior[c]
            # need not cauculate posterior
            prediction[i][np.argmax(posterior_numerator)] = 1
        return prediction


class Logistic_Regression:

    def __init__(self):
        self.w = None
        self.b = None

    def fit(self, data, labels, n_iterations, learning_rate):
        n_features = data.shape[1]
        shrink_labels = np.apply_along_axis(np.argmax, axis=1, arr=labels).reshape(data.shape[0], 1)
        self.w = np.random.randn(n_features).reshape(n_features, 1)
        self.b = np.random.randn(1).reshape(1, 1)
        for _ in range(n_iterations):
            self.w, self.b = self._gradient_descent(data, shrink_labels, self.w, self.b, learning_rate)

    def _gradient_descent(self, data, labels, current_weights, current_bias, learning_rate):
        # let w' = [b, w^T]^T, D' = [1, data]^T, x'_i = [1, sample^T]^T
        # loss function: 1 / n sum_i y_i * log(sigma(w'^T * x'_i)) + (1 - y_i) * log(1 - sigma(w'^T * x'_i)) 
        # gradient_k = c_i * x'_ik
        # gradient for single sample gradient = c_i * x'_i, where c_i =  sigma(w'^T * x'_i) - y_i
        # for mini-batch or full data, gradient = D' * C, where C = [c_1, c_2, ..., c_n]^T
        n = data.shape[0]
        D_ = np.concatenate((np.ones(n).reshape(n, 1), data), axis=1)  # [n, k + 1]
        w_ = np.concatenate((current_bias, current_weights))
        predict = sigmoid(np.matmul(D_, w_))
        C = predict - labels
        gradient = np.matmul(D_.T, C) / n
        new_weights = current_weights - learning_rate * gradient[1:] 
        new_bias = current_bias - learning_rate * gradient[0]
        return new_weights, new_bias

    def predict(self, samples):
        n_samples = samples.shape[0]
        prediction = np.zeros(shape=(n_samples, 2))
        probability = sigmoid(np.matmul(samples, self.w) + self.b)
        for i in range(n_samples):
            c = 1 if probability[i][0] > 0.5 else 0
            prediction[i][c] = 1
        return prediction


X, y = gen_class_data(2000, 2)
train_size = 1500
X_train = X[0: train_size]
y_train = y[0: train_size]
X_test = X[train_size:]
y_test = y[train_size:]

# print(X_train)
# print(y_train)

learning_rate = 0.01
n_iterations = 5000

lr = Logistic_Regression()
lr.fit(X_train, y_train, n_iterations, learning_rate)
y_predict = lr.predict(X_test)

# pgm = Probalistic_Generative_Model()
# pgm.fit(X_train, y_train)
# y_predict = pgm.predict(X_test)