truncatedGaussian.m

function [X meaneffective sigmaeffective] = truncatedGaussian(sigma, range, varargin)
% function X = TruncatedGaussian(sigma, range)
%          X = TruncatedGaussian(sigma, range, n)
%
% Purpose: generate a pseudo-random vector X of size n, X are drawn from
% the truncated Gaussian distribution in a RANGE braket; and satisfies
% std(X)=sigma.
% RANGE is of the form [left,right] defining the braket where X belongs.
% For a scalar input RANGE, the braket is [-RANGE,RANGE].
%
% X = TruncatedGaussian(..., 'double') or
% X = TruncatedGaussian(..., 'single') return an array X of of the
%    specified class.
%
% If input SIGMA is negative, X will be forced to have the same "shape" of
% distribution function than the unbounded Gaussian with standard deviation
% -SIGMA: N(0,-SIGMA). It is similar to calling RANDN and throw away values
% ouside RANGE. In this case, the standard deviation of the truncated
% Gaussian will be different than -SIGMA. The *effective* mean and
% the effective standard deviation can be obtained by calling:
%   [X meaneffective sigmaeffective] = TruncatedGaussian(...)
%
% Example:
% 
% sigma=2;
% range=[-3 5]
% 
% [X meaneff sigmaeff] = TruncatedGaussian(sigma, range, [1 1e6]);
% 
% stdX=std(X);
% fprintf('mean(X)=%g, estimated=%g\n',meaneff, mean(X))
% fprintf('sigma=%g, effective=%g, estimated=%g\n', sigma, sigmaeff, stdX)
% hist(X,64)
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Last update: 19/April/2009
%              12-Aug-2010, use asymptotic formula for unbalanced
%                           range to avoid round-off error issue

% We keep track this variables so as to avoid calling fzero if
% TruncatedGaussian is called succesively with the same sigma and range
persistent PREVSIGMA PREVRANGE PREVSIGMAC

% shape preserving?
shapeflag = (sigma<0);

% Force inputs to be double class
range = double(range);
if isscalar(range)
    % make sure it's positive
    range=abs(range);
    range=[-range range];
else
    range=sort(range); % right order
end
sigma = abs(double(sigma));

n=varargin;

if shapeflag
    % Prevent the same pdf as with the normal distribution N(0,sigma)
    sigmac = sigma;
else
    if diff(range)^2<12*sigma^2 % This imposes a limit of sigma wrt range
        warning('TruncatedGaussian:RangeSigmaIncompatible', ...
                'TruncatedGaussian: range and sigma are incompatible\n');
        sigmac = Inf;
    elseif isequal([sigma range], [PREVSIGMA PREVRANGE]) % See line #80
        sigmac = PREVSIGMAC; % Do not need to call fzero
    else
        % Search for "sigmac" such that the truncated Gaussian having
        % sigmac in the formula of its pdf gives a standard deviation
        % equal to sigma
        [sigmac res flag] = fzero(@scz,sigma,[],...
                                   sigma^2,range(1),range(2)); %#ok
        sigmac = abs(sigmac); % Force it to be positive
        if flag<0 % Someting is wrong
            error('TruncatedGaussian:fzerofailled', ...
                  'Could not estimate sigmac\n');
        end
        % Backup the solution
        [PREVSIGMA PREVRANGE PREVSIGMAC] = deal(sigma,range,sigmac);
    end
end

% Compute effective standard deviation
meaneffective=meantrunc(range(1), range(2), sigmac);
sigmaeffective=stdtrunc(range(1), range(2), sigmac);

% Inverse of the cdf functions
if isinf(sigmac)
    % Uniform distribution to maximize the standard deviation within the
    % range. It is like a Gaussian with infinity standard deviation
    if any(strcmpi(n,'single'))
        range = single(range);
    end
    cdfinv = @(y) range(1)+y*diff(range);
else
    c = sqrt(2)*sigmac;

    rn = range/c;
    asymthreshold = 4;
    
    if any(strcmpi(n,'single'))
        % cdfinv will be single class
        c = single(c);
        %e = single(e);
    end
    
    % Unbalanced range
    if prod(sign(rn))>0 && all(abs(rn)>=asymthreshold)
        % Use asymptotic expansion
        % based on a Sergei Winitzi's paper "A handly approximation for the
        % error function and its inverse", Feb 6, 2008.
        c = c*sign(rn(1));
        rn = abs(rn);
        left = min(rn);
        right = max(rn);
        
        a = 0.147;

        x2 = left*left;
        ax2 = a*x2;
        e1 = (4/pi+ax2) ./ (1+ax2);
        e1 = exp(-x2.*e1); % e1 < 3.0539e-008 for asymthreshold = 4

        x2 = right*right;
        ax2 = a*x2;
        e2 = (4/pi+ax2) ./ (1+ax2);
        e2 = exp(-x2.*e2); % e2 < 3.0539e-008 for asymthreshold = 4
        
        % Taylor series of erf(right)-erf(left) ~= sqrt(1-e2)-sqrt(1-e1)
        de =  -0.5*(e2-e1) -0.125*(e2-e1)*(e2+e1);
        
        % Taylor series of erf1 := erf(left)-1 ~= sqrt(1-e1)-1 
        erf1 = (-0.5*e1 - 0.125*e1^2);
        cdfinv = @(y) c*asymcdfinv(y, erf1, de, a);
        
    else
        e = erf(range/c);

        cdfinv = @(y) c*erfinv(e(1)+diff(e)*y);
    end
end

% Generate random variable
X = cdfinv(rand(n{:}));
% Clip to prevent some nasty numerical issues with of erfinv function
% when argument gets close to +/-1
X = max(min(X,range(2)),range(1));

return

end % TruncatedGaussian

%%
function x = asymcdfinv(y, erf1, de, a)
% z = erf(left) + de*y = 1 + erf1 + de*y, input argument of erfinv(.)
f = erf1 + de*y; % = z - 1; thus z = 1+f
% 1 - z^2 = -2f-f^2
l = log(-f.*(2 + f)); % log(-2f-f^2) = log(1-z.^2);
b = 2/(pi*a) + l/2;
x = sqrt(-b + sqrt(b.^2-l/a));
end % asymcdfinv

function m=meantrunc(lower, upper, s)
% Compute the mean of a trunctated gaussian distribution
if isinf(s)
    m = (upper+lower)/2;
else
    a = (lower/sqrt(2))./s;
    b = (upper/sqrt(2))./s;
    corr = sqrt(2/pi)*(-exp(-b.^2)+exp(-a.^2))./(erf(b)-erf(a));
    m = s.*corr;
end
end % vartrunc

function v=vartrunc(lower, upper, s)
% Compute the variance of a trunctated gaussian distribution
if isinf(s)
    v = (upper-lower)^2/12;
else
    a = (lower/sqrt(2))./s;
    b = (upper/sqrt(2))./s;
    if isinf(a)
        ea=0;
    else
        ea = a.*exp(-a.^2);
    end
    if isinf(b)
        eb = 0;
    else
        eb = b.*exp(-b.^2);
    end
    corr = 1 - (2/sqrt(pi))*(eb-ea)./(erf(b)-erf(a));
    v = s.^2.*corr;
end
end % vartrunc

function stdt=stdtrunc(lower, upper, s)
% Standard deviation of a trunctated gaussian distribution
arg = vartrunc(lower, upper, s)-meantrunc(lower, upper, s).^2;
%arg = max(arg,0);
stdt = sqrt(arg);
end % stdtrunc

function res=scz(sc, targetsigma2, lower, upper)
% Gateway for fzero, aim the standard deviation to a target value
res = vartrunc(lower, upper, sc) - targetsigma2 - ...
      meantrunc(lower, upper, sc).^2;
end % scz

% End of file TruncatedGaussian.m