Code covered by the BSD License

# New Regression Capabilities in r2012a

### Richard Willey (view profile)

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')

```