%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
U = A\b;
%get new z values
-z = reshape((M*U)',[13^2 1]);
-
-
+z = M*U;
%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);
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]);
%compute mean of autocovariance matrix
autocovar = autocovar(1:64,1:64,L:2*L-1)/N;
-
-
-
-
-
%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
%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'])
%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
%plot solution
figure
zs = reshape(Mp*U,[n n]);
-surf(xsp,xsp,zs)
\ No newline at end of file
+surf(xsp,xsp,zs)
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);
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;
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);