function theta = IRLS(X, y, link_fn, num_iter, theta)
% IRLS Iteratively Reweighted Least Squares algorithm for fitting a Generalized Linear Model.
% theta = IRLS(X, y, link_fn, num_iter, theta)
%
% link_fn can be 'logistic' or 'noisyor'.
% Currently we can only handle scalar output.
% See also glmtrain.
%
% In a GLIM, E[Y|X=x] = f(eta), where eta = theta'*x and f is the link fn.
%
% X(l,:) = l'th input
% y(l) = l'th output (target)

% If link_fn = logistic, f(eta) = 1/(1+exp(-eta)), f'(eta) = f(eta) * (1-f(eta))
% If link_fn = noisyor,  f(eta) = 1-exp(eta), f'(eta) = -exp(eta)

%nex = size(X,1);
%input_dim = size(X,2);
%theta = randn(input_dim,1);

y = y(:);
theta = theta(:);
for t=1:num_iter
  eta = X*theta;
  switch link_fn
    case 'logistic',
      f = 1 ./ (1 + exp(-eta));
      fp = f .* (1-f);
    case 'noisyor',
      fp = -exp(eta);
      f = 1 + fp;
    otherwise,
      error(['unrecognized link fn ' link_fn]);
  end
  W = diag(fp);
  e = (y - f) ./ (fp + eps);
  z = eta + e;
  theta = pinv(X'*W*X)*X'*W*z;
end
