www.gusucode.com > stats 源码程序 matlab案例代码 > stats/LogisticGrowthModelExample.m

    %% Nonlinear Mixed-Effects Model
% Enter and display data on the growth of five orange trees.

% Copyright 2015 The MathWorks, Inc.

CIRC = [30 58 87 115 120 142 145;
        33 69 111 156 172 203 203;
        30 51 75 108 115 139 140;
        32 62 112 167 179 209 214;
        30 49 81 125 142 174 177];
time = [118 484 664 1004 1231 1372 1582];

h = plot(time,CIRC','o','LineWidth',2);
xlabel('Time (days)')
ylabel('Circumference (mm)')
title('{\bf Orange Tree Growth}')
legend([repmat('Tree ',5,1),num2str((1:5)')],...
       'Location','NW')
grid on
hold on
%%
% Use an anonymous function to specify a logistic growth model.
model = @(PHI,t)(PHI(:,1))./(1+exp(-(t-PHI(:,2))./PHI(:,3)));
%%
% Fit the model using |nlmefit| with default settings (that is, assuming
% each parameter is the sum of a fixed and a random effect, with no
% correlation among the random effects):
TIME = repmat(time,5,1);
NUMS = repmat((1:5)',size(time));

beta0 = [100 100 100];
[beta1,PSI1,stats1] = nlmefit(TIME(:),CIRC(:),NUMS(:),...
                              [],model,beta0)
%%
% The negligible variance of the second random effect, |PSI1(2,2)|,
% suggests that it can be removed to simplify the model.
[beta2,PSI2,stats2,b2] = nlmefit(TIME(:),CIRC(:),...
    NUMS(:),[],model,beta0,'REParamsSelect',[1 3])
%%
% The loglikelihood |logl| is unaffected, and both the Akaike and Bayesian
% information criteria ( |aic| and |bic| ) are reduced, supporting the
% decision to drop the second random effect from the model.
%%
% Use the estimated fixed effects in |beta2| and the estimated random
% effects for each tree in |b2| to plot the model through the data.
PHI = repmat(beta2,1,5) + ...          % Fixed effects
      [b2(1,:);zeros(1,5);b2(2,:)];    % Random effects

tplot = 0:0.1:1600;
for I = 1:5
  fitted_model=@(t)(PHI(1,I))./(1+exp(-(t-PHI(2,I))./ ... 
       PHI(3,I)));
  plot(tplot,fitted_model(tplot),'Color',h(I).Color, ...
	   'LineWidth',2)
end