Download All

Code covered by the BSD License  

Highlights from
Statistical Learning Toolbox

from Statistical Learning Toolbox by Dahua Lin
Functions for statistical learning, pattern recognition and computer vision, covering many topics.

Description of sllogistreg
Home > sltoolbox > regression > sllogistreg.m

sllogistreg

PURPOSE ^

SLLOGISTREG Performs Multivariate Logistic Regression

SYNOPSIS ^

function [A, b, props, info] = sllogistreg(X, nums, varargin)

DESCRIPTION ^

SLLOGISTREG Performs Multivariate Logistic Regression

 $ Syntax $
   - [A, b] = sllogistreg(X, nums, ...)
   - [A, b, props] = sllogistreg(X, nums, ...)
   - [A, b, props, info] = sllogistreg(X, nums, ...)

 $ Arguments $
   - X:            The input sample matrix
   - nums:         The numbers of samples in different classes
   - A:            The solved logistic coefficients
   - b:            The solved logistic shift value
   - props:        The probability values under the resultant model
   - info:         The information of learning process
                   - exitflag: The exitflag given by fminunc
                   - numiters: The number of iterations
                   - fval:     The final objective value

 $ Description $
   - [A, b] = sllogistreg(X, nums, ...) solves the multivariate logistic 
     regression. Suppose there are n samples, from C classes, then C
     models are learned for the C classes. The k-th model is characterized
     by a coefficient vector a_k and a shift value b_k. Then a sample x,
     the probability that it belongs to the k-th class is given by the
     following formula:
           p(x, k) = P_k * h_k(x) / sum_l (P_l * h_l(x))
     here, we have
           h_k(x) = a_k^T * x + b_k
     The objective of logistic regression is to maximize the sum of 
     the probabilities that the samples belong to its own class:
           maximize sum_i p(x_i, c_i)
     here, c_i is the class label corresponding to x_i.
     In the input arguments, X should be a d x n matrix containing the
     x samples with the samples from the same class gathering together,
     nums should be a 1 x C row vector. 
     In the output arguments, A should be a d x C matrix, and b is a
     1 x C row vector. Each column in A together with the corresponding
     shift value in b describe a model for one class.

     You can specify the following parameter to control the numeric 
     optimization process:
       - 'weights':        The weights of the samples  (1 x n)
                           default = []: means all weights are 1
       - 'priors':         The priors of the classes (1 x C)
                           default = []: means equal priors
       - 'maxiter':        The maximum number of iterations
                           default = 300
       - 'tolF':           The tolerance of the change of objective func
                           default = 1e-6
       - 'tolX':           The tolerance of the change of parameters
                           default = 1e-6
       - 'display':        The level of display
                           default = 'off'
       - 'init':           A cell array as {A0, b0}
                           default = {}, using random initialization

   - [A, b, props] = sllogistreg(X, nums, ...) also output the 
     probability values that the input samples belong to the true class.
     It would be a 1 x n row vector.

   - [A, b, props, info] = sllogistreg(X, nums, ...) outputs the info
     about the optimization process
 
 $ Remarks $
   - This function is implemented based on the optimization function
     fminunc and fmincon, and thus requires the optimization toolbox 
     be installed.

 $ History $
   - Created by Dahua Lin, on Sep 16, 2006

CROSS-REFERENCE INFORMATION ^

This function calls:
  • sladdvec SLADDVEC adds a vector to columns or rows of a matrix
  • slmulvec SLMULVEC multiplies a vector to columns or rows of a matrix
  • slposteriori SLPOSTERIORI Computes the posterioris
  • slposterioritrue SLPOSTERIORITRUE Computes the posteriori that samples belong to true class
  • raise_lackinput RAISE_LACKINPUT Raises an error indicating lack of input argument
  • slexpand SLEXPAND Expand a set to multiple instance
  • slnums2bounds SLNUMS2BOUNDS Compute the index-boundaries from section sizes
  • slparseprops SLPARSEPROPS Parses input parameters
This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [A, b, props, info] = sllogistreg(X, nums, varargin)
0002 %SLLOGISTREG Performs Multivariate Logistic Regression
0003 %
0004 % $ Syntax $
0005 %   - [A, b] = sllogistreg(X, nums, ...)
0006 %   - [A, b, props] = sllogistreg(X, nums, ...)
0007 %   - [A, b, props, info] = sllogistreg(X, nums, ...)
0008 %
0009 % $ Arguments $
0010 %   - X:            The input sample matrix
0011 %   - nums:         The numbers of samples in different classes
0012 %   - A:            The solved logistic coefficients
0013 %   - b:            The solved logistic shift value
0014 %   - props:        The probability values under the resultant model
0015 %   - info:         The information of learning process
0016 %                   - exitflag: The exitflag given by fminunc
0017 %                   - numiters: The number of iterations
0018 %                   - fval:     The final objective value
0019 %
0020 % $ Description $
0021 %   - [A, b] = sllogistreg(X, nums, ...) solves the multivariate logistic
0022 %     regression. Suppose there are n samples, from C classes, then C
0023 %     models are learned for the C classes. The k-th model is characterized
0024 %     by a coefficient vector a_k and a shift value b_k. Then a sample x,
0025 %     the probability that it belongs to the k-th class is given by the
0026 %     following formula:
0027 %           p(x, k) = P_k * h_k(x) / sum_l (P_l * h_l(x))
0028 %     here, we have
0029 %           h_k(x) = a_k^T * x + b_k
0030 %     The objective of logistic regression is to maximize the sum of
0031 %     the probabilities that the samples belong to its own class:
0032 %           maximize sum_i p(x_i, c_i)
0033 %     here, c_i is the class label corresponding to x_i.
0034 %     In the input arguments, X should be a d x n matrix containing the
0035 %     x samples with the samples from the same class gathering together,
0036 %     nums should be a 1 x C row vector.
0037 %     In the output arguments, A should be a d x C matrix, and b is a
0038 %     1 x C row vector. Each column in A together with the corresponding
0039 %     shift value in b describe a model for one class.
0040 %
0041 %     You can specify the following parameter to control the numeric
0042 %     optimization process:
0043 %       - 'weights':        The weights of the samples  (1 x n)
0044 %                           default = []: means all weights are 1
0045 %       - 'priors':         The priors of the classes (1 x C)
0046 %                           default = []: means equal priors
0047 %       - 'maxiter':        The maximum number of iterations
0048 %                           default = 300
0049 %       - 'tolF':           The tolerance of the change of objective func
0050 %                           default = 1e-6
0051 %       - 'tolX':           The tolerance of the change of parameters
0052 %                           default = 1e-6
0053 %       - 'display':        The level of display
0054 %                           default = 'off'
0055 %       - 'init':           A cell array as {A0, b0}
0056 %                           default = {}, using random initialization
0057 %
0058 %   - [A, b, props] = sllogistreg(X, nums, ...) also output the
0059 %     probability values that the input samples belong to the true class.
0060 %     It would be a 1 x n row vector.
0061 %
0062 %   - [A, b, props, info] = sllogistreg(X, nums, ...) outputs the info
0063 %     about the optimization process
0064 %
0065 % $ Remarks $
0066 %   - This function is implemented based on the optimization function
0067 %     fminunc and fmincon, and thus requires the optimization toolbox
0068 %     be installed.
0069 %
0070 % $ History $
0071 %   - Created by Dahua Lin, on Sep 16, 2006
0072 %
0073 
0074 %% parse and verify input arguments
0075 
0076 if nargin < 2
0077     raise_lackinput('sllogistreg', 2);
0078 end
0079 
0080 [d, n] = size(X);
0081 if ~isvector(nums) || size(nums, 1) ~= 1
0082     error('sltoolbox:invalidarg', ...
0083         'nums should be a 1 x C row vector');
0084 end
0085 if sum(nums) ~= n
0086     error('sltoolbox:sizmismatch', ...
0087         'The numbers in nums are not consistent with the sample number');
0088 end
0089 
0090 C = length(nums);
0091 
0092 if C < 2
0093     error('sltoolbox:invalidarg', ...
0094         'There should be at least 2 classes.');
0095 end
0096 
0097 opts.weights = [];
0098 opts.priors = [];
0099 opts.maxiter = 300;
0100 opts.tolF = 1e-6;
0101 opts.tolX = 1e-6;
0102 opts.display = 'off';
0103 opts.init = {};
0104 opts = slparseprops(opts, varargin{:});
0105 
0106 w = opts.weights;
0107 if ~isempty(w)
0108     if ~isequal(size(w), [1, n])
0109         error('sltoolbox:sizmismatch', ...
0110             'w should be a 1 x n row vector');
0111     end
0112 end
0113 
0114 pri = opts.priors;
0115 if ~isempty(pri)
0116     if ~isequal(size(pri), [1, C])
0117         error('sltoolbox:sizmismatch', ...
0118             'pri should be a 1 x C row vector');
0119     end
0120 end
0121 
0122 if isempty(opts.init)
0123     is_inited = false;
0124 else
0125     A0 = opts.init{1};
0126     b0 = opts.init{2};
0127     if ~isequal(size(A0), [d, C])
0128         error('sltoolbox:sizmismatch', 'A0 should be a d x C matrix');
0129     end
0130     if ~isequal(size(b0), [1, C])
0131         error('sltoolbox:sizmismatch', 'b0 should be a 1 x C row vector');
0132     end
0133     is_inited = true;    
0134 end
0135 
0136 
0137 %% main
0138 
0139 %NOTE: we convert the problem of maximization to minimization of the
0140 %      negative object
0141 
0142 % augment problem
0143 
0144 Xa = [X; ones(1, n)];
0145 
0146 % prepare for optimization
0147 optimfunc = @(v) logistic_objfun(v, Xa, nums, d, n, C, pri, w);
0148 optimopts = optimset(optimset('fminunc'), ...
0149     'LargeScale', 'on', ...
0150     'GradObj', 'on', ...
0151     'MaxIter', opts.maxiter, ...
0152     'Display', opts.display, ...
0153     'TolFun', opts.tolF, ...
0154     'TolX', opts.tolX);
0155     
0156 % initialization
0157 if ~is_inited
0158     v0 = rand((d+1)*C, 1);
0159 else
0160     v0 = reshape([A0; b0], (d+1)*C, 1);
0161 end
0162 
0163 % do optimization
0164 [v, fval, exitflag, optimoutput] = fminunc(optimfunc, v0, optimopts);
0165 clear v0;
0166 clear optimfunc;
0167 clear Xa;
0168 
0169 % make output
0170 v = reshape(v, d+1, C);
0171 A = v(1:d, :);
0172 b = v(d+1, :);
0173 clear v;
0174 
0175 if nargout >= 3
0176     L = compute_logit(A, X) ;
0177     L = sladdvec(L, b', 1);
0178     props = slposterioritrue(L, nums, pri, 'log');            
0179 end
0180 
0181 if nargout >= 4
0182     info.exitflag = exitflag;
0183     info.numiters = optimoutput.iterations;
0184     info.fval = -fval; 
0185 end
0186 
0187 
0188 %% The core functions for objective and gradient evaluation
0189 
0190 function [f, g] = logistic_objfun(v, Xa, nums, d, n, C, pri, w)
0191 % v is a reshape version of Aa
0192 
0193 % get input
0194 Aa = reshape(v, d+1, C);
0195 L = compute_logit(Aa, Xa);
0196 clear Aa;
0197 
0198 % compute f
0199 P = slposteriori(L, pri, 'log');
0200 [sps, eps] = slnums2bounds(nums);
0201 pps = zeros(1, n);
0202 for k = 1 : C
0203     sk = sps(k); ek = eps(k);
0204     pps(sk:ek) = P(k, sk:ek);
0205 end
0206 
0207 if isempty(w)
0208     f = -sum(log(pps));
0209 else
0210     f = -sum(log(pps) .* w);
0211 end
0212 
0213 % compute g
0214 M = make_indicatormap(nums, C, n);
0215 M = M - P;
0216 clear P;
0217 if ~isempty(w)
0218     M = slmulvec(M, w, 2);
0219 end
0220 g = Xa * M';
0221 clear M;
0222 g = -g(:);
0223 
0224 
0225 function L = compute_logit(Aa, Xa)
0226 % L is the logit values (a_k' * x_i + b_k): C x n matrix
0227 L = Aa' * Xa;
0228 
0229 
0230 %% Auxiliary functions
0231 
0232 function M = make_indicatormap(nums, C, n)
0233 % C x n
0234 % C(i, j) = 1, when sample j belongs to class i
0235 
0236 M = zeros(C, n);
0237 I = slexpand(nums);
0238 J = 1:n;
0239 inds = sub2ind([C, n], I, J);
0240 clear I J;
0241 M(inds) = 1;
0242
Generated on Wed 20-Sep-2006 12:43:11 by m2html © 2003

Contact us at files@mathworks.com