广义线性模型

普通线性模型：

假设主要有以下几点：

1.响应变量Y和误差正态性：响应变量Y和误差项服从正态分布，并且误差是一个白噪声过程，因而具有零均值，同方差的特性。

2.预测量X和未知参数的非随机性：预测量X具有非随机性，可测不存在测量误差；未知参数认为是未知但不具有随机性的常数。

3.主要研究响应变量的均值E[Y]

4.连接方式：响应变量E[Y]与预测变量的线性组合。f(x)=x

广义线性模型是在普通线性模型的基础上，将上述四点模型假设进行推广而得出的应用范围更广，更实用的回归模型。

1.响应变量的分布推广至指数分布族：比如正态分布，泊松分布，二项分布，负二项分布，伽马分布，逆高斯分布

2.预测变量X和未知参数的非随机性

3.研究对象：广义线性模型的主要研究对象任然是响应变量的均值E[Y]

4.连接方式;理论上可以是任意的，而不再局限于f(x)=x.

因变量y不再只是正态分布，而是扩大为指数族中的任一分布。

解释变量x的线性组合不再直接用于解释因变量y的均值u，而是通过一个联系函数g来解释g(u).

线性回归是一种parametric learning algorithm,而局部加权线性回归是一种non-parametric learning algorithm。Parametric learning algorithm有固定的(指的是：值的大小是固定的)、有限的参数、通过训练样本，找到合适的参数后，对于之后未知的输入，我们可以直接利用这组参数得出其相应的预测输出。而non-parametric learning algorithm需要的计算量与输入的训练集大小成正比，对于每次新的输入，需要计算相应的参数后，才能求出相应的预测输出。

局部加权回归的缺点是：需要拿着训练数据去做预测。

当训练样本的features大于样本数时，那么矩阵X就不是满秩的。无法计算X.T*X。为了解决这个问题，引入了ridge regression的概念。ridge regession加入了lambda*I到矩阵X.T*X

当两个或更多的特征有关联时，如果用正规的最小二乘法回归，可能会有一个非常大的正的权值和一个负的非常大的权值。

自适应学习率的学习：

基于Armijo准则计算搜索方向上的最大步长，其基本思想是沿着搜索方向移动一个比较大的步长估计值，然后以迭代形式不断缩减步长，直到该步长使得函数值f(xk+a*dk)相对于当前函数值f(xk)的减小程度大于预设的期望值为止。

# -*- coding:utf8 -*-

import math
import matplotlib.pyplot as plt


def f(w, x):
    N = len(w)
    i = 0
    y = 0
    while i < N - 1:
        y += w[i] * x[i]
        i += 1
    y += w[N - 1]  # 常数项
    return y


def gradient(data, w, j):
    M = len(data)   # 样本数
    N = len(data[0])
    i = 0
    g = 0   # 当前维度的梯度
    while i < M:
        y = f(w, data[i])
        if (j != N - 1):
            g += (data[i][N - 1] - y) * data[i][j]
        else:
            g += data[i][N - 1] - y
        i += 1
    return g / M


def gradientStochastic(data, w, j):
    N = len(data)   # 维度
    y = data[N - 1] - f(w, data)
    if (j != N - 1):
        return y * data[j]
    return y    # 常数项


def isSame(a, b):
    n = len(a)
    i = 0
    while i < n:
        if abs(a[i] - b[i]) > 0.01:
            return False
        i += 1
    return True


def fw(w, data):
    M = len(data)   # 样本数
    N = len(data[0])
    i = 0
    s = 0
    while i < M:
        y = data[i][N - 1] - f(w, data[i])
        s += y ** 2
        i += 1
    return s / 2


def fwStochastic(w, data):
    y = data[len(data) - 1] - f(w, data)
    y **= 2
    return y / 2


def numberProduct(n, vec, w):
    N = len(vec)
    i = 0
    while i < N:
        w[i] += vec[i] * n
        i += 1


def assign(a):
    L = []
    for x in a:
        L.append(x)
    return L


# a = b
def assign2(a, b):
    i = 0
    while i < len(a):
        a[i] = b[i]
        i += 1


def dotProduct(a, b):
    N = len(a)
    i = 0
    dp = 0
    while i < N:
        dp += a[i] * b[i]
        i += 1
    return dp


# w当前值；g当前梯度方向；a当前学习率；data数据
def calcAlpha(w, g, a, data):
    c1 = 0.3
    now = fw(w, data)
    wNext = assign(w)
    numberProduct(a, g, wNext)
    next = fw(wNext, data)
    # 寻找足够大的a，使得h(a)>0
    count = 30
    while next < now:
        a *= 2
        wNext = assign(w)
        numberProduct(a, g, wNext)
        next = fw(wNext, data)
        count -= 1
        if count == 0:
            break

    # 寻找合适的学习率a
    count = 50
    while next > now - c1*a*dotProduct(g, g):
        a /= 2
        wNext = assign(w)
        numberProduct(a, g, wNext)
        next = fw(wNext, data)

        count -= 1
        if count == 0:
            break
    return a

# w当前值；g当前梯度方向；a当前学习率；data数据
def calcAlphaStochastic(w, g, a, data):
    c1 = 0.01   # 因为是每个样本都下降，所以参数运行度大些，即：激进一些
    now = fwStochastic(w, data)
    wNext = assign(w)
    numberProduct(a, g, wNext)
    next = fwStochastic(wNext, data)
    # 寻找足够大的a，使得h(a)>0
    count = 30
    while next < now:
        if a < 1e-10:
            a = 0.01
        else:
            a *= 2
        wNext = assign(w)
        numberProduct(a, g, wNext)
        next = fwStochastic(wNext, data)
        count -= 1
        if count == 0:
            break

    # 寻找合适的学习率a
    count = 50
    while next > now - c1*a*dotProduct(g, g):
        a /= 2
        wNext = assign(w)
        numberProduct(a, g, wNext)
        next = fwStochastic(wNext, data)

        count -= 1
        if count == 0:
            break
    return a


def normalize(g):
    s = 0
    for x in g:
        s += x * x
    s = math.sqrt(s)
    i = 0
    N = len(g)
    while i < N:
        g[i] /= s
        i += 1


def calcCoefficient(data, listA, listW, listLostFunction):
    M = len(data)       # 样本数目
    N = len(data[0])    # 维度
    w = [0 for i in range(N)]
    wNew = [0 for i in range(N)]
    g = [0 for i in range(N)]

    times = 0
    alpha = 100.0  # 学习率随意初始化
    same = False
    while times < 10000:
        i = 0
        while i < M:
            j = 0
            while j < N:
                g[j] = gradientStochastic(data[i], w, j)
                j += 1
            normalize(g)  # 正则化梯度
            alpha = calcAlphaStochastic(w, g, alpha, data[i])
            #alpha = 0.01
            numberProduct(alpha, g, wNew)

            print "times,i, alpha,fw,w,g:\t", times, i, alpha, fw(w, data), w, g
            if isSame(w, wNew):
                if times > 5:  #防止训练次数过少
                    same = True
                    break
            assign2(w, wNew)  # 更新权值

            listA.append(alpha)
            listW.append(assign(w))
            listLostFunction.append(fw(w, data))

            i += 1
        if same:
            break
        times += 1
    return w


if __name__ == "__main__":
    fileData = open("d8.txt")
    data = []
    for line in fileData:
        d = map(float, line.split('\t'))
        data.append(d)
    fileData.close()

    listA = []  # 每一步的学习率
    listW = []  # 每一步的权值
    listLostFunction = []  # 每一步的损失函数值
    w = calcCoefficient(data, listA, listW, listLostFunction)

    # 绘制学习率
    plt.plot(listA, 'r-', linewidth=2)
    plt.plot(listA, 'go')
    plt.xlabel('Times')
    plt.ylabel('Ratio/Step')
    plt.show()

    # 绘制损失
    listLostFunction = listLostFunction[0:100]
    listLostFunction[0] /= 2
    plt.plot(listLostFunction, 'r-', linewidth=2)
    plt.plot(listLostFunction, 'gv', alpha = 0.75)
    plt.xlabel('Times')
    plt.ylabel('Loss Value')
    plt.grid(True)
    plt.show()

    # 绘制权值
    X = []
    Y = []
    for d in data:
        X.append(d[0])
        Y.append(d[1])
    plt.plot(X, Y, 'cp', label=u'Original Data', alpha=0.75)
    x = [min(X), max(X)]
    y = [w[0] * x[0] + w[1], w[0] * x[1] + w[1]]
    plt.plot(x, y, 'r-', linewidth=3, label='Regression Curve')
    plt.legend(loc='upper left')
    plt.grid(True)
    plt.show()

以上代码为：分别用了梯度下降算法和随机梯度下降算法求线性回归的方程。学习率使用Armijo准则计算。

采用多项式插值法拟合简单函数，然后根据简单函数估计函数的极值点，这样选择合适的步长的效率会高很多。

现在拥有的数据为：Xk的函数值f(Xk)及其导数f'(Xk),再加上第一次尝试的步长a0。如果a0满足条件，显然算法退出；若a0不满足条件，则可以根据上述信息可以构造一个二次函数：

显然，导数为0处的最优值为：

若a1满足Armijo准则，则输出该学习率；否则继续迭代。

#include "stdafx.h"
#include <math.h>

double f(double x) {
	return x * x;
}

double g(double x) {//一阶导数
	return 2 * x;
}

double g2(double x) {//二阶导数
	return 2;
}

//x为当前值 d为x处的导数  a为输入的学习率
//使用Armijo求学习率
double GetA_Armijo(double x, double d, double a) {
	double c1 = 0.3;
	double now = f(x);
	double next = f(x-a*d);

	int count = 30;
	while (next < now) {
		a *= 2;
		next = f(x-a*d);

		count --;

		if (count == 0)
			break;
	}

	count = 50;

	while(next > now - c1*a*d*d) {
		a /= 2;
		next = f(x-a*d);
		count --;

		if(count == 0) {
			break;
		}
	}

	return a;
}
//使用二次插值法求学习率
// x为当前值 d为梯度的方向 a为前一步的学习率
double GetA_Quad(double x, double d, double a) {
	double c1 = 0.03;
	double now = f(x);
	double next = f(x-a*d);

	int count = 30;
	while (next < now) {
		a *= 2;
		next = f(x-a*d);

		count -= 1;
		if (count == 0) {
			break;
		}
	}

	while (next < now-c1*a*d*d) {
		a = d*a*a/(now + d*a-next);//二次插值求学习率

	}
}
int _tmain(int argc, _TCHAR* argv[])
{
	double x = 1.5;
	double d; //一阶导 梯度方向
	double a=0.01; //学习率

	for (int i = 0; i < 100; i++)
	{
		d=g(x);
		//a = GetA_Armijo(x,d,10);
	}
	return 0;
}

逻辑回归和广义线性模型有什么关系呢？

线性回归

Softmax Regression：