V_NORMCDFLOG calculates log of Normal Cumulative Distribution function p=(x,m,s) Inputs: X Input data (vector or matrix) M Mean of Normal distribution [default 0] S Std deviation of Normal distribution [default 1] Outputs: P P = log(normcdf(X)); same size as X The routine gives accurate values even if X is large and negative
0001 function p=v_normcdflog(x,m,s) 0002 %V_NORMCDFLOG calculates log of Normal Cumulative Distribution function p=(x,m,s) 0003 % 0004 % Inputs: 0005 % 0006 % X Input data (vector or matrix) 0007 % M Mean of Normal distribution [default 0] 0008 % S Std deviation of Normal distribution [default 1] 0009 % 0010 % Outputs: 0011 % 0012 % P P = log(normcdf(X)); same size as X 0013 % 0014 % The routine gives accurate values even if X is large and negative 0015 0016 % Copyright (C) Mike Brookes 2016 0017 % Version: $Id: v_normcdflog.m 10865 2018-09-21 17:22:45Z dmb $ 0018 % 0019 % VOICEBOX is a MATLAB toolbox for speech processing. 0020 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0021 % 0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0023 % This program is free software; you can redistribute it and/or modify 0024 % it under the terms of the GNU General Public License as published by 0025 % the Free Software Foundation; either version 2 of the License, or 0026 % (at your option) any later version. 0027 % 0028 % This program is distributed in the hope that it will be useful, 0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0031 % GNU General Public License for more details. 0032 % 0033 % You can obtain a copy of the GNU General Public License from 0034 % http://www.gnu.org/copyleft/gpl.html or by writing to 0035 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0037 persistent a b 0038 if isempty(a) 0039 a=0.996; 0040 % f=@(x) -0.5*(x.^2+log(2*pi))-real(log(-(a+x.^2)./x))-real(log(normcdf(x))); 0041 % b=fzero(f,-22); % cutoff value for conventional formula 0042 b=-22.2491306156561; % precalculated value 0043 end 0044 if nargin>2 0045 x=(x-m)/s; 0046 elseif nargin>1 0047 x=x-m; 0048 end 0049 t=x<b; % mask for large negative values 0050 p=zeros(size(x)); 0051 p(~t)=real(log(0.5*erfc(-x(~t).*sqrt(0.5)))); % use this formula normally 0052 p(t)=-0.5*(x(t).^2+log(2*pi))-real(log(-x(t)-a./x(t))); % use this approximation for large negative x 0053