function [condl, condu] = cond(X, theta, sigmal, sigmau)
% Lower and upper likelihoods for set of Gaussians (variance intervals)
% Author: Fabio Cozman 1/13/97
% 
% Calculates lower and upper bounds for a likelihood based 
% based on the Gaussian distribution (but considering an 
% interval of variance).
%
% The model assume a vector theta with the possible values
% for the parameter theta. The likelihoods are of the form:
%
%                                  1
%         p( X | theta ) = ---------------- Exp[ - (X - theta)^2/(2 sigma) ]
%                          Sqrt[2 Pi sigma]
%
% For a fixed theta and a measurement X, the function is:
%
%                        k1
%         f(sigma) = ----------- Exp[ - k2/sigma]
%                    Sqrt[sigma]
%
% This function has a single maximum at sigma = 2 k2 = (X - theta)^2
% (because f(sigma)->0 as sigma->0 and sigma->infinity, and
% the first derivative has a single zero at 2k2, and the second
% derivative is negative at 2k2).
% First derivative of f(sigma):
%                                k1 (2 k2 - s)
%                          ------------------------
%                          2 s^2 Sqrt[s] Exp[-k2/s]
% Second derivative of f(sigma):
%                                k1(4 k2^2 - 12 k2 s + 3 s^2)
%                                ----------------------------
%                                  4 s^4 Sqrt[s] Exp[-k2/s]
%
% Consequently, to find the maximum and minimum of f(sigma) in the 
% interval [ sigmal, sigmau ], it is necessary to look at the values
% of sigma: sigmal, 2k2, sigmau.
% The procedure is to check the values and fill the vectors condl and condu.

disp('Lower and upper likelihoods from Gaussians (by Fabio Cozman)');

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

k2 = 0.0;
size_theta = size(theta);
condl = zeros(size_theta);
condu = zeros(size_theta);

values = [ 0.0, 0.0 ];
vmax = [ 0.0, 0.0 ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Build the conditional bounds

for I = 1:size_theta(2)
	values(1) = normpdf(X, theta(I), sigmal);
	values(2) = normpdf(X, theta(I), sigmau);

	condl(I) = min(values);

	vmax(1) = max(values);
	k2 = (X - theta(I))^2/2;
	sigma_star = 2 * k2;
	if ( sigma_star < sigmal ) | ( sigma_star > sigmau )
		condu(I) = vmax(1);
	else
 		vmax(2) = normpdf(X, theta(I), sigma_star);
		condu(I) = max(vmax);
	end

end

