#147 bug fix + estimate log(tau) rather than tau

Error was in ExponentialPower.m. It was using the wrong subjective time function.

ddToolbox/models/parametric_models/exponential_power_models/DF_ExponentialPower.m

@@ -16,11 +16,12 @@
 	methods (Static, Access = protected)
 
 		function y = function_evaluation(x, theta)
+			tau = exp(theta.logtau);
 			if verLessThan('matlab','9.1')
-				y = exp( - bsxfun(@times, theta.k , bsxfun(@power, x, theta.tau)) );
+				y = exp( - bsxfun(@times, theta.k , bsxfun(@power, x, tau)) );
 			else
 				% use new array broadcasting in 2016b
-				y = exp( - theta.k .* x.^theta.tau );
+				y = exp( - theta.k .* x.^tau );
 			end
 		end
 

ddToolbox/models/parametric_models/exponential_power_models/ExponentialPower.m

@@ -6,15 +6,15 @@
 			obj = obj@SubjectiveTimeModel(data, varargin{:});
 
             obj.dfClass = @DF_ExponentialPower;
-            obj.subjectiveTimeFunctionFH = @SubjectiveTimePowerFunctionEP;
+            obj.subjectiveTimeFunctionFH = @SubjectiveTimePowerFunction;
 
 			% Create variables
-			obj.varList.participantLevel = {'k','tau','alpha','epsilon'};
-			obj.varList.monitored = {'k','tau','alpha','epsilon', 'Rpostpred', 'P', 'VA', 'VB'};
+			obj.varList.participantLevel = {'k','logtau','alpha','epsilon'};
+			obj.varList.monitored = {'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 = 'tau';
-            obj.varList.discountFunctionParams(2).label = 'tau';
+            obj.varList.discountFunctionParams(2).name = 'logtau';
+            obj.varList.discountFunctionParams(2).label = '$\log(\tau)$';
 
             obj.plotOptions.dataPlotType = '2D';
 		end

ddToolbox/models/parametric_models/exponential_power_models/ModelSeparateExpPower.m

@@ -20,7 +20,7 @@
 			nExperimentFiles = obj.data.getNExperimentFiles();
 			for chain = 1:nchains
 				initialParams(chain).k = unifrnd(-1, 1, [nExperimentFiles,1]);
-                initialParams(chain).tau = unifrnd(0.01, 1, [nExperimentFiles,1]);
+                initialParams(chain).logtau = normrnd(0, 5, [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/exponential_power_models/SubjectiveTimePowerFunction.m

@@ -39,13 +39,14 @@ function plot(obj, pointEstimateType)
 	methods (Static, Access = protected)
 
 		function y = function_evaluation(x, theta)
+			tau = exp(theta.logtau);
 			if verLessThan('matlab','9.1')
                 y = bsxfun(@times, ...
-                        bsxfun(@power, x, theta.tau),...
+                        bsxfun(@power, x, tau),...
                         theta.k);
 			else
 				% use new array broadcasting in 2016b
-				y = theta.k .* (x .^ theta.tau);
+				y = theta.k .* (x .^ tau);
 			end
 		end
 

ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.jags

@@ -1,4 +1,4 @@
-# RANDOM FACTORS:   k[p], tau[p], epsilon[p], alpha[p]
+# RANDOM FACTORS:   k[p], logtau[p], epsilon[p], alpha[p]
 # HYPER-PRIORS ON:  epsilon[p], alpha[p]
 
 # RANDOM (by participant) FACTORS APPEAR IN A LOOP OVER PARTICIPANTS
@@ -10,20 +10,20 @@ model{
 # DISCOUNT FUNCTION PARAMETERS =================================================
 # RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
 
-K_MEAN          <- 0                     # <---- currently guesstimating
-K_PRECISION     <- 1/(0.01)              # <---- currently guesstimating
-TAU_MEAN        <- 1                     # <---- currently guesstimating
-TAU_PRECISION   <- 1/1^2                 # <---- currently guesstimating
+K_MEAN            <- 0                     # <---- currently guesstimating
+K_PRECISION       <- 1/(0.01)              # <---- currently guesstimating
+LOG_TAU_MEAN      <- 0                     # <---- currently guesstimating
+LOG_TAU_PRECISION <- 1/2^2                 # <---- currently guesstimating
 
 for (p in 1:nRealExperimentFiles){       # no +1 because no shrinkage hyperprior
 	k[p]        ~ dnorm(K_MEAN, K_PRECISION)
-    tau[p]      ~ dnorm(TAU_MEAN, TAU_PRECISION) T(0,)
+    logtau[p]   ~ dnorm(LOG_TAU_MEAN, LOG_TAU_PRECISION)
 }
 
 # MODEL-SPECIFIC: CALCULATION OF PRESENT SUBJECTIVE VALUES
 for (t in 1:length(ID)) {
-    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]]) ) )
+    VA[t] <- A[t] * (exp( -k[ID[t]] * (DA[t]^exp(logtau[ID[t]])) ) )
+	VB[t] <- B[t] * (exp( -k[ID[t]] * (DB[t]^exp(logtau[ID[t]])) ) )
 }
 
 

ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.jags

@@ -1,4 +1,4 @@
-# RANDOM FACTORS:   k[p], tau[p], epsilon[p], alpha[p]
+# RANDOM FACTORS:   k[p], logtau[p], epsilon[p], alpha[p]
 # HYPER-PRIORS ON:  epsilon[p], alpha[p]
 
 # RANDOM (by participant) FACTORS APPEAR IN A LOOP OVER PARTICIPANTS
@@ -11,18 +11,18 @@ model{
 # RANDOM (BY PARTICIPANT) FACTORS; HYPER-PRIORS = NO 
 
 K_MEAN          <- 0                     # <---- currently guesstimating
-K_PRECISION     ~ dgamma(0.01, 0.01)     #sigma=0.1         # <---- currently guesstimating
+K_PRECISION     ~ dgamma(1, 1)     #sigma=0.1         # <---- currently guesstimating
 
 for (p in 1:nRealExperimentFiles){       # no +1 because no shrinkage hyperprior
 	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
+    logtau[p]   ~ dnorm(0, 1/0.5^2)
 }
 
 # 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]]) ) )
+    VA[t] <- A[t] * (exp( -k[ID[t]] * (DA[t]^exp(logtau[ID[t]])) ) )
+	VB[t] <- B[t] * (exp( -k[ID[t]] * (DB[t]^exp(logtau[ID[t]])) ) )
 }