function  [mu,sigma,table] = truncGauss(a,table)

%TRUNCGAUSS -- Mean and std.dev. of a truncated Gaussian variable.
%
%  [mu,sigma] = truncGauss(a)
%  [mu,sigma] = truncGauss(a,table)
%
%  Let X is a standard normal variable (cf. NORMPDF, NORMCDF, and ERF)
%  and construct Y = max(A,X), where A is a constant. The new variable
%  has mean MU greater than zero and std.dev. SIGMA less than 1.0.
%  The current function tabulates MU and SIGMA as a function of A.
%
%  The values are interpolated from a table built via Monte Carlo methods.
%  To generate the table, use the call  [mu,sigma]=truncGauss(a,N) ,
%  where A is a vector defining the grid (e.g. [-3:0.1:+3]) and
%  N is an integer specifying the number of Monte-Carlo samples.
%  The TABLE is a structure with three fields: 'grid','mu',and 'sigma'.
%  It is stored in a persistent variable within the local workspace
%  so subsequent calls to TRUNCGAUSS can re-use it.  Alternatively,
%  TABLE can be explicitly passed to the function or loaded from a
%  file. In this latter case, the second argument must be a string
%  indicating a filename. The default is 'truncGauss.mat'.
%
%  See also INTERP1, RANDN, NORMCDF, NORMPDF, NORMINV, ERF, ERFINV.

% The exact formula for MU is:   mu(a) = a.*normcdf(a) + normpdf(a)
% But beware: the implementations of NORMCDF, ERF, etc. involve
% approximations and are quite costly. Therefore I decided to tabulate.
%
% I didn't calculate a closed formula for SIGMA, but it seems that
% for A<0, the following is an excellent approximation
%   sigma(a) = 1 - 0.8322.*normcdf(a)
% The approximation was derived by trial and error. It overestimates
% SIGMA when A>0.  In particular, for A=+Inf the true SIGMA is 0,
% while the approximation gives 0.1678 = 1-0.8322.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2002/07/01 19:05 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

persistent  TRUNCGAUSS_TABLE

%-- Make sure there is a valid table, one way or another
if (nargin==1)    % nothing is supplied explicitly; the typical case
   if (isempty(TRUNCGAUSS_TABLE))
      load('truncGauss.mat') ;   % load TRUNCGAUSS_TABLE from disk
   %else the persistent variable is already set and everything is OK
   end
elseif (ischar(table))           % file name?
   load(table) ;
elseif (isstruct(table))         % explicit structure?
   TRUNCGAUSS_TABLE = table ;    % update the persistent variable
elseif  (isnumeric(table))       % integer N ==> request to (re)build
   TRUNCGAUSS_TABLE = build_table(a,table) ;
else
   error('TABLE must be a struct, a filename, or an integer.') ;
end
table = TRUNCGAUSS_TABLE ;

%-- Interpolate for the requested values of A
grid = table.grid ;
mu = interp1(grid,table.mu,a,'*linear') ;
sigma = interp1(grid,table.sigma,a,'*linear') ;

%-- Handle out-of-range points
idx = find(a<grid(1)) ;    % too small ==> no cutoff ==> mu=0, sigma=1
mu(idx) = 0 ; sigma(idx) = 1 ;

idx = find(a>grid(end)) ;  % too large ==> all cutoff ==> mu=A, sigma=0
mu(idx) = a(idx) ; sigma(idx) = 0 ;

%-- Return MU and SIGMA
%%%%% End of the main function

%%%%%%%%%  SUBFUNCTION  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  table = build_table(grid,N)

disp('Building TRUNCGAUSS table...') ;
grid = grid(:) ;
n = length(grid) ;
mu = zeros(n,1) ;
sigma = zeros(n,1) ;

X = randn(N,1) ;
for k = 1:n
   Y = max(grid(k),X) ;
   mu(k) = mean(Y) ;
   sigma(k) = std(Y,1) ;
end

table.grid = grid ;
table.mu = max(mu,0) ;        % mean(truncated) must be >= mean(untrunc)
table.sigma = min(sigma,1) ;  % similarly, sigma(trunc) <= sigma(untrunc)
table.N = N ;

%%%%% End of subfunction BUILD_TABLE

%%%%% End of file TRUNCGAUSS.M
