Gaussian pdf in Julia

Also switched names of parameters from mu to μ and std to σ.

Using math_const values will allow for creating a templated type that could do evaluations in BigFloat types, should someone want to do so.

src/univariate/normal.jl

@@ -1,47 +1,57 @@
 immutable Normal <: ContinuousUnivariateDistribution
-    mean::Float64
-    std::Float64
+    μ::Float64
+    σ::Float64
     function Normal(mu::Real, sd::Real)
     	sd > zero(sd) || error("std must be positive")
     	new(float64(mu), float64(sd))
     end
+    Normal(mu::Real) = new(float64(mu), 1.0)
+    Normal() = new(0.0, 1.0)
 end
 
-Normal(mu::Real) = Normal(mu, 1.0)
-Normal() = Normal(0.0, 1.0)
+@_jl_dist_2p Normal norm
 
-const Gaussian = Normal
+@Base.math_const twopi 6.283185307179586 big(2.)*pi
+@Base.math_const log2pi 1.8378770664093456 log(big(twopi))
+@Base.math_const sqrt2pi 2.5066282746310007 sqrt(big(twopi))
+@Base.math_const logsqrt2pi 0.9189385332046728 0.5*big(log2pi)
 
-@_jl_dist_2p Normal norm
+zval(d::Normal, x::Real) = (x - d.μ)/d.σ
 
-entropy(d::Normal) = 0.5 * log(2.0 * pi) + 0.5 + log(d.std)
+phi{T<:FloatingPoint}(z::T) = exp(-0.5*z*z)/sqrt2pi
 
-insupport(d::Normal, x::Number) = isreal(x) && isfinite(x)
+logphi{T<:FloatingPoint}(z::T) = -(0.5*z*z + logsqrt2pi)
+
+pdf(d::Normal, x::FloatingPoint) = phi(zval(d,x))/d.σ
+pdf(d::Normal, x::Integer) = pdf(d, float64(x))
+
+logpdf(d::Normal, x::FloatingPoint) = logphi(zval(d,x)) - log(d.σ)
+logpdf(d::Normal, x::Integer) = logphi(zval(d,float64(x))) - log(d.σ)
+
+const Gaussian = Normal
+
+entropy(d::Normal) = logsqrt2pi + 0.5 + log(d.σ)
+
+insupport(d::Normal, x::Real) = isfinite(x)
 
 kurtosis(d::Normal) = 0.0
 
-mean(d::Normal) = d.mean
+mean(d::Normal) = d.μ
 
-median(d::Normal) = d.mean
+median(d::Normal) = d.μ
 
-function mgf(d::Normal, t::Real)
-	m, s = d.mean, d.std
-	return exp(t * m + 0.5 * s^t * t^2)
-end
+mgf(d::Normal, t::Real) = exp(t * d.μ + 0.5 * d.σ^t * t^2)
 
-function cf(d::Normal, t::Real)
-	m, s = d.mean, d.std
-	return exp(im * t * m - 0.5 * s^t * t^2)
-end
+cf(d::Normal, t::Real) = exp(im * t * d.μ - 0.5 * d.σ^t * t^2)
 
-modes(d::Normal) = [d.mean]
+modes(d::Normal) = [d.μ]
 
-rand(d::Normal) = d.mean + d.std * randn()
+rand(d::Normal) = d.μ + d.σ * randn()
 
 skewness(d::Normal) = 0.0
 
-std(d::Normal) = d.std
+std(d::Normal) = d.σ
 
-var(d::Normal) = d.std^2
+var(d::Normal) = d.σ^2
 
 fit_mle{T <: Real}(::Type{Normal}, x::Array{T}) = Normal(mean(x), std(x))