广义线性模型

 指数分布族

前面学习了线性回归和logistic回归。我们知道对于(P(y|x; heta))

若y属于实数，满足高斯分布，得到基于最小二乘法的线性回归,若y取{0,1}，满足伯努利分布，得到Logistic回归。

这两个算法，其实都是广义线性模型的特例。

1) 伯努利分布

http://www.cnblogs.com/mikewolf2002/p/7667944.html

2) 高斯分布

http://www.cnblogs.com/mikewolf2002/p/7669977.html

指数分布族的定义：

     若一类概率分布可以写成如下形式，那么它就属于指数分布族：

[P(y;eta)=b(y)exp(eta^{T}T(y)-a(eta))]

     (eta)称作分布的自然参数，通常是一个实数。(T(y))称作充分统计量，通常(T(y)=y)，实际上是一个概率分布的充分统计量（统计学知识）。对于给定的(a，b，T)三个函数，上式定义了一个以(eta)为参数的概率分布集合，即改变(eta)可以得到不同的概率分布。

https://en.wikipedia.org/w/index.php?title=Exponential_family&oldid=648989632

证明伯努利分布是指数分布族：

伯努利分布式对0,1问题进行建模，它的分布律如下：

[Ber(phi ),p(y=1; varphi ) = varphi; p(y=0; varphi ) = 1- varphi]

伯努利分布是对0，1问题进行建模的，对于(Bernoulli（varphi）,yin {0,1})有(p(y=1;varphi)=varphi;p(y=0;varphi)=1−varphi)，下面将其推导成指数分布族形式：

egin{align*} L( heta) &=P(y;varphi)=varphi^y(1-varphi)^{1-y}=exp(log(varphi^y(1-varphi)^{1-y})) \ 
&=exp(ylog(varphi)+(1-y)log(1-varphi)) \ &= exp(ylog(frac{varphi}{1-varphi})+log(1-varphi))end{align*}

    将其与指数族分布形式对比，可以看出：

[b(y) = 1]

[T(y) = y]

[eta=logfrac{phi}{1-phi} => phi=frac{1}{1+e^{-eta}}]

[a(eta)=-log(1-phi)=log(1+e^{eta})]

    表明伯努利分布也是指数分布族的一种。从上述式子可以看到，(eta)的形式与logistic函数（sigmoid）一致，这是因为 logistic模型对问题的前置概率估计其实就是伯努利分布。

高斯分布

    下面对高斯分布进行推导，推导公式如下（为了方便计算，我们将方差(sigma^2)设置为1）：

[N(mu,1)=frac{1}{sqrt{2pi}}exp^{-frac{1}{2}(y-mu)^2}=frac{1}{sqrt{2pi}}exp^{-frac{1}{2}y^2+ymu-frac{1}{2}mu^2}=frac{1}{sqrt{2pi}}exp^{-frac{1}{2}y^2}exp^{ymu-frac{1}{2}mu^2}]

将上式与指数族分布形式比对，可知：

[b(y)=frac{1}{sqrt{2pi}}exp^{-frac{1}{2}y^2}]

[T(y)=y]

[eta=mu]

[a(eta)=frac{mu^2}{2}]

两个典型的指数分布族，伯努利和高斯分布。其实大多数概率分布都可以表示成指数分布族形式，如下所示：

• 伯努利分布（Bernoulli）：对 0、1 问题进行建模；
• 多项式分布（Multinomial）：对 K 个离散结果的事件建模；
• 泊松分布（Poisson）：对计数过程进行建模，比如网站访问量的计数问题，放射性衰变的数目，商店顾客数量等问题；
• 伽马分布（gamma）与指数分布（exponential）：对有间隔的正数进行建模，比如公交车的到站时间问题；
• β 分布：对小数建模；
• Dirichlet 分布：对概率分布进建模；
• Wishart 分布：协方差矩阵的分布；

广义线性模型GLM

选定了一个指数分布族后，怎样来推导出一个GLM呢？

假设：

（1）(y|x; heta~ExponentialFamily（eta）)，即假设试图预测的变量(y)在给定(x)，以( heta)作为参数的条件概率，属于以(eta)作为自然参数的指数分布族，例：若要统计网站点击量(y)，用泊松分布建模。

（2） 给定(x)，目标是求出以(x)为条件的(T(y))的期望(E[T(y)|x])，即让学习算法输出(h(x) = E[T(y)|x])

（3）(eta= heta^Tx),即自然参数和输入特征(x)之间线性相关，关系由θ决定。仅当(eta)是实数时才有意义,若(eta)是一个向量，则(eta_i= heta_i^Tx)

推导伯努利分布的GLM：

(y|x; heta~ExponentialFamily（eta）)，伯努利分布属于指数分布族,对给定的(x, heta)，学习算法进行一次预测的输出：

[h_ heta(x)=E[y|x; heta]=phi=frac{1}{1+e^{-eta}}=frac{1}{(1+e^{- heta^Tx})}]

得到logistic回归算法。

正则响应函数：(g(eta)=E[y;eta]),将自然参数(eta)和(y)的期望联系起来

正则关联函数:(g^{-1})

推导多项式分布的GLM：

多项式分布是在(k)个可能取值上的分布，即(yin {1,..,k}),如将收到的邮件分成(k)类，诊断某病人为(k)种病中的一种等问题。

（1）将多项式分布写成指数分布族的形式：

设多项式分布的参数：(phi_1,phi_2,…,phi_k) ，且(sumlimits_{i=1}^kphi_i=1) ，(phi_i)表示第(i)个取值的概率分布，最后一个参数可以由前(k-1)个推导出，所以只将前(k-1)个(phi_1,phi_2,…,phi_{k-1})视为参数。

多项式分布是少数几个(T(y)!=y)的分布，(T(1) sim T(k))都定义成一个(k-1)维的列向量，表示为：

[T(1)=egin{bmatrix}
1\
0\
0\
vdots\
0end{bmatrix}T(2)=egin{bmatrix}
0\
1\
0\
vdots\
0end{bmatrix}T(3)=egin{bmatrix}
0\
0\
1\
vdots\
0end{bmatrix}T(k-1)=egin{bmatrix}
0\
0\
0\
vdots\
1end{bmatrix}T(k)=egin{bmatrix}
0\
0\
0\
vdots\
0end{bmatrix}]

这样定义(T(y))是为了将多项式分布写成指数分布族形式。

*定义符号：指示函数，(1{.})

1{True} = 1, 1{False} = 0，即大括号内命题为真，值为1,否则为0。

例：1{2=3} = 0, 1{1+1=2} = 1

用(T(y)_i)表示(T(y))的第(i)个元素，则(T(y)_i=1{y=i})

根据参数(phi)的意义（(phi_i)表示第(i)个取值的概率分布），可推出：

egin{align*} p(y;phi) &=phi_1^{1{y=1}}phi_2^{1{y=2}}...phi_k^{1{y=k}} \ 
&=phi_1^{1{y=1}}phi_2^{1{y=2}}...phi_k^{1-sumlimits_{i=1}^{k-1}1{y=i}} \ &=phi_1^{(T(y))_1}phi_2^{(T(y))_2}...phi_k^{1-sumlimits_{i=1}^{k-1}(T(y))_i} \ &=>exp(T(y))_1log(phi_1)+ exp(T(y))_2log(phi_2)+...+(1-sumlimits_{i=1}^{k-1}(T(y))_i)log(phi_k) \ &=exp(T(y))_1log(phi_1/phi_k)+ exp(T(y))_2log(phi_2/phi_k)+...+exp(T(y))_{k-1}log(phi_{k-1}/phi_k)+log(phi_k)) \ &=b(y)exp(eta^{T}T(y)-a(eta))end{align*}

可得：

[eta=egin{bmatrix}
log(phi_1/phi_k)\
log(phi_2/phi_k)\
vdots\
log(phi_{k-1}/phi_k)
end{bmatrix}]

[a(eta)=-log(phi_k)]

[b(y)=1]

证明多项式分布是指数分布族。

再用(eta)表示(phi)：

[eta_i=logfrac{phi_i}{phi_k} \
e^{eta_i}=frac{phi_i}{phi_k} \
phi_ke^{eta_i}=phi_i \
phi_ksumlimits_{i=1}^ke^{eta_i}=sumlimits_{i=1}^kphi_i=1 \
phi_i=frac{e^{eta_i}}{sumlimits_{j=1}^ke^{eta_j}}]

（2）根据上述假设（3）中自然参数和输入(x)的线性关系，可求得：

egin{align*} p(y=i|x; heta) &=phi_i \ 
  &=frac{e^{eta_i}}{sumlimits_{j=1}^ke^{eta_j}} \ &=frac{e^{ heta_i^Tx}}{sumlimits_{j=1}^ke^{ heta_j^Tx}} \ end{align*}

（3）根据上述假设（2）中的输出(h(x) = E[T(y)|x])，可求得：

egin{align*} h_ heta(x) &=E[T(y)|x; heta] \ 
  &=egin{bmatrix}
1{y=1}|x; heta\
1{y=2}|x; heta\
vdots\
1{y=k-1}|x; heta
end{bmatrix} \ &= egin{bmatrix}
phi_1\
phi_2\
vdots\
phi_{k-1}
end{bmatrix}\ &=egin{bmatrix}
frac{e^{ heta_1^Tx}}{sumlimits_{j=1}^ke^{ heta_j^Tx}} \
frac{e^{ heta_2^Tx}}{sumlimits_{j=1}^ke^{ heta_j^Tx}} \
vdots\
frac{e^{ heta_{k-1}^Tx}}{sumlimits_{j=1}^ke^{ heta_j^Tx}}
end{bmatrix} end{align*}

称这种回归算法为softmax回归，是logistic回归的推广。

Softmax回归的训练方法和logistic相似，通过极大似然估计找到参数( heta),其对数似然性为：

egin{align*} l( heta) &=sumlimits_{i=1}^{m}logp(y^{(i)}|x^{(i)}; heta); \ 
  &=sumlimits_{i=1}^{m}logprodlimits_{l=1}^kigg(frac{e^{ heta_l^Tx^{(i)}}}{sumlimits_{j=1}^ke^{ heta_j^Tx^{(i)}}}igg)^{1{y^{(i)}=l}} = sumlimits_{i=1}^migg(sumlimits_{l=1}^kigg( 1{y^{(i)}=l}logigg( frac{e^{ heta_l^Tx^{(i)}}}{sumlimits_{j=1}^ke^{ heta_j^Tx^{(i)}} }igg)igg)igg)end{align*}

其中( heta)是一个(k imes(n+1))的矩阵，代表这(k)个类的所有训练参数，每个类的参数是一个(n+1)维的向量。

跟Logistic回归一样，softmax也可以用梯度下降法或者牛顿迭代法求解，但在softmax回归中直接用上述对数似然函数求导是不能更新参数的，因为它存在冗余的参数，通常用牛顿方法中的Hessian矩阵也不可逆，是一个非凸函数，那么可以通过添加一个权重衰减项来修改代价函数，使得代价函数是凸函数，并且得到的Hessian矩阵可逆

再通过梯度上升或牛顿方法找对数似然性的最大值，即可求出参数( heta)。

详细的公式推导见：http://ufldl.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92