Логистическая регрессия (пример)

b0 = log(mean(y)/(1-mean(y))); % 1st element, function of the mean value of y's
b = [b0 zeros(size(X,2)-1,1)]';  % column-vector of parameters

while 1==1
    z = X*b; % the logit^-1 variable is function of parameters
    p = 1./(1+exp(-z)); % recover the regression
    w = p.*(1-p); % calculate the weights of the samples
    u = z + (y-p)./w; % calculate the dependent variable for this step of least squares
    b_old = b; % store old parameters
    b = inv( X'*diag(w)*X ) * X' * diag(w) * u; % calculate new parameters with least squares
    % termanate the iterations if changes of the parameters are small
    if sumsqr(b - b_old) <= TolFun, break; end
end

%% Logistic regression example
% The Newton-Raphson algorithm is used to obtain the optimal parameters
% of the regression model
 
%% Create a demonstration data set
 
% independent variable, 20 samples
x = [[-8:1]'; [2:11]'];
% dependent variable of zeros and ones
y = [zeros(9,1); 1; 0; ones(9,1)];
% construct the matrix of independent samples
X = [ones(size(x,1),1) x];
 
%Plot the initial data
 
h = figure;
hold on
plot(X(:,2), y,'k*'); hold on
txtlegend = {'initial data'};
colors = {'r-','g-', 'b-', 'c-', 'm-', 'y-'};
ncolor  = 0;
 
%% Set the constant for iteration convergence
 
% the algorithm stops when the difference of the parameter is small
TolFun = 10^-3;
 
%% Define the initial value of the parameters
 
% 1st element, function of the mean value of y's
b0 = log(mean(y)/(1-mean(y)));
% column-vector of parameters
b = [b0 zeros(size(X,2)-1,1)]';
 
%%  The Newton-Raphson procedure
 
while 1==1
    % the logit^-1 variable is function of parameters
    z = X*b;
    % recover the regression
    p = 1./(1+exp(-z));
    % calculate the weights of the samples
    w = p.*(1-p);
    % calculate the dependent variable for this step of least squares
    u = z + (y-p)./w;
 
    % plot the results of this step
    plot(x, p, colors{mod(ncolor,length(colors))+1});
    ncolor = ncolor + 1; % change color
    txtlegend{end+1} = ['iteration ', num2str(ncolor)];
 
    % store old parameters
    b_old = b;
    % calculate new parameters with least squares
    b = inv( X'*diag(w)*X ) * X' * diag(w) * u;
 
    % termanate the iterations if changes of the parameters are small
    if sumsqr(b - b_old) <= TolFun, break; end
end
 
%%  Show the result
 
% claculate recovered dependent variable
p = 1./(1+exp(-X*b));
% plot the result
txtlegend{end+1} = 'recovered data';
plot(x, p,'rs');
legend(txtlegend);
axis tight
xlabel('x');
ylabel('p = (1+e^z)^{-1}, x = b_0+b_1x');
 
%% Notes
% In Matlab there are glmfit and glmval functuions. Use them for
% professional purposes.
 
return

Категории: Регрессионный анализ | Популярные и обзорные статьи | Библиотеки алгоритмов