load huntington.mat

n = size(hunt,1); % # of training points
y = hunt.disease; % binary: 1=disease

% linear model: t = \beta_0 + \beta_1*CAGs
X = [ones(n,1) hunt.CAGs];
beta = zeros([size(X,2), 1]);
n_param = length(beta);

alpha = 0.035; % step size
max_iter = 500000; % fixed step size requires many iterations
g = zeros(size(beta)); % pre-allocate the gradient
for i = 1:max_iter
    % model predictions using the current parameters;
    % used to compute the gradient
    t = X*beta;
    P = 1./(1 + exp(-t));
    % fill out each entry of the gradient
    for j = 1:n_param
        dLdP = (1-y)./(1-P) - y./P;
        dPdt = P .* (1-P);
        dtdB = X(:,j);
        g(j) = 1/n*sum(dLdP .* dPdt .* dtdB);
    end
    % update
    beta = beta - alpha*g;
end

beta