%
% Classical and Bayesian inference for the Behrens Fisher problem.
%
% X = N(mux,1/phix), Y = N(muy,1/phiy)
%
% 2 independent samples (nx,meanx,sdx), (ny,meany,sdy)
% We wish to carry out inference concerning mux-muy.
%
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
% PROGRAM
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
%
% Sample data.  We will suppose that the summary statistics for each sample
% are available.
%
nx=16;
meanx=1;
sdx=100;
ny=9;
meany=2;
sdy=1;
%
% Confidence level for interval. This should be fixed.
%
alpha=0.05;
%
% Combined variance estimator.
%
sdcomb=sqrt(sdx*sdx/nx+sdy*sdy/ny);
%
% Classical approximate 100(1-alpha)% interval
%
r=(sdcomb^4)/((sdx^4)/((nx^2)*(nx-1))+(sdy^4)/((ny^2)*(ny-1)));
tclass=icdf('t',(1-alpha*0.5),r);
intclass=[meanx-meany-sdcomb*tclass, meanx-meany+sdcomb*tclass]
%
% Bayesian interval given Jeffreys priors.
% Uses the Patil approximation to the Behrens Fisher distribution.
%
omega=atan((sdx/sdy)*sqrt(ny/nx));
f11=((nx-1)/(nx-3))*(sin(omega)^2);
f12=((ny-1)/(ny-3))*(cos(omega)^2);
f1=f11+f12;
f2=(f11^2)/(nx-5)+(f12^2)/(ny-5);
b=4+(f1^2)/f2;
a=sqrt(f1*(b-2)/b);
tbayes=icdf('t',(1-alpha*0.5),b)*a;
intbayes=[meanx-meany-sdcomb*tbayes, meanx-meany+sdcomb*tbayes]
%
% Monte Carlo version of Bayesian estimator.  This could be compared with
% Patil's interval.
%
montesize=100000;
mux=meanx+sdx*trnd((nx-1),1,montesize)/sqrt(nx);
muy=meany+sdy*trnd((ny-1),1,montesize)/sqrt(ny);
delta=sort(mux-muy);
intbayesmonte=[delta(montesize*alpha*0.5),delta(montesize*(1-0.5*alpha))]
