implment fit_mean_Σ for Normal and MvNormal (#35)

Author: bgctw

docs/src/index.md

@@ -84,6 +84,17 @@ fit(::Type{D}, ::QuantilePoint, ::QuantilePoint) where {D<:Distribution}
 ```@docs
 fit(::Type{D}, ::Any, ::QuantilePoint, ::Val{stats} = Val(:mean)) where {D<:Distribution, stats}
 ```
+
+## Fit to mean and uncertainty parameter
+For bayesian inversion it is often required to specify a distribution given
+the expected value (the predction of the population value) and a description of 
+uncertainty of an observation.
+
+```@docs
+fit_mean_Σ
+```
+
+
 ## Currently supported distributions
 Univariate continuous
 - Normal

docs/src/mvlognormal.md

@@ -7,10 +7,6 @@ CurrentModule = DistributionFits
 Can be fitted to a given mean, provided the Covariance of the underlying
 normal distribution.
 
-```@docs
-fit_mean_Σ(::Type{MvLogNormal}, mean::AbstractVector{T1}, Σ::AbstractMatrix{T2}) where {T1 <:Real,T2 <:Real}
-```
-
 ```jldoctest; output = false, setup = :(using DistributionFits)
 Σ = hcat([0.6,0.02],[0.02,0.7])
 μ = [1.2,1.3]

src/fitstats.jl

@@ -295,3 +295,15 @@ end
 # function fit_mean_quantile(d::Type{D}, mean::Real, qp::QuantilePoint) where D<:Distribution  
 #     error("fit_mean_quantile not yet implemented for Distribution of type: $D")
 # end
+
+"""
+    fit_mean_Σ(::Type{<:Distribution}, mean, Σ)
+
+Fit a Distribution to mean and uncertainty quantificator Σ. 
+
+The meaning of `Σ` depends on the type of distribution:
+- `MvLogNormal`, `MvNormal`: the Covariancematrix of the associated normal distribution
+- `LogNormal`, `Normal`: the scale parameter, i.e. the standard deviation at log-scale, `σ`
+"""
+function fit_mean_Σ end
+ 

src/multivariate/mvlognormal.jl

@@ -1,12 +1,3 @@
-"""
-    fit_mean_Σ(::Type{<:Distribution}, mean, Σ)
-
-Fit a Distribution to mean and uncertainty quantificator Σ. 
-
-The meaning of `Σ` depends on the type of distribution:
-- `MvLogNormal`: the Covariancematrix of the associated normal distribution
-- `LogNormal`: the scale parameter, i.e. the standard deviation at log-scale, `σ`
-"""
 function fit_mean_Σ(::Type{MvLogNormal}, mean::AbstractVector{T1}, Σ::AbstractMatrix{T2}) where {T1 <:Real,T2 <:Real}
     _T = promote_type(T1, T2)
     fit_mean_Σ(MvLogNormal{_T}, mean, Σ)

src/multivariate/mvnormal.jl

@@ -0,0 +1,15 @@
+function fit_mean_Σ(::Type{MvNormal}, mean::AbstractVector{T1}, Σ::AbstractMatrix{T2}) where {T1 <:Real,T2 <:Real}
+    _T = promote_type(T1, T2)
+    fit_mean_Σ(MvNormal{_T}, mean, Σ)
+end
+function fit_mean_Σ(::Type{MvNormal{T}}, mean::AbstractVector{T1}, Σ::AbstractMatrix{T2}) where {T, T1 <:Real,T2 <:Real}
+    meanT = T1 == T ? mean : begin
+        meanT = similar(mean, T)
+        meanT .= mean
+    end
+    ΣT = T2 == T ? Σ : begin
+        ΣT = similar(Σ, T)
+        ΣT .= Σ
+    end
+    MvNormal(meanT, ΣT)
+end

src/multivariates.jl

@@ -5,7 +5,7 @@ for fname in [
     # "multinomial.jl",
     # "dirichletmultinomial.jl",
     # "jointorderstatistics.jl",
-    # "mvnormal.jl",
+    "mvnormal.jl",
     # "mvnormalcanon.jl",
     # "mvlogitnormal.jl",
     "mvlognormal.jl",

src/univariate/continuous/normal.jl

@@ -35,3 +35,12 @@ end
 function fit_mode_quantile(D::Type{Normal{T}}, mode::Real, qp::QuantilePoint) where {T}
     fit(D, QuantilePoint(mode, 0.5), qp)
 end
+
+function fit_mean_Σ(::Type{Normal}, mean::T1, σ::T2) where {T1 <: Real, T2 <: Real}
+    _T = promote_type(T1,T2)
+    fit_mean_Σ(Normal{_T}, mean, σ)
+end
+function fit_mean_Σ(D::Type{Normal{T}}, mean::Real, σ::Real) where {T}
+    Normal{T}(mean, σ)
+end
+

test/multivariate/mvnormal.jl

@@ -0,0 +1,11 @@
+using DistributionFits
+using Test
+
+@testset "fit_mean_Σ" begin
+    m = [3.0f0, 4.0f0]
+    Σ = hcat([2.0f0,0.2],[0.2, 2.0f0])
+    d = fit_mean_Σ(MvNormal, m, Σ)
+    @test mean(d) == m
+    @test cov(d) == Σ
+end;
+

test/multivariate/test_multivariate.jl

@@ -123,6 +123,7 @@ end
 
 const tests = [
     "mvlognormal",
+    "mvnormal",
 ]
 #tests = ["mvlognormal"]
 

test/univariate/continuous/normal.jl

@@ -8,3 +8,11 @@ end;
     d = Normal(3.0f0, 2.0f0)
     test_univariate_fits(d)
 end;
+
+@testset "fit_mean_Σ" begin
+    m = 3.0f0
+    σ = 2.0f0
+    d = fit_mean_Σ(Normal, m, σ)
+    @test mean(d) == m
+    @test scale(d) == σ
+end;