【独家发布】Pattern Recognition and Machine Learning

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

eBook: 0.zip (8.18 MB)

2014-11-25 05:55:36 上传

Matlab Code:https://github.com/PRML/PRML

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏5 回帖

关键词：Recognition cognition Learning machine Pattern

1. function y = pdfDirichletLn(X, a)
   
2. % Compute log pdf of a Dirichlet distribution.
   
3. %   X: d x n data matrix satifying (sum(X,1)==ones(1,n) && X>=0)
   
4. %   a: d x k parameters
   
5. %   y: k x n probability density
   
6. % Written by Mo Chen (mochen80@gmail.com).
   
7. X = bsxfun(@times,X,1./sum(X,1));
   
8. if size(a,1) == 1
   
9.     a = repmat(a,size(X,1),1);
   
10. end
    
11. c = gammaln(sum(a,1))-sum(gammaln(a),1);
    
12. g = (a-1)'*log(X);
    
13. y = bsxfun(@plus,g,c');

复制代码

1. function y = pdfGaussLn(X, mu, sigma)
   
2. % Compute log pdf of a Gaussian distribution.
   
3. % Written by Mo Chen (mochen80@gmail.com).
   
4. 
   
5. [d,n] = size(X);
   
6. k = size(mu,2);
   
7. if n == k && size(sigma,1) == 1           
   
8.     X = bsxfun(@times,X-mu,1./sigma);
   
9.     q = dot(X,X,1);  % M distance
   
10.     c = d*log(2*pi)+2*log(sigma);          % normalization constant
    
11.     y = -0.5*(c+q);
    
12. elseif size(sigma,1)==d && size(sigma,2)==d && k==1   % one mu and one dxd sigma
    
13.     X = bsxfun(@minus,X,mu);
    
14.     [R,p]= chol(sigma);
    
15.     if p ~= 0
    
16.         error('ERROR: sigma is not PD.');
    
17.     end
    
18.     Q = R'\X;
    
19.     q = dot(Q,Q,1);  % quadratic term (M distance)
    
20.     c = d*log(2*pi)+2*sum(log(diag(R)));   % normalization constant
    
21.     y = -0.5*(c+q);
    
22. elseif size(sigma,1)==d && size(sigma,2)==k % k mu and k diagonal sigma
    
23.     lambda = 1./sigma;
    
24.     ml = mu.*lambda;
    
25.     q = bsxfun(@plus,X'.^2*lambda-2*X'*ml,dot(mu,ml,1)); % M distance
    
26.     c = d*log(2*pi)+2*sum(log(sigma),1); % normalization constant
    
27.     y = -0.5*bsxfun(@plus,q,c);
    
28. elseif size(sigma,1)==1 && (size(sigma,2)==k || size(sigma,2)==1) % k mu and (k or one) scalar sigma
    
29.     X2 = repmat(dot(X,X,1)',1,k);
    
30.     D = bsxfun(@plus,X2-2*X'*mu,dot(mu,mu,1));
    
31.     q = bsxfun(@times,D,1./sigma);  % M distance
    
32.     c = d*(log(2*pi)+2*log(sigma));          % normalization constant
    
33.     y = -0.5*bsxfun(@plus,q,c);
    
34. else
    
35.     error('Parameters mismatched.');
    
36. end

复制代码

1. function [y, sigma, p] = linInfer(X, model, t)
   
2. % Compute linear model reponse y = w'*x+b and likelihood
   
3. % X: d x n data
   
4. % t: 1 x n response
   
5. w = model.w;
   
6. b = model.w0;
   
7. y = w'*X+b;
   
8. if nargout > 1
   
9.     beta = model.beta;
   
10. 
    
11.     if isfield(model,'V')   % V*V'=inv(S) 3.54
    
12.         X = model.V'*bsxfun(@minus,X,model.xbar);
    
13.         sigma = sqrt(1/beta+dot(X,X,1));   % 3.59
    
14.     else
    
15.         sigma = sqrt(1/beta);
    
16.     end
    
17.    
    
18.     if nargin == 3 && nargout == 3
    
19. %         p = exp(logGaussPdf(t,y,sigma));
    
20.         p = exp(-0.5*(((t-y)./sigma).^2+log(2*pi))-log(sigma));
    
21.     end
    
22. end

复制代码

1. function [model, llh] = regressEbEm(X, t, alpha, beta)
   
2. % Fit empirical Bayesian linear model with EM
   
3. % X: d x n data
   
4. % t: 1 x n response
   
5. if nargin < 3
   
6.     alpha = 0.02;
   
7.     beta = 0.5;
   
8. end
   
9. [d,n] = size(X);
   
10. 
    
11. xbar = mean(X,2);
    
12. tbar = mean(t,2);
    
13. 
    
14. X = bsxfun(@minus,X,xbar);
    
15. t = bsxfun(@minus,t,tbar);
    
16. 
    
17. C = X*X';
    
18. Xt = X*t';
    
19. dg = sub2ind([d,d],1:d,1:d);
    
20. I = eye(d);
    
21. tol = 1e-4;
    
22. maxiter = 100;
    
23. llh = -inf(1,maxiter+1);
    
24. for iter = 2:maxiter
    
25.     A = beta*C;
    
26.     A(dg) = A(dg)+alpha;
    
27.     U = chol(A);
    
28.     V = U\I;
    
29. 
    
30.     w = beta*(V*(V'*Xt));
    
31.     w2 = dot(w,w);
    
32.     err = sum((t-w'*X).^2);
    
33.    
    
34.     logdetA = 2*sum(log(diag(U)));   
    
35.     llh(iter) = 0.5*(d*log(alpha)+n*log(beta)-alpha*w2-beta*err-logdetA-n*log(2*pi));
    
36.     if llh(iter)-llh(iter-1) < tol*abs(llh(iter-1)); break; end
    
37.    
    
38.     trS = dot(V(:),V(:));
    
39.     alpha = d/(w2+trS);   % 9.63
    
40.    
    
41.     gamma = d-alpha*trS;
    
42.     beta = n/(err+gamma/beta);
    
43. end
    
44. b = tbar-dot(w,xbar);
    
45. 
    
46. llh = llh(2:iter);
    
47. model.b = b;
    
48. model.w = w;
    
49. model.alpha = alpha;
    
50. model.beta = beta;
    
51. model.xbar = xbar;
    
52. model.V = V;

复制代码

1. function [model, llh] = regressEbFp(X, t, alpha, beta)
   
2. % Fit empirical Bayesian linear model with Mackay fixed point method
   
3. % X: d x n data
   
4. % t: 1 x n response
   
5. if nargin < 3
   
6.     alpha = 0.02;
   
7.     beta = 0.5;
   
8. end
   
9. [d,n] = size(X);
   
10. 
    
11. xbar = mean(X,2);
    
12. tbar = mean(t,2);
    
13. 
    
14. X = bsxfun(@minus,X,xbar);
    
15. t = bsxfun(@minus,t,tbar);
    
16. 
    
17. C = X*X';
    
18. Xt = X*t';
    
19. dg = sub2ind([d,d],1:d,1:d);
    
20. I = eye(d);
    
21. tol = 1e-4;
    
22. maxiter = 100;
    
23. llh = -inf(1,maxiter+1);
    
24. for iter = 2:maxiter
    
25.     A = beta*C;
    
26.     A(dg) = A(dg)+alpha;
    
27.     U = chol(A);
    
28.     V = U\I;
    
29. 
    
30.     w = beta*(V*(V'*Xt));
    
31.     w2 = dot(w,w);
    
32.     err = sum((t-w'*X).^2);
    
33.    
    
34.     logdetA = 2*sum(log(diag(U)));   
    
35.     llh(iter) = 0.5*(d*log(alpha)+n*log(beta)-alpha*w2-beta*err-logdetA-n*log(2*pi));
    
36.     if llh(iter)-llh(iter-1) < tol*abs(llh(iter-1)); break; end
    
37.    
    
38.     trS = dot(V(:),V(:));
    
39.     gamma = d-alpha*trS;
    
40.     alpha = gamma/w2;
    
41.     beta = (n-gamma)/err;
    
42. end
    
43. b = tbar-dot(w,xbar);
    
44. 
    
45. llh = llh(2:iter);
    
46. model.b = b;
    
47. model.w = w;
    
48. model.alpha = alpha;
    
49. model.beta = beta;
    
50. model.xbar = xbar;
    
51. model.V = V;

复制代码

1. function U = classFda(X, y, d)
   
2. % Fisher (linear) discriminant analysis
   
3. n = size(X,2);
   
4. k = max(y);
   
5. 
   
6. E = sparse(1:n,y,true,n,k,n);  % transform label into indicator matrix
   
7. nk = full(sum(E));
   
8. 
   
9. m = mean(X,2);
   
10. Xo = bsxfun(@minus,X,m);
    
11. St = (Xo*Xo')/n;
    
12. 
    
13. mk = bsxfun(@times,X*E,1./nk);
    
14. mo = bsxfun(@minus,mk,m);
    
15. mo = bsxfun(@times,mo,sqrt(nk/n));
    
16. Sb = mo*mo';
    
17. % Sw = St-Sb;
    
18. 
    
19. [U,A] = eig(Sb,St,'chol');
    
20. [~,idx] = sort(diag(A),'descend');
    
21. U = U(:,idx(1:d));

复制代码

1. function [model, llh] = classLogitBin(X, t, lambda)
   
2. % logistic regression for binary classification (Bernoulli likelihood)
   
3. if any(unique(t) ~= [0,1])
   
4.     error('t must be a 0/1 vector!');
   
5. end
   
6. if nargin < 3
   
7.     lambda = 1e-4;
   
8. end
   
9. [d,n] = size(X);
   
10. dg = sub2ind([d,d],1:d,1:d);
    
11. X = [X; ones(1,n)];
    
12. d = d+1;
    
13. 
    
14. tol = 1e-4;
    
15. maxiter = 100;
    
16. llh = -inf(1,maxiter);
    
17. converged = false;
    
18. iter = 1;
    
19. 
    
20. 
    
21. h = ones(1,n);
    
22. h(t==0) = -1;
    
23. w = zeros(d,1);
    
24. z = w'*X;
    
25. while ~converged && iter < maxiter
    
26.     iter = iter + 1;
    
27.     y = sigmoid(z);
    
28.    
    
29.     g = X*(y-t)'+lambda*w;
    
30.     r = y.*(1-y);
    
31.     R = spdiags(r(:),0,n,n);   
    
32.     H = X*R*X';
    
33.     H(dg) = H(dg)+lambda;
    
34.    
    
35.     p = -H\g;
    
36. 
    
37.     wo = w;
    
38.     w = wo+p;
    
39.     z = w'*X;   
    
40.     llh(iter) = -sum(log1pexp(-h.*z))-0.5*lambda*dot(w,w);
    
41.     converged = norm(p) < tol || abs(llh(iter)-llh(iter-1)) < tol;
    
42.     while ~converged && llh(iter) < llh(iter-1)
    
43.         p = 0.5*p;
    
44.         w = wo+p;
    
45.         z = w'*X;   
    
46.         llh(iter) = -sum(log1pexp(-h.*z))-0.5*lambda*dot(w,w);
    
47.         converged = norm(p) < tol || abs(llh(iter)-llh(iter-1)) < tol;
    
48.     end
    
49. end
    
50. llh = llh(2:iter);
    
51. model.w = w;

复制代码

1. function [model, llh] = classLogitMul(X, t, lambda, method)
   
2. % logistic regression for multiclass problem (Multinomial likelihood)
   
3. if nargin < 4
   
4.     method = 1;
   
5. end
   
6. 
   
7. if nargin < 3
   
8.     lambda = 1e-4;
   
9. end
   
10. 
    
11. X = [X; ones(1,size(X,2))];
    
12. 
    
13. if method == 1
    
14.     [W, llh] = NewtonSolver(X, t, lambda);
    
15. else
    
16.     [W, llh] = blockNewtonSolver(X, t, lambda);
    
17. end
    
18. model.W = W;
    
19. 
    
20. function [W, llh] = NewtonSolver(X, t, lambda)
    
21. 
    
22. [d,n] = size(X);
    
23. k = max(t);
    
24. 
    
25. tol = 1e-4;
    
26. maxiter = 100;
    
27. llh = -inf(1,maxiter);
    
28. converged = false;
    
29. iter = 1;
    
30. dk = d*k;
    
31. dg = sub2ind([dk,dk],1:dk,1:dk);
    
32. T = sparse(1:n,t,1,n,k,n);
    
33. W = zeros(d,k);
    
34. HT = zeros(d,k,d,k);
    
35. while ~converged && iter < maxiter
    
36.     iter = iter+1;
    
37.     Z = X'*W;
    
38.     logY = bsxfun(@minus,Z,logsumexp(Z,2));
    
39.     llh(iter) = dot(T(:),logY(:))-0.5*lambda*dot(W(:),W(:));
    
40.     converged = abs(llh(iter)-llh(iter-1)) < tol;
    
41.    
    
42.     Y = exp(logY);
    
43.     for i = 1:k
    
44.         for j = 1:k
    
45.             r = Y(:,i).*((i==j)-Y(:,j));
    
46.             HT(:,i,:,j) = bsxfun(@times,X,r')*X';
    
47.         end
    
48.     end
    
49.     G = X*(Y-T)+lambda*W;
    
50.     H = reshape(HT,dk,dk);
    
51.     H(dg) = H(dg)+lambda;
    
52.     W(:) = W(:)-H\G(:);
    
53. end
    
54. llh = llh(2:iter);
    
55. 
    
56. function [W, llh] = blockNewtonSolver(X, t, lambda)
    
57. 
    
58. [d,n] = size(X);
    
59. k = max(t);
    
60. 
    
61. dg = sub2ind([d,d],1:d,1:d);
    
62. tol = 1e-4;
    
63. maxiter = 100;
    
64. llh = -inf(1,maxiter);
    
65. converged = false;
    
66. iter = 1;
    
67. 
    
68. T = sparse(1:n,t,1,n,k,n);
    
69. W = zeros(d,k);
    
70. Z = X'*W;
    
71. logY = bsxfun(@minus,Z,logsumexp(Z,2));
    
72. Y = exp(logY);
    
73. while ~converged && iter < maxiter
    
74.     iter = iter+1;
    
75.     for j = 1:k
    
76.         r = Y(:,j).*(1-Y(:,j));
    
77.         H = bsxfun(@times,X,r')*X';
    
78.         H(dg) = H(dg)+lambda;
    
79. 
    
80.         g = X*(Y(:,j)-T(:,j))+lambda*W(:,j);
    
81.         W(:,j) = W(:,j)-H\g;
    
82.         Z(:,j) = X'*W(:,j);
    
83.         logY = bsxfun(@minus,Z,logsumexp(Z,2));
    
84.         Y = exp(logY);
    
85.     end
    
86.    
    
87.     llh(iter) = dot(T(:),logY(:))-0.5*lambda*dot(W(:),W(:));
    
88.     converged = abs(llh(iter)-llh(iter-1)) < tol;
    
89. end
    
90. llh = llh(2:iter);

复制代码

1. function [model, llh] = classRvmEbEm(X, t, alpha)
   
2. % Relevance Vector Machine classification training by empirical bayesian (ARD)
   
3. % using standard EM update
   
4. if nargin < 3
   
5.     alpha = 0.02;
   
6. end
   
7. n = size(X,2);
   
8. X = [X;ones(1,n)];
   
9. d = size(X,1);
   
10. alpha = alpha*ones(d,1);
    
11. 
    
12. tol = 1e-4;
    
13. maxiter = 100;
    
14. llh = -inf(1,maxiter);
    
15. infinity = 1e+10;
    
16. for iter = 2:maxiter
    
17.     used = alpha < infinity;
    
18.     alphaUsed = alpha(used);
    
19.     [w,V,lllh] = optLogitNewton(X(used,:),t,alphaUsed);
    
20.     w2 = w.^2;
    
21.    
    
22.     logdetS = -2*sum(log(diag(V)));   
    
23.     llh(iter) = lllh+0.5*(sum(log(alphaUsed))-logdetS-dot(alphaUsed,w2)-n*log(2*pi)); % 7.114
    
24.     if abs(llh(iter)-llh(iter-1)) < tol; break; end
    
25. 
    
26.     dgS = dot(V,V,2);
    
27.     alpha = 1./(w2+dgS);    % 9.67
    
28. end
    
29. llh = llh(2:iter);
    
30. 
    
31. model.used = used;
    
32. model.w = w;
    
33. model.alpha = alpha;

复制代码