% Nonlinear Model example April 5, 1999

global T r N

% Data generation

T = (0:1:50)'; T(1,1)=.01; % delay of test

% c = 2.046; d = .2698; % pareto II model parameters

% R = 1./(1 + c.*T).^d; % Proportion Recalled by Pareto-II model

c = 5.26; d = .289; % weibull model parameters

R = 1./exp((c.*T./1000).^d); % Proportion Recalled by Weibull

N = 250; % sample size of proportion

e = sqrt((R).*(1-R)./N); % error std

e = e.*randn(size(T,1),1); % error in est of proportion

r = R + e; % sample est of proportion

plot(T,r);

max = corrcoef([R,r]).^2;

'R^2 between true and sample est'

max(1,2)

u_chi = -sum(r.*log(r) + (1-r).*log(1-r));

parm = [1; 1]; % starting values

A = [ -1 0 ; 0 -1]; b = [ 0 ; 0 ]; % inequality constraints

'initial chi'

2.*N.*( f_eval(parm)- u_chi)

options=optimset('LargeScale','off');

[parm, chi] = fmincon('f_eval',parm,A,b,[],[],[],[],[],options);

'parm and chi-square'

[ parm ; 2.*N.*(chi - u_chi)]

% [parm, chi] = fminsearch('f_eval',parm)

function chi = f_eval(parm)

% evaluate model fit

% chi = f_eval(parm)

global T r N

c = parm(1,1);

d = parm(2,1);

p = 1./exp((c.*T./1000).^d);

chi = -sum( r.*log(p) + (1-r).*log(1-p) );

% chi = ((r-p)'*(r-p))./((r-mean(r))'*(r-mean(r)));