function [bounds, theta, pl, pu, alphal, alphau] = ex()
% Example for Quasi-Bayesian inferences with fractional programming.
% Author: Fabio Cozman 1/13/97
%
% Calculation of expected value interval for a parameter
% theta, with:
%
% prior defined by an epsilon-contaminated class
%    (1-epsilon) Beta(theta,15,11) + epsilon q(theta)
% for a Beta distribution translated to the interval 120-220.
%
% likelihood defined by the convex hull of a set of Gaussians
%    N(X; theta, sigma) where sigma in [sigmal, sigmau]

disp('Example: Quasi-Bayesian inference through fractional programming');
disp(' Author: Fabio Cozman');
disp(' Date: 1/13/97');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basic constants

thetaO = 120;
thetaF = 180;

epsilon = 0.15;
beta_a = 11;
beta_b = 15;

sigmal = 10;
sigmau = 11;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Measurements

X1 = 156.2684;
X2 = 150.7508;
X3 = 153.5161;
X4 = 143.0349;
X5 = 166.9614;
X6 = 150.5906;
X7 = 157.5622;
X8 = 154.0049;
X9 = 136.5862;
X10 = 153.7504;
X11 = 161.2516;
X12 = 157.2864;
X13 = 126.2255;
X14 = 147.2622;
X15 = 146.7706;
X16 = 153.1799;
X17 = 144.8883;
X18 = 149.9796;
X19 = 166.0651;
X20 = 158.4765;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construction of the vector theta

theta = thetaO:thetaF;
number_of_theta = thetaF - thetaO;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construction of the function that generates expected value

f = theta;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construction of the prior

A = - eye(number_of_theta + 1);
B = ( - (1-epsilon)/(thetaF - thetaO) * ...
	 betapdf( (theta - thetaO)/(thetaF - thetaO), beta_a, beta_b) )';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construction of the likelihood

[L1, U1] = cond(X1, theta, sigmal, sigmau);
[L2, U2] = cond(X2, theta, sigmal, sigmau);
[L3, U3] = cond(X3, theta, sigmal, sigmau);
[L4, U4] = cond(X4, theta, sigmal, sigmau);
[L5, U5] = cond(X5, theta, sigmal, sigmau);
[L6, U6] = cond(X6, theta, sigmal, sigmau);
[L7, U7] = cond(X7, theta, sigmal, sigmau);
[L8, U8] = cond(X8, theta, sigmal, sigmau);
[L9, U9] = cond(X9, theta, sigmal, sigmau);
[L10, U10] = cond(X10, theta, sigmal, sigmau);
[L11, U11] = cond(X11, theta, sigmal, sigmau);
[L12, U12] = cond(X12, theta, sigmal, sigmau);
[L13, U13] = cond(X13, theta, sigmal, sigmau);
[L14, U14] = cond(X14, theta, sigmal, sigmau);
[L15, U15] = cond(X15, theta, sigmal, sigmau);
[L16, U16] = cond(X16, theta, sigmal, sigmau);
[L17, U17] = cond(X17, theta, sigmal, sigmau);
[L18, U18] = cond(X18, theta, sigmal, sigmau);
[L19, U19] = cond(X19, theta, sigmal, sigmau);
[L20, U20] = cond(X20, theta, sigmal, sigmau);

L = L1 .* L2 .* L3 .* L4 .* L5 .* L6 .* L7 .* L8 .* L9 .* L10 .* ...
	L11 .* L12 .* L13 .* L14 .* L15 .* L16 .* L17 .* L18 .* L19 .* L20; 
U = U1 .* U2 .* U3 .* U4 .* U5 .* U6 .* U7 .* U8 .* U9 .* U10 .* ...
	U11 .* U12 .* U13 .* U14 .* U15 .* U16 .* U17 .* U18 .* U19 .* U20;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of the posterior

[bounds, alphal, alphau] = qb2(f, A, B, L, U);

pl =  L .* alphal(1:(number_of_theta+1))' + ...
       U .* alphal((number_of_theta+2):(number_of_theta+number_of_theta+2))' ;

pu =  L .* alphau(1:(number_of_theta+1))' + ...
       U .* alphau((number_of_theta+2):(number_of_theta+number_of_theta+2))' ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of the prior

ppl = (alphal(1:(number_of_theta+1)) + ...
	alphal((number_of_theta+2):(number_of_theta+number_of_theta+2)))/ ...
		sum(alphal);

ppu = (alphau(1:(number_of_theta+1)) + ...
	alphau((number_of_theta+2):(number_of_theta+number_of_theta+2)))/ ...
		sum(alphau);


