From d158298e56d30c4523d9ad5ce350050888fdd25d Mon Sep 17 00:00:00 2001 From: ilopezgp Date: Wed, 5 Aug 2020 16:26:12 -0700 Subject: [PATCH] Substitute unnecessary matrix inversion, use linear solve instead. Remove test for covariance inverse. --- src/MCMC.jl | 9 ++------- test/MCMC/runtests.jl | 1 - 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/src/MCMC.jl b/src/MCMC.jl index d9166daf3..98f3b87b0 100644 --- a/src/MCMC.jl +++ b/src/MCMC.jl @@ -33,8 +33,6 @@ struct MCMCObj{FT<:AbstractFloat, IT<:Int} obs_sample::Vector{FT} "covariance of the observational noise" obs_noise_cov::Array{FT, 2} - "inverse of obs_noise_cov" - obs_noise_covinv::Array{FT, 2} "array of length N_parameters with the parameters' prior distributions" prior::Array{Prior, 1} "MCMC step size" @@ -96,14 +94,12 @@ function MCMCObj( param = param_init_copy log_posterior = [nothing] iter = [1] - obs_noise_covinv = inv(obs_noise_cov) accept = [0] if algtype != "rwm" error("only random walk metropolis 'rwm' is implemented so far") end MCMCObj{FT,IT}(obs_sample, obs_noise_cov, - obs_noise_covinv, priors, [step], burnin, @@ -167,12 +163,11 @@ function log_likelihood(mcmc::MCMCObj{FT}, log_rho = FT[0] if gvar == nothing diff = g - mcmc.obs_sample - log_rho[1] = -FT(0.5) * diff' * mcmc.obs_noise_covinv * diff + log_rho[1] = -FT(0.5) * diff' * (mcmc.obs_noise_cov \ diff) else - gcov_inv = inv(Diagonal(gvar)) log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar)) diff = g - mcmc.obs_sample - log_rho[1] = -FT(0.5) * diff' * gcov_inv * diff + log_gpfidelity + log_rho[1] = -FT(0.5) * diff' * (Diagonal(gvar) \ diff) + log_gpfidelity end return log_rho[1] end diff --git a/test/MCMC/runtests.jl b/test/MCMC/runtests.jl index 4c05e460a..a630548b9 100644 --- a/test/MCMC/runtests.jl +++ b/test/MCMC/runtests.jl @@ -69,7 +69,6 @@ using CalibrateEmulateSample.GPEmulator # 1.0/sqrt(0.05) * obs_sample ≈ 4.472 @test mcmc.obs_sample ≈ [4.472] atol=1e-2 @test mcmc.obs_noise_cov == σ2_y - @test mcmc.obs_noise_covinv == inv(σ2_y) @test mcmc.burnin == burnin @test mcmc.algtype == mcmc_alg @test mcmc.iter[1] == max_iter + 1