function beta = mmqr (q, p, y, x, func, toler, beta)
% Use MM algorithm to minimize modified (smoothed) objective
% function for computing regression quantiles.
%
% Inputs:
% -- q is the quantile desired (number in (0,1)).
% -- p is the size of the parameter vector. This is actually
% redundant since it's just the length of beta, which is also
% input, but I didn't change the code to eliminate this
% redundancy. Note that p is not necessarily the number of columns
% of x, though this is usually the case.
% -- y is the response vector (n x 1).
% -- x is the predictor matrix (n x m).
% -- func is a character array containing the name of a function
% that evaluates f(x, beta). This MATLAB function takes (x, beta)
% as inputs where x is the n x m predictor matrix and beta is some
% p x 1 vector of parameter values. It returns an n x 1 vector.
% In the case of linear QR, m=p and the function (whatever it's
% called) returns x*beta.
% There must also be a second MATLAB function, with the same name
% as func except with the letter d appended to the front, that
% also takes (x,beta) as input but returns the $n x p$ matrix
% of partial derivatives of f(x,beta) evaluated at beta. In the
% linear QR example, m=p and this second function simply returns x.
% -- toler is a small number giving the minimum change in the value
% of the surrogate function before convergence is declared.
% -- beta is the starting value of beta (p x 1 vector).
flops(0); % reset flops count.
dfunc = ['d' func]; % name of differential function.
iteration=0;
change=realmax;
n=length(y);
tn=toler/n; e0=-tn/log(tn); epsilon=(e0-tn)/(1+log(e0));
% epsilon is the smoothing parameter. As suggested in the paper, it
% is chosen to roughly satisfy -epsilon(log epsilon)=toler/n using a
% Newton-Raphson step above.
res = y - feval(func,x,beta); % Compute residuals for current beta.
weights = 1./(epsilon+abs(res));
v = 1-2*q - res.*weights;
lastsurr = -sum(res.*v) - (1-2*q)*sum(res); % Last value of surrogate (up
% to a constant).
while change > toler
iteration = iteration + 1;
J=feval(dfunc,x,beta); % J is the first differential matrix.
step = -inv(J'*(weights(:,ones(p,1)).*J)) * (J' * v);
beta=beta+step;
res=y-feval(func,x,beta);
surr = sum(res.^2.*weights) + (4*q-2) * sum(res);
% Now do the step-halving procedure to ensure a decrease in the
% value of the surrogate function.
while surr > lastsurr
step=step/2;
beta=beta-step;
res=y-feval(func,x,beta);
surr = sum(res.^2.*weights) + (4*q-2) * sum(res);
end
change=lastsurr-surr;
weights = 1./(epsilon+abs(res));
v = 1-2*q - res.*weights;
lastsurr = -sum(res.*v) - (1-2*q)*sum(res);
end
lastobj=q*sum(res)-sum(res(res<0));
% At this point, beta is the estimate of regression coefficients;
% lastobj is the value of the objective function at beta;
% flops is the number of floating-point operations needed
% (according to the way MATLAB counts them, at least).