function p = gaussian_cdf(x)
% GAUSSIAN_CDF Cumulative distribution function of the standard Normal.
% p = gaussian_cdf(x)
%
% p = Phi(x) = Pr(Z <= x), where Z ~ N(0,1)
% 
% erf(x) = 2 Phi(x sqrt(2))-1 => Phi(x) = 0.5(erf(x/sqrt(2))+1)

% Q(x)    = Pr(X > x) = 1/sqrt(2) 1/sqrt(pi) int_x^infty exp(-z^2/2) dz
% erfc(x) = 2/sqrt(pi) int_x^infty exp(-t^2) dt
% Hence using the substitution z = sqrt(2) t (so dz = sqrt(2) dt) we get
% Q(x)    = 1/sqrt(2) 1/sqrt(pi) int_x^infty exp(-t^2) sqrt(2) dt
%         = 1/2 * erfc(x/sqrt(2))
%
% P(x) = 1 - Q(x) = 1 - 0.5*erfc(x/sqrt(2))
% but erfc(x) = 1 - erf(x) so P(x) = 0.5 + 0.5*erf(x/sqrt(2))

p = 0.5 + 0.5*erf(x/sqrt(2));
