function [beta, converged] = logisticReg(X, y, tol, beta)
  % Two class logistic regression
  % Input: X, data matrix  d+1 x N, 1st row of X is ones
  %        y, class indicator, 1 or 0, a 1 by N vector
  %        tol, (Optional) convergence tolerance, absolute step size 
  %             (Default: 1.0e-6)
  %        beta, (Optional) d+1 dimensional initial guess
  %             (Default: zeros(d+1,0), which is an excellent initial guess)
  % Output: beta = (d+1) vector, where
  %                p(y=1 | x) = exp(beta' * x)/(1 + exp(beta' * x))
  %                p(y=0 | x) = 1 - p(y=1 | x)
  %         converged - boolean, whether last update was less than tol.
  
  if ~exist('tol', 'var')
    tol = 1.0e-6;
  end
  if ~exist('beta', 'var')
    beta = zeros(size(X,1), 1);
  end
  its = 0;
  converged = false;
  while (its < 50)
    p = logistic(X, beta);
    W = p .* (1-p);
    
    H = (X .* repmat(W, size(X,1), 1)) * (X');
    
    betaStep = pinv(H) * ( X * (y - p)');
    
    if (max(abs(betaStep)) < tol)
      converged = true;
      break;
    end
    
    beta = beta + betaStep;
    
    its = its+1;
    
  end
  
  return