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