Merge pull request #176 from drbenvincent/dev

from Dev

ddToolbox/Data/Data.m

@@ -32,6 +32,7 @@
 			p.addParameter('files',[],@(x) iscellstr(x)|ischar(x));
 			p.parse(dataFolder, varargin{:});
 
+			assert(~isempty(p.Results.files), 'no filenames provided under ''files'' input argument')
 			try
 				table();
 			catch

ddToolbox/models/Model.m

@@ -92,6 +92,7 @@
 			obj = obj.postSamplingActivities();
 		end
 
+
 		function obj = postSamplingActivities(obj)
 
 			%% Post-sampling activities (for model sub-classes) -----------
@@ -101,26 +102,39 @@
 			obj = obj.calcDerivedMeasures();
 
 			%% Post-sampling activities (common to all models) ------------
+			% posterior prediction calculation
 			obj.postPred = PosteriorPrediction(obj.coda, obj.data, obj.observedData);
 
+			% calc and export convergence summary and parameter estimates
+			obj.export();
 
-			% TODO: This should be a method of CODA
- 			convergenceSummary(obj.coda.getStats('Rhat',[]), obj.plotOptions.savePath, obj.data.getIDnames('all'))
-
-			exporter = ResultsExporter(obj.coda, obj.data, obj.postPred.postPred, obj.varList, obj.plotOptions);
-			exporter.printToScreen();
-			exporter.export(obj.plotOptions.savePath, obj.plotOptions.pointEstimateType);
-			% TODO ^^^^ avoid this duplicate use of pointEstimateType
-
+			% plot or not
 			if ~strcmp(obj.plotOptions.shouldPlot,'no')
 				% TODO: Allow public calls of obj.plot to specify options.
 				% At the moment the options need to be provided on Model
 				% object construction
 				obj.plot()
 			end
-
+			
 			obj.tellUserAboutPublicMethods()
 		end
+
+		function export(obj)
+			% TODO: This should be a method of CODA
+			% TODO: better function name. What does it do? Calculate or
+			% export, or both?
+			convergenceSummary(obj.coda.getStats('Rhat',[]),...
+				obj.plotOptions.savePath,...
+				obj.data.getIDnames('all'))
+
+			exporter = ResultsExporter(obj.coda, obj.data,...
+				obj.postPred.postPred,...
+				obj.varList,...
+				obj.plotOptions);
+			exporter.printToScreen();
+			exporter.export(obj.plotOptions.savePath, obj.plotOptions.pointEstimateType);
+			% TODO ^^^^ avoid this duplicate use of pointEstimateType
+		end
 
 		%% Public MIDDLE-MAN METHODS
 

ddToolbox/models/nonparametric_models/DF_SLICE_PsychometricFunction.m

@@ -29,8 +29,9 @@
 
 		function plot(obj)
 			% 
-			x = [-2:0.01:2];
-
+			%x = [-2:0.01:2];
+			x = [0:0.01:2];
+
 			try
 				plot(x, obj.eval(x, 'nExamples', 100), '-', 'Color',[0.5 0.5 0.5 0.1])
 			catch

ddToolbox/models/nonparametric_models/separateNonParametric.jags

@@ -21,22 +21,41 @@
 model{
 
 groupALPHAmu <- 0           # TODO: NEEDS ATTENTION
-groupALPHAsigma <-0.01      # TODO: NEEDS ATTENTION
+groupALPHAsigma <-0.1      # TODO: NEEDS ATTENTION
 
 epsilon_alpha <- 1+1
-epsilon_beta <- 1+100
+epsilon_beta <- 1+200
+
+#SIGMA ~ dgamma(1,1)
+SIGMA <- 0.05
 
 # priors
+
+# precision of the indifference point (Rstar) is SIGMA multiplied by the number of time steps since the previous delay. This is how our uncertainty increase over time.
+prec[1] <- 1/((SIGMA * uniqueDelays[1])^2)
+for (d in 2:length(uniqueDelays)) {
+    prec[d] <- 1/((SIGMA * (uniqueDelays[d]-uniqueDelays[d-1]))^2)
+}
+
 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,)
 
     # PRIOR OVER DISCOUNT FRACTION
-    # prior over A/B should be centered on 1 if we are expecting discounting and anti-discounting, and can be truncates
-    # prior over A-B should be centered on 0 if we are expecting discounting and anti-discounting, and can be truncates
-	for (d in 1:length(uniqueDelays)) {
-        Rstar[p,d] ~ dnorm(1, 1/ 1^2) T(0,)   
+    # assume indifference point (RStar) at the first delay (after delay=0) is centered on 1
+    mu[p,1] <- 1
+    # assume indifference point at all subsequent delays are centered on the previous indifference point
+	for (d in 2:length(uniqueDelays)) {
+        mu[p,d] <- mu[p,d-1]
 	}
+
+    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  
+    }
+
 }
 
 for (t in 1:length(ID)) {

ddToolbox/models/parametric_models/bens_new_model/ModelSeparateExpPower.m

@@ -20,8 +20,8 @@
 			% Generate initial values of the root nodes
 			nExperimentFiles = obj.data.getNExperimentFiles();
 			for chain = 1:nchains
-				initialParams(chain).k = unifrnd(0.01, 0.5, [nExperimentFiles,1]);
-                initialParams(chain).tau = unifrnd(0.01, 2, [nExperimentFiles,1]);
+				initialParams(chain).k = unifrnd(-1, 1, [nExperimentFiles,1]);
+                initialParams(chain).tau = unifrnd(0.01, 1.2, [nExperimentFiles,1]);
 				initialParams(chain).epsilon = 0.01 + rand([nExperimentFiles,1])./10;
 				initialParams(chain).alpha = abs(normrnd(0.01,5,[nExperimentFiles,1]));
 			end

ddToolbox/models/parametric_models/bens_new_model/mixedExpPower.jags

@@ -33,6 +33,7 @@ for (t in 1:length(ID)) {
 # NOTE: because we have 2 discounting parameters, I've found that more restrictive priors need to be placed on the error parameters
 alpha_mu        ~ dnorm(0,1/(10^2)) T(0,)
 alpha_precision ~ dgamma(0.01,0.01)
+alpha_sigma     <- sqrt(1/alpha_precision)
 
 # error rates (epsilon) hyperprior
 groupW          ~ dbeta(1.1, 10.9)  # mode for lapse rate
@@ -43,8 +44,11 @@ epsilon_alpha <- groupW*(groupK-2)+1
 epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
-	alpha[p]    ~ dnorm(alpha_mu, alpha_precision) T(0,)
     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
 }
 
 

ddToolbox/models/parametric_models/bens_new_model/separateExpPower.jags

@@ -11,15 +11,16 @@ model{
 # RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
 
 K_MEAN          <- 0                     # <---- currently guesstimating
-K_PRECISION     <- 1/(0.1^2)     #sigma=0.1         # <---- currently guesstimating
+K_PRECISION     ~ dgamma(0.01, 0.01)     #sigma=0.1         # <---- currently guesstimating
 
 for (p in 1:nRealExperimentFiles){       # no +1 because no shrinkage hyperprior
-	k[p]        ~ dnorm(K_MEAN, K_PRECISION)
-    tau[p]      ~ dexp(1)
+	k[p]        ~ dt(K_MEAN, K_PRECISION, 1)
+    tau[p]      ~ dnorm(1, 1/0.2^2) T(0.0001, 1.5) #<------ check these truncation bounds
 }
 
 # MODEL-SPECIFIC: CALCULATION OF PRESENT SUBJECTIVE VALUES
 for (t in 1:length(ID)) {
+	#VA[t] <- A[t] # assuming DA=0
     VA[t] <- A[t] * (exp( -k[ID[t]] * (DA[t]^tau[ID[t]]) ) )
 	VB[t] <- B[t] * (exp( -k[ID[t]] * (DB[t]^tau[ID[t]]) ) )
 }

ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.jags

@@ -12,9 +12,12 @@ model{
 
 K_MEAN       ~ dnorm(0.01, 1/(0.5^2))
 K_PRECISION  ~ dgamma(0.001,0.001)
+K_SIGMA      <- sqrt(1/K_PRECISION)
 
 for (p in 1:(nRealExperimentFiles+1)){      # +1 for unobserved participant
-	k[p]     ~ dnorm(K_MEAN, K_PRECISION)
+    # using reparameterisation to avoid funnel of hell
+    k_offset[p] ~ dnorm(0,1)
+    k[p] <- K_MEAN + k_offset[p] * K_SIGMA
 }
 
 # MODEL-SPECIFIC: CALCULATION OF PRESENT SUBJECTIVE VALUES
@@ -39,8 +42,11 @@ epsilon_alpha <- groupW*(groupK-2)+1
 epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
-	alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
     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
 }
 
 

ddToolbox/models/parametric_models/exponential_models/mixedExp1.jags

@@ -10,7 +10,7 @@ model{
 # DISCOUNT FUNCTION PARAMETERS =================================================
 # RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
 
-K_MEAN          <- 0                     # <---- currently guesstimating
+K_MEAN          <- 0.01                  # <---- currently guesstimating
 K_PRECISION     <- 1/(0.5^2)             # <---- currently guesstimating
 
 for (p in 1:nRealExperimentFiles){       # no +1 because no shrinkage hyperprior
@@ -39,8 +39,11 @@ epsilon_alpha <- groupW*(groupK-2)+1
 epsilon_beta <- (1-groupW)*(groupK-2)+1
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
-	alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
     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
 }
 
 

ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.jags

@@ -29,9 +29,12 @@ groupCmu        ~ dnorm(0, 1/(10000^2))
 groupCsigma     ~ dunif(0, 10000)
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
-	m[p]        ~ dnorm(groupMmu, 1/(groupMsigma^2))
-	c[p]        ~ dnorm(groupCmu, 1/(groupCsigma^2))
-}
+    # using reparameterisation to avoid funnel of hell
+    m_offset[p] ~ dnorm(0,1)
+    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)) {
@@ -59,7 +62,10 @@ 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)
-	alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
+
+    # using reparameterisation to avoid funnel of hell
+    alpha_offset[p] ~ dnorm(0,1) T(0,)
+    alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
 }
 
 # MODEL IN-SPECIFIC CODE BELOW... SHOULD NOT CHANGE ACROSS MODELS ==============

ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.jags

@@ -29,8 +29,11 @@ groupCmu        ~ dnorm(0, 1/(100^2)) ## UPDATED SINCE PAPER
 groupCsigma     ~ dunif(0, 10) ## UPDATED SINCE PAPER
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
-	m[p]        ~ dnorm(groupMmu, 1/(groupMsigma^2))
-	c[p]        ~ dnorm(groupCmu, 1/(groupCsigma^2))
+    # using reparameterisation to avoid funnel of hell
+    m_offset[p] ~ dnorm(0,1)
+    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
@@ -59,7 +62,9 @@ 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)
-	alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
+    # using reparameterisation to avoid funnel of hell
+    alpha_offset[p] ~ dnorm(0,1) T(0,)
+    alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
 }
 
 # MODEL IN-SPECIFIC CODE BELOW... SHOULD NOT CHANGE ACROSS MODELS ==============

ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.jags

@@ -48,7 +48,10 @@ 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)
-	alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
+
+    # using reparameterisation to avoid funnel of hell
+    alpha_offset[p] ~ dnorm(0,1) T(0,)
+    alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
 }
 
 

ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.jags

@@ -7,8 +7,8 @@ model{
 # RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
 
 for (p in 1:nRealExperimentFiles){
-	m[p]       ~ dnorm(-0.243, 1/(100^2))  # dnorm(-0.243, 1/(0.072^2))
-	c[p]       ~ dnorm(0, 1/(1000^2))      # dnorm(0, 1/(1000^2))
+	m[p]       ~ dnorm(-0.243, 1/(0.5^2))  # dnorm(-0.243, 1/(0.072^2))
+	c[p]       ~ dnorm(0, 1/(10^2))      # dnorm(0, 1/(1000^2))
 }
 
 for (t in 1:length(ID)) {
@@ -22,8 +22,8 @@ for (t in 1:length(ID)) {
 }
 
 # RESPONSE ERROR PARAMETERS ====================================================
-epsilon_alpha <- 1.1
-epsilon_beta <- 10.9
+epsilon_alpha <- 1+1
+epsilon_beta <- 1+10
 for (p in 1:nRealExperimentFiles){
 	alpha[p]   ~ dexp(0.01)
 	epsilon[p] ~ dbeta(epsilon_alpha , epsilon_beta ) T(,0.5)

ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.stan

@@ -64,7 +64,7 @@ transformed parameters {
 model {
   m        ~ normal(-0.243, 0.5);
   c        ~ normal(0, 10);
-  epsilon  ~ beta(1.1, 10.9);
+  epsilon  ~ beta(1+1, 1+10);
   alpha    ~ exponential(0.01);
 
   R ~ bernoulli(P);

ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.jags

@@ -10,8 +10,9 @@ groupLogKmu       ~ dnorm(log(1/50),1/(2.5^2))
 groupLogKsigma    ~ dexp(0.5)
 
 for (p in 1:(nRealExperimentFiles+1)){ # +1 for unobserved participant
-	# small constant (0.0001) added below to avoid numerical problems.
-	logk[p]     ~ dnorm(groupLogKmu, 1/((groupLogKsigma+0.0001)^2))
+    # using reparameterisation to avoid funnel of hell
+    logk_offset[p] ~ dnorm(0,1)
+    logk[p] <- groupLogKmu + logk_offset[p] * groupLogKsigma
 }
 
 # MODEL-SPECIFIC: CALCULATION OF PRESENT SUBJECTIVE VALUES
@@ -36,7 +37,10 @@ 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)
-	alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
+
+    # using reparameterisation to avoid funnel of hell
+    alpha_offset[p] ~ dnorm(0,1) T(0,)
+    alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
 }
 
 

ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.jags

@@ -11,7 +11,9 @@ groupLogKmu <- log(1/50)
 groupLogKsigma <- 2.5
 
 for (p in 1:nRealExperimentFiles){
-	logk[p]     ~ dnorm(groupLogKmu, 1/(groupLogKsigma^2))
+    # using reparameterisation to avoid funnel of hell
+    logk_offset[p] ~ dnorm(0,1)
+    logk[p] <- groupLogKmu + logk_offset[p] * groupLogKsigma
 }
 
 
@@ -30,7 +32,10 @@ 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)
-	alpha[p]    ~ dnorm(groupALPHAmu, 1/(groupALPHAsigma^2)) T(0,)
+
+    # using reparameterisation to avoid funnel of hell
+    alpha_offset[p] ~ dnorm(0,1) T(0,)
+    alpha[p] <- groupALPHAmu + alpha_offset[p] * groupALPHAsigma
 }
 
 

ddToolbox/sampleWithMatlabStan.m

@@ -1,5 +1,5 @@
 function codaObject = sampleWithMatlabStan(...
-	modelFilename, observedData, mcmcparams, ~, ~)
+	modelFilename, observedData, mcmcparams, initialParameters, ~)
 
 assert(ischar(modelFilename))
 assert(isstruct(observedData))
@@ -17,11 +17,16 @@
 toc
 % ++++++++++++++++++++++++++++++++++++++++++++++++++
 
+% % convert initialParameters. % TODO: make this general
+% init.groupW				= [initialParameters(:).groupW]';
+% init.groupALPHAmu		= [initialParameters(:).groupALPHAmu]';
+% init.groupALPHAsigma	= [initialParameters(:).groupALPHAsigma]';
+
 %% Get our sampler to sample
 display('SAMPLING STAN MODEL...')
 tic
 obj.stanFit = stan_model.sampling(...
-	'data', observedData,...
+	'data', observedData,... %'init', init,...
 	'warmup', mcmcparams.nburnin,...	% warmup = burn-in
 	'iter', ceil( mcmcparams.nsamples / mcmcparams.nchains),...		% iter = number of MCMC samples
 	'chains', mcmcparams.nchains,...