机器学习-2.5-监督学习之生成型学习算法

生成型学习算法

1 生成型学习算法简述

我们之前讲到的算法。譬如逻辑回归，试图从\(x\)直接学习得到一组映射关系到\(y\)，即通过样本\((x,y)\)学习得到,使得其能够准确的预测\(y\)。这类算法叫做判别性学习算法。

本节我们将介绍生成学习型算法。举个例子，对于一个大象，我们可以为大象构造一个关于其特征的模型。对于一条狗，同样可以为狗构造一个关于关于其特征的模型。对于一个新动物，我们可以根据前述构建的模型，判断这个动物是像大象多一点还是像狗多一点，来确定这动物是大象还是狗。

我们可以根据样本得到大象的特征模型\(p(x|y=0)\)和狗的特征模型\(p(x|y=1)\)，同时也可以得到\(p(y)\)。根据贝叶斯公式，有如下：

$$\max (p(y|x)) = \max (\frac{p(x|y)p(y)}{p(x)}) = \max (\frac{p(x|y)p(y)}{p(x|y = 1)p(y = 1) + p(x|y = 0)p(y = 0)})$$

对于我们做预测的时候，可以不考虑分母。

$$\max (p(y|x)) = \max (p(x|y)p(y))$$

为什么不考虑分母?我们可以将\(p(x|y)p(y)\)简化为\(f(y)\),因此分母就是\(f(1)+f(0)\)。如果我们通过样本对完成了建模，因此分母就是一个常数了，因此没有计算的必要了。

2 高斯判别分析

2.1 多维高斯分布介绍

略，详见cs229-note2.pdf

2.2 高斯判别分析模型

对于一个分类问题,\(y \in (0,1)\)。假设\(x\)是连续的随机变量，服从高斯分布。具体描述如下：

$$y \sim Bernoulli(\phi)$$ $$x|y = 0 \sim N({\mu _0},\Sigma)$$ $$x|y = 1 \sim N({\mu _1},\Sigma)$$ $$p(x|y = 0) = \frac{1}{{{(2\pi )}^{\frac{n}{2}}}|\Sigma {|^{\frac{n}{2}}}}\exp ( - \frac{1}{2}{(x - {\mu _0})^T}{\Sigma ^{ - 1}}(x - {\mu _0}))$$ $$p(x|y = 1) = \frac{1}{{{(2\pi )}^{\frac{n}{2}}}|\Sigma {|^{\frac{n}{2}}}}\exp ( - \frac{1}{2}{(x - {\mu _1})^T}{\Sigma ^{ - 1}}(x - {\mu _1}))$$

2.3 公式推导

我们设置\(\phi = p(y = 0)\)。根据前面的分析，在给定样本的情况下，我们需要保证下式最大。

$$\prod\limits _{i = 1}^m {p({y^i}|{x^i})} = \prod\limits _{i = 1}^m {p({x^i}|{y^i})p({y^i})}$$

这里为了简写，我们设置\({z _0} = x - {\mu _0}\) ，\({z _1} = x - {\mu _1}\)。取自然对数后，得到最大释然函数:

$$\ell ({\mu _1},{\mu _1},\Sigma ,\phi ) = \sum\limits _{i = 1}^m {\ln \{ {{[\frac{1}{{{{(2\pi )}^{\frac{n}{2}}}|\Sigma {|^{\frac{n}{2}}}}}\exp ( - \frac{1}{2}z _0^T{\Sigma ^{ - 1}}{z _0})]}^{1\{ {y^i} = 0\} }}{{[\frac{1}{{{{(2\pi )}^{\frac{n}{2}}}|\Sigma {|^{\frac{n}{2}}}}}\exp ( - \frac{1}{2}z _1^T{\Sigma ^{ - 1}}{z _1})]}^{1\{ {y^i} = 1\} }}\} } $$ $$\ell ({\mu _0},{\mu _1},\Sigma ,\phi ) = \sum\limits _{i = 1}^m {\{ 1\{ {y^i} = 0\} [ - \frac{1}{2}z _0^T{\Sigma ^{ - 1}}{z _0} - \frac{n}{2}\ln 2\pi - \frac{n}{2}\ln |\Sigma | + \ln \phi ] + 1\{ {y^i} = 1\} } [ - \frac{1}{2}z _1^T{\Sigma ^{ - 1}}{z _1} - \frac{n}{2}\ln 2\pi - \frac{n}{2}\ln |\Sigma | + \ln (1 - \phi )]\}$$

然后我们分别对各个参数求偏导数:

$$\frac{{\partial \ell ({\mu _0},{\mu _1},\sum ,\phi )}}{{\partial \phi }} = \sum\limits _{i = 1}^m {(\frac{{1\{ {y^i} = 1\} }}{\phi } - \frac{{1\{ {y^i} = 0\} }}{{1 - \phi }})} = \sum\limits _{i = 1}^n {(\frac{{1\{ {y^i} = 1\} }}{\phi } - \frac{{1 - 1\{ {y^i} = 1\} }}{{1 - \phi }})} = 0$$

所以有:

$$\phi = \frac{1}{m}\sum\limits _{i = 1}^n {1\{ {y^i} = 1\} }$$

然后对\(\mu _0\)求导：

$${\nabla _{{\mu _0}}}\ell ({\mu _0},{\mu _1},\Sigma ,\phi ) = \sum\limits _{i = 1}^m {{\nabla _{{\mu _0}}}( - \frac{1}{2}{{({x^i} - {\mu _0})}^T}{\Sigma ^{ - 1}}({x^i} - {\mu _0}))1\{ {y^i} = 0\} }$$ $${\nabla _{{\mu _0}}}\ell ({\mu _0},{\mu _1},\Sigma ,\phi ) = - \frac{1}{2}\sum\limits _{i = 1}^m {{\nabla _{{\mu _0}}}( - \mu _0^T{\Sigma ^{ - 1}}{x^i} - {{({x^i})}^T}{\Sigma ^{ - 1}}{\mu _0} + \mu _0^T{\Sigma ^{ - 1}}{\mu _0})1\{ {y^i} = 0\} }$$

我们对上面的式子分头计算。

\(\Sigma\)为对角阵。由于求偏微分的数为一个常数，因此该值与其迹的相同。另外，这里使用了一些关于迹的公式。具体详见附录。

$${\nabla _{{\mu _0}}}\mu _0^T{\Sigma ^{ - 1}}x = {\nabla _{{\mu _0}}}tr(\mu _0^T{\Sigma ^{ - 1}}{x^i}) = {[{\nabla _{\mu _0^T}}tr(\mu _0^T{\Sigma ^{ - 1}}{x^i})]^T} = {[{({\Sigma ^{ - 1}}{x^i})^T}]^T} = {\Sigma ^{ - 1}}{x^i}$$ $${\nabla _{{\mu _0}}}{({x^i})^T}{\Sigma ^{ - 1}}{\mu _0} = {\nabla _{{\mu _0}}}tr({({x^i})^T}{\Sigma ^{ - 1}}{\mu _0}) = {\nabla _{{\mu _0}}}{\mu _0}{({x^i})^T}{\Sigma ^{ - 1}} = {({({x^i})^T}{\Sigma ^{ - 1}})^T} = {\Sigma ^{ - 1}}{x^i}$$ $${\nabla _{{\mu _0}}}\mu _0^T{\Sigma ^{ - 1}}{\mu _0} = {[{\nabla _{\mu _0^T}}\mu _0^T{\Sigma ^{ - 1}}{\mu _0}E]^T} = [E\mu _0^T{\Sigma ^{ - 1}} + {E^T}\mu _0^T{({\Sigma ^{ - 1}})^T}] = 2{\Sigma ^{ - 1}}{\mu _0}$$

所以有:

$${\nabla _{{\mu _0}}}\ell ({\mu _0},{\mu _1},\Sigma ,\phi ) = - \sum\limits _{i = 1}^m {{\nabla _{{\mu _0}}}({\Sigma ^{ - 1}}{\mu _0} - {\Sigma ^{ - 1}}{x^i})1\{ {y^i} = 0\} } = 0$$ $${\mu _0} = \frac{{\sum\limits _{i = 1}^m {1\{ {y^i} = 0\} {x^i}} }}{{\sum\limits _{i = 1}^m {1\{ {y^i} = 0\} } }}$$

同理有:

$${\mu _1} = \frac{{\sum\limits _{i = 1}^m {1\{ {y^i} = 1\} {x^i}} }}{{\sum\limits _{i = 1}^m {1\{ {y^i} = 1\} } }}$$

然后对\(\Sigma^{-1}\)求偏导数。

$${\nabla _{{\Sigma ^{ - 1}}}}\ell ({\mu _0},{\mu _1},\Sigma ,\phi ) = \sum\limits _{i = 1}^m {{\nabla _{{\Sigma ^{ - 1}}}}[( - \frac{1}{2}z _0^T{\Sigma ^{ - 1}}{z _0} + \frac{1}{2}\ln \frac{1}{{|\Sigma |}})1\{ {y^i} = 0\} + ( - \frac{1}{2}z _1^T{\Sigma ^{ - 1}}{z _1} + \frac{1}{2}\ln \frac{1}{{|\Sigma |}})1\{ {y^i} = 1\} ]}$$

这里分开计算:

$${\nabla _{\Sigma ^{ - 1}}}(z _0^T{\Sigma ^{ - 1}}{z _0}) = {\nabla _{\sum ^{ - 1}}}tr(z _0^T{\Sigma ^{ - 1}}{z _0}) = {\nabla _{\sum ^{ - 1}}}tr({\Sigma ^{ - 1}}{z _0}z _0^T) = {({z _0}z _0^T)^T} = {z _0}z _0^T$$ $${\nabla _{\Sigma ^{ - 1}}}\ln \frac{1}{|\Sigma |} = {\nabla _{\Sigma ^{ - 1}}}\ln |{\Sigma ^{ - 1}}| = \frac{1}{|{\Sigma ^{ - 1}}|}{\nabla _{\Sigma ^{ - 1}}}{\Sigma ^{ - 1}} = {({({\Sigma ^{ - 1}})^{ - 1}})^T}$$

所以有:

$${\nabla _{{\sum ^{ - 1}}}}\ell ({\mu _0},{\mu _1},\Sigma ,\phi ) = \sum\limits _{i = 1}^m {[( - \frac{1}{2}{z _0}z _0^T + \frac{1}{2}\Sigma )1\{ {y^i} = 0\} + ( - \frac{1}{2}{z _1}z _1^T + \frac{1}{2}\Sigma )1\{ {y^i} = 1\} ]} = 0$$ $$\Sigma = \frac{1}{m}\sum\limits _{i = 1}^m {[({x^i} - {\mu _0}){{({x^i} - {\mu _0})}^T}1\{ {y^i} = 0\} + ({x^i} - {\mu _1}){{({x^i} - {\mu _1})}^T}1\{ {y^i} = 1\} ]}$$ $$\Sigma = \frac{1}{m}\sum\limits _{i = 1}^m {[({x^i} - {\mu _{{y^i}}}){{({x^i} - {\mu _{{y^i}}})}^T}]}$$

2.3 高斯型算法实例

我们分别以 \((1,1)\)和\((2,2)\)为均值生成一组高斯分布。下面图是生成的样本。

前提假设是我们知道两组分类是符合高斯分布的，切假定协方差相同。但我们不知道两组数据的均值和协方差。根据前面的公式有如下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import multivariate_normal
import random
from matplotlib.lines import Line2D

## this is a program about GDA

global MAKE_DATA
global SHOW_PIC

def make_data(x):
  if MAKE_DATA == True:
    mean_1 = [1,1]
    mean_2 = [2,2]
    cov = [[0.1,0],[0,0.1]]
    arr1 = np.random.multivariate_normal(mean_1, cov, 100)
    arr2 = np.random.multivariate_normal(mean_2, cov, 50)
    print(arr1)
    print(arr2)
    x.append(arr1)
    x.append(arr2)
  if MAKE_DATA == False:
    arr1 = np.array([[ 0.88235916 , 1.01511634],[ 0.75243817 , 0.76520033],[ 0.95710848 , 1.41894337],[ 1.48682891 , 0.78885043],[ 1.24047011 , 0.71984948],[ 0.67611276 , 1.07909452],[ 1.03243669 , 1.08929695],[ 1.0296548  , 1.25023769],[ 1.54134008 , 0.39564824],[ 0.34645057 , 1.61499636],[ 0.77206174 , 1.23613698],[ 0.91446988 , 1.38537765],[ 0.99982962 , 1.34448471],[ 0.78745962 , 0.9046565 ],[ 0.74946602 , 1.07424473],[ 1.09294839 , 1.14711993],[ 0.39266844 , 0.78788004],[ 0.83112357 , 1.2762774 ],[ 1.05056188 , 1.13351562],[ 1.62101523 , 1.15035562],[ 0.70377517 , 1.1136416 ],[ 1.03715472 , 0.47905693],[ 0.94598381 , 0.8874837 ],[ 0.94447128 , 2.02796925],[ 0.72442242 , 1.09835206],[ 0.69046731 , 1.46232182],[ 1.20744606 , 1.10280041],[ 0.70665746 , 0.82139503],[ 1.08803887 , 1.4450361 ],[ 0.88530961 , 0.75727475],[ 0.98418545 , 0.80248161],[ 0.74970386 , 1.13205709],[ 0.72586454 , 1.06058385],[ 0.9071812  , 1.09975063],[ 0.75182835 , 0.93570147],[ 0.80052289 , 1.08168507],[ 0.40180652 , 0.9526211 ],[ 0.62312617 , 0.84385058],[ 0.68212516 , 1.25912717],[ 1.19773245 , 0.16399654],[ 0.96093132 , 0.43932091],[ 1.25471657 , 0.92371829],[ 1.12330272 , 1.26968747],[ 1.30361985 , 0.99862123],[ 1.23477665 , 1.1742804 ],[ 0.28471876 , 0.5806044 ],[ 1.89355099 , 1.19928671],[ 1.09081369 , 1.28467312],[ 1.40488635 , 0.90034427],[ 1.11672364 , 1.49070515],[ 1.35385212 , 1.35767891],[ 0.92746374 , 1.79096697],[ 1.89142562 , 0.98228303],[ 1.0555218  , 0.86070833],[ 0.69001255 , 1.12874741],[ 0.98137315 , 1.3398852 ],[ 1.02525371 , 0.77572865],[ 1.1354295  , 1.07098552],[ 1.50829164 , 1.43065998],[ 1.09928764 , 1.55540292],[ 0.64695084 , 0.79920395],[ 0.82059034 , 0.97533491],[ 0.56345455 , 1.08168272],[ 1.06673215 , 1.19448556],[ 0.96512548 , 1.5268577 ],[ 0.96914451 , 1.00902985],[ 0.72879413 , 0.92476415],[ 1.0931483  , 1.13572242],[ 1.34765121 , 0.83841006],[ 1.57813788 , 0.65915892],[ 0.59032608 , 0.82747946],[ 0.83838504 , 0.67588473],[ 1.35101322 , 1.21027851],[ 0.71762153 , 0.41839038],[ 0.61295604 , 0.66555018],[ 0.64379346 , 0.92925228],[ 1.1194968  , 0.65876736],[ 0.39495437 , 0.67246734],[ 1.05223282 , 0.17889116],[ 0.97810984 , 1.12794664],[ 0.98392719 , 0.73590255],[ 1.25587405 , 1.21853038],[ 1.01150226 , 1.01835571],[ 1.02251614 , 0.72704228],[ 1.00261519 , 0.95347185],[ 0.96362523 , 0.8607009 ],[ 0.88034659 , 1.2307104 ],[ 0.75907236 , 0.92799796],[ 0.54898709 , 1.69882285],[ 0.55032649 , 0.98831566],[ 1.33360789 , 1.19793298],[ 0.83231239 , 0.8946538 ],[ 1.05173094 , 1.26324289],[ 0.81482231 , 0.56198584],[ 1.03854797 , 1.0553811 ],[ 1.32669227 , 1.61115811],[ 1.13322152 , 1.68151695],[ 0.39754618 , 1.19392967],[ 0.61344185 , 1.05281434],[ 1.18415366 , 0.864884  ]])
    arr2 = np.array([[ 2.15366548 , 1.88035458],[ 2.36978774 , 1.76550283],[ 2.46261387 , 2.10568262],[ 1.90475526 , 1.95242885],[ 1.77712677 , 1.96004856],[ 1.5995514  , 2.1323943 ],[ 1.52727223 , 1.50295551],[ 1.80330407 , 1.57942301],[ 1.86487049 , 1.87234414],[ 1.9586354  , 1.96279729],[ 2.59668134 , 2.414423  ],[ 2.818419   , 1.76280366],[ 2.01511628 , 2.10858546],[ 2.15907962 , 1.81593012],[ 1.63966834 , 2.2209023 ],[ 2.47220599 , 1.70482956],[ 2.08760748 , 2.51601971],[ 1.50547722 , 1.8487145 ],[ 1.68125583 , 2.64968501],[ 2.01924282 , 2.0953572 ],[ 2.22563534 , 2.18266325],[ 2.2684291  , 2.23581599],[ 2.13787557 , 1.9999382 ],[ 1.02638695 , 1.68134967],[ 2.35614619 , 1.32072125],[ 2.20054871 , 1.47401445],[ 1.99454827 , 1.71658741],[ 1.83269065 , 2.47662909],[ 2.40097251 , 2.21823862],[ 2.54404652 , 1.85742018],[ 1.84150027 , 2.06350351],[ 1.69490855 , 1.70169334],[ 1.44745704 , 1.88295233],[ 2.24376639 , 1.67530495],[ 1.42911921 , 1.81854548],[ 1.33789289 , 2.27686128],[ 2.43509821 , 1.95032131],[ 1.9512447  , 1.4595415 ],[ 2.13041192 , 1.79372755],[ 2.2753866  , 2.23781951],[ 2.26753401 , 1.78149305],[ 2.06505449 , 2.01939606],[ 2.44426826 , 2.1437101 ],[ 2.16607141 , 2.31077167],[ 1.96097237 , 2.49100193],[ 1.37255424 , 1.60735016],[ 1.63947758 , 2.17852314],[ 2.13722666 , 2.00559707],[ 1.222696   , 1.67075059],[ 2.56982685 , 2.51218813]])
    x.append(arr1)
    x.append(arr2)
  if SHOW_PIC == 1:
    figure, ax = plt.subplots()
    ax.set_xlim(left=-1, right=4)
    ax.set_ylim(bottom=-1, top=4)
    for i in range(len(arr1)):
      plt.plot(arr1[i][0], arr1[i][1], 'b--', marker='+', color='g')
    for i in range(len(arr2)):
      plt.plot(arr2[i][0], arr2[i][1], 'b--', marker='o', color='b')
    plt.xlabel("x1")
    plt.ylabel("x2")
    plt.plot()
    plt.show()

def calcPhi(x):
  return len(x[1])/(len(x[0])+len(x[1]))

def calcMu(x,i):
  return x[i].mean(axis=0)

def calcSigma(x,mu_0,mu_1):
  sum=np.array([[0,0],[0,0]])
  x0=x[0]
  for i in range(len(x0)):
    z0=np.array([x0[i]-mu_0])
    z0T = np.array([x0[i]-mu_0]).transpose()
    sum = sum + np.dot(z0T, z0)
    print(np.dot(z0T, z0))
  x1=x[1]
  for i in range(len(x1)):
    z1=np.array([x1[i]-mu_1])
    z1T = np.array([x1[i]-mu_1]).transpose()
    sum = sum + np.dot(z1T, z1)
    print(np.dot(z1T, z1))
  return sum/(len(x0)+len(x1))

def classify(x,mu_0,mu_1,sigma,phi):
  var0 = multivariate_normal(mean=mu_0.tolist(), cov=sigma.tolist())
  var1 = multivariate_normal(mean=mu_1.tolist(), cov=sigma.tolist())
  # p(y=1|x) = p(x|y=1)/p(y=1) / p(x) = p(x|y=1)p(y=1) / (p(x|y=1)p(y=1)+p(x|y=0)p(y=0))
  return var1.pdf(x)*phi/(var1.pdf(x)*phi+var0.pdf(x)*(1-phi)) > 0.5

def testClassify(mu_0,mu_1,sigma,phi):
  if SHOW_PIC != 2:
    return
  figure, ax = plt.subplots()
  ax.set_xlim(left=-1, right=4)
  ax.set_ylim(bottom=-1, top=4)
  for i in range(100):
    x1 = random.randint(0,3000)/1000
    x2 = random.randint(0,3000)/1000
    if( classify([x1,x2],mu_0,mu_1,sigma,phi)) == True:
      plt.plot(x1, x2, 'b--', marker='+', color='g')
    else:
      plt.plot(x1, x2, 'b--', marker='o', color='b')
  # draw y = -x+3
  (line1_xs, line1_ys) = [(0, 3), (3, 0)]
  ax.add_line(Line2D(line1_xs, line1_ys, linewidth=1, color='blue'))
  plt.xlabel("x1")
  plt.ylabel("x2")
  plt.plot()
  plt.show()

if __name__ == "__main__":
  # 0 debug option
  MAKE_DATA = False
  SHOW_PIC = 1        # 0 - do not show   1 - show sample data    2 - show test data

  # 1 make the sample
  # the format of x is [ndarray0,ndarray1] , ndarray1 is the set of the first class, we set y=0.
  # ndarray1 is the set of the second class, we set y=1
  x = []
  make_data(x)

  if MAKE_DATA == True:
    exit(0)

  # 2 learn from the sample
  print("hello")
  phi = calcPhi(x)                # phi  = p(y=1)  means close (2,2)
  mu_0 = calcMu(x,0)
  mu_1 = calcMu(x,1)
  sigma = calcSigma(x,mu_0,mu_1)
  print(phi,", ",mu_0,", ",mu_1)
  print(sigma.tolist())

  print("---------")
  # 3 test classify
  testClassify(mu_0,mu_1,sigma,phi)

对于两个协方差相同的样本，我们知道两组数据有相同的概率分布曲线，具体都应该是一个个同心圆。因此，我们可以确定一条垂直于两点连线的曲线是对该样本的理论上正确的分割。因此，我们得到如下结果，并与理论分割的曲线，即\(y=-x+3\)比较，发现对我们随机产生的样本分类效果非常好。

附录 A 相关公式

Hessian矩阵定义:

$${\nabla _A}f(A) = \left( {\begin{array}{*{20}{c}} {\frac{{\partial f}}{{\partial {A _{11}}}}}& \ldots &{\frac{{\partial f}}{{\partial {A _{1n}}}}} \\\ \vdots & \ddots & \vdots \\\ {\frac{{\partial f}}{{\partial {A _{n1}}}}}& \cdots &{\frac{{\partial f}}{{\partial {A _{nn}}}}} \end{array}} \right)$$

矩阵的迹的定义:

$$trA = \sum\limits _{i = 1}^n {{A _{ii}}}$$

关于矩阵迹的公式:

$$tr(a) = a$$ $$tr(aA) = atr(A)$$ $$tr(aA) = atr(A)$$ $$tr(ABC) = tr(CAB) + tr(BCA)$$ $$trA = tr{A^T}$$ $$tr(A + B) = trA + trB$$ $${\nabla _A}tr(AB) = {B^T}$$ $${\nabla _{A^T}}f(A) = {({\nabla _A}f(A))^T}$$ $${\nabla _{A^T}}trAB{A^T}C = CAB + {C^T}A{B^T}$$ $${\nabla _{A^T}}trAB{A^T}C = CAB + {C^T}A{B^T}$$

对矩阵做展开能证明上述大部分公式，这里暂时略。

附录B

手写公式和word版博客地址:

1
https://github.com/zhengchenyu/MyBlog/tree/master/doc/mlearning/生成型学习算法

机器学习-2.4-监督学习之softmax

Softmax Regression

下面介绍一下指数族分布的另外一个例子。之前的逻辑回归中，可以用来解决二分类的问题。对于有k中可能结果的时候，问题就转化为多分类问题了，也就是接下来要说明的softmax regression问题。
s

1 函数定义

在进行下一步推导前，我们先定义部分辅助函数。我们这里事先定义一个k-1维向量\(T(y)\)，这里具体对应于指数族分布的\(T(y)\),具体定义如下:

$$T(1) = \left[ \begin{gathered} 1 \hfill \\\ 0 \hfill \\\ 0 \hfill \\\ ... \hfill \\\ 0 \hfill \\ \end{gathered} \right],T(2) = \left[ \begin{gathered} 0 \hfill \\\ 1 \hfill \\\ 0 \hfill \\\ ... \hfill \\\ 0 \hfill \\ \end{gathered} \right],T(3) = \left[ \begin{gathered} 0 \hfill \\\ 0 \hfill \\\ 1 \hfill \\\ ... \hfill \\\ 0 \hfill \\\ \end{gathered} \right],...,T(k - 1) = \left[ \begin{gathered} 0 \hfill \\\ 0 \hfill \\\ 0 \hfill \\\ ... \hfill \\\ 1 \hfill \\ \end{gathered} \right],T(k) = \left[ \begin{gathered} 0 \hfill \\\ 0 \hfill \\\ 0 \hfill \\\ ... \hfill \\\ 0 \hfill \\ \end{gathered} \right]$$

我们定义\(T(y) _1\)表示的\(T(y)\)第一个元素，其他依次类推。

然后再引入一个函数，具体定义如下:

$$1\{ true\} = 1$$ $$1\{ false\} = 0$$

可以通过带入值的方式验证上式。

2 模型推导

在softmax中,我们知道\(y \in \{ 1,2,…,k\}\)。我们设置\(y=1\)的概率为\(\phi _1\),
类似的有\(y=k\)的概率为\(\phi _k\)，并且我们轻易得到如下公式：

$$\sum\limits _{i = 1}^k {\phi _i} = 1$$ $${\phi _k} = \sum\limits _{i = 1}^{k - 1}{\phi _i}$$

由于我们已经知道\(p(y = i) = {\phi _i}\)，综上可以得到:

$$p(y) = \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 - \sum\limits _{i = 1}^{k - 1} {1\{ y = i\} } } = \phi _1^{T{(y)} _1}\phi _2^{T{{(y)} _2}}...\phi _k^{1 - \sum\limits _{i = 1}^{k - 1} {T{(y)} _i} }$$ $$p(y) = \exp (T{(y) _1}\ln {\phi _1} + T{(y) _2}\ln {\phi _2} + ... + (1 - \sum\limits _{i = 1}^{k - 1} {T{{(y)} _i}} )\ln {\phi _k})$$ $$p(y) = \exp (T{(y) _1}\ln \frac{\phi _1}{\phi _k} + T{(y) _2}\ln \frac{\phi _2}{\phi _k} + ... + T{(y) _{k - 1}}\ln \frac{\phi _{k - 1}}{\phi _k} + \ln {\phi _k})$$

我们对照指数族分布，其中\(T(y)\)我们已经定义，可以得到其他参数如下:

$$b(y) = 1$$ $$a(\eta ) = - \ln {\phi _k}$$ $$\eta = \left[ \begin{gathered} \ln \frac{\phi _1}{\phi _k} \hfill \\\ \ln \frac{\phi _2}{\phi _k} \hfill \\\ ... \hfill \\\ \ln \frac{\phi _{k - 1}}{\phi _k} \hfill \\\ \end{gathered} \right]$$

我们用\(\eta\)来表示\(\phi\)，目标是将我们的问题归整为一个变量的问题，从而更容易计算。对上面向量拆开计算得到:

$${\eta _i} = \ln \frac{\phi _i}{\phi _k}$$

继而可以转化为:

$${\eta _i} = \ln \frac{\phi _i}{\phi _k}$$

又由于我们已知\(\sum\limits _{j = 1}^k {\phi _j} = 1\)，所以有:

$$\sum\limits _{i = 1}^k {e^{\eta _i}{\phi _k}} = 1$$

继而有:

$${\phi _k} = \frac{1}{\sum\limits _{j = 1}^k {e^{\eta _j}}}$$

最后我们可以得到:

$${\phi _i} = \frac{e^{\eta _i}}{\sum\limits _{j = 1}^k {e^{\eta _j}}}$$

由于指数组分布假设是关于输入的线性函数，所以得到在已知\(x\)和\(\theta\)的情况下\(y=i\)的概率的公式，如下：

$$p(y = i|x,\theta ) = {\phi _i} = \frac{e^{\eta _i}}{\sum\limits _{j = 1}^k {e^{\eta _j}}} = \frac{e^{\theta _i^Tx}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}}$$

这里我们可以得到\({h_\theta }(x)\),如下:

$${h_\theta }(x) = E[T(y)|x,\theta ] = E[\left( \begin{gathered} 1\{ y = 1\} \hfill \\\ 1\{ y = 2\} \hfill \\\ ... \hfill \\\ 1\{ y = k - 1\} \hfill \\\ \end{gathered} \right)|x,\theta ]$$ $${h_ \theta }(x) = E[\left( \begin{gathered} {\phi _1} \hfill \\\ {\phi _2} \hfill \\\ ... \hfill \\\ {\phi _{k - 1}} \hfill \\\ \end{gathered} \right)|x,\theta ] = \left( \begin{gathered} \frac{e^{\theta _1^Tx}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}} \hfill \\\ \frac{e^{\theta _2^Tx}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}} \hfill \\\ ... \hfill \\\ \frac{e^{\theta _{k - 1}^Tx}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}} \hfill \\\ \end{gathered} \right)$$

其中,\(k\)为取值的可能集合的大小，\(\theta\)为一个\(k*n\)的矩阵。

到这里我们可以定义我们的损失函数了，如下:

$$J(x) = \frac{1}{2}\sum\limits _{i = 1}^m {{(T(y) - {h _\theta }(x))}^2}$$

上式是一个向量，我们可以对这个向量继续求平方和，来衡量准确度，如下：

$$J(x) = \frac{1}{2}\sum\limits _{i = 1}^m {|{{(T(y) - {h _\theta }(x))}^2}|}$$

损失函数已经定义完成，我们就有了算法的截止条件了。接下来就是找到算法的最优迭代方向，也就是计算其偏导数了。

我们使用最大释然性来计算迭代方向。建模的原则是：对于一个已知\(y=i\)的样本\((x,y)\)，我们需要找到一个方向去迭代\(\theta\)，使得\(\phi _i\)尽可能大，使得其他\(\phi\)尽可能小。当然由于所有\(\phi\)之和为1,所以我们只需要保证\(\phi _i\)尽可能大即可。因此,在已经的\(m\)个样本的情况下，我们只需要保证下面的式子最大:

$$\sum\limits _{i = 1}^m {p(y = {y^i}|x,\theta )} = \sum\limits _{i = 1}^m {\prod\limits _{l = 1}^k {p{(y = l|x,\theta )}^{1\{ {y^i} = l\}}} }$$

取自然对数后，我们定义最大似然函数:

$$l(\theta ) = \sum\limits _{i = 1}^m {\ln \prod\limits _{l = 1}^k {(\frac{e^{\theta _l^Tx}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}})^{1\{{y^i} = l\}}}}$$

当且仅当\({y^i} = l\)的时候, \({1{yi=l\} = 1})。其他值为0，我们可以将连乘转化，如下：

$$l(\theta ) = \sum\limits _{i = 1}^m {\ln (\frac{e^{\theta _{y^i}^Tx}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}})} = \sum\limits _{i = 1}^m {\ln [\theta _{y^i}^Tx - \ln \sum\limits _{j = 1}^k {e^{\theta _j^Tx}}]} $$

然后对其求偏导数,对\(y _i=f\)的时候后有：

$$\frac{\partial l(\theta )}{\partial {\theta _f}} = \sum\limits _{i = 1}^m {\ln [{x _f} - \frac{e^{\theta _f^Tx}{x _f}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}}]}$$

对\({y_i} \ne f\)的时候，有：

$$\frac{\partial l(\theta )}{\partial {\theta _f}} = \sum\limits _{i = 1}^m {\ln [0 - \frac{e^{\theta _f^Tx}{x _f}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}}]}$$

综上，可以得到:

$$\frac{\partial l(\theta )}{\partial {\theta _f}} = \sum\limits _{i = 1}^m {\ln [1\{ {y^i} = f\} - \frac{e^{\theta _f^Tx}}{\sum\limits _{j = 1}^k {e^{\theta _j^Tx}}}]{x _f}} = \sum\limits _{i = 1}^m {\ln [1\{{y^i} = f\} - {\phi _f}]{x _f}}$$

\(x _f\)是一个向量，因为公式一次迭代更新一个维度下的一组\(\theta\)。事实上，这里是一个向量求微分。如果我们对的某个元素\(\theta\)进行微分，我们依然能够得到这个公式。

3 实际问题的解决

我们随机制造一组数据，在\(2*2\)的空间内使用直线\({x _2} = 0.5*{x _1}\)和直线\({x _2} = -0.5*{x _1}+2\)具体的分类如下:

根据之前的推导，编写如下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
import random
import math
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.lines import Line2D

g_x_arr=[[1, 0.035, 1.344], [1, 0.662, 0.598], [1, 1.791, 1.889], [1, 0.158, 0.12], [1, 1.55, 1.835], [1, 2.0, 0.613], [1, 1.176, 0.368], [1, 0.564, 0.043], [1, 1.559, 1.507], [1, 1.998, 0.988], [1, 0.082, 0.941], [1, 0.542, 1.371], [1, 0.542, 0.671], [1, 1.34, 1.856], [1, 0.049, 0.089], [1, 1.933, 0.871], [1, 1.753, 1.024], [1, 0.315, 1.341], [1, 0.829, 1.26], [1, 0.686, 1.721], [1, 1.222, 1.129], [1, 0.55, 0.075], [1, 0.767, 0.346], [1, 1.516, 1.752], [1, 1.347, 0.905], [1, 0.127, 0.782], [1, 1.169, 1.272], [1, 1.301, 0.273], [1, 0.081, 0.739], [1, 0.203, 0.658], [1, 0.347, 1.064], [1, 0.793, 1.193], [1, 1.428, 0.326], [1, 0.509, 0.983], [1, 0.12, 0.884], [1, 0.251, 0.282], [1, 0.73, 0.445], [1, 1.889, 1.323], [1, 1.314, 1.795], [1, 1.297, 1.467], [1, 1.669, 0.613], [1, 0.753, 0.114], [1, 0.94, 1.972], [1, 0.738, 1.603], [1, 1.508, 1.237], [1, 0.979, 0.572], [1, 0.128, 1.254], [1, 0.569, 0.155], [1, 0.88, 0.211], [1, 0.405, 0.603], [1, 1.02, 1.9], [1, 0.438, 1.535], [1, 1.506, 1.638], [1, 1.712, 0.394], [1, 0.556, 0.124], [1, 0.444, 0.115], [1, 0.595, 1.009], [1, 0.165, 0.089], [1, 1.57, 0.634], [1, 1.429, 1.181], [1, 0.8, 0.671], [1, 1.914, 1.091], [1, 0.594, 0.569], [1, 0.935, 0.277], [1, 0.47, 0.522], [1, 0.94, 1.924], [1, 0.194, 1.933], [1, 0.612, 0.613], [1, 0.236, 0.894], [1, 1.888, 0.251], [1, 1.548, 0.191], [1, 1.543, 0.603], [1, 1.521, 0.02], [1, 0.923, 0.856], [1, 0.649, 1.31], [1, 0.379, 1.746], [1, 1.345, 0.902], [1, 0.937, 0.524], [1, 1.018, 0.68], [1, 1.738, 1.623], [1, 1.534, 1.9], [1, 0.139, 1.911], [1, 1.508, 1.173], [1, 0.798, 0.865], [1, 0.451, 1.186], [1, 1.63, 1.123], [1, 0.82, 0.848], [1, 1.213, 1.48], [1, 0.894, 0.664], [1, 1.456, 0.934], [1, 0.59, 1.525], [1, 0.522, 1.329], [1, 1.179, 1.396], [1, 0.527, 0.273], [1, 1.399, 1.215], [1, 0.966, 1.514], [1, 1.341, 0.028], [1, 0.479, 0.191], [1, 1.193, 0.724], [1, 0.714, 1.285]]
g_Ty_arr=[[0, 0], [0, 0], [0, 1], [0, 0], [0, 1], [1, 0], [1, 0], [1, 0], [0, 1], [1, 0], [0, 0], [0, 0], [0, 0], [0, 1], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 1], [0, 0], [1, 0], [1, 0], [0, 1], [0, 0], [0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 1], [0, 1], [0, 1], [1, 0], [1, 0], [0, 1], [0, 0], [0, 0], [0, 0], [0, 0], [1, 0], [1, 0], [0, 0], [0, 1], [0, 0], [0, 1], [1, 0], [1, 0], [1, 0], [0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 1], [0, 0], [1, 0], [0, 0], [0, 1], [0, 1], [0, 0], [0, 0], [1, 0], [1, 0], [1, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 1], [0, 1], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 1], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [1, 0], [1, 0], [0, 0], [0, 0]]
g_y_arr=[2, 2, 1, 2, 1, 0, 0, 0, 1, 0, 2, 2, 2, 1, 2, 0, 2, 2, 2, 1, 2, 0, 0, 1, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 1, 1, 1, 0, 0, 1, 2, 2, 2, 2, 0, 0, 2, 1, 2, 1, 0, 0, 0, 2, 2, 0, 2, 2, 1, 2, 0, 2, 1, 1, 2, 2, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 2, 2]

def phi(i,theta,x):
	theta_n = len(theta)
	numerator = math.exp(np.dot(theta[i], x))
	denominator = 0;
	for j in range(theta_n):                      # maybe need to optimize
		denominator = denominator + math.exp(np.dot(theta[j],x))
	return numerator/denominator

def h_theta(theta,x):
	theta_n = len(theta)
	ret = []
	for i in range(theta_n):
		ret.append(phi(i,theta,x))
	return ret

def one(y,i):
	if y == i:
		return 1
	else:
		return 0

def derl(f,x_arr,y_arr,theta):
	# 0.
	# f is the dimension to calc the partial derivative
	# 1. calc the size
	# theta_n is the size of theta, equals to the length of Ty' result set -1 .Here, equals lenght of {0,1,2} = 3-1 =2
	# m is the sum of samples
	# x_n is the dimension of the input variable 'x' + 1 (for x_0 =1). Here is 2 + 1 = 3
	theta_n = len(theta)
	m = len(x_arr)
	x_n = len(x_arr[0])
	# initial the output variable sum.Here is a vector, and the length of this vector is x_n
	sum = []
	for x_dim_index in range(x_n):
			sum.append(0)
	# 2. calc the partial derivative
	for i in range(x_n):
		sum[i] = 0
		for j in range(m):
			sum[i]=sum[i]+(one(y_arr[j],f)-phi(f,theta,x_arr[j]))*x_arr[j][i]
		sum[i]=sum[i]/m
	#sum=plus_vector(sum,theta[f],0.1)
	return sum

def plus_vector(arr1,arr2,a):
	return [x + a * y for x, y in zip(arr1, arr2)]

def update(theta,a,x_arr,y_arr):
	theta_n = len(theta)
	for i in range(theta_n):
		theta[i]=plus_vector(theta[i],derl(i, x_arr,y_arr,theta),a)
	return theta

def judge1(theta,x_arr,y_arr_vector,limit,debug):
	j_theta = J(theta,x_arr,y_arr_vector)
	if j_theta < limit:
		return True
	if debug:
		print("|J_theta(x)| = ", j_theta,"\n")
	return False

def J(theta,x_arr,y_arr_vector):
	sum=0
	m = len(x_arr)
	y_len = len(y_arr_vector[0])
	for i in range(m):
		for j in range(y_len):
			sum = sum + (phi(j,theta,x_arr[i]) - y_arr_vector[i][j])**2
	return sum

def calTheVaule():
	#theta=[[theta_1_1,theta_1_2],[theta_2_1,theta_2_2],[theta_3_1,theta_3_2]]
	theta = [[1,1,1],[1,1,1],[1,1,1]]
	a = 1
	count = 0
	while 1:
		if judge1(theta,g_x_arr,g_Ty_arr,1,True):
			break
		theta=update(theta,a,g_x_arr,g_y_arr)
		count = count + 1
		print("count=",count,"theta=",theta)

def calcIndex(theta,x):
  phi0 = phi(0,theta,x)
  phi1 = phi(1,theta,x)
  phi2 = phi(2,theta,x)
  if phi0 >= phi1 and phi0 >=phi1:
    return 0
  elif phi1>= phi0 and phi1 >= phi2:
    return 1
  else:
    return 2

def testValue():
  # limit =1， count = 33844
  theta = [[17.418683742474045, 13.989348587756815, -37.731307869980014],[-33.56035031194091, 1.3797244978249763, 34.095843086377855],[19.97548491920147, -11.65484983882828, 7.628008093823735]]
  samples = 100
  figure, ax = plt.subplots()
  # 设置x，y值域
  ax.set_xlim(left=0, right=2)
  ax.set_ylim(bottom=0, top=2)
  # 两条line的数据
  (line1_xs, line1_ys) = [(0, 2), (0, 1)]
  (line2_xs, line2_ys) = [(0, 2), (2, 1)]
  # 创建两条线，并添加
  ax.add_line(Line2D(line1_xs, line1_ys, linewidth=1, color='blue'))
  ax.add_line(Line2D(line2_xs, line2_ys, linewidth=1, color='blue'))
  for i in range(samples):
    x_0 = random.randint(0,2000)/1000
    x_1 = random.randint(0,2000)/1000
    index = calcIndex(theta,[1,x_0,x_1])
    if index == 0:
      plt.plot(x_0, x_1, 'b--', marker='x', color='r')
    elif index == 1:
      plt.plot(x_0, x_1, 'b--', marker='+', color='g')
    else:
      plt.plot(x_0, x_1, 'b--', marker='o', color='b')
  plt.xlabel("x1")
  plt.ylabel("x2")
  plt.plot()
  plt.show()

if __name__== "__main__":
  #calTheVaule()
  testValue()

经过33844次迭代之后，我们得到\(\theta\) = [[17.418683742474045, 13.989348587756815, -37.731307869980014],[-33.56035031194091, 1.3797244978249763, 34.095843086377855],[19.97548491920147, -11.65484983882828, 7.628008093823735]]

迭代的次数越多数据会越准确，这里迭代了上三万次，实际上可以通过牛顿法来减少迭代的次数，这里就不再重新代码了，牛顿法可以参见前面的文章。

然后，我们随机制造一组值，看看分类效果，具体如下：

附录 公式推导手写版

机器学习-2.3-监督学习之通用线性模型

Generalized Linear Models

回归分析中我们使用了高斯分布,分类问题服从伯努利分布。这两种方法都是第一个更通用的分布的特例，我们称其Generalized Linear Models (GLMs)

1 指数族分布

我们定义如下形式为指数族分布：

$$p(y;\eta ) = b(y)exp({\eta ^T}T(y) - a(\eta ))$$

其中,\(\eta \)是自然参数(natural parameter )，\(T(y) \)叫做充分统计量(sufficient statistic),一般\(T(y)=y\)。\(a(\eta) \)是log partition function 。\({e^{ - a\left( \eta \right)}}\)在归整化很重要，这个参数确保\(p(y;\eta )\)归整为1。

2 指数族分布的特例

其中伯努利分布和高斯分布都是指数族分布的一个特例。

2.1 伯努利分布的指数族表达

伯努利分布的指数表示:

$$p(y,\varphi ) = {\varphi ^y}{(1 - \varphi )^{1 - y}} = \exp (y\log \frac{\varphi }{1 - \varphi } + \log (1 - \varphi ))$$

对照指数函数族的表达式，有\(b(y) = 1\),\(T(y) = y\),\(\eta = \log \frac{\varphi }{1 - \varphi}\)等价于\(\varphi=\frac{e\eta}{e\eta + 1}\),\(a(y) = \log \frac{1}{1 - \varphi } = \log ({e^\eta } + 1)\)

2.2 高斯分布的指数族表达

高斯分布的指数表示:

$$p(y;\mu ) = \frac{1}{{\sqrt {2\pi } }}\exp ( - \frac{1}{2}{(y - \mu )^2}) = \frac{1}{{\sqrt {2\pi } }}\exp ( - \frac{1}{2}{y^2})\exp (\mu y - \frac{1}{2}{\mu ^2})$$

对照指数函数族的表达式，有\(b(y) = \frac{1}{\sqrt {2\pi }}\exp ( - \frac{1}{2}{y^2})\),\( \eta = \mu \),\(T(y) = y\),\(a(\eta ) = \frac{1}{2}{u^2} = \frac{1}{2}{\eta ^2}\)。

3 指数函数组的三个假设

三个假设分别为:

  • \(y|x;\theta\)服从指数分布
  • 给定\(x\)的情况下，我们的目标是预测\(T(y)\)。在大多数的例子中，\(T(y)=y\),因此就意味着我们通过学习得到函数h, 使得\(h(x)=E(y|x)\), 意思就是给定x，得到函数h使得其与训练样本的数学期望相等。
  • 自然参数\(\theta\)与输入\(x\)是线性关系。

对照上面的过程，我们可知伯努利分布和高斯分布都符合这三个假设

4 GLMs公式推导

还是之前的问题，已知有\(m\)个样本，通过学习模拟正确的\(h(x)\),使得其能够有效的推测\(y\)。(这里我们假定了\(T(y)=y\))

得到最大释然函数:

$$l(\eta ) = \sum\limits_{i = 1}^m {b({y^i}){\text{ }}exp({\eta ^T}T({y^i}) - a(\eta ))} $$

然后为了求其最大值，我们对其取自然对数后求倒数。

$$\frac{{\partial \ln (l(\eta ))}}{{\partial \eta }} = \sum\limits _{i = 1}^m {\frac{{\partial (\ln b({y^i}) + {\eta ^T}T({y^i}) - a(\eta ))}}{{\partial \eta }}} = \sum\limits _{i = 1}^m {(T({y^i})} - \frac{{\partial a(\eta )}}{{\partial \eta }})$$

在我们之前指数族分布的定义3中，\(\eta\)是各个维度\(x\)的线性表达。有\(\eta = {\theta ^T}x\)。所以，我们对每一个维度求偏导数，有：

$$\frac{{\partial \ln (l(\eta ))}}{{\partial {\theta _j}}} = \frac{{\partial \ln (l(\eta ))}}{{\partial \eta }}\frac{{\partial \eta }}{{\partial {\theta _j}}} = \sum\limits _{i = 1}^m {(T({y^i})} - \frac{{\partial a(\eta )}}{{\partial \eta }})x _j^i$$

如果指数分布为伯努利分布，我们将式子带入其中，能够轻易的关于\(\theta _j\)的偏导数，且与之前推导的结果一致，具体过程不详细说了。高斯分布同理。

机器学习 2.2 监督学习之逻辑回归

监督学习

本章继续研究监督学习的问题。

1 回归问题

我们随机制造一个样本空间\({x _2} \geqslant 2{x _1} + 1\)，对于对空间进行分类。具体例子如下：

对于这些样本，假定可以通过二维线性曲线进行分类。然后模拟出某条曲线，使样本得到一个很好的分类。

这个分类问题的样本的y值只有0或1(0-1分布或伯努利分布)，即表示在直线上面或直线下面。如果用之前的线性回归结果将非常糟糕。因为对于分类问题，线性回归得到的线性曲线，得到的\(y > 1\)或\(y < 0\)这些是没有意义的，肯定是错误的。因此引入了logistic函数，转换为Logistic regression，logistic函数如下：

$$h(z) = \frac{1}{{1 + {e^{ - z}}}}$$

函数曲线如下：

如果我们设置\(z = {\theta _1}{x _1} + {\theta _0} - {x _2}\)。因此，对于分类问题就可以转化为：对于\({x _2} \geqslant {\theta _1}{x _1} + {\theta _0}\)就可以转化为h(z)趋向于0的程度，对于\({x _2} \leqslant {\theta _1}{x _1} + {\theta _0}\)就可以转化为h(z)趋向于1的程度。我们简化为如下公式:

$$z = {\theta _1}{x _1} + {\theta _0}{x _0} + {\theta _2}{x _2} = {\theta ^T}x$$

其中,\({x _0} = 1\)，\({\theta _2} = - 1\)。

2 回归问题推导

我们求辅助函数的导数，以备后来使用：

$$h'(z) = (\frac{1}{{1 + {e^{ - z}}}})' = - \frac{{ - {e^{ - z}}}}{{{{(1 + {e^{ - z}})}^2}}} = \frac{1}{{1 + {e^{ - z}}}}(1 - 1 + {e^{ - z}}) = h(z)(1 - h(z))$$

接下来我们从公式角度对问题进行建模分析与推导。

这里我们将我们的函数\({h_\theta }(x)\)转化为如下形式:

$${h_\theta }(x) = \frac{1}{{1 + {e^{ - {\theta ^T}x}}}}$$

对于0-1分布，由于\({h_\theta }(x)\)只能为0或1。所以，在已知\(x\)和\({\theta}\)的情况下，我们可以认为有\({h_\theta }(x)\)的概率得到\(y=1\)。因此有如下公式：

$$p(y = 1|x,\theta ) = {h_\theta }(x)$$ $$p(y = 0|x,\theta ) = 1 - {h_\theta }(x)$$

综合两个公式有如下公式：

$$p(y|x,\theta ) = {h_\theta }{(x)^y}{(1 - {h _\theta }(x))^{(1 - y)}}$$

即在已经样本的情况下得到准确的\(y\)，即要使用极大释然方法，得到样本概率乘积的最大值,即得到下式的最大值：

$$L(\theta ) = \prod\limits _{i = 1}^m {{h _\theta }{{(x{}^i)}^{{y^i}}}{{(1 - {h _\theta }({x^i}))}^{(1 - {y^i})}}} $$

取自然对数，对其中一个维度\({\theta _j}\)求偏导数，如下：

$$\ell ({\theta _j}) = \frac{{\partial \ln (L(\theta ))}}{{\partial {\theta _j}}} = \frac{\partial } {{\partial {\theta _j}}}(\sum\limits _{i = 1}^m {({y^i}\ln ({h _\theta }(x{}^i)) + (1 - {y^i})} \ln (1 - {h _\theta }(x{}^i))))$$

我们使用一个辅助函数\(z = {\theta ^T}x\), 并且\(\ln x\)的倒数是\(\frac{1}{x}\)。
又由于复合函数有如下求导法则：

$$\frac{{df(g(x))}}{{dx}} = \frac{{df(z)}}{{dz}}g(x)'$$

其中\(z = g(x)\)。。(这里再次强调下标识维度，上表是样本编号)。所以有：

$$\ell ({\theta _j}) = \frac{{\partial L(\theta )}}{{\partial {\theta _j}}} = \frac{\partial }{{\partial {\theta _j}}}(\sum\limits _{i = 1}^m {({y^i}\ln ({h _\theta }(x{}^i)) + (1 - {y^i})} \ln (1 - {h _\theta }(x{}^i))))$$

这里设\(z = {\theta ^T}x\)，有：

$$\ell ({\theta _j}) = \sum\limits _{i = 1}^m {({y^i}\frac{1}{{h(z)}}h'(z) + (1 - {y^i})} \frac{1}{{1 - h(z)}}( - h'(z)))$$

将之前的倒数公式带入其中,有:

$$\ell ({\theta _j}) = \sum\limits _{i = 1}^m {(x _j^i(y - {h _\theta }({x^i})))} $$

因此这样就可以使用如下的式子学习样本值：

$${\theta _j} = {\theta _j} + \alpha \sum\limits _{i = 1}^m {(x _j^i(y - {h _\theta }({x^i})))} $$

3 代码

从这一章开始讲代码从matlab转化更为通用的python。

我们选取合适的步长，当\(J({\theta})\)小于合适的阈值(这里配置为4)，就认为迭代完成。分类代码如下:

import random
import math
import numpy as np
from matplotlib import pyplot as plt
k = 2
b = 1
x1_res = [-0.711, -0.2798, 0.4258, 0.4108, 0.0492, -0.842, -0.6036, -0.8394, 0.395, -0.8612, 0.806, 0.1476, 0.5464,
          0.3992,
          0.8326, 0.76, -0.7288, 0.0718, -0.4436, 0.8448, -0.017, -0.108, -0.4032, -0.9678, -0.4322, -0.9024, -0.0146,
          -0.5346, -0.3808, 0.0498, 0.6376, -0.3946, 0.5896, -0.9504, 0.2046, -0.0994, 0.468, -0.9236, 0.0996, 0.4366,
          -0.167, 0.9618, 0.5796, -0.383, 0.0254, -0.367, 0.0734, -0.8434, -0.9948, -0.0188]
x2_res = [-0.5108, -0.1254, 0.6166, 0.5296, 0.5758, -0.5972, 0.3994, 0.9808, 0.7134, -0.9414, -0.098, -0.2778, 0.667,
          0.1916, -0.0298, -0.5124, -0.509, -0.9754, 0.261, 0.4256, -0.7104, 0.8398, 0.7094, -0.5558, -0.799, 0.666,
          0.392,
          0.6688, 0.7832, 0.6372, -0.788, -0.1138, 0.0922, -0.6254, -0.1456, 0.3608, 0.4694, -0.8946, 0.9708, 0.8736,
          -0.5254, 0.0252, 0.6554, -0.2486, 0.9064, -0.35, -0.7724, 0.131, 0.4446, 0.2468]
          
def make_samples(x0_arr,x1_arr,x2_arr,y_arr):
  # classified by x_2 = 2*x_1 + 1 . (x_1,x_2) is the point of 2-D plane
  # y is 0 or 1, means 'x_2 >= 2*x_1 + 1 '  or 'x_2 < 2*x_1 + 1'
  for i in range(len(x1_res)):
    x0_arr.append(1)
    x1_arr.append(x1_res[i])
    x2_arr.append(x2_res[i])
    y_arr.append(0 if (x2_arr[i] >= k * x1_arr[i] + b ) else 1)

def h(z):
  return 1/(1 + math.exp(-1*z))

def h_xita(x0,x1,x2,xita):
  return h(transvection([x0,x1,x2], xita))

def transvection(a,b):
  return np.dot(a, b)

def J(x0_arr,x1_arr,x2_arr,y_arr,xita):
  sum = 0;
  length = len(x1_arr)
  for i in range(length):
    sum = sum + ( h_xita(x0_arr[i],x1_arr[i],x2_arr[i],xita) - y_arr[0] )**2
  return sum/2

def updateXita(x0_arr,x1_arr,x2_arr,y_arr,xita,a,update_x):
  sum = 0;
  length = len(x0_arr)
  for i in range(length):
    sum = sum + ( y_arr[i] - h_xita(x0_arr[i],x1_arr[i],x2_arr[i],xita) ) * update_x[i]
  return sum*a

def showPic(k,b,point_x1_arr,point_x2_arr):
  ## show line
  X = np.arange(-1,1,0.001)
  Y=[]
  for i in range(X.size):
    Y.append(k*X[i]+b)
  plt.plot(X, Y, 'b--', label="logistic")
  ## show points
  for i in range(len(point_x1_arr)):
    if(point_x2_arr[i]>2*point_x1_arr[i]+1):
      plt.plot(point_x1_arr[i],point_x2_arr[i],'b--',marker = 'x', color = 'g')
    else:
      plt.plot(point_x1_arr[i],point_x2_arr[i],'r--',marker = '+', color = 'r')
  plt.xlabel("x1")
  plt.ylabel("x2")
  plt.figure(figsize=(8, 4))
  plt.show()

if __name__== "__main__":
  x0_arr = []
  x1_arr = []
  x2_arr = []
  y_arr = []
  make_samples(x0_arr,x1_arr,x2_arr,y_arr)
  xita = [1,20,-1]
  a = 0.01
  limit = 4
  count = 0
  while 1:
    j = J(x0_arr, x1_arr, x2_arr, y_arr, xita)
    print("j=",j)
    if j < limit:
      break;
    xita[0] = xita[0] + updateXita(x0_arr, x1_arr, x2_arr, y_arr, xita, a, x0_arr)
    xita[1] = xita[1] + updateXita(x0_arr, x1_arr, x2_arr, y_arr, xita, a, x1_arr)
    count = count + 1
    print("cout=",count)
    print("xita0=",xita[0]," xita1=",xita[1])
  showPic(xita[1],xita[0],x1_arr,x2_arr)
  print(xita)

经过数论迭代分析，得到\( \theta _0 = 2.7841420917875297 \)，\( \theta _1 = 4.850051144306185 \)。

得到如下结果:

4 代码的改进

经过试验发现上面的算法的问题是如果选取合适的阈值。这里将算法修改为两次的\(J({\theta})\)的差值小于0.00001，这样就任务函数已经几近收敛,然后结束迭代。另外，如果\(J({\theta})\)的值某一次迭代增大了,说明算法异常或者迭代参数设置过大。修改后的代码如下:

  while 1:                                                                             
    j = J(x0_arr, x1_arr, x2_arr, y_arr, xita)                                         
    print("j=",j)                                                                      
    if j > j_last:                                                                     
      print("算法异常或选取迭代系数过大")                                                           
      break;                                                                           
    if j_last - j < 0.00001:                                                           
      print("out ",j_last - j)                                                         
      break;                                                                           
    j_last = j                                                                         
    xita[0] = xita[0] + updateXita(x0_arr, x1_arr, x2_arr, y_arr, xita, a, x0_arr)     
    xita[1] = xita[1] + updateXita(x0_arr, x1_arr, x2_arr, y_arr, xita, a, x1_arr)

得到如下结果:

5 牛顿法改进

这里说明另外一种快速求值的方法，牛顿法。

如上所示为牛顿法的简易过程。我们初始值为\(x _0\)，然后在\((x _0,f(x _0))\)做切线得到\(x _0\)，依次类推会接近于得到\(x _n\),使得\(f(x _n)\)接近于0。具体算法的证明请见文献2。

根据上面的算法，我们可以轻易的推导出

$${x _{i + 1}} = {x _i} - \frac{{f({x _i})}}{{f'({x _i})}}$$

联系之前的问题，我们可以转化快速得到\({\theta _j}\)使得\(\ell ({\theta _j}) = 0\)问题。（注: 之前的算法使找到最优下降方向使得\(J({\theta})\)最小，而牛顿法是直接得到最优的\({\theta _j}\)）

本问题的公式如下:

$$\ell ({\theta _j}) = \sum\limits _{i = 1}^m {(x _j^i(y - {h _\theta }({x^i})))} $$ $$\ell '({\theta _j}) = - \sum\limits _{i = 1}^m {{{(x _j^i)}^2}} $$ $${\theta _j}: = {\theta _j} - \frac{{\ell ({\theta _j})}}{{\ell '({\theta _j})}} = {\theta _j} - \frac{{\sum\limits _{i = 1}^m {(x _j^i(y - {h _\theta }({x^i})))} }}{{ - \sum\limits _{i = 1}^m {{{(x _j^i)}^2}} }}$$

代码如下:

import sys
import math
import numpy as np
from matplotlib import pyplot as plt
k = 2
b = 1
x1_res = [-0.711, -0.2798, 0.4258, 0.4108, 0.0492, -0.842, -0.6036, -0.8394, 0.395, -0.8612, 0.806, 0.1476, 0.5464,
          0.3992,
          0.8326, 0.76, -0.7288, 0.0718, -0.4436, 0.8448, -0.017, -0.108, -0.4032, -0.9678, -0.4322, -0.9024, -0.0146,
          -0.5346, -0.3808, 0.0498, 0.6376, -0.3946, 0.5896, -0.9504, 0.2046, -0.0994, 0.468, -0.9236, 0.0996, 0.4366,
          -0.167, 0.9618, 0.5796, -0.383, 0.0254, -0.367, 0.0734, -0.8434, -0.9948, -0.0188]
x2_res = [-0.5108, -0.1254, 0.6166, 0.5296, 0.5758, -0.5972, 0.3994, 0.9808, 0.7134, -0.9414, -0.098, -0.2778, 0.667,
          0.1916, -0.0298, -0.5124, -0.509, -0.9754, 0.261, 0.4256, -0.7104, 0.8398, 0.7094, -0.5558, -0.799, 0.666,
          0.392,
          0.6688, 0.7832, 0.6372, -0.788, -0.1138, 0.0922, -0.6254, -0.1456, 0.3608, 0.4694, -0.8946, 0.9708, 0.8736,
          -0.5254, 0.0252, 0.6554, -0.2486, 0.9064, -0.35, -0.7724, 0.131, 0.4446, 0.2468]

def make_samples(x0_arr,x1_arr,x2_arr,y_arr):
  # classified by x_2 = 2*x_1 + 1 . (x_1,x_2) is the point of 2-D plane
  # y is 0 or 1, means 'x_2 >= 2*x_1 + 1 '  or 'x_2 < 2*x_1 + 1'
  for i in range(len(x1_res)):
    x0_arr.append(1)
    x1_arr.append(x1_res[i])
    x2_arr.append(x2_res[i])
    y_arr.append(0 if (x2_arr[i] >= k * x1_arr[i] + b ) else 1)

def h(z):
  return 1/(1 + math.exp(-1*z))

def h_xita(x0,x1,x2,xita):
  return h(transvection([x0,x1,x2], xita))

def transvection(a,b):
  return np.dot(a, b)

def l(x0_arr,x1_arr,x2_arr,y_arr,xita,a,update_x):
  sum = 0;
  for i in range(len(x0_arr)):
    sum = sum + ( y_arr[i] - h_xita(x0_arr[i],x1_arr[i],x2_arr[i],xita) ) * update_x[i]
  return sum*a

def derl(update_x):
  sum = 0;
  for i in range(len(update_x)):
    sum = sum + update_x[i]*update_x[i];
  return -1 * sum

def showPic(k,b,point_x1_arr,point_x2_arr):
  ## show line
  X = np.arange(-1,1,0.001)
  Y=[]
  for i in range(X.size):
    Y.append(k*X[i]+b)
  plt.plot(X, Y, 'b--', label="logistic")
  ## show points
  for i in range(len(point_x1_arr)):
    if(point_x2_arr[i]>2*point_x1_arr[i]+1):
      plt.plot(point_x1_arr[i],point_x2_arr[i],'b--',marker = 'x', color = 'g')
    else:
      plt.plot(point_x1_arr[i],point_x2_arr[i],'r--',marker = '+', color = 'r')
  plt.xlabel("x1")
  plt.ylabel("x2")
  plt.figure(figsize=(8, 4))
  plt.show()

if __name__== "__main__":
  x0_arr = []
  x1_arr = []
  x2_arr = []
  y_arr = []
  make_samples(x0_arr,x1_arr,x2_arr,y_arr)
  xita = [1,20,-1]
  xita_bak = [10000,10000]
  count = 0
  while 1:
    xita[0] = xita[0] - l(x0_arr, x1_arr, x2_arr, y_arr, xita, 1,x0_arr)/derl(x0_arr)
    xita[1] = xita[1] - l(x0_arr, x1_arr, x2_arr, y_arr, xita, 1, x1_arr) / derl(x1_arr)
    if abs(l(x0_arr, x1_arr, x2_arr, y_arr, xita, 1,x0_arr)) <= 0.001 :
      print("xita = ",xita)
      break;
    xita_bak[1]=xita[1]
    xita_bak[0]=xita[0]
    count = count + 1
    print("count=", count, " xita = ", xita)

  showPic(xita[1],xita[0],x1_arr,x2_arr)
  print(xita)

得到\({\theta _0}\)和\({\theta _1}\)分别为
2.7762287194293407, 4.835356147838368。得到的曲线如下：

6 参考文献

  • 1 << cs229-notes1 >>
  • 2 << 计算机科学计算 >> 施吉林,张宏伟，金光日编

未来代码将会维护在https://github.com/zhengchenyu/mlearning

机器学习 2.1 监督学习之线性回归

监督学习

监督学习说明

如上图所示，监督学习实质就是: 在给定训练集合，使用某种学习算法得到学习函数h，能够较为准确的预测y。

1 线性回归

下面有一组房价关于房子面积以及卧室数目的样本数据。试着从这些样本数据中构造线性模型，以预测房价。

房间数目 房间大小 房价
2 88.0 1760288
2 88.0 1762136
… … …

1.1 梯度下降法

1.1.1 梯度下降法说明

假设房价与房间数目、房屋带下呈线性关系，有如下关系:

$${h_\theta }(x) = {\theta _0} + {\theta _1}{x _1} + {\theta _2}{x _2}$$

假设\({x_0} = 1\),则有如下公式:

$${h_\theta }(x) = {\theta _0}{x _0} + {\theta _1}{x _1} + {\theta _2}{x _2} = \sum\limits _{i = 0}^n{\theta _i}{x _i}$$

假设y为样本实际的值，因此有如下损失函数可以定义：

$$J(\theta ) = \frac{1}{2}\sum\limits _{i = 0}^m {{{({h _\theta }({x^i}) - {y^i})}^2}} $$

我们只要保证上面的公式趋于0就是合适的，因此这个问题就转化为得到合适的使得\(J(\theta )\)最小。实际也就转化为利用最小二乘法进行回归分析。

做这样一个想想，对于二维的情况。\(J(\theta )\)是一个关于的二维函数的话，他一定类似于锅状的凸函数。根据梯度下降算法(参考<<最优化方法>>)，沿着某一个方向的负梯度方向就是在该方向上下降最大的方向。因此，我们可以选取各个维度的负梯度作为下降的方向，为了方便各个梯度方向选取同样的步长。因此，我们对任意维度采用如下公式进行递归(为步长)：

$${\theta _j}: = {\theta _j} - \alpha \frac{{\partial J(\theta )}}{{\partial {\theta _j}}} = {\theta _j} - \alpha \sum\limits _{i = 0}^m{({h _\theta }({x^i}) - {y^i}) {x _j}}$$

我们简化，我们仅仅对房间大小与房价的一维关系进行回归分析。

注：对于每一次迭代运算都需要对所有全量信息进行计算。这样带来了大量的计算。工程上可以采用采样的方式计算新的，这里数据较少，就不用该方法了。

1.1.2 梯度下降法案例

我们有一组数据，为房间大小与房价的关系。我们的目标是通过使用梯度下降算法，模拟得到满足学习函数，以预测房价。可以看如下点图:

房价点图

根据上一节的公式编写如下代码进行回归分析:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
% 以下为主算法流程
% y = xita0 +xita1*x
mfilename='/Users/zcy/Desktop/mlearning/prices'
[v,x,y]=textread(mfilename,'%n%n%n');
 
limit = 1000;
y=y./10000;
% step
a = 0.0000001;
 
xita0 = 0.01;
xita1 = 1.6;
err = JFun(x,y,xita0,xita1);
disp(['err = ',num2str(err)])
count = 0; 
while err>limit 
    xita0 = xita0 - a*gradient(x,y,xita0,xita1,true);
    xita1 = xita1 - a*gradient(x,y,xita0,xita1,false)    
    err = JFun(x,y,xita0,xita1);
    count = count + 1;
    disp(['err = ',num2str(err),', count = ',num2str(count)])
end
 
scatter(x,y,'k')
hold on
x1 = 70:0.01:100;
y1 = xita0+x1.*xita1;
plot(x1,y1)

JFun.m文件流程:

1
2
3
4
5
function [s]  = JFun(x,y,xita0,xita1)
s = 0;
for i = 1:size(x)
   s = s + 0.5*((xita0*1+xita1*x(i)-y(i))^2);
end

gradient.m文件:

1
2
3
4
5
6
7
8
9
10
11
function [s]  = gradient(x,y,xita0,xita1,b)
s = 0;
if b 
    for i = 1:size(x)
        s = s + (xita0*1+xita1*x(i)-y(i))*1;
    end
else 
    for i = 1:size(x)
        s = s + (xita0*1+xita1*x(i)-y(i))*x(i);
    end
end

这里为了示意，代码仅考虑两维，实际不建议代码这样写。另外，选取合适的limit值和步长需要经过打印调试。选取不当，容易造成函数无法收敛。

回归得到如下曲线：

回归曲线1

修改limit为500之后，回归更加准确。

回归曲线2

1.2 线性回归的解析解

略，详见斯坦福大学机器学习讲义。

1.3 线性回归的概率解释

本节从概率的角度进行线性回归分析。对我们的任意一组变量有如下公式：

$${y^i} = {\theta ^T}{x^i} + {\varepsilon ^i}$$

其中\(\varepsilon \)为误差，假设我们的误差服从正太分布
\({\varepsilon ^i} \sim {\rm N}(0,{\sigma ^2})\)。即有如下公式：

$$p({\varepsilon ^i}) = \frac{1}{{\sqrt {2\pi } \sigma }}{e^{( - \frac{{{{({\varepsilon ^i})}^2}}}{{2{\sigma ^2}}})}}$$

代入上面的公式有：

$$p({y^i}|{x^i},\theta ) = \frac{1}{{\sqrt {2\pi } \sigma }}{e^{( - \frac{{{{({y^i} - {\theta ^T}{x^i})}^2}}}{{2{\sigma ^2}}})}}$$

该公式的含义为在已知\(x ^i \)和\(\theta \)的情况下，得到准确的\(y ^i \)概率。我们可以利用最大似然估计，最大释然估计实质就是对所有采样值调整参数得到最准确y值的方法。(详细请看概率论)。有如下公式:

$$L(\theta ) = \prod\limits _{i = 1}^m {p({y^i}|{x^i},\theta )} = \prod\limits _{i = 1}^m {\frac{1}{{\sqrt {2\pi } \sigma }}{e^{( - \frac{{{{({y^i} - {\theta ^T}{x^i})}^2}}}{{2{\sigma ^2}}})}}} $$

上面的公式是对各个数据采样得到准确y的概率的乘积，问题也就转化为调整使的上式最大，进一步转化为求的上式的倒数为0的情况。对上式去log，如下:

$$[\log (L(\theta )) = \log (\prod\limits _{i = 1}^m {\frac{1}{{\sqrt {2\pi } \sigma }}{e^{( - \frac{{{{({y^i} - {\theta ^T}{x^i})}^2}}}{{2{\sigma ^2}}})}}} )$$ $$\log (L(\theta )) = \sum\limits _{i = 1}^m {\log \frac{1}{{\sqrt {2\pi } \sigma }}{e^{( - \frac{{{{({y^i} - {\theta ^T}{x^i})}^2}}}{{2{\sigma ^2}}})}}} = m\log \frac{1}{{\sqrt {2\pi } \sigma }} - \frac{1}{{{\sigma ^2}}} \cdot \frac{1}{2}\sum\limits _{i = 1}^m {{{({y^i} - {\theta ^T}{x^i})}^2}} $$

看上式最后一项，前面部分为常数，因此问题有转化为求最小值的问题了。

因此从概率角度上考虑误差的情况下，我们前面基于最小二乘的算法也是正确的。

1.4 局部加权回归(Loess)

1.4.1 局部加权回归说明

有一组非线性去先，譬如\(y = {x^2}\)。如果使用之前的方法拟合(即假设其为y=kx+b型曲线)，必然不会得到理想的结果。对于这样的问题，我们想求出\(f(20)\),加入我们只对\(x\)在20附近的样本进行一维拟合，我们就可以得到一个精确值。因此，我们引入如下公式：

$$J(\theta ) = \frac{1}{2}\sum\limits_ {i = 0}^m {{\omega ^i}{{({h_ \theta }({x^i}) - {y^i})}^2}} $$ $${\omega ^i} = {e^{ - \frac{{{{({x^i} - x)}^2}}}{{2{\tau ^2}}}}}$$

其中，在原来的价值函数中引入权重\(\omega\)。\(\tau\)是波长，该值越小，越趋于向20附近的值进行回归计算。对于求\(f(20)\)的情况，在的公式中\(x\)恒为20,可以知道该值在\(x ^i\)趋近于20的时候趋近于1，趋近于\( \pm \infty \)的时候趋近于0。也就是相当于仅仅是在20附近进行拟合。由于新引入的权重函数与\(\theta\)无关。所以得到如下公式:

$${\theta _j}: = {\theta _j} - \alpha \frac{{\partial J(\theta )}}{{\partial {\theta _j}}} = {\theta _j} - \alpha \sum\limits _{i = 0}^m {{\omega ^i}({h _\theta }({x^i}) - {y^i})} {x _j}$$

我们这里是一个一维的问题，可以得到如下公式:

$$J(\theta ) = \frac{1}{2}\sum\limits _{i = 0}^m {{\omega ^i}{{({\theta _1}{x^i} + {\theta _0} - {y^i})}^2}} $$ $${\omega ^i} = {e^{ - \frac{{{{({x^i} - 20)}^2}}}{{2{\tau ^2}}}}}$$

1.4.2 局部加权回归案例

主流程代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
% y = xita0 +xita1*x
x = -100:1:100;
y = x.*x;
% x= 20, y = 40x - 400
x0 = 20;       
tao = 5;
a = 0.0003;
 
xita0 = -1;
xita1 = 1;
limit = 7850;
err = JFun_Loess(x,y,xita0,xita1,x0,tao);
count = 0; 
while err>limit 
    xita0 = xita0 - a*gradient_Loess(x,y,xita0,xita1,true,x0,tao);
    xita1 = xita1 - a*gradient_Loess(x,y,xita0,xita1,false,x0,tao);    
    err = JFun_Loess(x,y,xita0,xita1,x0,tao);
    count = count + 1;
end
 
scatter(x,y,'k')
hold on
x1 = -100:1:100;
y1 = xita0+x1.*xita1;
disp(['xita0 = ',num2str(xita0),', xita1 = ',num2str(xita1)])
plot(x1,y1)

JFun_loess:

1
2
3
4
5
function [s]  = JFun_Loess(x,y,xita0,xita1,x0,tao)
s = 0;
for i = 1:size(x')
    s = s + 0.5*((xita0*1+xita1*x(i)-y(i))^2)*exp(-0.5*(x(i)-x0)*(x(i)-x0)/tao/tao);
end

gradient_Loess:

1
2
3
4
5
6
7
8
9
10
11
function [s]  = gradient_Loess(x,y,xita0,xita1,b,x0,tao)
s = 0;
if b 
    for i = 1:size(x')
        s = s + (xita0*1+xita1*x(i)-y(i))*1*exp(-0.5*(x(i)-x0)*(x(i)-x0)/tao/tao);
    end
else
    for i = 1:size(x')
        s = s + (xita0*1+xita1*x(i)-y(i))*x(i)*exp(-0.5*(x(i)-x0)*(x(i)-x0)/tao/tao);
    end
end

在x=20处分析，理论值为y = 40x – 400。经过模拟得到的值为y=39.6825x-368.253。
下图为得到的结果，实际上已经非常接近这个二次曲线在20这个位置的切线了。

局部加权分析结果图

2 回归分析相关概念

另外引申一个概念:

之前我们把房价与房屋大小认为是一个一维曲线。如上图，可以看出很多点被反映到曲线中。我们可以称其为欠拟合过程。与之对应的，如果使用高阶函数进行拟合，即假如样本有5个点，我们可以通过一个4阶函数的曲线进行完整拟合，但是这样的曲线往往并不是一个良好的房价与房屋大小的反映，这被称为过拟合。(该部分内容可以想见矩阵与数值分析)

机器学习 1 绪论

机器学习的几种方法

  • 监督学习

    根据训练资料中建立一个模式，并依此模式推测新的实例。一般来讲监督学习的任务一般是分类和回归分析，往往需要人工预先准备好范例(训练资料)。

  • 非监督学习

不同于监督学习，非监督学习仅对网络提供输入范例，会从这些范例中找到潜在类别规则，主要目的是对原始资料进行分类。

  • 强化学习

    强化学习强调如何基于环境而行动，以取得最大的预期利益。其灵感来源于心理学中的行为主义理论，即有机体如何在环境给予的奖励或惩罚的刺激下，逐步形成对刺激的预期，产生能获得最大利益的习惯性行为。

一些规范

n表示样本的维度，下标。m 表示样本数目，上标。