From: Wolfgang Bangerth Date: Wed, 23 Feb 2022 05:03:05 +0000 (-0700) Subject: Import the updated version from David's repository. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d47f42e13781208902315a69a5d9557e07e1dc42;p=code-gallery.git Import the updated version from David's repository. --- diff --git a/MCMC-Laplace/Matlab/forward_solver.m b/MCMC-Laplace/Matlab/forward_solver.m index 833532b..33e0a1d 100644 --- a/MCMC-Laplace/Matlab/forward_solver.m +++ b/MCMC-Laplace/Matlab/forward_solver.m @@ -21,14 +21,14 @@ A = zeros(33^2,33^2); %loop over cells to build A for i=0:31 - for j=0:31 %build A by summing over contribution from each cell + for j=0:31 %build A by summing over contribution from each cell %find local coefficient in 8x8 grid - thet = theta(floor(i/4)+1,floor(j/4)+1); + theta_loc = theta(floor(i/4)+1,floor(j/4)+1); %update A by including contribution from cell (i,j) dof = [lbl(i,j),lbl(i,j+1),lbl(i+1,j+1),lbl(i+1,j)]; - A(dof,dof) = A(dof,dof) + thet*A_loc; + A(dof,dof) = A(dof,dof) + theta_loc*A_loc; end end @@ -43,6 +43,4 @@ A = sparse(A); U = A\b; %get new z values -z = reshape((M*U)',[13^2 1]); - - +z = M*U; diff --git a/MCMC-Laplace/Matlab/get_statistics.m b/MCMC-Laplace/Matlab/get_statistics.m index ef3c48a..c99353a 100644 --- a/MCMC-Laplace/Matlab/get_statistics.m +++ b/MCMC-Laplace/Matlab/get_statistics.m @@ -8,7 +8,7 @@ %OUTPUTS: %theta_mean = overall mean of chains %covars = covariance matrices of each independent chain -%autocov = mean of autocovariance matrix over all the chains +%autocovar = mean of autocovariance matrix over all the chains %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [theta_mean,covars,autocovar] = get_statistics(data,theta_means); @@ -22,7 +22,7 @@ covars = zeros(64,64,N); autocovar = zeros(64,64,2*L-1); %compute covariance matrices and mean autocovariance matrix -for n=1:N %loop over independent Markov chains +for n=1:N %loop over independent Markov chains %get data from chain n data_ = reshape(permute(data(:,:,:,n),[3 2 1]),[L 64]); @@ -39,8 +39,3 @@ end %compute mean of autocovariance matrix autocovar = autocovar(1:64,1:64,L:2*L-1)/N; - - - - - diff --git a/MCMC-Laplace/Matlab/main.m b/MCMC-Laplace/Matlab/main.m index f4c4786..45f86a4 100644 --- a/MCMC-Laplace/Matlab/main.m +++ b/MCMC-Laplace/Matlab/main.m @@ -17,7 +17,7 @@ load precomputations.mat %define lag time and data matrix data = zeros(8,8,L,N); %data matrix of samples at lag times -theta_means = zeros(8,8,N); %overall mean of theta +theta_means = zeros(8,8,N); %overall mean of theta tic @@ -63,12 +63,13 @@ parfor n=1:N %update theta means theta_means(:,:,n) = theta_mean/N_L; + end toc %compute statistics on data set -[theta_mean,covars,autocovars] = get_statistics(data,theta_means); +[theta_mean,covars,autocovar] = get_statistics(data,theta_means); %save data to Matlab workspace, labeled by N and N_L save (['data_N_' num2str(N) '_N_L_ ' num2str(N_L) '.mat']) diff --git a/MCMC-Laplace/Matlab/plot_solution.m b/MCMC-Laplace/Matlab/plot_solution.m index fb353f4..39ddb2b 100644 --- a/MCMC-Laplace/Matlab/plot_solution.m +++ b/MCMC-Laplace/Matlab/plot_solution.m @@ -23,14 +23,14 @@ Mp = reshape(Mp,[n^2 33^2]); %run forward solver on mean of theta A = zeros(33^2,33^2); for i=0:31 - for j=0:31 %build A by summing over contribution from each cell + for j=0:31 %build A by summing over contribution from each cell %find local coefficient in 8x8 grid - thet = theta(floor(i/4)+1,floor(j/4)+1); + theta_loc = theta(floor(i/4)+1,floor(j/4)+1); %update A by including contribution from cell (i,j) dof = [lbl(i,j),lbl(i,j+1),lbl(i+1,j+1),lbl(i+1,j)]; - A(dof,dof) = A(dof,dof) + thet*A_loc; + A(dof,dof) = A(dof,dof) + theta_loc*A_loc; end end @@ -49,4 +49,4 @@ close all %plot solution figure zs = reshape(Mp*U,[n n]); -surf(xsp,xsp,zs) \ No newline at end of file +surf(xsp,xsp,zs) diff --git a/MCMC-Laplace/Matlab/precomputations.m b/MCMC-Laplace/Matlab/precomputations.m index 76cec1b..4d0426e 100644 --- a/MCMC-Laplace/Matlab/precomputations.m +++ b/MCMC-Laplace/Matlab/precomputations.m @@ -18,7 +18,7 @@ lbl = @(i,j) 33*j+i+1; inv_lbl = @(k) [k-1-33*floor((k-1)/33),floor((k-1)/33)]; %construct measurement matrix, M -xs = 1/14:1/14:13/14; %measurement points +xs = 1/14:1/14:13/14; %measurement points M = zeros(13,13,33^2); for k=1:33^2 c = inv_lbl(k); @@ -30,7 +30,7 @@ for k=1:33^2 end M = reshape(M,[13^2 33^2]); -%construct exact coefficient matrix, thetas0 +%construct exact coefficient matrix, theta_hat theta_hat = ones(8,8); theta_hat(2:3,2:3) = 0.1; theta_hat(6:7,6:7) = 10; @@ -53,10 +53,10 @@ b(boundaries) = zeros(128,1); %enforce boundary conditions on b exact_values %set global parameters and functions for simulation -sig = 0.05; %likelihood standard deviation -sig_pr = 2; %prior (log) standard deviation -sig_prop = 0.0725; %proposal (log) standard deviation -theta0 = ones(8,8); %initial theta values +sig = 0.05; %likelihood standard deviation +sig_pr = 2; %prior (log) standard deviation +sig_prop = 0.0725; %proposal (log) standard deviation +theta0 = ones(8,8); %initial theta values forward_solver_ = @(theta) ... forward_solver(theta,lbl,A_loc,Id,boundaries,b,M); log_probability_ = @(theta,z) log_probability(theta,z,z_hat,sig,sig_pr);