]> https://gitweb.dealii.org/ - code-gallery.git/commitdiff
Import the updated version from David's repository.
authorWolfgang Bangerth <bangerth@colostate.edu>
Wed, 23 Feb 2022 05:03:05 +0000 (22:03 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 23 Feb 2022 05:03:05 +0000 (22:03 -0700)
MCMC-Laplace/Matlab/forward_solver.m
MCMC-Laplace/Matlab/get_statistics.m
MCMC-Laplace/Matlab/main.m
MCMC-Laplace/Matlab/plot_solution.m
MCMC-Laplace/Matlab/precomputations.m

index 833532b10f001ee813b9115ecb225b15821b0ebb..33e0a1db38bcf0e1f24a2adf287f07d004988cf0 100644 (file)
@@ -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;
index ef3c48a78c5d4c34aa8b3ad7ce61b566b06a2360..c99353a616e3b84cbd63b723c20189cf5b1feb29 100644 (file)
@@ -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;
-
-
-
-
-
index f4c4786f741b4ac353148c57c2e756d8b4481afd..45f86a4380df510fe29c0b3205863f5b2bd54d0e 100644 (file)
@@ -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'])
index fb353f424b14e31fffe627a293d672e21f15d9a9..39ddb2bbc26328af333f95f02828f16173fffdf5 100644 (file)
@@ -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)
index 76cec1b684af4b4a07d9a2d73572e9c60fd9bad7..4d0426e558c8c8ecbde442dc4c4a27735490336a 100644 (file)
@@ -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);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.