%...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW %----------------------------------------------------------------------------------85 % ~ ~ ~ MCMC_TOY.M ~ ~ ~ % %-----------------------------------------------------------------------------------% %... % %...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. % %... % %...Please be advised: % %... NB0: This program calls mcmc_fwd.m, which in turn may call other functions % %... please make sure all files are present before running. % %... NB1: parameter ranges may have to be tuned - in the absence of an intelligent % %... way to do this, try a first guess parameter range, then modify depending % %... on how much of the parameter space is taken up by non-negligible PDF % %... values % %... NB2: depending on the choice of model, observations and parameter ranges, the % %... sampler may not "find" the parameter regions of high likelihood. This % %... can result in poor sampler performance. This is especially true for % %... strongly nonlinear, multimodal or curved (degenerate) PDFs. % %... % %... Poly2 description: (NOT FUNCTIONAL) % %... 2nd order polynomial of the form: % %... f(x) = c*x^2 + b*x + a % %... % %... Poly3 description: (NOT FUNCTIONAL) % %... 3rd order polynomial of the form: % %... f(x) = d*x^3 + c*x^2 + b*x + a % %... % %... Lotka description: (NOT FUNCTIONAL) % %... Nonlinear predator-prey model as included in Matlab % %... yp = diag([1- .01*y(2), -1 + .02*y(1)])*y % %... % %... 1DADV description: % %... Nonlinear advection (diffusion?) of one-dimensional fluid, started from an % %... initial sinusoidal state. Parameters to be inverted are amplitude and freq- % %... uency of initial sinusoids as well as background amplitude. % %... % %... Dampharm description: % %... Damped harmonic osciallator with parameters of spring constant, damping % %... coefficient and initial amplitude to invert. % %... % %... Gaussian description: % %... A 'hack' to sample from a multivariate Gaussian distribution. Some linear % %... combination of proposal parameter values is returned as the observation. % %... With assumed Gaussian observation error, this results in a multivatiate % %... Gaussian posterior. Useful as a sanity check. % %... % %... % %...INPUTS: % %... model: Model whose parameter are inverted by MCMC % %... 'sine' ..... Sinusoid (semi-operational) % %... 'poly2' ..... 2nd order polynomial (not operational) % %... 'poly3' ..... 3rd order polynomial (not operational) % %... 'lotka' ..... Lotka-Volterra predator/prey model (not operational)% %... 'lorenz' ..... Lorenz 1963 3-variable model (not operational) % %... '1dadv' ..... 1D advection (Tomi's code ) % %... 'dampharm'..... damped harmonic osciallator % %... 'gaussian'..... point value with gaussian obs error = gaussian % %... % %... nmc: Number of MCMC cycles % %... rand_seed: Random number seeder, determines first and subsequent proposals. % %... % %...IMPORTANT!!! % %... -> when cov_prop=true, p_sigma is not a standard deviation, but is instead % %... a covariance matrix! This is reflected in the diagnostic output. % %... % %...OUTPUTS: % %... % %...UPDATES: % %... October 13 2010 ... work started % %... November 02 2010 ... main bits (sine) coded. % %... August 25 2011 ... covariant proposal added % %... August 31 2011 ... begin coding on advanced samplers % %... September 7 2011 ... begin coding lorenz 3-component model % %... % %-----------------------------------------------------------------------------------% %...Created by Marcus van Lier-Walqui, supported by NSF grant AGS-1019184. % %...Code based on work by Derek Posselt, supported by NASA MAP grants NNX09AJ43G % %...and NNX09AJ46G as well as Office of Naval Research grant N00173-10-1-G035. % %... % %...2011, RSMAS, University of Miami, Miami FL. % %... % %...COME ON YOU GUNNERS!!! % %... % %-----------------------------------------------------------------------------------% %...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW function mcmc_toy(model,nmc,rand_seed) % |\ /| \ / | \ /\ / % % | \/ | \/ |__ \/ \/ % %-----------------------------------------------------------------------------------% % ~ ~ ~ MCMC inversion setup (general, then for each fcn) ~ ~ ~ % %-----------------------------------------------------------------------------------% max_burn = 10000; %...Maximum number of cycles for burn-in. This or nmc/2 is chosen set_guess = 0; %...set initial parameter values to guess (1) or random (0) cov_check = 200; %...number of MCMC iterations between proposal variance updates cov_mem = 200; %...number of MCMC iterations to store for computing proposal var cov_prop = true; %...whether proposal search co-varies (true) or is independent %rand_seed = 9; %...Seed for random number generator (integer > 0 ) 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) %write = 1; %...Set level of file output (0=none,1=some,2=tons) 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 iburn_spin = cov_check*10; uniform_dist= false;%...Whether to sample from a uniform parameter dist (true) %...Create random number stream rstrm = RandStream('mt19937ar','seed',rand_seed); %---------------------------------- Sinusoid ---------------------------------------% if strcmp(model,'sine') %------------------ Initialize model parameters ---------------------- %...these are somewhat "fixed" - the average user needs not mess %...with these parameters. nparams = 3; %...Number of parameters gridsize = 20; %...Size of spatial domain xax = 1:gridsize; %... %... amplitude frequency phase % | | | p_true = ([ 5. 3. 0. ]); %..."Truth" parameter values p_min = ([ 2. 1. -1. ]); %...minimum allowed parameter values p_max = ([ 8. 5. 1. ]); %...maximum allowed parameter values p_guess= ([ 3. 4. 0.2 ]); %...first guess (if set_guess=1) p_flag = ([ 1 1 1 ]); %...are parameters inverted? (1=yes) p_sigma = (p_max-p_min) * 0.1; %...proposal sigma set to 0.1*range nobs = 2.; %...Note, this is the array index at which observations obs_loc_arr(1) = 1; %...occur, not the array value. Add or comment out obs_loc_arr(2) = 7; %...observations as neeeded. %obs_loc_arr(3) = %...Initialize observation sigma - this part isn't quite clear to me... needs to %...be tested! obs_sigma = 0.05 * (p_max(1)-p_min(1)); %...observation sigma scaled by amplitude %..."Allocate" save arrays for output of data after model run param_save = zeros(nparams,nmc); fwd_obs_save = zeros(nobs,nmc); % phi_save = zeros(nobs,nobs,nmc); phi_tot_save = zeros(1,nmc); obs_arr = zeros(1,nobs); fwd_obs_arr = zeros(1,nobs); %--------------------------------- 1D Advection ------------------------------------% elseif strcmp(model, '1dadv') % %------------------ Initialize model parameters ---------------------- %...Following the format in Tomi's matlab code (for the most part) %...NOTE: these parameters are specified in the forward code mcmc_fwd.m %ntime = 500; %...Number of time steps %dt = 10.; %...Time step %nspace = 100; %...Number of spatial grid intervals %L = 36000.; %...Domain size in some distance units (e.g. km) nparams = 5; %...Number of control parameters considered %...Description of parameters........................................ %... wl1 ... wavelength #1 as fraction of domain size %... wl2 ... wavelength #2 as fraction of domain size %... amp0 ... base amplitude (mean amplitude) %... amp1 ... wave #1 amplitude %... amp2 ... wave #2 amplitude %.................................................................... % wl1 wl2 amp0 amp1 amp2 % | | | | | p_true = ([ 0.5 0.125 10. 5. 2. ]); %...True parameter values p_min = ([ 0.4 0.08 7. 0.1 0.1 ]); %...minimum parameter range p_max = ([ 0.7 0.2 11. 9. 6.0 ]); %...maximum parameter range p_guess = ([ 0.7 0.2 6. 3. 0.9 ]); %... p_flag = ([ 1 1 1 1 1 ]); %...are params inverted? %...Determine inital proposal variance or 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.1)^2; p_cond(i,i) = p_sigma(i,i)*0.0001; end else p_sigma = (p_max-p_min)*0.05; %...proposal sigma end p_step = (p_max-p_min)*0.01; %...Histogram bin steps for plotting routine paramshort = {'\lambda_1' '\lambda_2' 'A_0' 'A_1' 'A_2'}; %...Specification of observation settings........................... %nobs = 3.; %...Number of obs %obs_loc_arr(1,:) = ([200 30]); %...([itime jspace]) of observation #1 %obs_loc_arr(2,:) = ([200 60]); %...([itime jspace]) of observation #2 %obs_loc_arr(3,:) = ([200 90]); %...([itime jspace]) of observation #3 %nobs = 5; %...Number of obs %obs_loc_arr(1,:) = ([100 60]); %...([itime jspace]) of observation #1 %obs_loc_arr(2,:) = ([200 60]); %...([itime jspace]) of observation #2 %obs_loc_arr(3,:) = ([300 60]); %...([itime jspace]) of observation #3 %obs_loc_arr(4,:) = ([400 60]); %...([itime jspace]) of observation #4 %obs_loc_arr(5,:) = ([500 60]); %...([itime jspace]) of observation #5 nobs = 8; %...Number of obs obs_loc_arr(1,:) = ([30 20]); %...([itime jspace]) of observation #1 obs_loc_arr(2,:) = ([30 40]); %...([itime jspace]) of observation #2 obs_loc_arr(3,:) = ([30 60]); %...([itime jspace]) of observation #3 obs_loc_arr(4,:) = ([30 80]); %...([itime jspace]) of observation #4 obs_loc_arr(5,:) = ([100 20]); %...([itime jspace]) of observation #5 obs_loc_arr(6,:) = ([100 40]); %...([itime jspace]) of observation #5 obs_loc_arr(7,:) = ([100 60]); %...([itime jspace]) of observation #5 obs_loc_arr(8,:) = ([100 80]); %...([itime jspace]) of observation #5 %-------------- Specify observational error (covariance) ----------------- %obs_cov = ([ 0.01 0.0 0.0 ;... %...Obs #1 % 0.0 0.01 0.0 ;... %...Obs #2 % 0.0 0.0 0.01 ]); %...Obs #3 var_factor = 0.100; obs_cov = eye(nobs)*var_factor; %...Independent obs %...Scale observational error by (avg of amplitude 1 & 2 ranges)^2 obs_cov = obs_cov*(((p_max(4)-p_min(4))+(p_max(5)-p_max(5)))*.5)^2; %...Invert the covariance matrix inv_cov = inv(obs_cov); %------------------------ 1D Advection - State Estimation --------------------------% %...Unlike the above experiment, here we attempt to invert the intial state using % %...MCMC. In order to make this problem tractable, we will allow for the reduction %... elseif strcmp(model,'1dadv_state') %--------------------------- 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 elseif 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 range p_max = ([ 12. 12. 12. ]);% 12. 12. ]); %...maximum parameter range p_guess = ([ 9. 9. 9. ]);% 9. 9. ]); %... p_flag = ([ 1 1 1 ]);% 1 1 ]); %...are params inverted? if (cov_prop) p_sigma = eye(nparams); %...THIS IS COVARIANCE! for i=1:nparams p_sigma(i,i) = sqrt((p_max(i)-p_min(i))*0.1); 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 = ([ 10. 10. 10. 10. (8./3.) 28. ]); p_min = ([ 3.5 1.5 7.5 8. 2.0 26. ]); p_max = ([ 18.5 19.5 18.5 12. 3.5 30. ]); p_guess = ([ 11. 9. 10.8 9. 3. 20. ]); p_flag = ([ 1 1 1 1 1 1 ]); 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 = 9; % time index dimension obs_loc_arr(1,:) = ([ 100 1 ]); obs_loc_arr(2,:) = ([ 100 2 ]); obs_loc_arr(3,:) = ([ 100 3 ]); obs_loc_arr(4,:) = ([ 200 1 ]); obs_loc_arr(5,:) = ([ 200 2 ]); obs_loc_arr(6,:) = ([ 200 3 ]); obs_loc_arr(7,:) = ([ 150 1 ]); obs_loc_arr(8,:) = ([ 150 2 ]); obs_loc_arr(9,:) = ([ 150 3 ]); %-------------- Specify observational error (covariance) ----------------- %...For lack of better assumptions, use independent obs. Otherwise, things %...could get messsy.... var_factor = 0.100; 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); %---------------------------- 2nd Order Polynomial ---------------------------------% elseif (strcmp(model,'poly2') .or. strcmp(model,'poly3')) % %------------------ Initialize model parameters ---------------------- %...these are somewhat "fixed" - the average user needs not mess %...with these parameters. %gridsize = 20; %...Size of domain %xax = 1:gridsize; %... %yax = 1:gridsize; %... %a_true = 3.; %... %b_true = 5.; %c_true = 2.; % %%---------------- Initialize user-defined parameters ----------------- %%...These parameters should be defined by the user to test and "tune" %%...the behavior of the inversion. % %nobs = 2; %...number of observation %obs_point = (5,15); % %%------------------ Generate "truth" measurement --------------------- %obs_arr= 0. end %-----------------------------------------------------------------------------------% % ~ ~ ~ END USER INPUT, BEGIN MCMC ALGORITH ~ ~ ~ % %-----------------------------------------------------------------------------------% %-------------------------- 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'); %-----------------------------------------------------------------------------------% % ~ ~ ~ Run First Guess ~ ~ ~ % %-----------------------------------------------------------------------------------% %...Set up diagnostic output file if (diag_to_file) diid = fopen([save_name '_diagnostics.txt'],'w'); else diid = 1; end %...Save model setting 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'); %------------------------ Generate 1st guess forward obs ---------------------------% if (set_guess == 1) params = p_guess; else params = p_true; %...these will be overwritten unless they are not 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 end [fwd_obs_arr] = mcmc_fwd(nparams,nobs,params,obs_loc_arr,model); %------------------------- Compute cost for 1st guess ------------------------------% %...calculate cost function - using matlab-specific shortcuts %...If uniform_dist = true, then the prior distribution (bounded unifrom) will be %...sampled, rather than the likelihood p(y|x). phi = zeros(nobs,nobs); phi_tot = 0.; if (uniform_dist) phi = -.05; phi_tot = -.05*nobs*nobs; else for i=1:nobs for j=1:nobs phi(i,j) = -0.5*abs((fwd_obs_arr(i)-obs_arr(i))*inv_cov(i,j)*... (fwd_obs_arr(j)-obs_arr(j))); %phi_tot = phi_tot+phi(i,j); end end phi_tot = sum(sum(phi)); end %phi = -0.5 * (obs_arr - fwd_obs_arr).^2 / (obs_sigma.^2); %phi_tot = sum(phi); %...set initial Maximum Likelihood Estimate phi_mle = phi; phi_tot_mle = phi_tot; params_mle = params; fwd_obs_arr_mle = fwd_obs_arr; if (diagnostic>0) fprintf(diid,'Printing Output for first guess parameters......................\n'); for i=1:nparams fprintf(diid,'parameter no.%u = %6.4g , ',i,params(i)); end fprintf(diid,'\n'); for i=1:nobs fprintf(diid,'forward obs no.%u = %6.4g, ',i,fwd_obs_arr(i)); end fprintf(diid,'\n'); %for i=1:nobs % fprintf(diid,'phi for obs no.%u = %6.4g, ',i,phi(i)) %end fprintf(diid,'\n'); fprintf(diid,'Total phi = %6.4g \n',phi_tot); fprintf(diid,'----------------------------------------------------------------\n'); end %-----------------------------------------------------------------------------------% % ~ ~ ~ MCMC Main loop (Burn-in, then sampling) ~ ~ ~ % %-----------------------------------------------------------------------------------% %...Initialize counters icount = 1; bcount = 1; cov_count = 1; n_accept = 0; burn_in = true; acc_rate = 0.; running = true; %..."Allocate" storage/test arrays p_store = zeros(nparams,cov_mem); b_test = zeros(1,cov_check); tic %...Calculate initial value of Cholesky decomp. covariance matrix if (cov_prop) %...Chop down parameter vectors and covariance matrix to prevent problems with %...unperturbed elements. nprtb = 0; for i = 1:nparams if (p_flag(i) == 1) nprtb = nprtb + 1; params_cut(nprtb) = params(i); 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 while ( running ) found_pert = true; accepted = false; %------------------------------ perturb parameters -------------------------------% p_params = p_true; %...initially, set parameter values to truth %...Perturb parameters to be inverted, %...Check to make sure they are within range %...For covarying parameter perturbation if (cov_prop) %size(R'*randn(nprtb,1)) p_params_cut = params_cut + (R'*randn(rstrm,nprtb,1))'; %...Generate perturbation vector %...Where appropriate, fill in perturbed parameters, check for out-of-range nprtb_tmp = 0; for i =1:nparams if (p_flag(i) == 1) nprtb_tmp = nprtb_tmp + 1; p_params(i) = p_params_cut(nprtb_tmp); if ((p_params(i) > p_max(i)) || (p_params(i) < p_min(i))) found_pert = false; if(diagnostic>=1) fprintf(diid,'Parameter out of range: p_param(%u) = %6.4g,%6.4g,%6.4g\n',... i,p_params(i),p_min(i),p_max(i)); end end end end %...For independent parameter perturbation else for i=1:nparams if (p_flag(i) == 1) p_params(i) = params(i) + p_sigma(i)*randn(rstrm); if ((p_params(i) > p_max(i)) || (p_params(i) < p_min(i))) found_pert = false; if(diagnostic>=1) fprintf(diid,'Parameter out of range: p_param(%u) = %6.4g,%6.4g,%6.4g\n',... i,p_params(i),p_min(i),p_max(i)); end end end end end if (found_pert == false) %...if any parameters are out of range, reject parameter set and cycle the loop if (diagnostic>1) fprintf(diid,... 'Rejecting parameter set (out of range), iteration %u...... \n',icount); for i=1:nparams fprintf(diid,'Proposed param(%u) = %6.4g, Previous param(%u) = %6.4g \n',... i,p_params(i),i,params(i)); end fprintf(diid,'------------------------------------------------------------\n'); end else %...otherwise, if params are within range, run the forward model %----------------------------- Run forward model ---------------------------------% %...print diagnostics if (diagnostic>1) fprintf(diid,'Proposed parameter set:'); for i=1:nparams fprintf(diid,'param(%u) = %6.4g, ',i,p_params(i)); end fprintf(diid,'\n'); end %...Save previous values of observations fwd_obs_arr_save = fwd_obs_arr; %...Run forward model [fwd_obs_arr] = mcmc_fwd(nparams,nobs,p_params,obs_loc_arr,model); %---------------------------- Compute Cost Function ------------------------------% %...Note that the code has been changed to consider covarying observational %...error. p_phi = zeros(nobs,nobs); p_phi_tot = 0.; if (uniform_dist) p_phi = -0.05; p_phi_tot = -0.05*nobs*nobs; else for i=1:nobs for j=1:nobs p_phi(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 p_phi_tot=sum(sum(p_phi)); end %p_phi = -0.5 * (obs_arr - fwd_obs_arr).^2 / (obs_sigma.^2); %p_phi_tot = sum(p_phi); %------------------------------- Test for Acceptance -----------------------------% %...If new liklihood is greater than old, then accept if (p_phi_tot > phi_tot) %...Set accepted=true and increment accepted counter accepted = true; if (burn_in == 0),n_accept = n_accept + 1;end %...Diagnostic output if (diagnostic>1) fprintf(diid,'Accepted parameter set: better fit!, iteration %u....... \n',icount); for i=1:nparams fprintf(diid,'Proposal param(%u) = %6.4g, Previous param(%u) = %6.4g \n',... i,p_params(i),i,params(i)); end for i=1:nobs fprintf(diid,'Forward obs(%u) = %6.4g, Obs(%u) = %6.4g \n',i,fwd_obs_arr(i),... i,obs_arr(i)); end %for i=1:nobs % fprintf(diid,'New phi for obs(%u) = %6.4g, Old phi for obs(%u) = %6.4g \n',... % i,p_phi(i),i,phi(i)) %end fprintf(diid,'New total phi = %6.4g, Old total phi = %6.4g \n',p_phi_tot,phi_tot); fprintf(diid,'----------------------------------------------------------------\n'); end %...Replace former likelihood and parameters with new phi = p_phi; phi_tot = p_phi_tot; params = p_params; if (cov_prop),params_cut = p_params_cut;end %...Check to see if we've found a new MLE if (phi_tot > phi_tot_mle) phi_mle = p_phi; phi_tot_mle = p_phi_tot; params_mle = p_params; fwd_obs_arr_mle = fwd_obs_arr; fprintf(diid,'----------------------------------------------------------------\n'); fprintf(diid,'Found new MLE!\n'); fprintf(diid,'MLE Likelihood = %6.4g \n',exp(phi_tot_mle)); for i=1:nparams fprintf(diid,'MLE param(%u) = %6.4g, ',i,p_params(i)); end fprintf(diid,'\n'); for i=1:nobs fprintf(diid,'MLE fwd obs(%u) = %6.4g, ',i,fwd_obs_arr(i)); end fprintf(diid,'\n'); fprintf(diid,'----------------------------------------------------------------\n'); end %...If not, accept probabilistically (is that a word?) elseif ( exp( p_phi_tot - phi_tot) >= rand(rstrm) ) %...Set accepted=true and increment accepted counter accepted = true; if (not(burn_in)),n_accept = n_accept + 1;end %...Diagnostic output if (diagnostic>1) fprintf(diid,'Accepted parameter set: probabilistic, iteration %u..... \n',icount); for i=1:nparams fprintf(diid,' Proposal param(%u) = %6.4g, Previous param(%u) = %6.4g \n',... i,p_params(i),i,params(i)); end for i=1:nobs fprintf(diid,'Forward obs(%u) = %6.4g, Obs(%u) = %6.4g \n',i,fwd_obs_arr(i),... i,obs_arr(i)); end %for i=1:nobs % fprintf(diid,' New phi for obs(%u) = %6.4g, Previous phi, obs(%u) = %6.4g \n',... % i,p_phi(i),i,phi(i)) %end fprintf(diid,' New total phi = %6.4g, Previous total phi = %6.4g \n',... p_phi_tot,phi_tot); fprintf(diid,'----------------------------------------------------------------\n'); end %...Replace former liklihood and parameters with new phi = p_phi; phi_tot = p_phi_tot; params = p_params; if (cov_prop),params_cut = p_params_cut;end else if (diagnostic>1) fprintf(diid,'Rejecting parameter set, iteration %u................... \n',icount); for i=1:nparams fprintf(diid,' Proposal param(%u) = %6.4g, Previous param(%u) = %6.4g \n',... i,p_params(i),i,params(i)); end for i=1:nobs fprintf(diid,'Forward obs(%u) = %6.4g, Obs(%u) = %6.4g \n',i,fwd_obs_arr(i),... i,obs_arr(i)); end %for i=1:nobs % fprintf(diid,' Proposal phi for obs(%u) = %6.4g, Old phi, obs(%u) = %6.4g \n',... % i,p_phi(i),i,phi(i)) %end fprintf(diid,' Proposal total phi = %6.4g, Old total phi = %6.4g \n',... p_phi_tot,phi_tot); fprintf(diid,'----------------------------------------------------------------\n'); end %...Set value of forward observations back to previous fwd_obs_arr = fwd_obs_arr_save; end %------------------------- End Acceptance Test ---------------------------------% %...If we are in burn-in, fill parameter storage array (?). This feels clunky and %...perhaps should be re-written when time permits. The advantage of this appears %...to be that the newest parameter set will always be at the end. The disadvan- %...tage is that it requires shifting all values for every model run. if (burn_in) 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 end %...Fill burn-in test array and compute new acceptance rate. Once again, this %...feels clunky and should perhaps be re-written. See notes for filling of the %...parameter storage array. Same conclusions apply. 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; %...Diagnostic output if (diagnostic>0) fprintf(diid,'Acceptance rate: %6.4g \n',acc_rate); end %---------------------------- Proposal Sigma Update ----------------------------% %...Check for whether we are at an update iteration. This is if (burn_in && (bcount >= cov_mem) && (cov_count >= cov_check)) %...Diagnostic output if (diagnostic>0) fprintf(diid,['Reached covariance update: bcount = %6.4g, cov_count = %6.4g,' ... 'cov_check = %6.4g \n'],bcount,cov_count,cov_check) ; end %...Check if burn-in is complete, if so, do not update proposal variance and %...set burn_in to false, print diagnostics. if ( (icount > iburn_spin) && ( abs(acc_rate-target_rate) <= acc_thresh) ) burn_in = false; if (cov_prop) fprintf(diid,'-----------------------------------------------------\n'); fprintf(diid,'Finished Burn-In!\n'); fprintf(diid,'Acceptance rate = %6.4g \n',acc_rate); fprintf(diid,'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 if (diagnostic>0) fprintf(diid,'-----------------------------------------------------\n'); fprintf(diid,'Finished Burn-In!\n'); fprintf(diid,'Acceptance rate = %6.4g \n',acc_rate); fprintf(diid,'Parameter sigma = %6.4g \n',p_sigma); fprintf(diid,'-----------------------------------------------------\n'); end end %...Also make sure burn-in isn't taking to long. If it takes longer than nmc/2 %...end it and move on to main sample loop elseif ( (bcount*2 >= nmc) || (bcount >= max_burn) ) burn_in = false; if (diagnostic>0) fprintf(diid,'-------------------------------------------------------\n'); fprintf(diid,'Number of burn iterations exceeded nmc/2 or max_burn!\n'); fprintf(diid,'bcount = %6.4g, max_burn = %6.4g, nmc = %6.4g \n',... bcount,max_burn,nmc); fprintf(diid,'Acceptance rate = %6.4g \n',acc_rate); if (cov_prop) fprintf(diid,'Parameter covariance = \n'); for i=1:nparams for j=1:nparams fprintf(diid,' %6.4g ',p_sigma(i,j)); end fprintf(diid,'\n'); end else fprintf(diid,'Parameter sigma = %6.4g \n',p_sigma); end fprintf(diid,'-------------------------------------------------------\n'); end %...Otherwise, update the proposal sigma values else %...reset counter cov_count = 1; %...Calculate mean p_mean = mean(p_store,2); %...Test for parameter values equal to zero - for certain experiments it may %...be useful to check for negative values as well. if (p_mean == 0.0) fprintf(diid,'Found parameter mean = 0.0!! STOPPING!!!\n'); return end %...Store old psigma values p_sigma_old = p_sigma; if (cov_prop) p_sigma_cut_old = p_sigma_cut; end %...Update proposal covariance or standard deviation...... if (cov_prop) p_sigma = cov(p_store') + p_cond; %...Chop down parameter vectors and covariance matrix to prevent problems with %...unperturbed elements. nprtb = 0; for i = 1:nparams if (p_flag(i) == 1) nprtb = nprtb + 1; params_cut(nprtb) = params(i); 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 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 %...burn-in check procedure and parameter sigma update %...If we are not at an update iteration, increment the counter else cov_count = cov_count +1; end %...if statement testing if we're in burn in and at a cov_check point end %...if statement testing if parameters are within range %--------------------------- SAVE DATA FOR OUTPUT ------------------------------% %...Increment our counters, why don't you. icount = icount + 1; if (burn_in) bcount = bcount+1; end if (mod(icount,1000)==0) fprintf(2,'Iteration # %u \n',icount); toc tic end fwrite(fid_par,params,'double'); fwrite(fid_fwd,fwd_obs_arr,'double'); fwrite(fid_phi_tot,phi_tot,'double'); %param_save(:,icount) = params; %fwd_obs_save(:,icount) = fwd_obs_arr; %phi_save(:,:,icount) = phi; %phi_tot_save(1,icount) = phi_tot; %-------------------------- TEST FOR FINISHED SAMPLE ---------------------------% if ( (icount - bcount + 1) >= nmc ) %...Diagnostic output if (diagnostic>0) fprintf(diid,'----------------------------------------------------------------\n'); fprintf(diid,'FINISHED SAMPLING!!! \n'); fprintf(diid,'Number Accepted = %u \n',n_accept); fprintf(diid,'Number of samples = %u \n',(icount-bcount)); fprintf(diid,'Total iterations = %u \n',icount); for i=1:nparams fprintf(diid,'MLE param(%u) = %6.4g, ',i,params_mle(i)); end fprintf(diid,'\n'); for i=1:nobs fprintf(diid,'Forward Obs(%u) = %6.4g, MLE fwd Obs(%u) = %6.4g, '... ,i,fwd_obs_arr(i),i,fwd_obs_arr_mle(i)); %fprintf(diid,'MLE cost fcn(%u) = %6.4g \n',i,phi_mle(i)) end fprintf(diid,'MLE cost function (total) = %6.4g \n',phi_tot_mle); fprintf(diid,'----------------------------------------------------------------\n'); end running = false; end %...If loop testing stage of MCMC algorithm end %...While loop - main MCMC cycle toc fclose(fid_par); fclose(fid_fwd); fclose(fid_phi_tot); fclose(fid_obs); pct_accept = n_accept/nmc; fprintf(diid,'TOTAL ACCEPTANCE PERCENTAGE = %6.4g \n',pct_accept); fprintf(diid,'Total iterations = %u \n',icount); %---------------------------------- SAVE DATA --------------------------------------% %save([save_name '.mat'], 'param_save', 'fwd_obs_save', 'phi_tot_save',... % 'params_mle', 'p_true'); %----------------------------- 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'); icount_save = icount-1; if (disc_burn) arr_temp = param_save; clear param_save; param_save = arr_temp(:,max_burn:icount_save); icount = 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'); 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_save],'double','l'); if (disc_burn) arr_temp = fwdobs_save; clear fwdobs_save; fwdobs_save = arr_temp(:,max_burn:icount_save); 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'); 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 %------------------------------------ EXIT -----------------------------------------%