%% load in the data from workplacesmoking.xlsx %% column 1 is the dependent variable while %% columns 2-(k+1) are independent variables %% column 2 contains the constant [w,varlist]=xlsread('workplacesmoking.xlsx'); dv=varlist{1,1}; %% get dimension of w nk1=size(w); % number of observations n=nk1(1); % number of independent variables % k+1 is the no of columns of w. there are k % covariates (including the constant) k=nk1(2)-1; %% extract y which is in the 1st column of w y=w(:,1); %% extract x which is in columns 2 through k x=w(:,2:(k+1)); %% generate starting values--which are from %% an OLS model xpxi=inv(x'*x); beta_start=xpxi*x'*y; %% start quasi-newton maximum likelihood %% set up initial values cc=1e-8; %%convergence criteria maxit=200; %%maximum number of interations beta=beta_start; logli=logitloglike(y,x,beta); grad=logitgrad(y,x,beta); hess=logithess(y,x,beta); cov=inv(-hess); %% open file to print out results file1=fopen('logit_mle_output2.txt','w'); fprintf(file1,'----------------------------------------------------------------\n'); %% now iterate to obtain solution for it=1:maxit; betalast=beta; loglast=logli; fprintf(file1,'iteration number %4.f, loglikelihood = %9.4f \n' ,it,logli); step=1; beta=betalast+step*cov*grad; disp(it) logli=logitloglike(y,x,beta); j=0; display(logli) display(loglast) %% if the log likelihood has not improved %% keep cutting step size in half if logli=loglast, break, end; end; end; display(j) if j>=10, break, end; grad=logitgrad(y,x,beta); hess=logithess(y,x,beta); cov=inv(-hess); criti=grad'*cov*grad if criti=10; fprintf(file1,'----------------------------------------------------------------\n'); fprintf(file1,'warning, after 10 halvings of correction vector, loglikelihood \n'); fprintf(file1,'has not improved. the results are not reliable \n'); fprintf(file1,'----------------------------------------------------------------\n'); end; if it>=maxit && criti>cc; fprintf(file1,'----------------------------------------------------------------\n'); fprintf(file1,'warning, maximum iterations reached and program has not \n'); fprintf(file1,'converged. the results are not reliable \n'); fprintf(file1,'----------------------------------------------------------------\n'); end; %% print out final results ym=sum(y)/n; dof=n-k; stderr=sqrt(diag(cov)); zscore=beta./stderr; pvalue=2*(1-normcdf(zscore,0,1)); fprintf(file1,'----------------------------------------------------------------\n'); fprintf(file1,'MLE estimates of Poisson model \n'); fprintf(file1,'----------------------------------------------------------------\n'); fprintf(file1,'Dependent variable = %s \n' ,dv); fprintf(file1,'Mean of dep variable = %9.4f \n' ,ym); fprintf(file1,'----------------------------------------------------------------\n'); fprintf(file1,'analysis of MLE procedure \n'); fprintf(file1,'\n'); fprintf(file1,'observations = %9.f \n' ,n); fprintf(file1,'parameters = %9.f \n' ,k); fprintf(file1,'degrees of freedom = %9.f \n' ,dof); fprintf(file1,'number of iterations = %9.f \n' ,it); fprintf(file1,'final convergence crit = %e \n' ,criti); fprintf(file1,'value of log likeliood = %9.4f \n' ,logli); c1='Covariate'; c2='Beta'; c3='Std Error'; c4='zscore'; c5='P-value: B=0'; fprintf(file1,'----------------------------------------------------------------\n'); fprintf(file1,'%12s %12s %12s %12s %12s \n', c1,c2,c3,c4,c5); fprintf(file1,'----------------------------------------------------------------\n'); for i=1:k; rowname=varlist{1,i+1}; fprintf(file1,'%12s %12.6f %12.6f %12.6f %12.6f \n', rowname,beta(i,1),stderr(i,1),zscore(i,1),pvalue(i,1)); end; fprintf(file1,'----------------------------------------------------------------\n'); fclose(file1);