Code covered by the BSD License  

Highlights from
New Regression Capabilities in r2012a

image thumbnail

New Regression Capabilities in r2012a

by

 

Example code to accompany the "New Regression Capabilities" presentation

GLMs.m
%% Using GeneralizedLinearModel for Logistic and Poisson Regression
% Copyright (c) 2012, The MathWorks, Inc.

%%  Generalized Linear Model Demo

Weight = [2100 2300 2500 2700 2900 3100 ...
     3300 3500 3700 3900 4100 4300]';
Tested = [48 42 31 34 31 21 23 23 21 16 17 21]';
Failed = [1 2 0 3 8 8 14 17 19 15 17 21]';

myFit = GeneralizedLinearModel.fit(Weight, [Failed Tested], ...
    'y ~ x1', 'distr', 'binomial')

%% Show the complete set of methods

methods(myFit)

%%

clear all
clc

% Create a set of 100 X values between 0 and 15

X = linspace(0,15,100);
X = X';

% Define a function that describes the average relationship between X and
% Y.  

mu = @(x) exp(-1.5 +.25*x);
Y = mu(X);

plot(X,Y); 
xlabel('Time');
ylabel('Count');

%%  Use the Y = f(X) to generate a dataset with known characteristics

%  At each "X", the observed Y is represented as a draw from a Poisson
%  distribution with mean (lambda) = Y

Noisy_Y = poissrnd(Y,100,1);

% Create figure
figure1 = figure;

% Create axes
axes1 = axes('Parent',figure1,'YTick',[0:max(Noisy_Y)], 'YGrid','on');
hold(axes1,'all');

% Create scatter
scatter(X,Noisy_Y, 'b');

xlabel('Time')
ylabel('Count')

% Note the following:

% The mean of a poisson process is also the variance of a poisson process.
% The blue curve doesn't exhibit constant variance (this violates one of
% the key assumptions underlying Ordinary Least Squares)

% The output from a Poisson process is always a non-negative integer.  If
% lambda is close to zero AND observations can never going
% non-negative, the shape of the distribution is bounded in one direction.
% As lambda gets larger, this boundary has less and less impact on the
% distribution of the errors.

hold on
plot(X,Y); 

%% Generate a Poisson Regression Model

myFit = GeneralizedLinearModel.fit(X, Noisy_Y,'log(y) ~ x1','distr','poisson');
disp(myFit)

%% Plot the residuals

plotResiduals(myFit, 'Fitted')

Contact us