function [Mb,accept]=b_logit(data,m,scale);



%B_LOGIT Produces a simulated sample from a logistic regression model

%        using the Metropolis algorithm.

%

%   [MB, ACCEPT]=b_logit(DATA,M,SCALE) returns a matrix MB of simulated parameter 

%   values and ACCEPT, the acceptance proportion of the algorithm, where DATA 

%   is the data matrix (first column the observed proportions, second column the

%   sample sizes, and remaining columns the design matrix), M is the number of 

%   iterations of the algorithm, and SCALE is the scale factor used in the

%   Metropolis-Hastings algorithm.



[beta,var]=logit(data);        % Finds maximum likelihood fit -- use as starting value



a=chol(var);                   % Takes cholesky decomposition



b=beta; p=length(beta);

Mb=zeros(m,p);                 % Sets up storage

lb=bl_like(data,b');           % Evaluates log likelihood at current value



if nargin<3

  scale=1;                     % Sets scaling factor of metropolis algorithm equal to 1

end



h=waitbar(0,'Simulation in progress');

accept=0;



for i=1:m

   

  bc=b+scale*a'*randn(p,1);    % Candidate value selected

  lbc=bl_like(data,bc');       % Evaluates log likelihood at candidate value

  prob=exp(lbc-lb);            % Computes acceptance probability

  if rand<prob                 % If random uniform is smaller tha probability, accept

    lb=lbc; b=bc; accept=accept+1; % candidate value

  end

 

  Mb(i,:)=b';                  % Store simulated value

 

  if i/20==fix(i/20)           % Shows wait bar

    figure(1)

    waitbar(i/m)

  end



end



accept=accept/m;               % Computes acceptance proportion

close(h)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function [beta,var]=logit(data)

%

% data matrix of form [y n x] - y is observed proportion

% sample size vector n 



s=size(data);

k=s(2);

y=data(:,1);

n=data(:,2);

x=data(:,3:k);

p=(y+.5)/2;

eta=log(p./(1-p));

w=diag(n.*p.*(1-p));

z=eta+(y-p)./(p.*(1-p));

b1=(x'*w*x)\(x'*w*z);

b0=zeros(size(b1));

%

% loop until convergence

%

it=0;

while (it<=10)*(norm(b1-b0)>.001*norm(b0))

	b0=b1;

	eta=x*b0;

	p=exp(eta)./(1+exp(eta));

	z=eta+(y-p)./(p.*(1-p));

   w=diag(n.*p.*(1-p));

	b1=(x'*w*x)\(x'*w*z);

	it=it+1;

end

beta=b1;

var=inv(x'*w*x);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function val=bl_like(data,beta)

% Binomial logistic loglikelihood



[npt,p]=size(beta);

y=data(:,1); x=data(:,3:(p+2)); n=data(:,2);

lp=x*beta';

val=sum(n.*y.*lp-n.*log(1+exp(lp)));

