diff --git a/ddToolbox/WAIC.m b/ddToolbox/WAIC.m
new file mode 100644
index 00000000..85e3c8e3
--- /dev/null
+++ b/ddToolbox/WAIC.m
@@ -0,0 +1,138 @@
+classdef WAIC
+    %WAIC WAIC object
+    %   Extended description here
+    %
+	% References
+	% Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., 
+	% & Rubin, D. B. (2013). Bayesian Data Analysis, Third Edition. 
+	% CRC Press.
+    properties (SetAccess = protected)
+		lppd, pWAIC, WAIC_value, WAIC_standard_error
+		nSamples, nCases
+	end
+	
+	properties (Hidden = true)
+		log_lik
+		lppd_vec, pWAIC_vec, WAIC_vec
+		modelName
+	end
+    
+    methods
+        
+        function obj = WAIC(log_lik)
+			
+			[obj.nCases, obj.nSamples] = size(log_lik);
+			
+			obj.log_lik = log_lik;
+			clear log_lik
+			
+			% Calculate lppd
+			% Equation 7.5 from Gelman et al (2013) 
+			obj.lppd_vec = log( mean( exp(obj.log_lik) , 2) );
+			obj.lppd = sum(obj.lppd_vec);
+			
+			% Calculate effective number of samples, pWAIC
+			% Equation 7.12 from Gelman et al (2013) 
+			obj.pWAIC_vec = var(obj.log_lik,0,2);
+			obj.pWAIC = sum(obj.pWAIC_vec);
+			
+			% Calculate WAIC
+			obj.WAIC_value = -2 * obj.lppd + 2 * obj.pWAIC;
+			
+			% Calculate WAIC standard error
+			obj.WAIC_vec = -2 * obj.lppd_vec + 2 * obj.pWAIC_vec;
+			obj.WAIC_standard_error = sqrt(obj.nCases)*std(obj.WAIC_vec);
+			
+		end
+		
+		function comparisonTable = compare(obj)
+			% Compare WAIC info from mulitple models
+			assert(numel(obj)>1, 'expecting an array of >1 WAIC object')
+			
+			% Build a table of values
+			model = {obj.modelName}';
+			WAIC = [obj.WAIC_value]';
+			pWAIC = [obj.pWAIC]';
+			lppd = [obj.lppd]';
+			SE = [obj.WAIC_standard_error]';
+			dWAIC = WAIC - min(WAIC);
+			weight = exp(-0.5.*dWAIC) ./ sum(exp(-0.5.*dWAIC));
+			
+			% dSE is the SE of the difference in WAIC (not SE!) between 
+			% each model and the top ranked model
+			[~, i_best_model] = min([obj.WAIC_value]);
+			for m = 1:numel(obj)
+				if m == i_best_model
+					dSE(m,1) = NaN;
+				else
+					% Calculate SE of difference (of WAIC values) between
+					% model m and i_best_model
+					WAIC_diff = obj(i_best_model).WAIC_vec - obj(m).WAIC_vec;
+					dSE(m,1) = sqrt(obj(m).nCases)*std(WAIC_diff);
+				end
+			end
+			% create table
+			comparisonTable = table(model, WAIC, pWAIC, dWAIC, weight, SE, dSE, lppd);
+			% sort so best models (lowest WAIC) values are at top of table
+			comparisonTable = sortrows(comparisonTable,{'WAIC'},{'ascend'});
+		end
+		
+		function plot(obj)
+			% produce a WAIC comparison plot
+			
+			comparisonTable = obj.compare();
+			
+			% define y-value positions for each model
+			y = [1:1:size(comparisonTable,1)];
+
+			ms = 6;
+			clf
+			hold on
+			
+			% in-sample deviance as solid circles
+			in_sample_deviance = -2*comparisonTable.lppd;
+			isd = plot(in_sample_deviance, y, 'ko',...
+				'MarkerFaceColor','k',...
+				'MarkerSize', ms);
+			
+			% WAIC as empty cirlcles, with SE errorbars
+			%waic = plot(comparisonTable.WAIC, y, 'ko');
+			waic_eb = errorbar(comparisonTable.WAIC,y,comparisonTable.SE,...
+				'horizontal',...
+				'o',...
+				'LineStyle', 'none',...
+				'Color', 'k',...
+				'MarkerFaceColor','w',...
+				'MarkerSize', ms);
+			
+			% plot dSE models
+			waic_diff = errorbar(comparisonTable.dWAIC([2:end])+min(comparisonTable.WAIC),...
+				y([2:end])-0.2, comparisonTable.dSE([2:end]),...
+				'horizontal',...
+				'^',...
+				'LineStyle', 'none',...
+				'Color', [0.5 0.5 0.5],...
+				'MarkerFaceColor','w',...
+				'MarkerSize', ms);
+			
+			% formatting
+			xlabel('deviance');
+			set(gca,...
+				'YTick', y,...
+				'YTickLabel', comparisonTable.model,...
+				'YDir','reverse');
+			ylim([min(y)-1, max(y)+1]);
+			
+			vline(min(comparisonTable.WAIC), 'Color',[0.5 0.5 0.5]);
+			
+			legend([isd, waic_eb, waic_diff],...
+				{'in-sample deviance', 'WAIC (+/- SE)', 'SE of WAIC difference (+/- SE)'},...
+				'location', 'eastoutside');
+			
+			title('WAIC Model Comparison')
+
+		end
+        
+	end
+    
+end
\ No newline at end of file
diff --git a/ddToolbox/models/Model.m b/ddToolbox/models/Model.m
index f3152ea3..f9bd94e6 100644
--- a/ddToolbox/models/Model.m
+++ b/ddToolbox/models/Model.m
@@ -7,6 +7,10 @@
 		coda % handle to coda object
 		data % handle to Data class
 	end
+	
+	properties (Hidden = false, SetAccess = protected, GetAccess = public)
+		WAIC_stats
+	end
 
 	%% Private properties
 	properties (SetAccess = protected, GetAccess = protected)
@@ -255,7 +259,18 @@ function export(obj)
 			%% Calculate AUC
 			MAX_DELAY = 365;
 			obj.auc = obj.calcAreaUnderCurveForAll(MAX_DELAY);
-
+			
+			%% Calculate WAIC
+			% first, prepare a table of log likelihood values. Rows are
+			% observations, columns are MCMC samples.
+			samples = obj.coda.getSamples({'log_lik'});
+			% collapse over chains
+			[chains, samples_per_chain, N] = size(samples.log_lik);
+			log_lik = reshape(samples.log_lik, chains*samples_per_chain, N)';
+			%  second, create WAIC object
+			obj.WAIC_stats = WAIC(log_lik);
+			
+			obj.WAIC_stats.modelName = class(obj);
 		end
 
 		function auc = calcAreaUnderCurveForAll(obj, MAX_DELAY)
diff --git a/ddToolbox/models/nonparametric_models/NonParametric.m b/ddToolbox/models/nonparametric_models/NonParametric.m
index df75eff5..801f01c5 100644
--- a/ddToolbox/models/nonparametric_models/NonParametric.m
+++ b/ddToolbox/models/nonparametric_models/NonParametric.m
@@ -1,56 +1,56 @@
 classdef (Abstract) NonParametric < Model
 	%NonParametric  NonParametric is a subclass of Model for examining models that do NOT make parametric assumptions about the discount function.
-	
+
 	properties (Access = private)
 		%AUC_DATA
 	end
-	
+
 	methods (Access = public)
-		
+
 		function obj = NonParametric(data, varargin)
 			obj = obj@Model(data, varargin{:});
 			obj.dfClass = @DF_NonParametric;
 			% Create variables
 			obj.varList.participantLevel = {'Rstar'};
-			obj.varList.monitored = {'Rstar', 'alpha', 'epsilon', 'Rpostpred', 'P'};
-			
+			obj.varList.monitored = {'log_lik', 'Rstar', 'alpha', 'epsilon', 'Rpostpred', 'P'};
+
 			obj.varList.discountFunctionParams(1).name = 'Rstar';
 			obj.varList.discountFunctionParams(1).label = 'Rstar';
-			
+
 			obj.plotOptions.dataPlotType = '2D';
 		end
-		
-		
+
+
 		% % TODO: do this by injecting new AUC values into CODA?
 		% function [auc] = getAUC(obj)
 		% 	% return AUC measurements.
 		% 	% This will return an object array of stocastic objects
 		% 	names = obj.data.getIDnames('all');
-		% 	
+		%
 		% 	for ind = 1:numel(names)
 		% 		personInfo = obj.getExperimentData(ind);
 		% 		discountFunction = obj.dfClass('delays',personInfo.delays,...
 		% 			'theta', personInfo.dfSamples);
-		% 		
+		%
 		% 		% append to object array
 		% 		auc(ind) = discountFunction.AUC;
 		% 	end
 		% end
-		
+
 	end
-	
+
     methods (Hidden = true)
         function dispModelInfo(obj)
             display('Discount function: fits indiferrence points to each delay independently')
         end
     end
-	
+
 	methods (Static, Access = protected)
-		
+
 		function observedData = additional_model_specific_ObservedData(observedData)
 			observedData.uniqueDelays = sort(unique(observedData.DB))';
 			observedData.delayLookUp = calcDelayLookup();
-			
+
 			function delayLookUp = calcDelayLookup()
 				delayLookUp = observedData.DB;
 				for n=1: numel(observedData.uniqueDelays)
@@ -59,16 +59,16 @@ function dispModelInfo(obj)
 				end
 			end
 		end
-		
+
 	end
-	
+
 	methods (Access = private)
-		
+
 		% TODO: All this will be removed once refactoring of
 		function personStruct = getExperimentData(obj, p)
-			
-			
-			
+
+
+
 			% Create a structure with all the useful info about a person
 			% p = person number
 			participantName = obj.data.getIDnames(p);
@@ -81,29 +81,29 @@ function dispModelInfo(obj)
 			personStruct.delays = obj.observedData.uniqueDelays;
 			personStruct.dfSamples = obj.extractDiscountFunctionSamples(p);
 			personStruct.data = obj.data.getExperimentData(p);
-			
+
 			% TODO: THIS IS A BOTCH ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 			% the model currently assumes all participants have been
 			% tested on the same set of delays. But this is not
 			% necessarily true. So here we need to exclude inferences
 			% about discount fractions for delays that this person was
 			% never tested on.
-			
+
 			% want to remove:
 			% - columns of personStruct.dfSamples
 			% - personStruct.delays
 			% where this person was not tested on this delay
-			
+
 			temp = obj.data.getRawDataTableForParticipant(p);
 			delays_tested = unique(temp.DB);
-			
+
 			keep = ismember(personStruct.delays, delays_tested);
-			
+
 			personStruct.delays = personStruct.delays(keep);
 			personStruct.dfSamples = personStruct.dfSamples(:,keep);
 			% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 		end
-		
+
 		function dfSamples = extractDiscountFunctionSamples(obj, personNumber)
 			samples = obj.coda.getSamples({'Rstar'});
 			[chains, nSamples, participants, nDelays] = size(samples.Rstar);
@@ -113,21 +113,21 @@ function dispModelInfo(obj)
 				dfSamples(:,d) = vec(personSamples(:,:,d));
 			end
 		end
-		
+
 	end
-	
-	
-	
-	
-	
+
+
+
+
+
 	% ==========================================================================
 	% ==========================================================================
 	% PLOTTING
 	% ==========================================================================
 	% ==========================================================================
-	
+
 	methods (Access = public)
-		
+
 		function plot(obj, varargin) % overriding from Model base class
             %plot Plots EVERYTHING
             %   PLOT(model) or model.PLOT will call all plot functions.
@@ -136,16 +136,16 @@ function plot(obj, varargin) % overriding from Model base class
             %   [...] = model.PLOT(PARAM1,VAL1,PARAM2,VAL2,...) specifies one
             %   or more of the following name/value pairs:
             %
-            %      'shouldExportPlots' Either true or false. Default is true. 
-            
+            %      'shouldExportPlots' Either true or false. Default is true.
+
             close all
-			
+
 			% parse inputs
 			p = inputParser;
 			p.FunctionName = mfilename;
 			p.addParameter('shouldExportPlots', true, @islogical);
 			p.parse(varargin{:});
-			
+
 			obj.plotDiscountFunctionGrid();
 			% Export
 			if obj.plotOptions.shouldExportPlots
@@ -154,7 +154,7 @@ function plot(obj, varargin) % overriding from Model base class
 					'suffix', obj.modelFilename,...
 					'formats', obj.plotOptions.exportFormats);
 			end
-			
+
 			obj.plotDiscountFunctionsOverlaid();
 			% Export
 			if obj.plotOptions.shouldExportPlots
@@ -163,65 +163,65 @@ function plot(obj, varargin) % overriding from Model base class
 					'suffix', obj.modelFilename,...
 					'formats', obj.plotOptions.exportFormats);
 			end
-			
+
 			% EXPERIMENT PLOT ==================================================
 			obj.psychometric_plots();
 			obj.plotAllExperimentFigures();
-			
+
 			% Posterior prediction plot
             dfPlotFunc = @(n, fh) obj.plotDiscountFunction(n, 'axisHandle',fh);
             obj.postPred.plot(obj.plotOptions, obj.modelFilename, dfPlotFunc)
-			
-			
+
+
 			%% TODO...
 			% FOREST PLOT OF AUC VALUES ========================================
 			% TODO: Think about plotting this with GRAMM
 			% https://github.com/piermorel/gramm
 			%
 			%figUnivariateSummary(alldata)
-			
+
 		end
-		
-		
+
+
 		function plotExperimentOverviewFigure(obj, ind)
             %model.plotExperimentOverviewFigure(N) Creates a multi-panel figure
             %   model.plotExperimentOverviewFigure(N) creates a multi-panel figure
             %   corresponding to experiment N, where N is an integer.
-            
+
 			latex_fig(12, 14, 3)
 			h = layout([1 2 3]);
-			
+
 			% opts.pointEstimateType	= obj.plotOptions.pointEstimateType;
 			% opts.timeUnits			= obj.timeUnits;
-			
+
 % 			% create cell arrays of relevant variables
 % 			discountFunctionVariables = {obj.varList.discountFunctionParams.name};
 % 			responseErrorVariables    = {obj.varList.responseErrorParams.name};
-			
+
 			obj.plotPosteriorErrorParams(ind, 'axisHandle', h(1))
 			obj.plotPosteriorAUC(ind, 'axisHandle', h(2))
 			obj.plotDiscountFunction(ind, 'axisHandle', h(3))
-			
+
 		end
-		
+
 	end
-	
+
 	methods (Access = protected)
-		
+
 		function psychometric_plots(obj)
 			% TODO: plot data on these figures
-			
+
 			names = obj.data.getIDnames('all');
 			for ind = 1:numel(names) % loop over files
 				fh = figure('Name', ['participant: ' names{ind}]);
 				latex_fig(12,10, 8)
-				
+
 				personStruct = getExperimentData(obj, ind);
-				
+
 				% work out a good subplot arrangement
 				nSubplots = numel(personStruct.delays);
 				subplot_handles = create_subplots(nSubplots, 'square');
-				
+
 				% plot a set of psychometric functions, one for each delay tested
 				for d = 1:nSubplots
 					subplot(subplot_handles(d))
@@ -230,7 +230,7 @@ function psychometric_plots(obj)
 					samples = obj.coda.getSamplesAtIndex_asStochastic(ind,{'alpha','epsilon'});
 					samples.indifference = Stochastic('rstar');
 					samples.indifference.addSamples(personStruct.dfSamples(:,d));
-					
+
 					psycho = DF_SLICE_PsychometricFunction('samples', samples);
 					psycho.plot();
 					%% plot response data TODO: move this to Data ~~~~~~~~~
@@ -256,6 +256,6 @@ function psychometric_plots(obj)
 				close(fh)
 			end
 		end
-		
+
 	end
 end
diff --git a/ddToolbox/models/nonparametric_models/separateNonParametric.jags b/ddToolbox/models/nonparametric_models/separateNonParametric.jags
index f25f900c..08268a4c 100644
--- a/ddToolbox/models/nonparametric_models/separateNonParametric.jags
+++ b/ddToolbox/models/nonparametric_models/separateNonParametric.jags
@@ -14,7 +14,7 @@
 # Parameters
 # - alpha
 # - epsilon
-# - Rstar: 
+# - Rstar:
 
 
 
@@ -39,7 +39,7 @@ for (d in 2:length(uniqueDelays)) {
 
 for (p in participantIndexList){
     epsilon[p]  ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-    
+
     # note reparameterisation is not working in this case. As in, it causes an error.
     alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
 
@@ -53,21 +53,22 @@ for (p in participantIndexList){
 
     for (d in 1:length(uniqueDelays)) {
         Rstar[p,d] ~ dnorm(mu[p,d], prec[d]) T(0,)
-        #Rstar[p,d] ~ dt(mu[p,d], prec[d], 1) T(0,) # equals Cauchy when k=1  
+        #Rstar[p,d] ~ dt(mu[p,d], prec[d], 1) T(0,) # equals Cauchy when k=1
     }
-    
+
 }
 
 for (t in 1:length(ID)) {
 	rewardratio[t] <- A[t]/B[t]
-    
+
     # if we are dealing with negative rewards, we want alpha to flip to negative, in order to 'horizontally flip' the psychometric function
     signflip[t] <- ifelse( B[t]>0, 1, -1 )
-    
+
 	df[t] <- Rstar[ID[t], delayLookUp[t]]
 	P[t] <- epsilon[ID[t]] + (1-2*epsilon[ID[t]]) * (1 - phi( (rewardratio[t]-df[t]) / alpha[ID[t]]*signflip[t]  ))
 
 	R[t]         ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 	Rpostpred[t] ~ dbern(P[t]) # posterior predicted response
 }
 
diff --git a/ddToolbox/models/parametric_models/beta_delta_models/BetaDelta.m b/ddToolbox/models/parametric_models/beta_delta_models/BetaDelta.m
index 2fc9cefc..dbefe8b7 100644
--- a/ddToolbox/models/parametric_models/beta_delta_models/BetaDelta.m
+++ b/ddToolbox/models/parametric_models/beta_delta_models/BetaDelta.m
@@ -10,7 +10,7 @@
 
 			% Create variables
 			obj.varList.participantLevel = {'beta','delta','epsilon'};
-			obj.varList.monitored = {'beta','delta', 'alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'beta','delta', 'alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
     		obj.varList.discountFunctionParams(1).name = 'beta';
 			obj.varList.discountFunctionParams(1).label = '$\beta$';
 			obj.varList.discountFunctionParams(2).name = 'delta';
diff --git a/ddToolbox/models/parametric_models/beta_delta_models/hierarchicalBetaDelta.jags b/ddToolbox/models/parametric_models/beta_delta_models/hierarchicalBetaDelta.jags
index 076e2348..bf4970ee 100644
--- a/ddToolbox/models/parametric_models/beta_delta_models/hierarchicalBetaDelta.jags
+++ b/ddToolbox/models/parametric_models/beta_delta_models/hierarchicalBetaDelta.jags
@@ -58,6 +58,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/beta_delta_models/mixedBetaDelta.jags b/ddToolbox/models/parametric_models/beta_delta_models/mixedBetaDelta.jags
index a9f01ab0..7ae892f3 100644
--- a/ddToolbox/models/parametric_models/beta_delta_models/mixedBetaDelta.jags
+++ b/ddToolbox/models/parametric_models/beta_delta_models/mixedBetaDelta.jags
@@ -60,6 +60,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/beta_delta_models/separateBetaDelta.jags b/ddToolbox/models/parametric_models/beta_delta_models/separateBetaDelta.jags
index fc110974..1eb7aa5b 100644
--- a/ddToolbox/models/parametric_models/beta_delta_models/separateBetaDelta.jags
+++ b/ddToolbox/models/parametric_models/beta_delta_models/separateBetaDelta.jags
@@ -47,6 +47,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/EbertPrelec.m b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/EbertPrelec.m
index c26d25c6..8b87a651 100644
--- a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/EbertPrelec.m
+++ b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/EbertPrelec.m
@@ -4,28 +4,28 @@
 
 		function obj = EbertPrelec(data, varargin)
 			obj = obj@SubjectiveTimeModel(data, varargin{:});
-            
+
             obj.dfClass = @DF_EbertPrelec;
             obj.subjectiveTimeFunctionFH = @SubjectiveTimePowerFunctionEP;
 
 			% Create variables
 			obj.varList.participantLevel = {'k','tau','alpha','epsilon'};
-			obj.varList.monitored = {'k','tau','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'k','tau','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
             obj.varList.discountFunctionParams(1).name = 'k';
             obj.varList.discountFunctionParams(1).label = 'discount rate, $k$';
             obj.varList.discountFunctionParams(2).name = 'tau';
             obj.varList.discountFunctionParams(2).label = 'tau';
-            
+
             obj.plotOptions.dataPlotType = '2D';
 		end
 
     end
-    
+
     methods (Hidden = true)
         function dispModelInfo(obj)
             display('Discount function: V = reward * exp(-(k*delay)^tau)')
             display('Ebert, J. E. J., & Prelec, D. (2007). The Fragility of Time: Time-Insensitivity and Valuation of the Near and Far Future. Management Science, 53(9), 1423-1438.')
         end
     end
-    
+
 end
diff --git a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/hierarchicalEbertPrelec.jags b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/hierarchicalEbertPrelec.jags
index 4e483326..82bfab52 100644
--- a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/hierarchicalEbertPrelec.jags
+++ b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/hierarchicalEbertPrelec.jags
@@ -4,23 +4,23 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES
 
 groupKmu       ~ dnorm(0,1/(0.5^2))
 groupKsigma    ~ dexp(0.5)
 
-groupTAUmu          ~ dnorm(1,1/(1^2))       
-groupTAUsigma     ~ dexp(0.5) 
+groupTAUmu          ~ dnorm(1,1/(1^2))
+groupTAUsigma     ~ dexp(0.5)
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     # using reparameterisation to avoid funnel of hell
     k_offset[p]     ~ dnorm(0,1)
     k[p]            <- groupKmu + k_offset[p] * groupKsigma
-    
+
     # TODO: use reparameterisation to avoid funnel of hell
     #TAU_offset[p]  ~ dnorm(0,1)
     #tau[p]         <- groupTAUmu + TAU_offset[p] * groupTAUsigma
-    
+
     tau[p]          ~ dnorm(groupTAUmu, groupTAUsigma) T(0,)
 }
 
@@ -63,6 +63,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/mixedEbertPrelec.jags b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/mixedEbertPrelec.jags
index af376568..5c8cce8b 100644
--- a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/mixedEbertPrelec.jags
+++ b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/mixedEbertPrelec.jags
@@ -8,12 +8,12 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
-K_MEAN          <- 0       
+K_MEAN          <- 0
 K_PRECISION     <- 1/(0.5^2)
-TAU_MEAN        <- 1       
-TAU_PRECISION   <- 1/(1^2)   
+TAU_MEAN        <- 1
+TAU_PRECISION   <- 1/(1^2)
 
 for (p in 1:nRealExperimentFiles){       # no +1 because no shrinkage hyperprior
 	k[p]        ~ dnorm(K_MEAN, K_PRECISION)
@@ -45,7 +45,7 @@ epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     epsilon[p]  ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-    
+
     # using reparameterisation to avoid funnel of hell
     alpha_offset[p] ~ dnorm(0,1) T(0,)
     alpha[p] <- alpha_mu + alpha_offset[p] * alpha_sigma
@@ -62,6 +62,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/separateEbertPrelec.jags b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/separateEbertPrelec.jags
index 94c19b20..6d9dc5d9 100644
--- a/ddToolbox/models/parametric_models/ebert_and_prelec_2007/separateEbertPrelec.jags
+++ b/ddToolbox/models/parametric_models/ebert_and_prelec_2007/separateEbertPrelec.jags
@@ -8,12 +8,12 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
-K_MEAN          <- 0       
+K_MEAN          <- 0
 K_PRECISION     <- 1/(0.5^2)
-TAU_MEAN        <- 1       
-TAU_PRECISION   <- 1/(1^2)   
+TAU_MEAN        <- 1
+TAU_PRECISION   <- 1/(1^2)
 
 for (p in 1:nRealExperimentFiles){       # no +1 because no shrinkage hyperprior
 	k[p]        ~ dnorm(K_MEAN, K_PRECISION)
@@ -56,6 +56,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/exponential_log_models/ExponentialLog.m b/ddToolbox/models/parametric_models/exponential_log_models/ExponentialLog.m
index 112c37a0..bb05eedc 100644
--- a/ddToolbox/models/parametric_models/exponential_log_models/ExponentialLog.m
+++ b/ddToolbox/models/parametric_models/exponential_log_models/ExponentialLog.m
@@ -4,23 +4,23 @@
 
 		function obj = ExponentialLog(data, varargin)
 			obj = obj@SubjectiveTimeModel(data, varargin{:});
-            
+
             obj.dfClass = @DF_ExponentialLog;
             obj.subjectiveTimeFunctionFH = @SubjectiveTimeLogFunction;
-            
+
 			% Create variables
 			obj.varList.participantLevel = {'k','logtau','alpha','epsilon'};
-			obj.varList.monitored = {'k','logtau','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'k','logtau','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
             obj.varList.discountFunctionParams(1).name = 'k';
             obj.varList.discountFunctionParams(1).label = 'discount rate, $k$';
             obj.varList.discountFunctionParams(2).name = 'logtau';
             obj.varList.discountFunctionParams(2).label = '$\log(\tau)$';
-            
+
             obj.plotOptions.dataPlotType = '2D';
 		end
-    
+
     end
-    
+
     methods (Hidden = true)
         function dispModelInfo(obj)
             display('Discount function: V = reward * exp(-k*log(1+delay*tau))')
diff --git a/ddToolbox/models/parametric_models/exponential_log_models/mixedExpLog.jags b/ddToolbox/models/parametric_models/exponential_log_models/mixedExpLog.jags
index 18f51060..ff6ce6d5 100644
--- a/ddToolbox/models/parametric_models/exponential_log_models/mixedExpLog.jags
+++ b/ddToolbox/models/parametric_models/exponential_log_models/mixedExpLog.jags
@@ -8,7 +8,7 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 K_MEAN          <- 0                     # <---- currently guesstimating
 K_PRECISION     <- 1/(0.01)              # <---- currently guesstimating
@@ -45,7 +45,7 @@ epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     epsilon[p]  ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-    
+
     # using reparameterisation to avoid funnel of hell
     alpha_offset[p] ~ dnorm(0,1) T(0,)
     alpha[p] <- alpha_mu + alpha_offset[p] * alpha_sigma
@@ -62,6 +62,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/exponential_log_models/separateExpLog.jags b/ddToolbox/models/parametric_models/exponential_log_models/separateExpLog.jags
index 456c8280..23f74d26 100644
--- a/ddToolbox/models/parametric_models/exponential_log_models/separateExpLog.jags
+++ b/ddToolbox/models/parametric_models/exponential_log_models/separateExpLog.jags
@@ -8,7 +8,7 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 K_MEAN          <- 0                     # <---- currently guesstimating
 K_PRECISION     ~ dgamma(0.001, 0.001)     #sigma=0.1         # <---- currently guesstimating
@@ -54,6 +54,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/exponential_models/Exponential1.m b/ddToolbox/models/parametric_models/exponential_models/Exponential1.m
index 887d2577..29b3e3b0 100644
--- a/ddToolbox/models/parametric_models/exponential_models/Exponential1.m
+++ b/ddToolbox/models/parametric_models/exponential_models/Exponential1.m
@@ -4,18 +4,18 @@
 
 		function obj = Exponential1(data, varargin)
 			obj = obj@Parametric(data, varargin{:});
-            
+
             obj.dfClass = @DF_Exponential1;
 
 			% Create variables
 			obj.varList.participantLevel = {'k','alpha','epsilon'};
-			obj.varList.monitored = {'k','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'k','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
             obj.varList.discountFunctionParams(1).name = 'k';
             obj.varList.discountFunctionParams(1).label = 'discount rate, $k$';
-            
+
             obj.plotOptions.dataPlotType = '2D';
 		end
-		
+
 	end
 
     methods (Hidden = true)
@@ -23,5 +23,5 @@ function dispModelInfo(obj)
             display('Discount function: V = exp(-k.delay)')
         end
     end
-    
+
 end
diff --git a/ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.jags b/ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.jags
index 6ee27d23..5dbac523 100644
--- a/ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.jags
+++ b/ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.jags
@@ -8,7 +8,7 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES
 
 K_MEAN       ~ dnorm(0.01, 1/(0.5^2))
 K_PRECISION  ~ dgamma(0.001,0.001)
@@ -31,7 +31,7 @@ for (t in 1:length(ID)) {
 
 # comparison acuity (alpha) hyperprior
 groupALPHAmu        ~ dnorm(0,1/(100^2)) T(0,)
-groupALPHAsigma     ~ dunif(0,100) 
+groupALPHAsigma     ~ dunif(0,100)
 
 # error rates (epsilon) hyperprior
 groupW          ~ dbeta(1.1, 10.9)  # mode for lapse rate
@@ -43,7 +43,7 @@ epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     epsilon[p]  ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-    
+
     # using reparameterisation to avoid funnel of hell
     alpha_offset[p] ~ dnorm(0,1) T(0,)
     alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
@@ -60,6 +60,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t])
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/exponential_models/mixedExp1.jags b/ddToolbox/models/parametric_models/exponential_models/mixedExp1.jags
index a22d4fff..0096a827 100644
--- a/ddToolbox/models/parametric_models/exponential_models/mixedExp1.jags
+++ b/ddToolbox/models/parametric_models/exponential_models/mixedExp1.jags
@@ -8,7 +8,7 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 K_MEAN          <- 0.01                  # <---- currently guesstimating
 K_PRECISION     <- 1/(0.5^2)             # <---- currently guesstimating
@@ -28,7 +28,7 @@ for (t in 1:length(ID)) {
 
 # comparison acuity (alpha) hyperprior
 groupALPHAmu        ~ dnorm(0,1/(100^2)) T(0,)
-groupALPHAsigma     ~ dunif(0,100) 
+groupALPHAsigma     ~ dunif(0,100)
 
 # error rates (epsilon) hyperprior
 groupW          ~ dbeta(1.1, 10.9)  # mode for lapse rate
@@ -40,7 +40,7 @@ epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     epsilon[p]  ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-    
+
     # using reparameterisation to avoid funnel of hell
     alpha_offset[p] ~ dnorm(0,1) T(0,)
     alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
@@ -57,6 +57,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t])
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/exponential_models/separateExp1.jags b/ddToolbox/models/parametric_models/exponential_models/separateExp1.jags
index 0299887c..a52ed640 100644
--- a/ddToolbox/models/parametric_models/exponential_models/separateExp1.jags
+++ b/ddToolbox/models/parametric_models/exponential_models/separateExp1.jags
@@ -8,7 +8,7 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 K_MEAN          <- 0.01             # <---- currently guesstimating
 K_PRECISION     <- 1/(0.5^2)        # <---- currently guesstimating
@@ -52,6 +52,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t])
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/exponential_power_models/ExponentialPower.m b/ddToolbox/models/parametric_models/exponential_power_models/ExponentialPower.m
index 2b5cd659..1d482264 100644
--- a/ddToolbox/models/parametric_models/exponential_power_models/ExponentialPower.m
+++ b/ddToolbox/models/parametric_models/exponential_power_models/ExponentialPower.m
@@ -4,23 +4,23 @@
 
 		function obj = ExponentialPower(data, varargin)
 			obj = obj@SubjectiveTimeModel(data, varargin{:});
-            
+
             obj.dfClass = @DF_ExponentialPower;
             obj.subjectiveTimeFunctionFH = @SubjectiveTimePowerFunction;
-            
+
 			% Create variables
 			obj.varList.participantLevel = {'k','logtau','alpha','epsilon'};
-			obj.varList.monitored = {'k','logtau','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'k','logtau','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
             obj.varList.discountFunctionParams(1).name = 'k';
             obj.varList.discountFunctionParams(1).label = 'discount rate, $k$';
             obj.varList.discountFunctionParams(2).name = 'logtau';
             obj.varList.discountFunctionParams(2).label = '$\log(\tau)$';
-            
+
             obj.plotOptions.dataPlotType = '2D';
 		end
-    
+
     end
-    
+
     methods (Hidden = true)
         function dispModelInfo(obj)
             display('Discount function: V = reward * exp(-k*(delay^tau))')
diff --git a/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.jags b/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.jags
index e7b11697..27c66f3f 100644
--- a/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.jags
+++ b/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.jags
@@ -8,7 +8,7 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 K_MEAN            <- 0                     # <---- currently guesstimating
 K_PRECISION       <- 1/(0.01)              # <---- currently guesstimating
@@ -45,7 +45,7 @@ epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     epsilon[p]  ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-    
+
     # using reparameterisation to avoid funnel of hell
     alpha_offset[p] ~ dnorm(0,1) T(0,)
     alpha[p] <- alpha_mu + alpha_offset[p] * alpha_sigma
@@ -62,6 +62,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.jags b/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.jags
index cb7bc1b3..f13f2ef4 100644
--- a/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.jags
+++ b/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.jags
@@ -8,7 +8,7 @@ model{
 
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 K_MEAN          <- 0                     # <---- currently guesstimating
 K_PRECISION     ~ dgamma(1, 1)     #sigma=0.1         # <---- currently guesstimating
@@ -54,6 +54,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m
index 563f3fd1..564808ec 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m
+++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m
@@ -1,28 +1,28 @@
 classdef (Abstract) Hyperbolic1MagEffect < Parametric
 	%Hyperbolic1MagEffect  Hyperbolic1MagEffect is a subclass of Model for examining the 1-parameter hyperbolic discounting function.
-	
+
 	methods (Access = public)
-		
+
 		function obj = Hyperbolic1MagEffect(data, varargin)
 			obj = obj@Parametric(data, varargin{:});
-			
+
 			obj.dfClass = @DF_HyperbolicMagnitudeEffect;
-			
+
 			% Create variables
 			obj.varList.participantLevel = {'m', 'c','alpha','epsilon'};
-			obj.varList.monitored = {'m', 'c','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'm', 'c','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
 			obj.varList.discountFunctionParams(1).name = 'm';
 			obj.varList.discountFunctionParams(1).label = 'slope, $m$';
 			obj.varList.discountFunctionParams(2).name = 'c';
 			obj.varList.discountFunctionParams(2).label = 'intercept, $c$';
-			
+
 			obj.plotOptions.dataPlotType = '3D';
 		end
-        
+
     end
-    
+
     methods (Access = private)
-        
+
         % MIDDLE-MAN METHOD
         function logk = getLogDiscountRate(obj, reward, index, varargin)
             %% Set up discount function
@@ -31,44 +31,44 @@
             discountFunction = obj.dfClass('samples', dfSamples);
             %             % inject a DataFile object into the discount function
             %             discountFunction.data = obj.data.getExperimentObject(index);
-            
+
             %% Do the actual call
             logk = discountFunction.getLogDiscountRate(reward, index, varargin{:});
         end
-        
+
     end
-        
+
     methods (Hidden = true)
         function dispModelInfo(obj)
             % TODO: Display the discount function
             %display('Discount function: V = reward * exp(-k*(delay^tau))')
         end
     end
-        
-        
+
+
     % ==========================================================================
     % ==========================================================================
     % PLOTTING
     % ==========================================================================
     % ==========================================================================
-        
+
 	methods (Access = public)
-    
+
 		% OVERRIDDING THIS METHOD FROM A SUPERCLASS
 		function plotExperimentOverviewFigure(obj, ind)
             %model.plotExperimentOverviewFigure(N) Creates a multi-panel figure
             %   model.plotExperimentOverviewFigure(N) creates a multi-panel figure
             %   corresponding to experiment N, where N is an integer.
-            
+
 			latex_fig(12, 14, 3);
 			h = layout([1 2 3 4 5 6]);
 			opts.pointEstimateType = obj.plotOptions.pointEstimateType;
-			
+
             % % create cell array
 			% discountFunctionVariables = {obj.varList.discountFunctionParams.name};
 			% responseErrorVariables = {obj.varList.responseErrorParams.name};
-			
-            
+
+
 			% %% Plot 1: density plot of (alpha, epsilon)
 			% obj.coda.plot_bivariate_distribution(h(1),...
 			% 	responseErrorVariables(1),...
@@ -76,7 +76,7 @@ function plotExperimentOverviewFigure(obj, ind)
 			% 	ind,...
 			% 	opts);
 			obj.plotPosteriorErrorParams(ind, 'axisHandle', h(1))
-			
+
 			% % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 			% % TODO #166 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 			% %%  Set up psychometric function
@@ -93,8 +93,8 @@ function plotExperimentOverviewFigure(obj, ind)
 			% % TODO #166 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 			% % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             obj.plotPsychometricFunction(ind, 'axisHandle', h(2))
-            
-			
+
+
 			% % TODO: this checking needs to be implemented in a smoother, more robust way
 			% if ~isempty(dfSamples) || ~any(isnan(dfSamples))
 			% 	%% Bivariate density plot of discounting parameters
@@ -105,20 +105,20 @@ function plotExperimentOverviewFigure(obj, ind)
 			% 		opts);
 			% end
 			obj.plotPosteriorDiscountFunctionParams(ind, 'axisHandle', h(3))
-			
-            
+
+
 			% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 			% TODO #166 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            
+
             % a list of reward values we are interested in
             rewards = [10, 100, 500]; % <----  TODO: inject this
             % create colours for colour coding of conditional k plotStyle
             col = linspace(0.1, 0.9, numel(rewards));
-            
+
 			%% Set up magnitude effect function -----------------------
 			discountFunctionVariables = obj.getDiscountFunctionVariables();
 			dfSamples = obj.coda.getSamplesAtIndex_asStochastic(ind, discountFunctionVariables);
-			
+
 			me = MagnitudeEffectFunction('samples', dfSamples);
 			me.maxRewardValue = obj.data.getMaxRewardValue(ind);
 			% plot magnitude effect
@@ -129,36 +129,36 @@ function plotExperimentOverviewFigure(obj, ind)
 			for n=1:numel(rewards)
 				vline(rewards(n), 'Color', [col(n) col(n) col(n)], 'LineWidth', 2);
 			end
-			
-			
+
+
 			subplot(h(5)) % -----------------------------------------------
 			%title('P(log(k) | reward)')
 			discountFunction = obj.dfClass('samples', dfSamples);
-			
+
 			discountFunction.getLogDiscountRate(rewards, ind ,...
 				'plot', true,...
 				'plot_mode', 'conditional_only');
-			
+
 			% TODO: this checking needs to be implemented in a
 			% smoother, more robust way
 			if ~isempty(dfSamples) || ~any(isnan(dfSamples))
 				subplot(h(6)) % -------------------------------------------
-				
+
 				plotOptions.pointEstimateType = obj.plotOptions.pointEstimateType;
 				plotOptions.dataPlotType = obj.plotOptions.dataPlotType;
 				plotOptions.timeUnits = obj.timeUnits;
 				%plotOptions.plotMode = p.Results.plot_mode;
 				plotOptions.maxRewardValue = obj.data.getMaxRewardValue(ind);
 				plotOptions.maxDelayValue = obj.data.getMaxDelayValue(ind);
-				
+
 				discountFunction.plot(plotOptions);
 			end
 			% TODO #166 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 			% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-			
-			
+
+
 		end
-		
+
 	end
-	
+
 end
diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.jags b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.jags
index e8792ec0..de46073f 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES
 
 # priors over group M mean
 groupMmu_MEAN          <- -0.243
@@ -34,7 +34,7 @@ for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     m[p] <- groupMmu + m_offset[p] * groupMsigma
     c_offset[p] ~ dnorm(0,1)
     c[p] <- groupCmu + c_offset[p] * groupCsigma
-} 
+}
 
 # MODEL-SPECIFIC: CALCULATION OF PRESENT SUBJECTIVE VALUES
 for (t in 1:length(ID)) {
@@ -62,7 +62,7 @@ epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
 	epsilon[p] ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-	
+
     # using reparameterisation to avoid funnel of hell
     alpha_offset[p] ~ dnorm(0,1) T(0,)
     alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
@@ -78,6 +78,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEmvnorm.jags b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEmvnorm.jags
index 8f51d2fd..9c101dcb 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEmvnorm.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEmvnorm.jags
@@ -55,6 +55,7 @@ for (t in 1:length(ID)) {
 
 	# response likelihood
 	R[t]         ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 	Rpostpred[t] ~ dbern(P[t]) # posterior predicted response
 }
 
diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.jags b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.jags
index f48b5776..df01a9f6 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES
 
 # priors over group M mean
 groupMmu_MEAN          <- -0.243
@@ -77,6 +77,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.jags b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.jags
index 93129c1d..bd230c36 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 # slope (m)
 groupMmu        <- -0.243
@@ -48,7 +48,7 @@ epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
 	epsilon[p]  ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)
-    
+
     # using reparameterisation to avoid funnel of hell
     alpha_offset[p] ~ dnorm(0,1) T(0,)
     alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
@@ -65,6 +65,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.jags b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.jags
index a1ef589d..ebb494cc 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 for (p in 1:nRealExperimentFiles){
 	m[p]       ~ dnorm(-0.243, 1/(0.5^2))  # dnorm(-0.243, 1/(0.072^2))
@@ -40,6 +40,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m b/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m
index d572acf5..67d797f7 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m
+++ b/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m
@@ -10,7 +10,7 @@
 
 			% Create variables
 			obj.varList.participantLevel = {'logk','alpha','epsilon'};
-			obj.varList.monitored = {'logk','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'logk','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
     	obj.varList.discountFunctionParams(1).name = 'logk';
 			obj.varList.discountFunctionParams(1).label = '$\log(k)$';
 
@@ -24,5 +24,5 @@ function dispModelInfo(obj)
             display('Discount function: V = 1 / (1+k*delay)')
         end
     end
-    
+
 end
diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.jags b/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.jags
index fe00cc1e..625329fe 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES
 
 groupLogKmu       ~ dnorm(log(1/50),1/(2.5^2))
 groupLogKsigma    ~ dexp(0.5)
@@ -54,6 +54,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.jags b/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.jags
index 7bf3c767..8e1905ed 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 # mean half-life of 50 days from a sample from my lab, with a std of ~2.5. Note k = 1/halflife.
 groupLogKmu <- log(1/50)
@@ -55,6 +55,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK.jags b/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK.jags
index 00b4d257..c7460dfe 100644
--- a/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK.jags
+++ b/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 # mean half-life of 50 days from a sample from my lab, with a std of ~2.5. Note k = 1/halflife.
 logk_MEAN      <- log(1/50)
@@ -40,6 +40,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperboloid_models/Hyperboloid.m b/ddToolbox/models/parametric_models/hyperboloid_models/Hyperboloid.m
index 96812504..291b0596 100644
--- a/ddToolbox/models/parametric_models/hyperboloid_models/Hyperboloid.m
+++ b/ddToolbox/models/parametric_models/hyperboloid_models/Hyperboloid.m
@@ -8,10 +8,10 @@
 
 			obj.dfClass = @DF_Hyperboloid;
             obj.subjectiveTimeFunctionFH = @SubjectiveTimeWebber;
-            
+
 			% Create variables
 			obj.varList.participantLevel = {'logk','S','alpha','epsilon'};
-			obj.varList.monitored = {'logk','S','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.monitored = {'log_lik', 'logk','S','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
             obj.varList.discountFunctionParams(1).name = 'logk';
 			obj.varList.discountFunctionParams(1).label = '$\log(k)$';
             obj.varList.discountFunctionParams(2).name = 'S';
@@ -27,5 +27,5 @@ function dispModelInfo(obj)
             display('Discount function: V = 1 / (1+k*delay)^S')
         end
     end
-    
+
 end
diff --git a/ddToolbox/models/parametric_models/hyperboloid_models/hierarchicalHyperboloid.jags b/ddToolbox/models/parametric_models/hyperboloid_models/hierarchicalHyperboloid.jags
index 24e2bd57..c75de388 100644
--- a/ddToolbox/models/parametric_models/hyperboloid_models/hierarchicalHyperboloid.jags
+++ b/ddToolbox/models/parametric_models/hyperboloid_models/hierarchicalHyperboloid.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = YES
 
 groupLogKmu       ~ dnorm(log(1/50),1/(2.5^2))
 groupLogKsigma    ~ dexp(0.5)
@@ -16,7 +16,7 @@ for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
     # using reparameterisation to avoid funnel of hell
     logk_offset[p] ~ dnorm(0,1)
     logk[p] <- groupLogKmu + logk_offset[p] * groupLogKsigma
-    
+
     S[p]       ~ dnorm(S_MEAN, S_PRECISION) T(-1,)
 }
 
@@ -59,6 +59,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperboloid_models/mixedHyperboloid.jags b/ddToolbox/models/parametric_models/hyperboloid_models/mixedHyperboloid.jags
index 95380d89..1fc1755d 100644
--- a/ddToolbox/models/parametric_models/hyperboloid_models/mixedHyperboloid.jags
+++ b/ddToolbox/models/parametric_models/hyperboloid_models/mixedHyperboloid.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 # mean half-life of 50 days from a sample from my lab, with a std of ~2.5. Note k = 1/halflife.
 groupLogKmu <- log(1/50)
@@ -16,7 +16,7 @@ for (p in 1:nRealExperimentFiles){
     # using reparameterisation to avoid funnel of hell
     logk_offset[p] ~ dnorm(0,1)
     logk[p] <- groupLogKmu + logk_offset[p] * groupLogKsigma
-    
+
     S[p]       ~ dnorm(S_MEAN, S_PRECISION) T(-1,)
 }
 
@@ -59,6 +59,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+  log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/ddToolbox/models/parametric_models/hyperboloid_models/separateHyperboloid.jags b/ddToolbox/models/parametric_models/hyperboloid_models/separateHyperboloid.jags
index 74e7cd21..edaf50e4 100644
--- a/ddToolbox/models/parametric_models/hyperboloid_models/separateHyperboloid.jags
+++ b/ddToolbox/models/parametric_models/hyperboloid_models/separateHyperboloid.jags
@@ -4,7 +4,7 @@
 model{
 
 # DISCOUNT FUNCTION PARAMETERS =================================================
-# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
+# RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO
 
 # mean half-life of 50 days from a sample from my lab, with a std of ~2.5. Note k = 1/halflife.
 logk_MEAN      <- log(1/50)
@@ -43,6 +43,7 @@ for (t in 1:length(ID)) {
 # response likelihood
 for (t in 1:length(ID)) {
 	R[t]  ~ dbern(P[t]) # likelihood of actual response
+	log_lik[t] <- logdensity.bern(R[t], P[t])
 }
 
 # POSTERIOR PREDICTION
diff --git a/demo/demo_WAIC.m b/demo/demo_WAIC.m
new file mode 100644
index 00000000..b70e9059
--- /dev/null
+++ b/demo/demo_WAIC.m
@@ -0,0 +1,65 @@
+% demo_WAIC
+
+%% setup
+toolbox_path = '~/git-local/delay-discounting-analysis/ddToolbox';
+addpath(toolbox_path)
+datapath = '~/git-local/delay-discounting-analysis/demo/datasets/non_parametric';
+
+addpath(toolbox_path)
+ddAnalysisSetUp();
+
+% Running this multiple times will result in slightly different model
+% comparison results. For real model comparison contexts, eg for research
+% papers, then you'd want to ensire your posteriors are good
+% approximations. Best way to start is by increasing the number of MCMC
+% samples.
+mcmc_params = struct('nsamples', 10000,...
+		'nchains', 4,...
+		'nburnin', 2000);
+
+%% Set up the data to be analysed
+myData = Data(datapath, 'files', allFilesInFolder(datapath, 'txt'));
+	
+%% Fit multiple models to the dataset
+
+modelA = ModelHierarchicalLogK(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelA'),...
+	'mcmcParams', mcmc_params);
+
+modelB = ModelHierarchicalExp1(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelB'),...
+	'mcmcParams', mcmc_params);
+	
+modelC = ModelHierarchicalHyperboloid(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelC'),...
+	'mcmcParams', mcmc_params);
+
+modelD = ModelHierarchicalMEUpdated(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelD'),...
+	'mcmcParams', mcmc_params);
+
+modelE = ModelHierarchicalBetaDelta(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelE'),...
+	'mcmcParams', mcmc_params);
+
+modelF = ModelHierarchicalEbertPrelec(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelF'),...
+	'mcmcParams', mcmc_params);
+
+%% Examine WAIC stats for models
+waic = [modelA.WAIC_stats,...
+	modelB.WAIC_stats,...
+	modelC.WAIC_stats,...
+	modelD.WAIC_stats,...
+	modelE.WAIC_stats,...
+	modelF.WAIC_stats];
+
+waic_comparison_table = waic.compare()
+
+waic.plot()
\ No newline at end of file
diff --git a/docs/index.md b/docs/index.md
index 391b6ffa..1c3c929e 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -11,6 +11,7 @@ Vincent, B. T. (2016) **[Hierarchical Bayesian estimation and hypothesis testing
 
 ## Tutorials
 - [Run the demo code](tutorial/run_demo_code)
+- [Model comparison](tutorial/model_comparison)
 
 ## How to
 - [Install the toolbox](howto/install.md)
diff --git a/docs/tutorial/WAICplot.png b/docs/tutorial/WAICplot.png
new file mode 100644
index 00000000..505af424
Binary files /dev/null and b/docs/tutorial/WAICplot.png differ
diff --git a/docs/tutorial/WAICtable.png b/docs/tutorial/WAICtable.png
new file mode 100644
index 00000000..7aeefedc
Binary files /dev/null and b/docs/tutorial/WAICtable.png differ
diff --git a/docs/tutorial/model_comparison.md b/docs/tutorial/model_comparison.md
new file mode 100644
index 00000000..f0779b41
--- /dev/null
+++ b/docs/tutorial/model_comparison.md
@@ -0,0 +1,112 @@
+# Tutorial: model comparison with WAIC
+
+Model comparison is tricky, so care needs to be taken here.
+
+Below we describe an example where you might be interested in working out which discount function is the best model for a given dataset.  That is, to evaluate how well a number discount functions can account for all data in a dataset - which is different from asking which is the best discount function for each individual participant in a dataset.
+
+Note that the WAIC model comparison functionality was implemented with reference to both Gelman et al (2013) and McElreath (2016).
+
+The code is available in `demo_WAIC.m`, but let's walk through it. First some setup code
+
+```matlab
+toolbox_path = '~/git-local/delay-discounting-analysis/ddToolbox';
+addpath(toolbox_path)
+datapath = '~/git-local/delay-discounting-analysis/demo/datasets/non_parametric';
+
+addpath(toolbox_path)
+ddAnalysisSetUp();
+```
+
+Now we'll set up MCMC parameters
+
+```matlab
+mcmc_params = struct('nsamples', 10000,...
+		'nchains', 4,...
+		'nburnin', 2000);
+```
+
+Then create a `Data` object that we'll analyse with multiple models.
+
+```matlab
+myData = Data(datapath, 'files', allFilesInFolder(datapath, 'txt'));
+```
+
+And conduct parameter estimation on multiple models. In a real context you'd perhaps want to use mode descriptive names for the the models and output folders, but this will be enough to get the right idea. This will take a little while to compute.
+
+```matlab
+modelA = ModelHierarchicalLogK(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelA'),...
+	'mcmcParams', mcmc_params);
+
+modelB = ModelHierarchicalExp1(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelB'),...
+	'mcmcParams', mcmc_params);
+
+modelC = ModelHierarchicalHyperboloid(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelC'),...
+	'mcmcParams', mcmc_params);
+
+modelD = ModelHierarchicalMEUpdated(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelD'),...
+	'mcmcParams', mcmc_params);
+
+modelE = ModelHierarchicalBetaDelta(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelE'),...
+	'mcmcParams', mcmc_params);
+
+modelF = ModelHierarchicalEbertPrelec(...
+	myData,...
+	'savePath', fullfile(pwd,'output','modelF'),...
+	'mcmcParams', mcmc_params);
+
+```
+
+At this point, you'd want to do your due diligence to ensure that the parameter estimates are kosher approximations of the true posterior. So, you could look at the Rhat chain convergence diagnostics, visually examine the MCMC chains, and also the model 'fits' and posterior predictive plots. Without this step, you cannot be really sure that the model comparison provides meaningful information.
+
+But assuming we've done that, each model contains a `WAIC` object, which we can access with the `model.WAIC_stats` property. We'll create an array of these WAIC objects from each of the models.
+
+```matlab
+waic = [modelA.WAIC_stats,...
+	modelB.WAIC_stats,...
+	modelC.WAIC_stats,...
+	modelD.WAIC_stats,...
+	modelE.WAIC_stats,...
+	modelF.WAIC_stats];
+```
+
+We can now compare the models with:
+
+```matlab
+waic_comparison_table = waic.compare()
+```
+
+![](WAICtable.png)
+
+with the following columns. Readers are referred to McElreath (2016) for more details
+- **WAIC:** Smaller values are better, lower out of sample deviance.
+- **pWAIC:** Estimated effective number of parameters.
+- **dWAIC:** Difference between each model's WAIC and the model with lowest WAIC.
+- **weight:** is the Akaike weight
+- **SE:** is the standard error of the WAIC value over data points (trials)
+- **dSE:** standard error of the difference in WAIC values and the top-ranked model.
+
+We can also produce a nice WAIC plot (suitable for dropping into a paper) with
+
+```matlab
+waic.plot()
+```
+
+![](WAICplot.png)
+
+Again, readers are referred to McElreath (2016) for a thorough explanation of this plot, the style of which is directly copied from that book.
+
+## Reference
+
+Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis, Third Edition. CRC Press.
+
+McElreath, R. (2016). Statistical Rethinking: A Bayesian Course with Examples in R and Stan. CRC Press.