%% this subroutine calculates the hessian
%% of the poisson log likelikood
%% inputs are y,x,beta,n,k and the function
%% returns a (kxk) matrix of 2nd derivatives
function [hess]=calchess(y,x,beta)
nk1=size(x);   %% get dimensions of x
k=nk1(2);      %% no. of parameters
n=nk1(1);      %% no. of observations
xbeta=x*beta;
exp_xbeta=exp(xbeta);
prob=exp_xbeta./(1+exp_xbeta);
hess=zeros(k,k);
for i=1:n;
xi=x(i,:);
probi=prob(i,:);
hess=hess-xi'*xi.*probi.*(1-probi);
end;
end