WIP: Refactoring Multivariate Normal distributions

After refactoring, the class hierarchy would look like:
```julia
abstract AbstractMvNormal <: ContinuousMultivariateDistribution

# Here, ZeroVector is an empty type that indicates a vector is zero. 
# Some operations are defined to make the interactions with zero vectors easy:
#
#  + (x::Vector, z::ZeroVector) = x
#
# This type subsumes both zero-mean and non-zero-mean cases, while makes it possible
# to seamlessly leverage the zero-mean-case to simplify computation
#
immutable MvNormal{Cov<:AbstractPDMat,Mean<:Union(Vector{Float64},ZeroVector{Float64})} <: AbstractMvNormal
    μ::Mean
    Σ::Cov
end

immutable MvNormalCanon{Prec<:AbstractPDMat,HVec<:Union(Vector{Float64},ZeroVector{Float64})} <: AbstractMvNormal
    h::HVec
    J::Prec
end
```
Convenient constructors are provided, so that end users don't need to worry about the type hierarchy:
```julia
function MvNormal{Cov<:AbstractPDMat}(μ::Vector{Float64}, Σ::Cov)
    dim(Σ) == length(μ) || throw(DimensionMismatch("The dimensions of μ and Σ are inconsistent."))
    MvNormal{Cov,Vector{Float64}}(μ, Σ)
end

MvNormal{Cov<:AbstractPDMat}(Σ::Cov) = MvNormal{Cov,ZeroVector{Float64}}(ZeroVector(Float64, dim(Σ)), Σ)

MvNormal(μ::Vector{Float64}, Σ::Matrix{Float64}) = MvNormal(μ, PDMat(Σ))
MvNormal(μ::Vector{Float64}, σ::Vector{Float64}) = MvNormal(μ, PDiagMat(abs2(σ)))
MvNormal(μ::Vector{Float64}, σ::Real) = MvNormal(μ, ScalMat(length(μ), abs2(float64(σ))))

MvNormal(Σ::Matrix{Float64}) = MvNormal(PDMat(Σ))
MvNormal(σ::Vector{Float64}) = MvNormal(PDiagMat(abs2(σ)))
MvNormal(d::Int, σ::Real) = MvNormal(ScalMat(d, abs2(float64(σ))))
```
The constructor will automatically choose the most economical form.

Not finished yet, still ongoing work.