% EXAMPLE2.m

% Summarizing posterior density for two-parameter problem

% [Albert (1996), Chapter 11]

%

% Have 2 x 2 contingency table with two unknown proportions p1 and p2.

% We observe the table [1 10;3 12]. If the proportions are transformed 

% to the log odds ratio W and the log odds ratio U, then the log posterior

% density with a noninformative prior is given by

%

% h(W,U)=T1-11 log(1+exp(T1))+3 T2-15 log(1+exp(T2))

%

% where T1=(W+U)/2, T2=(U-W)/2.  This log posterior is programmed in

% the MATLAB function logpost2.m



% Laplace method



data = [1 11 3 15];    % data consists of number of successes and sample size

                       % for each sample

mode = [1 1];          % starting guess at posterior mode

numiter = 10;          % number of iterations of Newton-Raphson



[mode,var,lint]=laplace('logpost2',mode,numiter,data)

                      % outputs posterior mode, estimate at posterior variance,

                      % and Laplace estimate at log integral



% Adaptive quadrature



mom=[0 1 0 1 0];      % estimate at E(W), SD(W), E(U), SD(U), corr(W,U)

numiter = 10;         % number of iterations of algorithm



[lint,mom]=ad_quad2('logpost2',mom,numiter,data)

                     % outputs posterior mean and variance,

                     % and value of log integral



% Gibbs sampling



start=[0 0];      % starting value for simulation

scale=[1 1];      % scale factors for random walk increment density

m=1000;           % number of iterations of algorithm



sim_sample=gibbs('logpost2',start,m,scale,data); 

                  % outputs simulated sample matrix from posterior density



