Comparing different models of rate-choose computation

FlemingDaw_PsyRev_figure567.m

@@ -49,7 +49,7 @@
 
     sigma_conf = sigma_conf_vec(sigmai);
 
-    % Pre-action confidence is equivalent to max of first-order conf
+    % Compute first-order model for comparison
     mean_cor_pre(sigmai) = max([computeFirstOrderConf(xp, -1, sigma_act), computeFirstOrderConf(xp, 1, sigma_act)]);
     pD_pre(sigmai) = computeFirstOrderConf(xp, 1, sigma_act);
 
@@ -109,7 +109,7 @@
 
     rho = rho_vec(sigmai);
 
-    % Pre-action confidence is max of first-order conf
+    % Simulate first-order confidence for comparison
     mean_cor_pre(sigmai) = max([computeFirstOrderConf(xp, -1, sigma_act), computeFirstOrderConf(xp, 1, sigma_act)]);
     pD_pre(sigmai) = computeFirstOrderConf(xp, 1, sigma_act);
 

model development/compare_rateChoose_max_jags.m

@@ -0,0 +1,111 @@
+% Check correspondence between rate-choose max and sampling
+
+clear all
+close all
+
+addpath(genpath('~/Dropbox/Utils/matjags'));
+
+% Set d
+d = 1;
+% Evaluate comparision for each parameter
+sigma = 1;
+sigma_a = 1.5;
+sigma_p = 1;
+xp_space = linspace(-2,2,20);
+
+for xpi = 1:length(xp_space)
+
+    % JAGS
+    jags_conf(xpi) = sampleMetaConf_singleTrial_noAction(xp_space(xpi), sigma, sigma_a, sigma_p);
+
+    % Analytic
+    ana_conf(xpi) = max([computeMetaConf_old(xp_space(xpi), -1, sigma, sigma_a, sigma_p), computeMetaConf_old(xp_space(xpi), 1, sigma, sigma_a, sigma_p)]);
+
+end
+
+h = figure;
+set(gcf, 'Position', [200 200 800 300]);
+subplot(1,3,1)
+plot(xp_space, jags_conf', 'LineWidth', 2);
+xlabel('xp', 'FontSize', 14);
+ylabel('confidence', 'FontSize', 14);
+axis square
+set(gca, 'FontSize', 12);
+
+subplot(1,3,2)
+plot(xp_space, ana_conf', 'LineWidth', 2);
+xlabel('xp', 'FontSize', 14);
+ylabel('confidence', 'FontSize', 14);
+axis square
+set(gca, 'FontSize', 12);
+
+subplot(1,3,3)
+plot(jags_conf(:), ana_conf(:), 'bx ', 'MarkerSize', 7, 'LineWidth', 1.5);
+line([0 1], [0 1], 'LineStyle', '--', 'Color', 'k');
+xlabel('Sampling', 'FontSize', 14);
+ylabel('Max', 'FontSize', 14);
+axis square
+set(gca, 'FontSize', 12);
+
+%% Simulate rate-choose experimental predictions under alternative definition without conditioning on action
+theta = [0 0.032 0.064 0.128 0.256 0.512 1];
+sigma = 1;
+sigma_act = 1;
+sigma_conf = 1;
+N = 1000;
+
+% Simulate x's, decisions and confidence for first order + second order
+meta_cor_rateChoose = nan(length(theta), N);
+meta_cor_chooseRate = nan(length(theta), N);
+d = nan(length(theta), N);
+a = nan(length(theta), N);
+
+for s = 1:length(theta)
+
+    for i = 1:N
+
+        if rand > 0.5
+            d(s,i) = 1;
+        else
+            d(s,i) = -1;
+        end
+
+        x = normrnd(d(s,i)*theta(s), sigma);
+        xa = normrnd(x, sigma_act);
+        xp = normrnd(x, sigma_conf);
+
+        if xa > 0
+            a(s,i) = 1;
+        else
+            a(s,i) = -1;
+        end
+
+        % Compute rate-choose confidence, second-order
+        meta_cor_rateChoose(s,i) = sampleMetaConf_singleTrial_noAction(xp, sigma, sigma_act, sigma_conf);
+
+        % Compute choose-rate confidence, second-order
+        meta_cor_chooseRate(s,i) = computeMetaConf_old(xp, a(s,i), sigma, sigma_act, sigma_conf);
+
+    end
+end
+
+for s = 1:length(theta)
+    meta_conf_rateChoose_err(s) = mean(meta_cor_rateChoose(s, d(s,:) ~= a(s,:)));
+    meta_conf_rateChoose_cor(s) = mean(meta_cor_rateChoose(s, d(s,:) == a(s,:)));
+    meta_conf_chooseRate_err(s) = mean(meta_cor_chooseRate(s, d(s,:) ~= a(s,:)));
+    meta_conf_chooseRate_cor(s) = mean(meta_cor_chooseRate(s, d(s,:) == a(s,:)));
+end
+
+h4a = figure;
+set(gcf, 'Position', [200 200 400 400])
+plot(theta, meta_conf_chooseRate_cor, 'g', 'LineWidth', 2);
+hold on
+plot(theta, meta_conf_chooseRate_err, 'r', 'LineWidth', 2);
+plot(theta, meta_conf_rateChoose_cor, 'g--', 'LineWidth', 2);
+plot(theta, meta_conf_rateChoose_err, 'r--', 'LineWidth', 2);
+xlabel('Stimulus strength (\theta)', 'FontSize', 20);
+ylabel('Confidence','FontSize', 20);
+set(gca, 'YLim', [0.5 1], 'FontSize', 20);
+legend({'Choose-rate cor', 'Choose-rate err', 'Rate-choose cor', 'Rate-choose err'}, 'Location', 'SouthWest');
+legend boxoff
+box off

model development/sampleMetaConf_singleTrial_noAction.m

@@ -0,0 +1,65 @@
+function conf = sampleMetaConf_singleTrial_noAction(xp, sigma, sigma_a, sigma_p)
+% function conf = selfSelf_singleTrial_noAction(xp, sigma, sigma_a, sigma_p)
+% 
+% Runs single-trial inference (via JAGS) of meta-inference model
+% xp - perceptual sample
+% aciton - assumed that we don't condition on action here
+%
+% sigma = shared noise (SD)
+% sigma_a = action noise (SD)
+% sigma_p = perception noise (SD)
+% 
+% SF 2014
+
+if nargin < 1
+    xp = 0.5;
+    sigma = 0.7;
+    sigma_a = 1;
+    sigma_p = 1;
+elseif nargin > 0 & nargin < 4
+    error('Data and precisions not specified')
+end
+
+%% Sampling
+% MCMC Parameters
+mcmc_params.nchains = 1; % How Many Chains?
+mcmc_params.nburnin = 10; % How Many Burn-in Samples?
+mcmc_params.nsamples = 1000;  %How Many Recorded Samples?
+mcmc_params.nthin = 1; % How Often is a Sample Recorded?
+mcmc_params.doparallel = 0; % Parallel Option
+mcmc_params.dic = 1;
+% Initialize Unobserved Variables
+for i=1:mcmc_params.nchains
+    S.x = 0;
+    mcmc_params.init0(i) = S;
+end
+
+% Assign variables to the observed nodes
+datastruct = struct('x2', xp, 'sigma', sigma, 'sigma1', sigma_a, 'sigma2', sigma_p);
+
+% Select model file and parameters to monitor
+model_file = 'sampleMetaConf.txt';
+monitorparams = {'d','a'};
+
+% Use JAGS to Sample
+[samples, stats] = matjags( ...
+    datastruct, ...
+    fullfile(pwd, model_file), ...
+    mcmc_params.init0, ...
+    'doparallel' , mcmc_params.doparallel, ...
+    'nchains', mcmc_params.nchains,...
+    'nburnin', mcmc_params.nburnin,...
+    'nsamples', mcmc_params.nsamples, ...
+    'thin', mcmc_params.nthin, ...
+    'dic', mcmc_params.dic,...
+    'monitorparams', monitorparams, ...
+    'savejagsoutput' , 0 , ...
+    'verbosity' , 0 , ...
+    'cleanup' , 1 , ...
+    'workingdir' , 'tmpjags' );
+act = samples.a(:);
+act(act == 0) = -1;
+
+conf = sum(samples.d(:) == act(:))./length(samples.d(:));
+
+