花了一些时间刷了一遍李航博士的《统计学习方法》,于是就有了一篇课后习题解答。由于自己的智商又低数学又差,很多难度较大的习题实在没有能力完成了。

第1章 统计学习方法概论

  • 1.1 说明伯努利模型的极大似然估计以及贝叶斯估计中的统计学习方法三要素。伯努利模型是定义在取值为0与1的随机变昼上的概率分布。假设观测到伯努利模型n次独立的数据生成结果,其中k次的结果为1,这时可以用极大似然估计或贝叶斯估计来估计结果为1的概率。

伯努利模型的极大似然估计:

以定义概率为 。得出似然函数:

需要使用微分法求最值,将方程两边同时对取微分:

其解为 ,显然

伯努利模型的贝叶斯估计:

定义概率为 之间的取值是等可能的,因此先验概率密度函数 。得出似然函数:

根据似然函数和先验概率密度函数,可以求解 的条件概率密度函数:

最终 的期望为:

统计学习方法三要素:

  极大似然估计 贝叶斯估计
模型
策略 最大化似然函数 求参数期望
算法
  • 1.2 通过经验风险最小化推导极大似然估计。证明模型是条件概率分布,当损失函数是对数损失函数时,经验风险最小化等价于极大似然估计。

极大似然估计的似然函数为:

两边取对数:

因此,当损失函数是对数损失函数时,经验风险最小化等价于极大似然估计。

第2章 感知机

  • 2.1 Minsky 与 Papert 指出:感知机因为是线性模型,所以不能表示复杂的函数,如异或 (XOR)。验证感知机为什么不能表示异或。

对于异或函数XOR,全部的输入与对应的输出如下:

1 1 -1
1 -1 1
-1 1 1
-1 -1 -1

如果我们将异或函数的样本展示出来应该是这个样子的:

我们显然无法使用一条直线将两类样本划分,异或问题是线性不可分的。

  • 2.2 模仿例题 2.1,构建从训练数据求解感知机模型的例子。

例题 2.1 的数据集如下:

x = [3 3; 4 3; 1 1];
y = [1; 1; -1];

感知机的模型为:

function y = hypothesis(x, w, b)
	y = ((x * w + b) >= 0) * 2 - 1;
end

感知机的训练过程为:

(1)选取初值

(2)从训练集中选取数据

(3)如果

(4)转至(2),直到训练集中没有误分类点

function [w b] = traing(x, y, eta)
	sample_count = size(x, 1);
	feature_count = size(x, 2);
	w = zeros(feature_count, 1);
	b = 0;
	miss = sum(hypothesis(x, w, b) != y);
	while (miss > 0)
		i = ceil(rand() * sample_count);
		xi = x(i, :);
		yi = y(i);
		if ((xi * w + b) * yi <= 0)
			w = w + eta * yi * transpose(xi);
			b = b + eta * yi;
		end
		miss = sum(hypothesis(x, w, b) != y);
	end
end

利用数据训练感知机:

[w b] = traing(x, y, 0.1)
hypothesis(x, w, b)

感知机训练结果如下:

w =

   0.20000
   0.10000

b = -0.50000
ans =

   1
   1
  -1
  • 2.3 证明以下定理:样本集线性可分的充分必要条件是正实例点所构成的凸壳与负实例点所构成的凸壳互不相交。

第3章 k邻近算法

  • 3.1 参照图 3.1,在 二维空间中给出实例点,画出 k 为 1 和 2 时的 K 近邻法构成的空间划分,并对其进行比较,体会 K 值选择与模型复杂度及预测准确率的关系。

  • 3.2 利用例题 3.2 构造的屈树求点 x=(3,4.5) 的最近邻点。

    1. 首先找到包含(3,4.5)的叶节点(4,7),将叶节点作为“当前最近点”,距离约2.69;
    2. 回退到(4,7)的父节点(5,4),距离(3,4.5)约2.06,将(5,4)作为“当前最近点”。以(3,4.5)为圆心,到“当前最近点”(5,4)距离为半径的圆显然和(5,4)的另一个子节点(2,3)区域相交,因此移动到(2,3);
    3. 移动到(2,3)后发现,距离(3,4.5)约1.81,因此将(2,3)作为“当前最近点”,由于(2,3)是叶节点,因此直接回退;
    4. 回到(5,4)的根节点(7,2),到(7,2)距离大于到“当前最近点”距离,同时(3,4.5)和“当前最近点”距离构成的圆和根节点的另一个子节点的区域不想交,所以搜索结束,得到最近点(3,4.5)。

  • 3.3 参照算法 3.3,写出输出为 x 的 K 近邻的算法。

在寻找最近邻节点的时候需要维护一个”当前最近点“,而寻找 K 近邻的时候,就需要维护一个”当前 K 近邻点集“。首先定义一个”当前 K 近邻点集“插入新点操作:如果”当前 K 近邻点集“元素数量小于K,那么直接将新点插入集合;如果”当前 K 近邻点集“元素数量等于K,那么将新节点替换原来集合中最远的节点。

(1)在 kd 树中找出包含目标点 x 的叶结点:从根结点出发,递归地向下访问树。若目标点 x 当前维的坐标小于切分点的坐标,则移动到左子结点,否则移动到右子结点,直到子结点为叶结点为止;

(2)如果”当前 K 近邻点集“元素数量小于K或者叶节点距离小于”当前 K 近邻点集“中最远点距离,那么将叶节点插入”当前 K 近邻点集“;

(3)递归地向上回退,在每个结点进行以下操作:

​ (a)如果”当前 K 近邻点集“元素数量小于K或者当前节点距离小于”当前 K 近邻点集“中最远点距离,那么将该节点插入”当前 K 近邻点集“,

​ (b)检查另一子结点对应的区域是否与以目标点为球心、以目标点与于”当前 K 近邻点集“中最远点间的距离为半径的超球体相交。如果相交,可能在另一个子结点对应的区域内存在距目标点更近的点,移动到另一个子结点 . 接着,递归地进行最近邻搜索;如果不相交,向上回退;

(4)当回退到根结点时,搜索结束,最后的”当前 K 近邻点集“即为 x 的 K 近邻点集。

KD树Octave实现如下:

count = 100;
dimen = 100;
test = 100;

function ret = distance(a, b)
	ret = sqrt(sum((a.-b).^2,2));
end

function ret = find_k_nearest_naive(x, p, k)
	[val, r] = sort(distance(x,p));
	ret = r(1:k);
end

function ret = kdtree_build_recur(x, r, d)
	count = length(r);
	dimen = size(x, 2);
	if (count == 1)
		ret = struct('point', r(1), 'dimen', d);
	else
		mid = ceil(count / 2);
		ret = struct('point', r(mid), 'dimen', d);
		d = mod(d,dimen)+1;
		% Build left sub tree
		if (mid > 1)
			left = r(1:mid-1);
			leftxd = x(left,d);
			[val, leftrr] = sort(leftxd);
			leftr = left(leftrr);
			ret.left = kdtree_build_recur(x, leftr, d);
		end
		% Build right sub tree
		if (count > mid)
			right = r(mid+1:count);
			rightxd = x(right,d);
			[val, rightrr] = sort(rightxd);
			rightr = right(rightrr);
			ret.right = kdtree_build_recur(x, rightr, d);
		end
	end
end

function ret = kdtree_build(x)
	[val, r] = sort(x(:,1));
	ret = struct('data',x,'root', kdtree_build_recur(x,r,1));
end

function ret = kdtree_cand_fartest(x, p, cand)
	[val, index] = max(distance(x, p)(cand));
	ret = cand(index);
end

function ret = kdtree_cand_insert(x, p, cand, k, point)
	if (length(cand) < k)
		ret = [cand; point];
	else
		fartest = kdtree_cand_fartest(x, p, cand);
		cand(find(cand == fartest)) = point;
		ret = cand;
	end 
end

function ret = kdtree_find_recur(x, node, p, ret, k)
	point = node.point;
	d = node.dimen;
	if (x(point,d) > p(d))
		% Search in left sub tree
		if (isfield(node, 'left'))
			ret = kdtree_find_recur(x, node.left, p, ret, k);
		end
		% Add current point if neccessary
		fartest = kdtree_cand_fartest(x, p, ret);
		if (length(ret) < k || distance(x(point,:),p) < distance(x(fartest,:),p))
			ret = kdtree_cand_insert(x, p, ret, k, point);
		end
		% Search in right sub tree if neccessary
		fartest = kdtree_cand_fartest(x, p, ret);
		radius = distance(x(fartest,:),p);
		if (isfield(node, 'right') && (length(ret) < k || p(d) + radius > x(point,d)))
			ret = kdtree_find_recur(x, node.right, p, ret, k);
		end
	else
		% Search in right sub tree
		if (isfield(node, 'right'))
			ret = kdtree_find_recur(x, node.right, p, ret, k);
		end
		% Add current point if neccessary
		fartest = kdtree_cand_fartest(x, p, ret);
		if (length(ret) < k || distance(x(point,:),p) < distance(x(fartest,:),p))
			ret = kdtree_cand_insert(x, p, ret, k, point);
		end
		% Search in left sub tree if neccessary
		fartest = kdtree_cand_fartest(x, p, ret);
		radius = distance(x(fartest,:),p);
		if (isfield(node, 'left') && (length(ret) < k || p(d) - radius <= x(point,d)))
			ret = kdtree_find_recur(x, node.left, p, ret, k);
		end
	end
end

function ret = kdtree_find(tree, p, k)
	x = tree.data;
	root = tree.root;
	ret = kdtree_find_recur(x, root, p, [], k);
end

correct = 0;
for i = 1:test
	x = rand(count, dimen);
	y = rand(1, dimen);
	k = ceil(rand(1,1)*count);
	ret = kdtree_build(x);
	a = sort(kdtree_find(ret, y, k));
	b = sort(find_k_nearest_naive(x, y, k));
	correct += sum(a == b) == k;
end
printf('correct:%d/%d\n', correct, test);

第4章 朴素贝叶斯法

  • 4.1 用极大似然估计法推出朴素贝叶斯法中的概率估计公式(4.8)及公式 (4.9)。

证明公式(4.8):

由于朴素贝叶斯法假设Y是定义在输出空间上的随机变量,因此可以定义概率为 。令 ,得出似然函数:

需要使用微分法求最值,将方程两边同时对取微分:

其解为 ,显然

证明公式(4.9):

,令 ,似然函数:

同样使用微分法求最值可以得到:

  • 4.2 用贝叶斯估计法推出朴素贝叶斯法中的慨率估计公式(4.10)及公式(4.11)。

第5章 决策树

  • 5.1 根据表 5.1 所给的训练数据集,利用信息增益比(C4.5 算法)生成决策树。

数值比字符串更容易处理,所以在学习过程中需要将字符串映射到数值:

特征(编号) 0 1 2
年龄(1) 青年 中年 老年
有工作(2)  
有自己的房子(3)  
信贷情况(4) 一般 非常好
类别  

训练数据集:

x = [
	0 0 0 0;
	0 0 0 1;
	0 1 0 1;
	0 1 1 0;
	0 0 0 0;
	1 0 0 0;
	1 0 0 1;
	1 1 1 1;
	1 0 1 2;
	1 0 1 2;
	2 0 1 2;
	2 0 1 1;
	2 1 0 1;
	2 1 0 2;
	2 0 0 0;
];
y = [0; 0; 1; 1; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 0];

计算数据集D的经验熵方法为:

function ans = entropy(y)
	ans = 0;
	output_values = unique(y);
	for i = (1:length(output_values))
		p = sum(y == output_values(i)) / length(y);
		ans -= p * log2(p);
	end
end 

计算特征A对数据集D的经验条件熵H(D|A):

function ans = entropycond(x, y, feature_index)
	ans = 0;
	sample_count = length(y);
	feature_values = unique(x(:,feature_index));
	for i = (1:length(feature_values))
		feature_value = feature_values(i);
		subset = find(x(:,feature_index) == feature_value);
		ans += entropy(y(subset)) * length(subset) / sample_count;
	end
end

计算信息增益比:

function ans = entropyincr(x, y, feature_index)
	ans = 1 - entropycond(x, y, feature_index) / entropy(y);
end

C4.5 生成算法如下:

(1)如果 D 中所有实例属于同一类,则置 T 为单结点树,井将 作为该结点的类,返回 T;

(2)如果 A= ,则登 T 为单结点树,并将 D 中实例数最大的类 C食作为该结点的类,返回 T;

(3)否则,按式 (5.10)计算 A 中各特征对 D 的信息增益比,选择信息增益比最大的特征

(4)如果 的信息增益比小千阔值 , 则置 T 为单结点树,井将 D 中实例数最大的类 作为该结点的类,返回 T;

(5)否则,对 的每一可能值, 依 将 D 分割为子集若干非空 ,将 中实例数最大的类作为标记,构建子结点,由结点及其子结点构成树 T, 返回 T;

(6)对结点 i,以 为训练集,以为特征集,递归地调用步(1)~步(5),得到子树 ,返回

function ans = generate(x, y, feature_indices, epsilon)
	if (length(feature_indices) == 0)
		ans = struct('leaf', true, 'output', mode(y));
	elseif (length(unique(y)) == 1)
		ans = struct('leaf', true, 'output', y(1));
	else 
		[rate, index] = max(arrayfun(@(i) entropyincr(x, y, i), feature_indices))
		if (rate < epsilon)
			ans = struct('leaf', true, 'output', mode(y));
		else
			feature_values = unique(x(:,index))
			feature_indices(find(feature_indices == index)) = [];
			child = struct();
			for j = (1:length(feature_values))
				subset = find(x(:,index) == feature_values(j));
				child.(mat2str(feature_values(j))) = generate(x(subset,:), y(subset), feature_indices, epsilon);
			end
			ans = struct('leaf', false, 'feature', index, 'child', child);
		end
	end
end

训练模型:

generate(x, y, (1:1:4), 0)

得到结果:

ans =
  scalar structure containing the fields:
    leaf = 0
    feature =  3
    child =
      scalar structure containing the fields:
        0 =
          scalar structure containing the fields:
            leaf = 0
            feature =  2
            child =
              scalar structure containing the fields:
                0 =
                  scalar structure containing the fields:
                    leaf =  1
                    output = 0
                1 =
                  scalar structure containing the fields:
                    leaf =  1
                    output =  1
        1 =
          scalar structure containing the fields:
            leaf =  1
            output =  1

  • 5.2 已知如表 5.2 所示的训练数据,试用平方误差损失准则生成一个二叉回归树。

训练数据集:

x = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10];
y = [4.50; 4.75; 4.91; 5.34; 5.80; 7.05; 7.90; 8.23; 8.70; 9.00];

寻找最优切分变量j和最优切分点s的方法为:

其中,

function [j, s, minval, c1, c2] = divide(x, y, feature_count)
	cost = zeros(feature_count, length(x)-1);
	for i = 1:feature_count
		for k = 1:(length(x)-1)
			value = x(k,i);
			y1 = y(find(x(:,i) <= value));
			c1 = mean(y1);
			y2 = y(find(x(:,i) > value));
			c2 = mean(y2);
			cost(i,k) = sum((y1 - c1).^2) + sum((y2 - c2).^2);
		end
	end
	[minval, row] = min(min(cost,[],2));
	[minval, col] = min(min(cost,[],1));
	j = row;
	s = x(col,j);
	c1 = mean(y(find(x(:,j) <= s)));
	c2 = mean(y(find(x(:,j) > s)));
end

最小二乘回归树生成算法过程如下:

在训练数据集所在的输入空间中,递归地将每个区域划分为两个子区域并决定每个子区域上的输出值,构建二叉决策树:

(1)选择最优切分变量 J 与切分点 s;

(2)用选定的对 (J,S) 划分区域并决定相应的输出值;

(3)继续对两个子区域调用步骤(1)(2)直至满足停止条件;

(4)将输入空间划分为 M 个区域,生成决策树:

在具体实现中,停止条件可以是区域不可以再分,或者小于某个阈值:

function tree = fittree(x, y, feature_count, epsilon)
	[j, s, minval, c1, c2] = divide(x, y, feature_count);
	tree = struct('feature', j, 'value', s);
	if (minval < epsilon || length(find(x(:,j) <= s)) <= 1)
		tree.left = c1;
	else 
		tree.left = fittree(x(find(x(:,j) <= s)), y(find(x(:,j) <= s)), feature_count, epsilon);
	end
	if (minval < epsilon || length(find(x(:,j) > s)) <= 1)
		tree.right = c2;
	else 
		tree.right = fittree(x(find(x(:,j) > s)), y(find(x(:,j) > s)), feature_count, epsilon);
	end
end

如果将阈值设为0.1:tree = fittree(x, y, 1, 0.1)那么可以得到:

tree =
  scalar structure containing the fields:
    feature =  1
    value =  5
    left =
      scalar structure containing the fields:
        feature =  1
        value =  3
        left =  4.7200
        right =  5.5700
    right =
      scalar structure containing the fields:
        feature =  1
        value =  7
        left =
          scalar structure containing the fields:
            feature =  1
            value =  6
            left =  7.0500
            right =  7.9000
        right =
          scalar structure containing the fields:
            feature =  1
            value =  8
            left =  8.2300
            right =  8.8500

也就是

  • 5.3 证明 CART 剪枝算法中,当确定时,存在唯一的最小子树 使损失函数 最小。

  • 5.4 证明 CART 剪枝算法中求出的子树序列分别是区间 的最优子树 ,这里

第6章 逻辑斯谛回归与最大熵模型

  • 6.1 确认逻辑斯谛分布属于指数分布族。

根据维基百科对于指数分布族的定义,逻辑斯谛分布不属于指数分布族。

In general, distributions that result from a finite or infinite mixture of other distributions, e.g. mixture model densities and compound probability distributions, are not exponential families. Examples are typical Gaussian mixture models as well as many heavy-tailed distributions that result from compounding (i.e. infinitely mixing) a distribution with a prior distribution over one of its parameters, e.g. the Student’s t-distribution (compounding a normal distribution over a gamma-distributed precision prior), and the beta-binomial and Dirichlet-multinomial distributions. Other examples of distributions that are not exponential families are the F-distribution, Cauchy distribution, hypergeometric distribution and logistic distribution.

  • 6.2 写出逻辑斯谛回归模型学习的梯度下降算法。

对于逻辑斯谛模型:

对数似然函数:

似然函数偏导:

梯度函数为:

根据梯度下降算法(针对求最大值做了修改):

(1)取初始值 , 置

(2)计算

(3)计算梯度 ,当 时,停止迭代,令 ;否则 ,求 ,使

(4)置 ,计算

时,停止迭代,令

(5)否则,置 ,转(3)。

编写程序,对数据进行学习:

% Answers for reference: w = [-16.38; 0.1483; 0.1589]

x = load('data6-2-x.txt');
y = load('data6-2-y.txt');

x = [ones(length(x),1) x];

function ans = like(x,y,w)
	sum(x*w.*y-log(1+exp(x*w)));
end

function ans = grad(x,y,w)
	transpose(sum(x.*y-(x.*exp(x*w))./(1+exp(x*w))));
end

function w0 = logistic_gradient_descent(x,y,trial,epsilon)
	w0 = zeros(size(x,2),1);
	l0 = like(x,y,w0);
	while (true)
		g = grad(x,y,w0);
		if (sum(abs(g))<epsilon)
			break;
		end
		w1 = w0;
		l1 = l0;
		temp = ones(size(x,2),1)/(10^trial);
		lambda = temp;
		for i = 1:length(temp)
			for j = 1:trial
				temp(i) = temp(i)*10;
				if (like(x,y,w0+g.*temp) > l1)
					lambda = temp;
					w1 = w0+g.*temp;
					l1 = like(x,y,w1);
				end
			end
		end
		if (sum(abs(w0-w1))<epsilon)
			break;
		end
		if (abs(l0-l1)<epsilon)
			break;
		end
		w0 = w1;
		l0 = l1;
	end
end

w = logistic_gradient_descent(x,y,5,0.00001)

得到结果:

w =

  -16.19103
    0.14732
    0.15672
  • 6.3 写出最大熵模型学习的DFP算法。(关于一般的DFP算法参见附录B)

假设有一组各种天气情况下是否打球的数据(揭开机器学习的面纱:最大熵模型100行代码实现[Python版]):

play outlook temperature humidity windy
no sunny hot high FALSE
no sunny hot high TRUE
yes overcast hot high FALSE
yes rainy mild high FALSE
yes rainy cool normal FALSE
no rainy cool normal TRUE
yes overcast cool normal TRUE
no sunny mild high FALSE
yes sunny cool normal FALSE
yes rainy mild normal FALSE
yes sunny mild normal TRUE
yes overcast mild high TRUE
yes overcast hot normal FALSE
no rainy mild high TRUE

数值比字符串更容易处理,所以在学习过程中需要将字符串映射到数值:

特征(编号) 0 1 2
outlook(1) sunny overcast rainy
temperature(2) cool mild hot
humidity(3) normal high  
windy(4) FALSE TRUE  
play no yes  

得到训练数据:

X = [0 2 1 0; 0 2 1 1; 1 2 1 0; 2 1 1 0; 2 0 0 0; 2 0 0 1; 1 0 0 1; 
	 0 1 1 0; 0 0 0 0; 2 1 0 0; 0 1 0 1; 1 1 1 1; 1 2 0 0; 2 1 1 1];
Y = [0 0 1 1 1 0 1 0 1 1 1 1 1 0];

经验分布的计算方法为:

function p = px(X, x)
	p = sum(prod(X==x, 2))/length(X);
end

function p = pxy(X, Y, x, y)
	X = [X transpose(Y)];
	x = [x y];
	p = px(X, x);
end

特征函数为:

function f = feature(X, Y, i, x, y)
	f = sum(X(i,:)!=x)==0 && Y(i)==y;
end

最大熵模型为

function [p output] = maxent(X, Y, w, x)
	count = length(X);
	output = unique(Y);
	z = 0;
	p = [];
	for y = output
		indicators = zeros(1, count);
		for i = (1:count)
			indicators(i) = feature(X, Y, i, x, y);
		end
		p = [p exp(indicators*w)];
		z += exp(indicators*w);
	end
	p = p/z;
end

对数似然函数(也就是优化的目标函数)为

function l = like(X, Y, w)
	count = length(w);
	in = unique(X, 'rows');
	out = unique(Y);
	l = 0;
	for x = transpose(in)
		for y = out
			for i = (1:count)
				l += pxy(X,Y,transpose(x),y)*w(i)*feature(X,Y,i,transpose(x),y);
			end
		end
	end
	for x = transpose(in)
		z = 0;
		for y = out
			indicators = zeros(1, count);
			for i = (1:count)
				indicators(i) = feature(X, Y, i, transpose(x), y);
			end
			z += indicators*w;
		end
		l -= px(X,transpose(x))*log(z);
	end
end

似然函数的梯度为:

function g = grad(X, Y, w)
	count = length(w);
	g = zeros(count, 1);
	in = unique(X, 'rows');
	out = unique(Y);
	for i = (1:count)
		for x = transpose(in)
			p = maxent(X,Y,w,transpose(x));
			for j = (1:length(out))
				y = out(j);
				g(i) -= px(X,transpose(x))*p(j)*feature(X,Y,i,transpose(x),y);
				g(i) += pxy(X,Y,transpose(x),y)*feature(X,Y,i,transpose(x),y);
			end
		end
	end
end

DFP算法:

(1)选定起始点,取为正定对称矩阵,精度要求

(2)计算。若,则停止计算,得近似解;否则转(3)

(3)置

(4)一维搜索:求使得

(5)置

(6)计算,若,则停止计算,得近似解;否则,计算

(7)置,转(3)

最优化算法通常是最小化目标函数,应用到最大化的时候需要修改算法或者目标函数:

function w = train(X, Y, e, trial=6)
	count = length(X);
	% (1) 
	w = rand(count,1);
	G = eye(count);
	% (2)
	g = grad(X,Y,w);
	if (sum(g) < e)
		return;
	end
	while (true)
		% (3)
		p = G*g;
		% (4)
		nextws = w+p*10.^((1:trial)-trial/2);
		likes = [];
		for nextw = nextws 
			likes = [likes like(X,Y,nextw)];
		end
		[v, i] = max(likes);
		% (5)
		d = p*10.^(i-trial/2);
		w += d;
		% (6)
		g2 = grad(X,Y,w);
		if (sum(g2) < e)
			return;
		end
		y = g2-g;
		G2 = G+d*d'/(transpose(d)*y)-G*y*transpose(y)*G/(transpose(y)*G*y);
		% (7)
		G = G2;
		g = g2;
	end
end

使用算法对数据进行学习预测:

w = train(X, Y, 0.00001);
maxent(X, Y, w, [0 2 1 0])
maxent(X, Y, w, [1 2 1 0])
maxent(X, Y, w, [0 0 1 1])

结果为:

ans = 

   1.0000e+00   4.0417e-10

ans =

   8.3811e-10   1.0000e+00

ans =

   0.50000   0.50000

不难发现,特征函数只是简单地判断输入输出对是否存在于样本中,导致模型的泛化能力几乎是没有的。

第7章 支持向量机

  • 7.1 比较感知机的对偶形式与线性可分支持向景机的对偶形式。

  • 7.2 已知正例点,负例点,试求最大间隔分离超平面和分类决策函数,并在图上画出分离超平面、间隔边界及支持向量。

相对于其他的算法,支持向量机的实现实在是太过于复杂,因此此处需要借助LIBSVM

% Data
x = [1 2; 2 3; 3 3; 2 1; 3 2];
y = [1; 1; 1; -1; -1];

% Training
model = svmtrain(y,x,'-s 0 -t 0 -c 100');
w = model.SVs' * model.sv_coef;
b = -model.rho;

% Draw
hf = figure();
hold on
% Draw pointsw = model.SVs' * model.sv_coef;
for i = (1:length(y))
	if (y(i) == 1)
		if (find(model.sv_indices==i))
			plot(x(i,1), x(i,2), 'xb', 'markersize', 16, 'linewidth', 2);
		else 
			plot(x(i,1), x(i,2), 'xb', 'markersize', 16, 'linewidth', 1);
		end
	elseif (y(i) == -1)
		if (find(model.sv_indices==i))
			plot(x(i,1), x(i,2), 'or', 'markersize', 16, 'linewidth', 2);
		else
			plot(x(i,1), x(i,2), 'or', 'markersize', 16, 'linewidth', 1);
		end
	end
end
% Draw margin boundary
for i = model.sv_indices'
	m = (x(i,:)*w+b);
	fplot(@(x) (-w(1)*x-b)/w(2)+m/sqrt(sum(w.^2))*w(1), [0, 4], '--k');
end
% Draw hyperplane
fplot(@(x) (-w(1)*x-b)/w(2), [0, 4], '-k');
% Set axis
set(gca, 'linewidth', 2, 'fontsize', 16);
axis([0, 4, 0, 4]);
% Export to SVG
print(hf, 'svm.svg', '-dsvg');

LIBSVM是不能直接得到的,经过计算,得到:

w =

  -1
   2

b = -2

因此分离超平面为: 分类决策函数为: 所有的支持向量下标在model.sv_indices中,通过计算几何间隔,可以绘制出间隔边界。

  • 7.3 线性支持向蜇机还可以定义为以下形式:

试求其对偶形式。

  • 7.4 证明内积的正整数幕函数:

是正定核函数,这里是正整数,

第8章 提升方法

  • 8.1 某公司招聘职员考查身体、业务能力、发展潜力这 3 项。身体分为合格1、不合格0两级,业务能力和发展潜力分为上1、中2、下3三级分类为合格1 、不合格-1两类。已知10个人的数据,如下表所示。假设弱分类器为决策树桩。.试用AdaBoost算法学习一个强分类器。
  I 2 3 4 5 6 7 8 9 I0
身体 0 0 1 1 1 0 1 1 1 0
业务 1 3 2 1 2 1 1 1 3 2
潜力 3 1 2 3 3 2 2 1 1 1
分类 -1 -1 -1 -1 -1 -1 1 1 -1 -1

为了适配AdaBoost算法,在建立决策树的时候需要加入权重,因此:

  • 计算数据集D的经验熵方法为:
  • 计算特征A对数据集D的经验条件熵H(D|A):
  • 数据集D中实例权重和最大的类作为该结点的类
function ans = entropycond(x, y, w, feature_index)
	ans = 0;
	sample_count = length(y);
	feature_values = unique(x(:,feature_index));
	for i = (1:length(feature_values))
		feature_value = feature_values(i);
		subindex = find(x(:,feature_index)==feature_value);
		suby = y(subindex);
		subw = w(subindex);
		mat = (transpose(unique(suby))==suby).*subw;
		p = sum(mat)/sum(sum(mat));
		ans -= sum(p.*log2(p)).*sum(subw);
	end
end

function feature_index = decision_root_training(x, y, w)
	fcount = size(x,2);
	dcount = length(y);
	feature_index = 1;
	max_entropy = 1;
	for i = 1:fcount
		ent = entropycond(x,y,w,i);
		if (ent<max_entropy)
			feature_index = i;
			max_entropy = ent;
		end
	end
end

function y1 = decision_root_pred(x, y, w, feature_index, x1)
	subx = x(:,feature_index);
	subx1 = x1(:,feature_index);
	feature_values = unique(subx);
	feature_pred = [];
	y1 = [];
	for i = 1:length(feature_values)
		feature_value = feature_values(i);
		subindex = find(x(:,feature_index)==feature_value);
		suby = (y.*w)(subindex);
		feature_pred = [feature_pred;sign(sum(suby))];
	end
	for i = 1:length(subx1)
		y1 = [y1;feature_pred(find(feature_values==subx1(i)))];
	end
end

根据AdaBoost算法:

(1)初始化训练数据的权值分布

(2)对

(a)使用具有权值分布 的训练数据集学习,得到基本分类器

(b)计算 在训练数据集上的分类误差率

(c)计算 的系数

(d)更新训练数据集的权值分布

(3)构建基本分类器的线性组合

function [a b w] = adaboost_training(x, y, m)
	a = b = [];
	scount = length(y);
	w = ones(1,scount)/scount;
	for i = 1:m
		b = [b; decision_root_training(x, y, transpose(w(i,:)))];
		y1 = decision_root_pred(x, y, transpose(w(i,:)), b(i), x);
		e = sum(w(i,:)*(y!=y1));
		a = [a; 1/2*log((1-e)/e)];
		vec = transpose(w(i,:)).*exp(-a(i).*y.*y1);
		w = [w; transpose(vec/sum(vec))];
	end
end

function y1 = adaboost_pred(x, y, a, b, w, x1)
	y1 = zeros(length(x),1);
	for i = 1:(size(w,1)-1)
		y1 += a(i)*decision_root_pred(x,y,transpose(w(i,:)),b(i),x1);
	end
	y1 = sign(y1);
end

训练数据集,当迭代次数为5次以上后,模型完全匹配训练数据集:

x = [0 1 3; 0 3 1; 1 2 2; 1 1 3; 1 2 3; 0 1 2; 1 1 2; 1 1 1; 1 3 1; 0 2 1];
y = [-1; -1; -1; -1; -1; -1; 1; 1; -1; -1];
[a, b, w] = adaboost_training(x, y, 5);
adaboost_pred(x, y, a, b, w, x)

结果为:

ans =

  -1
  -1
  -1
  -1
  -1
  -1
   1
   1
  -1
  -1
  • 8.2 比较支持向昼机、 AdaBoost 、逻辑斯谛回归模型的学习策略与算法。
方法 学习策略 学习算法
支持向量机 极小化正则化合页损失,软间隔最大化 序列最小最优化算法(SMO)
AdaBoost 极小化加法模型指数损失 前向分步加法算法
逻辑斯谛回归 极大似然估计,正则化的极大似然估计 改进的迭代尺度算法,梯度下降,拟牛顿法

第9章 EM算法及其推广

  • 9.1 如例9.1的三硬币模型。假设观洌数据不变,试选择不问的初值,例如, ,求模型参数的极大似然估计。

E步:

M步:

y = [1; 1; 0; 1; 0; 0; 1; 0; 1; 1];
pi = 0.48;
p = 0.55;
q = 0.67;
theta = 1;
while (theta > 0.00001)
	pio = pi;
	po = p;
	qo = q;
	mu = pi*p.^y.*(1-p).^y./(pi*p.^y.*(1-p).^y+(1-pi)*p.^y.*(1-p).^y);
	pi = mean(mu);
	p = sum(mu.*y)/sum(mu);
	q = sum((1-mu).*y)/sum(1-mu);
	theta = abs(pio-pi) + abs(po-p) + abs(qo-q);
end

得到:

pi =  0.48000
p =  0.60000
q =  0.60000
  • 9.2 证明引理9.2。

  • 9.3 已知观测数据

-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49 , 56, 60, 75

试估计两个分囊的高斯混合模型的5个参数。

首先定义高斯分布密度函数:

EM算法步骤为:

(1)取参数的初始值开始迭代

(2)E步:依据当前模型参数,计算分模型k对应的观测数据的响应度

(3)M步:计算新一轮迭代的模型参数

(4)重复第(2)步和第(3)步,直到收敛。

% Data
y = [-67 -48 6 8 14 16 23 24 28 29 41 49 56 60 75];

% EM algorithm
K = 2;
N = length(y);
% (1) Initialize parameters
mu = rand(1,K);
al = rand(1,K);
sg = rand(1,K)*10;
while (true)
	% (2) E step
	gm = zeros(N,K);
	for j = (1:N) 
		for k = (1:K)
			gm(j,k) = al(k)*normpdf(y(j),mu(k),sg(k));
		end
		gm(j,:) /= sum(gm(j,:));
	end
	% (3) M step
	mu2 = y*gm./sum(gm);
	sg2 = sum(gm.*((y-mu').^2)')./sum(gm);
	al2 = sum(gm)/N;
	% (4) Check convergence
	di = sqrt(sum((mu2-mu).^2)+sum((sg2-sg).^2)+sum((al2-al).^2));
	if di < 0.0001
		break
	end
	mu = mu2;
	sg = sg2;
	al = al2;
end

EM算法求解得到的是局部最优解,所以每次运行程序得到的结果可能会不一样。

  • 9.4 EM 算法可以用到朴素贝叶斯法的非监督学习,试写出其算法。

第10章 隐马尔科夫模型

  • 10.1 给定盒子和球组成的隐马尔可夫模型 ,其中,

,试用后向算法计算

根据后向算法:

(1)

(2)对于

(3)

计算:

A = [0.5 0.2 0.3; 0.3 0.5 0.2; 0.2 0.3 0.5];
B = [0.5 0.5; 0.4 0.6; 0.7 0.3];
p = [0.2; 0.4; 0.4];
T = 4;
O = [1 2 1 2];
b = [1; 1; 1];
for t = (T-1:-1:1)
	b = A * B(:,O(t+1)) .* b;
end
sum(p .* B(:,O(1)) .* b)

得到:

ans =  0.061328

给定盒子和球组成的隐马尔可夫模型 ,其中,

,用前后向概率计算

计算:

A = [0.5 0.1 0.4; 0.3 0.5 0.2; 0.2 0.2 0.6];
B = [0.5 0.5; 0.4 0.6; 0.7 0.3];
p = [0.2; 0.3; 0.5];
T = 8;
O = [1; 2; 1; 1; 2; 1; 2; 2];
a4 = p .* B(O(1));
for t = (1:3)
	a4 = A * a4 .* B(:,O(t+1));
end
b4 = [1; 1; 1];
for t = (T-1:-1:4)
	b4 = A * b4 .* B(:,O(t+1));
end 
a4(3) * b4(3) / sum(a4 .* b4)

得到:

ans =  0.33856
  • 10.2 在习题10.1中,试用维特比算法求最优路径

根据维特比算法:

(1)初始化

(2)递推。对

(3)终止

(4)最优路径回溯。对 编写程序:

% Example 10.1
A = [0 1 0 0; 0.4 0 0.6 0; 0 0.4 0 0.6; 0 0 0.5 0.5];
B = [0.5 0.5; 0.3 0.7; 0.6 0.4; 0.8 0.2];
pi = [0.25; 0.25; 0.25; 0.25];
O = [1; 1; 2; 2; 1];

% Example 10.3
% A = [0.5 0.2 0.3; 0.3 0.5 0.2; 0.2 0.3 0.5];
% B = [0.5 0.5; 0.4 0.6; 0.7 0.3];
% pi = [0.2; 0.4; 0.4];
% O = [1; 2; 1]

function path = viterbi(A, B, pi, O)
	T = length(O);
	delta = pi.*B(:,O(1));
	psi = [];
	for t = 2:T
		[value, index] = max(delta.*A);
		psi = [psi;index];
		delta = transpose(value).*B(:,O(t));
	end
	[p, it] = max(delta);
	path = [it];
	for t = T-1:-1:1
		it = psi(t,it);
		path = [it; path];
	end
end

viterbi(A, B, pi, O)

得到:

ans =

   4
   3
   2
   3
   4
  • 10.3 试用前向概率和后向概率推导
  • 10.4 比较维特比算法中变量的计算和前向算法中变量的计算的主要区别。

计算变量的时候直接对上个的结果进行数值计算,而计算变量需要在上个结果计算的基础上选择最大值。

第11章 条件随机场

  • 11.1 写出图11.3中无向图描述的概率图模型的因子分解式。
  • 11.2 证明,其中1是元素均为1的m维列向量。

    (1)证明

根据

因此

(2) 证明

根据

因此

  • 11.3 写出条件随机场模型学习的梯度下降法。

条件随机场的对数极大似然函数为:

因此梯度下降算法的目标函数就是:

目标函数的梯度函数为:

按照梯度下降算法:

(1) 取初始值,置t=0

(2) 计算

(3) 计算梯度,当,停止迭代,令;否则,令,求,使

(4) 置,计算

时,停止迭代,令

(5) 否则,置,转(3)

  • 11.4 参考图11.6的状态路径图,假设随机矩阵 分别是

    求以为起点为终点的所有路径的状态序列 y 的概率及概率最大的状态序列。

使用Octave计算:

M = zeros(2,2,4);
M(:,:,1) = [0 0; 0.5 0.5];
M(:,:,2) = [0.3 0.7; 0.7 0.3];
M(:,:,3) = [0.5 0.5; 0.6 0.4];
M(:,:,4) = [0 1; 0 1];

% Generate paths
pa = [2];
for i = (1:3)
	np = [];
	for row = transpose(pa)
		np = [np ; [transpose(row) 1]];
		np = [np ; [transpose(row) 2]];
	end
	pa = np;
end
pa = [pa ones(length(pa),1)*2];

% Compute probability
z = eye(2);
for i = (1:length(M))
	z = z * M(:,:,i);
end
z = z(2,2);
pr = [];
for row = transpose(pa)
	p = 1;
	for i = (1:length(row)-1)
		a = row(i);
		b = row(i+1);
		p = p * M(a,b,i);
	end
	pr = [pr; p];
end 

% Find the path with maximum probability
[value index] = max(pr);
mp = pa(index,:);

对于每一个状态序列:

pa =

   2   1   1   1   2
   2   1   1   2   2
   2   1   2   1   2
   2   1   2   2   2
   2   2   1   1   2
   2   2   1   2   2
   2   2   2   1   2
   2   2   2   2   2

对应的概率为:

pr =

   0.075000
   0.075000
   0.210000
   0.140000
   0.175000
   0.175000
   0.090000
   0.060000

其中概率最大的状态序列为:

mp =

   2   1   2   1   2