function [bounds, alphal, alphau] = qb2(f, A, B, L, U)
% Linear fractional programming for Quasi-Bayesian inferences 
% Author: Fabio Cozman 1/12/97
%
% Calculates lower and upper bounds for the expected value of 
% the function f given linear constraints on probability.
%
% A finite universe is assumed; a vector theta is assumed to
% contain all the possible values in this universe. The function 
% f is specified by a vector f; values for the function are obtained 
% by the dot product f . theta.
%
% The vector of probability values for the values in theta is
% denoted by ptheta.
%
% The constraints in prior probability are defined as:
%
%                         A ptheta <= B
%
% Additional constraints are imposed by the axioms of probability:
%
%                         \sum ptheta_i = 1
%                         ptheta_i >= 0
%
% The likelihood is given by the vectors L (lower likelihood)
% and U (upper likelihood). 
%
% The lower and upper probability envelopes are produced in
% the vectors alphal and alphau respectively.
%
% The optimization problem that is solved is:
%
%                          \sum_i  f_i alpha_i [L|U]_i
%            max/min  =  ------------------------------
%                          \sum_j   alpha_j [L|U]_j
% 
% The method of solution is through the Charnes-Cooper transform
% and then application of linear programming.

disp('Quasi-Bayesian inferences (by Fabio Cozman)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization

% Check the size of the inputs

[ widthF, size_of_thetaF ] = size(f);
[ number_of_constraintsA, size_of_thetaA ] = size(A);
[ number_of_constraintsB, widthB ] = size(B);
[ widthL, size_of_thetaL ] = size(L);
[ widthU, size_of_thetaU ] = size(U);

if (widthF ~= 1)
	disp('Function f must be defined by a row vector!');
	return;
end
if (widthL ~= 1) 
	disp('Lower likelihood must be defined by a row vector!');
	return;
end
if (widthU ~= 1)
	disp('Upper likelihood must be defined by a row vetor!');
	return;
end
if (number_of_constraintsA ~= number_of_constraintsB)
	disp('Matrices A and B must have the same number of rows!');
	return;
else
	number_of_constraints = number_of_constraintsA;
end
if (size_of_thetaA ~= size_of_thetaF) |  ...
   (size_of_thetaA ~= size_of_thetaL) |  ...
   (size_of_thetaA ~= size_of_thetaU) 
	disp('The vectors f, L and U, and matrix A must use the same');
	disp('	number of values for theta!');
	return
else
	size_of_theta = size_of_thetaA;
end

% Obtain the matrix C = A - B 1
% The constrains are reduced to C ptheta <= 0
%	(given the unitary constraint)

C = A - B * ones(1, size_of_theta);

% Adjust the unitary constrain (Charnes-Cooper transform)
%	and enforce the constraint through matrix C

Cunity(2:(number_of_constrains+1),:) =[ C, C ];
Cunity(1,:) = [ L, U ];

Bunity = zeros(number_of_constraints+1, 1);
Bunity(1) = 1;

% Adjust the objective function (Charnes-Cooper transform)

objective = [ f .* L, f .* U ];

% Generate the lower bounds for the linear programming routine

vlb = zeros(size(objective));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of lower bound
alphal = lp(objective, Cunity, Bunity, vlb, [], [], 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of upper bound
alphau = lp(-objective, Cunity, Bunity, vlb, [], [], 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generation of expected value interval
bounds = [ objective * alphal , objective * alphau ]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%










