% This is an approximation to what I did in the lab on October 8.

% Problem 1 - I-75 problem

% first I define a new function, called i75.m, which computes the
% log posterior -- here is the function:

type i75

function val=i75(m,data)

y1=data(1); y2=data(2); s=5;
z=(70-m)/s;
val=y1*log(phi(z))+y2*log(1-phi(z));

% define data

data=[2 5];

i75(70,data)
ans =
   -4.8520

% use adaptive quadrature - input to ad_quad1 are
% name of log posterior function, guesses at mean and standard deviation,
% number of iterations, and data

[a,b,c]=ad_quad1('i75',[70 2],10,data)
a =
   -2.3398
b =
   73.0039    2.5563
c =
  Columns 1 through 7 
   60.5814   63.8477   66.6531   69.2564   71.7642   74.2435   76.7513
  Columns 8 through 10 
   79.3546   82.1600   85.4263

% output is estimate at log integral, posterior mean and stand deviation,
% and grid where density is computed

% try laplace approximation

% input name of function, guess at mode, number of iterations, and data

[a,b,c]=laplace('i75',70,10,data)
a =
   72.8297
b =
    6.3085
c =
   -2.3480

% here output is posterior mean, posterior variance, and estimate at log integral

% by graphing posterior, can see that normal approximation is pretty good

x=linspace(60,80,100); f=exp(i75(x,data)); plot(x,f)

% Problem 3

% log posterior for analysis using a t(4) sampling density is stored in 
% Matlab function tlike2.m

% enter data into Matlab

data=[49 -67 8 16 6 23 28 41 14 29 56 24 75 60 -48];

% I guess that posterior mode of (m, log s) is (24, 2)

[m,v]=laplace('tlike2',[24 2],10,data)
m =
   26.4619    3.2642
v =
   61.8676   -0.1503
   -0.1503    0.0648

% contrast posterior mean with sample mean, which is posterior mean with normal errors

mean(data)
ans =
   20.9333

diary off