%...DESCRIPTION: % %... Matlab program to run a Monte Carlo-Markov to invert the prameter values of a % %... simple function. The user can choose from a number of simple functions below % %... to test performance. Much of the code here has been 'inspired' by (read: % %... shamelessly derived from) Derek Posselt's MCMC code. % %... Please do browse through the code before using. I've tried to lay out the % %... code as clearly and logically as possible. At the very least, read to where % %... the MCMC algorithm begins. And of course, read all these initial remarks. % %... DELAYED REJECTION: option max_delay_lev activates delayed rejection in the % %... Metropolis-Hastings sampler for values greater than 1. Currenlty, the maximum % %... value for max_delay_lev is 2, corresponding to a single delayed rejection. % %... The delayed rejection sampling algorithm is described in Green & Mira 2001. % %... Briefly, when a proposal is rejected (whether on the basis of prior or like- % %... lihood probability), a second proposal takes place with reduced variance. This% %... second proposal is accepted with a modified likelihood. The second propasal % %... at least twice as computationally expensive as the first, and so the increase % %... in sampler efficiency may be offset by an increase in computer time. % %... ADAPTIVE METROPOLIS: option adaptive=true replaces traditional Metropolis % %... sampling with adaptive Metropolis sampling (Haario et al 2001). This non- % %... Markovian sampler continuously adapts the proposal covariance based on the % %... full chain (after a 'traditional' burn-in period). Though all chain samples % %... are used, adaption only occurs ever n_adapt cycles. % function dram_schematic(model,nmc,rand_seed) %...September 23, 2011 --- DR coded, rangecheck not working properly, all else ok %...September 24, 2011 --- AM coded, not tested. %...October 10, 2011 --- Proposal scaling parameter added (Gelman et al 1996) %-----------------------------------------------------------------------------------% % ~ ~ ~ MCMC inversion setup (general, then for each fcn) ~ ~ ~ % %-----------------------------------------------------------------------------------% global rstrm global nprtb max_burn = 5000; %...Max number of cycles for burn-in. This or nmc/2 is chosen set_guess = false; %...set initial parameter values to guess (1) or random (0) cov_check = 100; %...num of MCMC iterations between proposal variance updates cov_mem = 100; %...num of MCMC iterations to store for computing proposal var cov_prop = true; %...whether proposal search co-varies (true) or is independent diagnostic = 2; %...Set level of diagnostic output (0=none,1=some,2=tons) diag_to_file= true; %...Save diagnostic output to file (true) or std output (false) target_rate = 0.25; %...Target acceptance rate acc_thresh = 0.10; %...plus/minus for acceptance rate plot_figs = true; %...Whether to plot PDF disc_burn = true; %...Whether to discard burn-in samples before plotting uniform_dist= false; %...Whether to sample from a uniform parameter dist (true) max_delay_lev= 2; %...Maximum (for now) is two adaptive = true; %...Whether to do MH (false) or Adaptive Metropolis (true) n_adapt = 100; %...How often proposal is adapted gelman_scale= false; %...Whether to use the Gelman (1996) scaling parameter %...Create random number stream rstrm = RandStream('mt19937ar','seed',rand_seed); iburn_spin = cov_check*10; %...set ranges, truth and inital proposal, generate intial chol. decomp. matrix %--------------------------- Damped Harmonic Oscillator ----------------------------% %...As the title suggests - we will allow for the inversion of initial amplitude as % %...well as model coefficients (damping coefficient, spring constant). The ODE that %...is calculated is y" + by' + c = 0, with y(t=0)=a, y'(t=0) = 0. This is solved %...using the Matlab ODE solver ode45. The function is defined in mcmc_dampharm.m if strcmp(model,'dampharm') %tspan = 0:1:50; nparams = 3; % amp dmp spr % a b c p_true = ([ 50. 0.2 3. ]); p_min = ([ 30. 0.05 2.8 ]); p_max = ([100. 0.45 3.3 ]); p_guess = ([ 45. 0.25 2. ]); p_flag = ([ 1 1 1 ]); p_step = (p_max-p_min)*0.01; paramshort = {'amp' 'dmp' 'spr'}; if (cov_prop) p_sigma = eye(nparams); %...THIS IS COVARIANCE! p_cond = eye(nparams); for i=1:nparams p_sigma(i,i) = ((p_max(i)-p_min(i))*0.1)^2; p_cond(i,i) = p_sigma(i,i)*0.0001; end else p_sigma = (p_max-p_min)*0.1; %...THIS IS STDDEV! end %------------- Specify number and location of observations --------------- %...Note that we assume the grid defined by tspan in mcmc_fwd.m - this is %...duplicated above for the sake of convenience. nobs = 3; % obs_loc_arr(1) = 5; obs_loc_arr(2) = 15; obs_loc_arr(3) = 30; %-------------- Specify observational error (covariance) ----------------- var_factor = 0.100; %...Independent obs, scaled by amplitude obs_cov = eye(nobs)*sqrt(var_factor*(p_true(1))); %...Covariance from file %load dampharm_cov_1.mat %obs_cov = cov_arr; %...Invert the covariance matrix inv_cov = inv(obs_cov); %---------------------------- Multivariate Gaussian --------------------------------% elseif strcmp(model,'gaussian') % nparams = 3; %5; % V1 V2 V3 V4 V5 p_true = ([ 10. 10. 10. ]);% 10. 10. ]); %...True parameter values p_min = ([ 8. 8. 8. ]);% 8. 8. ]); %...minimum parameter val p_max = ([ 12. 12. 12. ]);% 12. 12. ]); %...maximum parameter val p_guess = ([ 9. 9. 9. ]);% 9. 9. ]); %...parameter first guess p_flag = ([ 1 1 1 ]);% 1 1 ]); %...are params inverted? if (cov_prop) p_sigma = eye(nparams); %...THIS IS COVARIANCE! p_cond = eye(nparams); for i=1:nparams p_sigma(i,i) = sqrt((p_max(i)-p_min(i))*0.1); p_cond(i,i) = p_sigma(i,i)*0.0001; end else p_sigma = (p_max-p_min)*0.1; %...proposal sigma end p_step = (p_max-p_min)*0.01; %...Histogram bin steps for plotting routine paramshort = {'var1' 'var2' 'var3'};% 'var4' 'var5'}; %------------- Specify number and location of observations --------------- nobs = 3; % 5; % V1 V2 V3 V4 V5 obs_loc_arr(1,:) = ([1.0 0.5 0.0 ]); % 0.0 0.0]); %...Obs 1 obs_loc_arr(2,:) = ([0.5 1.0 0.0 ]); % 0.0 0.0]); %...Obs 2 obs_loc_arr(3,:) = ([0.0 0.0 1.0 ]); % 0.0 0.0]); %...Obs 3 %obs_loc_arr(4,:) = ([0.0 0.0 0.0 ]); % 1.0 0.0]); %...Obs 4 %obs_loc_arr(5,:) = ([0.0 0.0 0.0 ]); % 0.0 1.0]); %...Obs 5 %--------------- Specify observation error (covariance) ----------------- var_factor = 0.200; %...Independent observational error obs_cov = eye(nobs)*var_factor; %...Specified error covariance matrix % Ob1 Ob2 Ob3 Ob4 Ob5 %obs_cov = ([ 1.0 0.38 0.0 0.0 0.0;... %...Obs 1 % 0.38 1.0 0.0 0.0 0.0;... %...Obs 2 % 0.0 0.0 1.0 0.0 0.0;... %...Obs 3 % 0.0 0.0 0.0 1.0 0.0;... %...Obs 4 % 0.0 0.0 0.0 0.0 1.0])... %...Obs 5 % *var_factor; % Ob1 Ob2 Ob3 %obs_cov = ([ 0.7 0.3 0.0;... % 0.3 0.7 0.0;... % 0.0 0.0 1.0])... % *var_factor; %...Covariance from file %load('cov_gaussian_1.mat','corr_arr'); %obs_cov = corr_arr*var_factor; %obs_cov = cov_arr*var_factor; %...Invert the covariance matrix inv_cov = inv(obs_cov); %---------------------------- Lorenz 3-Component (1963) ----------------------------% elseif (strcmp(model,'lorenz')) % nparams = 6; % X0 Y0 Z0 sig rho beta p_true = ([ 1. 1. 1. 10. (8./3.) 28. ]); p_min = ([ -2. -5.5 -15.5 2. 1.0 26. ]); p_max = ([ 18. 10.5 10.5 18. 5.5 30. ]); p_guess = ([ 2. 9. 10.8 9. 3. 20. ]); p_flag = ([ 1 0 1 1 1 0 ]); p_step = (p_max-p_min)*0.01; paramshort = {'X_0' 'Y_0' 'Z_0' '\sigma' '\rho' '\beta'}; %...Set proposal sigma/covariance if (cov_prop) p_sigma = eye(nparams); %...THIS IS COVARIANCE! p_cond = eye(nparams); for i=1:nparams p_sigma(i,i) = ((p_max(i)-p_min(i))*0.10)^2; p_cond(i,i) = p_sigma(i,i)*0.0001; end else p_sigma = (p_max-p_min)*0.1; %...THIS IS STDDEV! end %------------- Specify number and location of observations --------------- nobs = 4; % time index dimension obs_loc_arr(1,:) = ([ 21 1 ]); obs_loc_arr(2,:) = ([ 21 3 ]); obs_loc_arr(3,:) = ([ 31 1 ]); obs_loc_arr(4,:) = ([ 31 3 ]); %-------------- Specify observational error (covariance) ----------------- %...For lack of better assumptions, use independent obs. Otherwise, things %...could get messsy.... var_factor = 0.504; obs_cov = eye(nobs)*var_factor; %...Scale (somewhat arbitrarily) by the range of intital value parameters obs_cov = obs_cov*(2.); %...Invert the covariance matrix inv_cov = inv(obs_cov); end nprtb = sum(p_flag); %-------------------------- Open Files for Saving Data ----------------------------- save_name = ['mcmc_' model '_' int2str(var_factor*1000) '_' int2str(rand_seed)]; fid_par = fopen([save_name '_params.dat'],'w'); fid_fwd = fopen([save_name '_fwd_obs.dat'],'w'); fid_phi_tot = fopen([save_name '_phi_tot.dat'],'w'); fid_obs = fopen([save_name '_obs.dat'],'w'); %...Set up diagnostic output file if (diag_to_file) diid = fopen([save_name '_diagnostics.txt'],'w'); else diid = 1; end %...Save model settings save([save_name '_setting.mat'],'model','nparams','nobs','paramshort','obs_cov',... 'var_factor','p_true','obs_loc_arr','nmc','cov_check','cov_mem','cov_prop',... 'diagnostic','diag_to_file','target_rate','acc_thresh','uniform_dist'); %------------------------ Generate Simulated Observations --------------------------% [obs_arr] = mcmc_fwd(nparams,nobs,p_true,obs_loc_arr,model); fwrite(fid_obs,obs_arr,'double'); %-----------------------------------------------------------------------------------% % ~ ~ ~ Run First Guess ~ ~ ~ % %-----------------------------------------------------------------------------------% if not(set_guess) params = p_true; %...these will be overwritten if they are perturbed for i=1:nparams if (p_flag(i) == 1) params(i) = p_min(i) + (p_max(i) - p_min(i))*rand(rstrm); end end else params = p_guess end %...Run forward model and calculate likelihood, save output [fwd_obs_arr] = mcmc_fwd(nparams,nobs,params,obs_loc_arr,model); fwd_obs_save = fwd_obs_arr; like = likelihood(nobs,fwd_obs_arr,obs_arr,inv_cov); fwrite(fid_fwd,obs_arr,'double'); fwrite(fid_par,params,'double'); fwrite(fid_phi_tot,like,'double'); p_like = 0.; %...Set counters, etc main_chain = true; delay_lev = 1; icount = 1; bcount = 1; cov_count = 1; n_accept = 0; acc_rate = 0.; delay_count = 0; delaying = false; burn_in = true; %..."Allocate" storage & test arrays p_store = zeros(nparams,cov_mem); b_test = zeros(1,cov_check); R = zeros(nprtb); %...Generate scale parameter (Gelman 1996) based on number of perturbed parameters if (gelman_scale) s_n = 5.76/(nprtb); else s_n = .9; end tic %...Prepare initial proposal covariance for use in drawing proposals nprtb = 0; for i = 1:nparams if (p_flag(i) == 1) nprtb = nprtb + 1; nprtb_tmp = 0; for j = 1:nparams if (p_flag(j) == 1) nprtb_tmp = nprtb_tmp + 1; p_sigma_cut(nprtb,nprtb_tmp) = s_n*p_sigma(i,j); end end end end R = chol(p_sigma_cut); %...Cholesky decomposition of cov %...Simulated Annealing pre-sampler.................................................. if ( sim_annealing ) %...Open files for saving simulated annealing chain and accept data fid_anchain = fopen([save_name '_anchain.dat'],'w'); fid_anaccept = fopen([save_name '_anaccept.dat'],'w'); %...Set proposal variance for simulated annealing p_sigma_anneal = 0.1*(p_max-p_min); %...Simulated annealing sampler parameters k = 1; n_accept = 0; kmax = 5000; i_invT = log(k+1)*.005; while ( k < kmax ) invT = log(k+1)*.005; %...Generate independent proposal perturbation if ( false ) [pert_vec] = perturb(true,p_flag,p_sigma_anneal,nparams,1,R); else [pert_vec] = perturb(false,p_flag,p_sigma_anneal,nparams,1,R); end p_params = params + pert_vec; %...Check if proposal is in range. If so, run fwd model, calculate likelihood [found_pert] = rangecheck(p_min,p_max,p_params,p_flag,nparams); if (not(found_pert)) p_like = 0. else [fwd_obs_arr] = mcmc_fwd(nparams,nobs,p_params,obs_loc_arr,model); p_like = likelihood(nobs,fwd_obs_arr,obs_arr,inv_cov); end %...Accept proposal depending on ratio of likelihoods and 'temperature' if (rand(rstrm) < ((p_like/like)^invT)) accepted = true; n_accept = n_accept + 1; params = p_params; like = p_like; fwrite(fid_anaccept,params,'double'); p_saacc_mean = (k-1)*p_saacc_mean/k + params/k; else accepted = false; end fwrite(fid_anchain,params,'double'); %...Calculate simulated annealing statistics p_sachain_mean = (k-1)*p_sachain_mean/k + params/k; end end %...Main loop........................................................................ while(main_chain) %...draw proposal perturbation found_pert = false; accepted = false; [pert_vec]=perturb(cov_prop,p_flag,p_sigma,nparams,delay_lev,R); p_params = params + pert_vec; %...test for in-range [found_pert] = rangecheck(p_min,p_max,p_params,p_flag,nparams); %...Out of range, reject if (not(found_pert) && (delay_lev >= max_delay_lev)) accepted = false; delaying = false; %...Out of range, delay rejection elseif (not(found_pert) && (delay_lev < max_delay_lev)) accepted = false; delaying = true; p_like_save(delay_lev) = 0.; pert_vec_save(:,delay_lev) = pert_vec; %...Perturb found! run fwd model, calculate likelihood, accept probabilistically else if (not(found_pert)), 'there is something wrong',end [fwd_obs_arr] = mcmc_fwd(nparams,nobs,p_params,obs_loc_arr,model); p_like = likelihood(nobs,fwd_obs_arr,obs_arr,inv_cov); %...Test for acc/rej. Make special consideration for delayed rejection if (delay_lev == 1) if (rand(rstrm) < (p_like/like)) accepted = true; delaying = false; %...Reject probabilistically, reject elseif (delay_lev >= max_delay_lev) accepted = false; delaying = false; %...Reject probabilistically, delay rejection else accepted = false; delaying = true; p_like_save(delay_lev) = p_like; pert_vec_save(:,delay_lev) = pert_vec; end elseif (delay_lev == 2) delaying = false; p_star = p_params - squeeze(pert_vec_save(:,delay_lev-1))'; inrange = rangecheck(p_min,p_max,p_star,p_flag,nparams); if (inrange) [fwd_star] = mcmc_fwd(nparams,nobs,p_star,obs_loc_arr,model); [like_star] = likelihood(nobs,fwd_star,obs_arr,inv_cov); alpha_star = min(1,like_star/p_like); if (rand(rstrm) < ((p_like/like)*(1-alpha_star)/... (1-min(1,p_like_save(delay_lev-1)/like)))) accepted = true; else accepted = false; end else if (rand(rstrm) < (p_like/like)/(1-min(1,p_like_save(delay_lev-1)/like))) accepted = true; else accepted = false; end end end %%...Accept probabilistically %if ((p_like/like) > rand(rstream)) % accepted = true; %%...Reject probabilistically, reject %elseif (delay_lev >= max_delay_lev) % accepted = false; %%...Reject probabilistically, delay rejection %else % accepted = false; % delaying = true; % p_like_save(delay_lev) = p_like; % pert_vec_save(:,delay_lev) = pert_vec; %end end %...Skip increment, writing to file, burn-in..................................... if (delaying) delay_lev = delay_lev + 1; delay_count = delay_count + 1; %...Increment counters, save output.............................................. else icount = icount + 1; if (accepted) n_accept = n_accept+1; params = p_params; like = p_like; fwd_obs_save = fwd_obs_arr; %...Save previous fwd obs value else end %...Write to file fwrite(fid_par,params,'double'); fwrite(fid_fwd,fwd_obs_save,'double'); fwrite(fid_phi_tot,like,'double'); delay_lev = 1; if (mod(icount,1000)==0) fprintf(2,'Iteration # %u \n',icount); toc tic end %...Adapt proposal variance, if at adaption interval (AM or burn-in for MH) %...In the case of burn-in: if (burn_in) %...Fill storage for burn-in calculation for i=1:nparams for j=1:cov_mem-1 p_store(i,j) = p_store(i,j+1); end p_store(i,cov_mem) = params(i); end b_test(1:cov_check-1) = b_test(2:cov_check); if (accepted) b_test(cov_check) = 1; else b_test(cov_check) = 0; end %...Compute new acceptance rate. Note that acc_rate for the first 99 model runs %...will be inaccurate (before the array is filled). acc_rate = sum(b_test)/cov_check; %...Check if we are at an update iteration if ((bcount >= cov_mem) && (cov_count >= cov_check)) %...Check if burn-in is complete. if so, set burn-in to false if ((icount > iburn_spin) && ( abs(acc_rate - target_rate) <= acc_thresh)) burn_in = false; n_accept = 0; if (adaptive) p_sigma_running=p_sigma; adapt_count = cov_mem; end %...Make sure burn-in isn't taking too long elseif ((bcount*2 >= nmc) || (bcount >= max_burn)) burn_in = false; n_accept = 0; if (adaptive) p_sigma_running=p_sigma; adapt_count = cov_mem; end %...Otherwise, update proposal variance else cov_count = 1; p_mean = (mean(p_store,2))'; %...store old psigma values p_sigma_old = p_sigma; if (cov_prop),p_sigma_cut_old = p_sigma_cut;end %...update proposal covariance and R-matrix for cov_prop if (cov_prop) p_sigma = cov(p_store') + p_cond; nprtb = 0; for i = 1:nparams if (p_flag(i) == 1) nprtb = nprtb + 1; nprtb_tmp = 0; for j = 1:nparams if (p_flag(j) == 1) nprtb_tmp = nprtb_tmp + 1; p_sigma_cut(nprtb,nprtb_tmp) = s_n*p_sigma(i,j); end end end end R = chol(p_sigma_cut); %...Cholesky decomposition of cov else p_sigma = std(p_store,0,2); end %...Diagnostic output if (diagnostic>0) if (cov_prop) for i=1:nparams fprintf(diid,'Param(%u) mean = %6.4g \n',i,p_mean(i)); end fprintf(diid,'Old parameter covariance = \n'); for i=1:nparams for j=1:nparams fprintf(diid,' %6.4g ',p_sigma_old(i,j)); end fprintf(diid,'\n'); end fprintf(diid,'New parameter covariance = \n'); for i=1:nparams for j=1:nparams fprintf(diid,' %6.4g ',p_sigma(i,j)); end fprintf(diid,'\n'); end fprintf(diid,'----------------------------------------------------\n'); else for i=1:nparams fprintf(diid,['Param(%u) mean = %6.4g, Old psigma(%u) = %6.4g,'... ' New psigma(%u) = %6.4g\n'],i,p_mean(i),i,p_sigma_old(i)... ,i,p_sigma(i)); end end end end %...is burn-in complete? else cov_count = cov_count + 1; end %...are we at an update? bcount = bcount + 1; elseif (adaptive) %...Iterate adaptive counter adapt_count = adapt_count + 1; inv_adapt_count = 1./adapt_count; %...Update running mean p_mean_new = (adapt_count-1)*inv_adapt_count*p_mean + inv_adapt_count*params; %...Update running covariance p_sigma_running = (adapt_count-1)*inv_adapt_count*p_sigma_running + ... s_n*inv_adapt_count*(adapt_count*(p_mean'*p_mean) - ... (adapt_count+1)*p_mean_new'*p_mean_new + ... params'*params); %...Set running mean if (mod(adapt_count,n_adapt) == 0), fprintf(diid,'Old mean, New mean = %6.4g , %6.4g \n',p_mean,p_mean_new); end p_mean = p_mean_new; %...Set p_sigma = p_sigma_running if at adapt cycle if (mod(adapt_count,n_adapt) == 0), fprintf(diid,' covariance in adaption loop = \n'); for i=1:nparams for j=1:nparams fprintf(diid,' %6.4g ',p_sigma(i,j)); end fprintf(diid,'\n'); end fprintf(diid,'-----------------------------------------------------\n'); p_sigma = p_sigma_running + p_cond; nprtb = 0; for i = 1:nparams if (p_flag(i) == 1) nprtb = nprtb + 1; nprtb_tmp = 0; for j = 1:nparams if (p_flag(j) == 1) nprtb_tmp = nprtb_tmp + 1; p_sigma_cut(nprtb,nprtb_tmp) = p_sigma(i,j); end end end end R = chol(p_sigma_cut); %...Cholesky decomposition of cov end end %...TEST FOR FINISHED SAMPLING......................................... if( (icount - bcount +1) >= nmc ),main_chain = false;end end %...are we delaying? end %...while main_chain icount toc fclose(fid_par); fclose(fid_fwd); fclose(fid_phi_tot); fclose(fid_obs); pct_accept = n_accept/nmc %----------------------------- PLOT RELEVANT PLOTS ---------------------------------% if (plot_figs) load('~/Matlab_Scripts/colormap4.mat','mycmap'); h=figure('position',[0 0 100*nparams 100*nparams],'paperpositionmode','auto',... 'color','none','InvertHardcopy','off'); fid = fopen([save_name '_params.dat'],'r','l'); param_save = fread(fid,[nparams icount],'double','l'); if (disc_burn) arr_temp = param_save; clear param_save; param_save = arr_temp(:,max_burn:icount-1); icount_new = icount-max_burn; clear arr_temp; end for i=1:(nparams-1) %...Plot 1D marginals of posterior parameter PDF subplot(nparams,nparams,((i-1)*nparams+1)),... hist(param_save((i+1),:),(p_min((i+1)):p_step((i+1)):... p_max((i+1)))); hh = findobj(gca,'Type','patch'); set(hh,'FaceColor','black','EdgeColor','none'); axis tight; ymax=ylim; view(270,270); set(gca,'xdir','reverse','xtick',[],'ytick',[]); axis off hold on text(((p_max(i+1)+p_min(i+1))*0.5),ymax(2),... char(paramshort(i+1)),... 'fontsize',14,'fontweight','bold','horizontalalignment','right',... 'fontname','avantgarde','verticalalignment','middle'); axis off hold off xlim([p_min((i+1)) p_max((i+1))]); %------------------------- subplot(nparams,nparams,((nparams*(nparams-1)+i+1))); hist(param_save((i),:),(p_min(i):p_step(i):p_max(i))); hh = findobj(gca,'Type','patch'); set(hh,'FaceColor','black','EdgeColor','none'); axis tight ymax=ylim; set(gca,'ydir','reverse','xtick',[],'ytick',[]) axis off hold on text(((p_max(i)+p_min(i))*0.5),ymax(2),... char(paramshort(i)),... 'fontsize',14,'fontweight','bold','horizontalalignment','center',... 'verticalalignment','top','fontname','avantgarde'); axis off hold off xlim([p_min((i)) p_max((i))]) for j=1:(nparams-1) %...Plot 2D (joint) marginals of posterior parameter PDF if (i <= j) figure(h); %...Generate truth "crosshairs" x = p_min((j+1)):p_step((j+1)):p_max((j+1)); y = p_true((i))+0*x; yy = p_min((i)):p_step((i)):p_max((i)); xx = p_true((j+1))+0*yy; %...Process and plot the data dum = [param_save((j+1),:)' param_save((i),:)']; dumm = hist3(dum,{p_min((j+1)):p_step((j+1)):p_max((j+1)) ... p_min((i)):p_step((i)):p_max((i)) }); subplot(nparams,nparams,(1+i+nparams*(j-1))),... contourf((p_min((i)):p_step((i)):p_max((i))),... (p_min((j+1)):p_step((j+1)):p_max((j+1))),dumm,... 30,'edgecolor','none') set(gcf,'Colormap',mycmap) set(gca,'fontsize',5) %...plot crosshairs hold on plot(y,x,'white','linewidth',0.5) plot(yy,xx,'white','linewidth',0.5) plot(y,x,'red','linewidth',0.5) plot(yy,xx,'red','linewidth',0.5) hold off end end end print(h,'-depsc2','-painters',['Plot_' save_name '_parpar.eps']) %...Plot joint observation-parameter PDF......................................... %load('~/Matlab_Scripts/MC_colormap.mat','mycmap'); load('~/Matlab_Scripts/Newgray_inv.mat','mycmap'); hh=figure('position',[0 0 100*(nparams+1) 100*(nobs+1)],'paperpositionmode',... 'auto','color','none','InvertHardcopy','off'); fod = fopen([save_name '_fwd_obs.dat'],'r','l'); fwdobs_save = fread(fod,[nobs icount],'double','l'); if (disc_burn) arr_temp = fwdobs_save; clear fwdobs_save; fwdobs_save = arr_temp(:,max_burn:icount-1); clear arr_temp; end %...Figure out range of observation values in chain obs_min=zeros(1,nobs); obs_max=zeros(1,nobs); for j=1:nobs obs_min(j)=max(min(fwdobs_save(j,:)),0.); obs_max(j)=max(max(fwdobs_save(j,:)),1.e-10); end obs_step = (obs_max-obs_min)*0.01; %...Create 1D marginals for parameters for i=1:nparams subplot(nobs+1,nparams+1,nobs*(nparams+1)+i+1),... hist(param_save(i,:),(p_min(i):p_step(i):p_max(i))); axis tight; ymax=ylim; set(gca,'ydir','reverse','xtick',[],'ytick',[]) axis off hold on text(((p_max(i)+p_min(i))*0.5),ymax(2),char(paramshort(i)),'fontsize',14,... 'fontweight','bold','horizontalalignment','center','verticalalignment',... 'top','fontname','avantgarde'); axis off hold off xlim([p_min(i) p_max(i)]) end %...Create 1D marginals for observations for j=1:nobs subplot(nobs+1,nparams+1,1+(j-1)*(nparams+1)),... hist(fwdobs_save(j,:),(obs_min(j):obs_step(j):obs_max(j))); axis tight ymax=ylim; view(270,270) set(gca,'xdir','reverse','xtick',[],'ytick',[]) axis off hold on text(obs_min(j),ymax(2),['Obs' num2str(j)],'fontsize',16,'fontweight','bold',... 'horizontalalignment','left','verticalalignment','bottom','rotation',90,... 'fontname','avantgarde'); axis off hold off xlim([obs_min(j) obs_max(j)]) end %...Create 2D joint marginals for i=1:nparams for j=1:nobs dum = [fwdobs_save(j,:)' param_save(i,:)']; dumm = hist3(dum,{obs_min(j):obs_step(j):obs_max(j) ... p_min(i):p_step(i):p_max(i)}); subplot(nobs+1,nparams+1,1+i+(nparams+1)*(j-1)),... contourf((p_min(i):p_step(i):p_max(i)),(obs_min(j):obs_step(j):obs_max(j)),... dumm,30,'edgecolor','none') set(gcf,'Colormap',mycmap) set(gca,'fontsize',5) end end %...Save output print(hh,'-depsc2','-painters',['Plot_' save_name '_obspar.eps']) %...Plot observation PDFs......................................................... %load('~/Matlab_Scripts/MC_colormap.mat','mycmap'); load('~/Matlab_Scripts/Newgray_inv.mat','mycmap'); h=figure('position',[0 0 100*(nobs) 100*(nobs)],'paperpositionmode',... 'auto','color','none','InvertHardcopy','off'); %...Plot 1D marginals of posterior fwd observation PDF for i=1:(nobs-1) subplot(nobs,nobs,((i-1)*nobs+1));... hist(fwdobs_save(i+1,:),(obs_min(i+1):obs_step(i+1):obs_max(i+1))); hh = findobj(gca,'Type','patch'); set(hh,'Facecolor','black','edgecolor','none'); axis tight; ymax = ylim; view(270,270); set(gca,'xdir','reverse','xtick',[],'ytick',[]); axis off hold on text((obs_max(i+1)+obs_min(i+1))*0.5,ymax(2),['Obs' int2str(i+1)],'fontsize',... 14,'fontweight','bold','horizontalalignment','right','fontname',... 'avantgarde','verticalalignment','middle'); axis off hold off xlim([obs_min(i+1) obs_max(i+1)]); %------------------------------------------------- subplot(nobs,nobs,(nobs*(nobs-1)+i+1)); hist(fwdobs_save(i,:),(obs_min(i):obs_step(i):obs_max(i))); hh = findobj(gca,'Type','patch'); set(hh,'facecolor','black','edgecolor','none'); axis tight ymax=ylim; set(gca,'ydir','reverse','xtick',[],'ytick',[]) axis off hold on text((obs_max(i)+obs_min(i))*0.5,ymax(2),['Obs' int2str(i)],'fontsize',14,... 'fontweight','bold','horizontalalignment','center','verticalalignment',... 'top','fontname','avantgarde'); axis off hold off xlim([obs_min(i) obs_max(i)]); %...Plot 2D (joint) marginals of fwd obs PDF for j=1:(nobs-1) if (i <= j) figure(h) %...Process and plot the data dum = [fwdobs_save(j+1,:)' fwdobs_save(i,:)']; dumm = hist3(dum,{obs_min(j+1):obs_step(j+1):obs_max(j+1) ... obs_min(i):obs_step(i):obs_max(i) }); subplot(nobs,nobs,(1+i+nobs*(j-1))),... contourf((obs_min(i):obs_step(i):obs_max(i)),... (obs_min(j+1):obs_step(j+1):obs_max(j+1)),dumm,30,'edgecolor','none') set(gcf,'Colormap',mycmap) set(gca,'fontsize',5) end end end print(h,'-depsc2','-painters',['Plot_' save_name '_obsobs.eps']) end %...Subfunction PERTURB.............................................................. %function [pert_vec]=perturb(cov_prop,p_flag,p_sigma,nparams,delay_lev,varargin) function [pert_vec]=perturb(cov_prop,p_flag,p_sigma,nparams,delay_lev,R) global rstrm global nprtb pert_vec = zeros(1,nparams); if (cov_prop) %R = varargin{1}; p_params_pert = (R'*randn(rstrm,nprtb,1))'/sqrt(delay_lev); counter = 1; for i =1:nparams if (p_flag(i) == 1) pert_vec(i) = p_params_pert(counter); counter = counter + 1; else pert_vec(i) = 0; end end else for i=1:nparams if (p_flag(i) == 1) pert_vec(i) = p_sigma(i)*randn(rstrm)/sqrt(delay_lev); end end end %...Subfunction RANGECHECK........................................................... function [inrange]=rangecheck(p_min,p_max,param,p_flag,nparams) inrange = true; for i = 1:nparams if (p_flag(i) == 1) if ((param(i) > p_max(i)) || (param(i) < p_min(i))) inrange = false; end end end %...Subfunction LIKELIHOOD........................................................... function [like]=likelihood(nobs,fwd_obs_arr,obs_arr,inv_cov) like = 0; for i = 1:nobs for j =1:nobs %...Gaussian error like_res(i,j) = -0.5*abs((fwd_obs_arr(i)-obs_arr(i))*inv_cov(i,j)*... (fwd_obs_arr(j)-obs_arr(j))); end end like = exp(sum(sum(like_res))); %...Subfunction ACCEPT_REJ........................................................... %function [accepted]=accept_rej(like,p_like,delay_lev,varargin) %%function [accepted]=accept_rej(like,p_like,delay_lev,p_like_save,pert_vec_save) %if (delay_lev == 1) % if (rand(rstrm) < (p_like/like)) % accepted = true; % else % accepted = false; % end %elseif (delay_lev == 2) % p_like_save = varargin{1}; % pert_vec_save = varargin{2}; % p_star = p_params - p_vec_save(delay_lev-1); % inrange = rangecheck(p_min,p_max,p_star,nparams); % if (inrange) % [fwd_star] = mcmc_fwd(nparams,nobs,p_star,obs_loc_arr,model); % [like_star] = likelihood(nobs,fwd_star,obs_arr,inv_obs_cov); % alpha_star = min(1,like_star/p_like); % if (rand(rstrm) < ((p_like/like)*(1-alpha_star)/(1-min(1,p_like_save(delay_lev-1)/like)))) % accepted = true; % else % accepted = false; % end % else % if (rand(rstrm) < (p_like/like)/(1-min(1,p_like_save(delay_lev-1)/like))) % accepted = true; % else % accepted = false; % end % end % %end