diff --git a/src/distributions/normal_family/normal_family.jl b/src/distributions/normal_family/normal_family.jl
index 467bfa03..ee872bbc 100644
--- a/src/distributions/normal_family/normal_family.jl
+++ b/src/distributions/normal_family/normal_family.jl
@@ -667,9 +667,7 @@ getsufficientstatistics(::Type{MvNormalMeanCovariance}) = (identity, (x) -> x *
 getlogpartition(::NaturalParametersSpace, ::Type{MvNormalMeanCovariance}) = (η) -> begin
     (η₁, η₂) = unpack_parameters(MvNormalMeanCovariance, η)
     k = length(η₁)
-    C = fastcholesky(-η₂)
-    l = logdet(C)
-    Cinv = LinearAlgebra.inv!(C)
+    Cinv, l = cholinv_logdet(-η₂)
     return (dot(η₁, Cinv, η₁) / 2 - (k * log(2) + l)) / 2
 end
 
diff --git a/test/distributions/normal_family/normal_family_tests.jl b/test/distributions/normal_family/normal_family_tests.jl
index a92ff400..c46b1538 100644
--- a/test/distributions/normal_family/normal_family_tests.jl
+++ b/test/distributions/normal_family/normal_family_tests.jl
@@ -335,3 +335,34 @@ end
         @test sort(svd(fi_dist).S) ≈ sort(svd(approxFisherInformation).S) rtol = 1e-1
     end
 end
+
+@testitem "Diffrentiabilty of ExponentialFamily(ExponentialFamily.MvNormalMeanCovariance) logpdf" begin
+    include("./normal_family_setuptests.jl")
+    for i in 1:5, d in 2:3
+        rng = StableRNG(d * i)
+        μ = 10randn(rng, d)
+        L = LowerTriangular(randn(rng, d, d) + d * I)
+        Σ = L * L'
+        n_samples = 1
+        dist = MvNormalMeanCovariance(μ, Σ)
+
+        samples = rand(rng, dist, n_samples)
+
+        θ = pack_parameters(MvNormalMeanCovariance, (μ, Σ))
+        ef = convert(ExponentialFamilyDistribution, MvNormalMeanCovariance(μ, Σ))
+
+        nat_space2mean_space = (η) -> begin
+            dist = convert(Distribution, ExponentialFamilyDistribution(MvNormalMeanCovariance, η))
+            μ, Σ = mean(dist), cov(dist)
+            pack_parameters(MvNormalMeanCovariance, (μ, Σ))
+        end
+
+        for sample in eachcol(samples)
+            mean_gradient = ForwardDiff.gradient(Base.Fix2(gaussianlpdffortest, sample), θ)
+            nat_gradient = ForwardDiff.gradient((η) -> logpdf(ExponentialFamilyDistribution(MvNormalMeanCovariance, η), sample), getnaturalparameters(ef))
+            jacobian = ForwardDiff.jacobian(nat_space2mean_space, getnaturalparameters(ef))
+            #autograd failing to compute jacobian of matrix part correclty. Comparing only vector (mean) part.
+            @test nat_gradient[1:d] ≈ (jacobian'*mean_gradient)[1:d]
+        end
+    end
+end