mvnormal.jl

srand(12345)

function makeLtri(x)
    l = length(x)
    p = round(Int, (sqrt(1 + 8l) - 1) / 2)
    L = zeros(eltype(x), (p, p))
    k = 1
    for j in 1:p
        for i in j:p
            L[i, j] = x[k]
            k += 1
        end
    end
    L
end

function getLtri(L)
    p, q = size(L)
    x = zeros(eltype(L), (round(Int, p * (p + 1) / 2),))
    k = 1
    for j in 1:p
        for i in j:p
            x[k] = L[i, j]
            k += 1
        end
    end
    x
end

# each of these functions computes
# x' ⋅ Σ^-1 x = ||L'^-1 x||^2
p = 5
L2 = makeLtri(rand(round(Int, p * (p + 1) / 2)))
# Σ2 = PDMats.PDMat(L2 * L2')
μ3 = rand(p)
function f1(x)
    L = makeLtri(x[p+1:end])
    d = MvNormal(x[1:p], L * L')
    sqmahal(d, zeros(p))
end
function f2(x)
    Σ2 = L2 * L2'
    d = MvNormal(x, Σ2)
    sqmahal(d, zeros(p))
end
function f3(x)
    L = makeLtri(x)
    d = MvNormal(μ3, L * L')
    sqmahal(d, zeros(p))
end

function gradμ(μ, l)
    L = makeLtri(l)
    Σ = L * L'
    inv(Σ) * 2μ
end
function gradΣ(μ, l)
    L = makeLtri(l)
    Σ = L * L'
    X = μ * μ'
    getLtri(-2 * inv(Σ) * X * inv(L'))
end

function grad1(x)
    μ = x[1:p]
    l = x[p+1:end]
    vcat(gradμ(μ, l), gradΣ(μ, l))
end
grad2(x) = gradμ(x, getLtri(L2))
grad3(x) = gradΣ(μ3, x)

agrad1 = ForwardDiff.gradient(f1)
agrad2 = ForwardDiff.gradient(f2)
agrad3 = ForwardDiff.gradient(f3)

Ntests = 10
for _ in 1:Ntests
    m = rand(p)
    l = rand(round(Int, p * (p + 1) / 2))
    @assert grad1(vcat(m, l)) ≈ agrad1(vcat(m, l))
    @assert grad2(m) ≈ agrad2(m)
    @assert grad3(l) ≈ agrad3(l)
end

normal.jl

function f1(x)
    d = Normal(x[1], x[2])
    mean(d) / std(d)
end
function f2(x)
    d = Normal(x[1], 3.5)
    mean(d) / std(d)
end
function f3(x)
    d = Normal(1.1, x[1])
    mean(d) / std(d)
end

grad1(x) = [1/x[2], -x[1]/x[2]^2]
grad2(x) = [1/3.5]
grad3(x) = [-1.1/x^2]

agrad1 = ForwardDiff.gradient(f1)
agrad2 = ForwardDiff.gradient(f2)
agrad3 = ForwardDiff.gradient(f3)

pts = Vector{Float64}[
    [1., 2.],
    [1.15, 0.75],
    [0.5, 3.75]
]
for p in pts
    @assert grad1(p) ≈ agrad1(p)
    @assert grad2(p[1]) ≈ agrad2([p[1]])
    @assert grad3(p[2]) ≈ agrad3([p[2]])
end

runtests.jl

using Distributions
using ForwardDiff
using Base.Test

tests = [
    "normal",
    "mvnormal"
    ]

print_with_color(:blue, "Running tests:\n")

srand(12345)

for t in tests
    test_fn = "$t.jl"
    print_with_color(:green, "* $test_fn\n")
    include(test_fn)
end