Distilling noise characteristics and prior expectations in multisensory causal inference.

Code ↔ Paper

18 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 18 matches

[1] § Methods › Model fitting ↔ utils/cmaes.m, lines 1–60 · score 0.74 · CMA ES, Covariance Matrix, Evolution Strategy, Adaptation, minimization
[2] § Methods › Model fitting ↔ utils/cmaes_modded.m, lines 1–60 · score 0.74 · CMA ES, Covariance Matrix, Evolution Strategy, Adaptation, minimization
[3] § Methods › Models and free parameters › Parametric models. ↔ manuscript_allplots.m, lines 238–249 · score 0.70 · TwoGaussians, Const GaussianLaplace, SingleGaussian, Exp GaussianLaplace, parametric models, causal
[4] § Methods › Models and free parameters › Parametric models. ↔ manuscript_allplots_old.m, lines 236–252 · score 0.70 · TwoGaussians, Const GaussianLaplace, SingleGaussian, Exp GaussianLaplace, parametric models, causal
[5] § Methods › Model recovery analysis ↔ manuscript_allplots.m, lines 373–496 · score 0.68 · ground truth generative, model recovery, fitted model, Rows, matrix, NLL
[6] § Results › Distilled parametric shapes for priors and noise › All-tasks fits. ↔ manuscript_allplots.m, lines 238–249 · score 0.65 · TwoGaussians, Const GaussianLaplace, SingleGaussian, Exp GaussianLaplace, Model comparison, parametric model
[7] § Results › Distilled parametric shapes for priors and noise › All-tasks fits. ↔ manuscript_allplots_old.m, lines 236–252 · score 0.64 · TwoGaussians, Const GaussianLaplace, SingleGaussian, Exp GaussianLaplace, Model comparison, parametric model
[8] § Methods › Models and free parameters › Parametric model names. ↔ manuscript_allplots_old.m, lines 941–1038 · score 0.64 · Gaussian Laplace prior, Exponential sensory noise, auditory stimulus
[9] § Results › Distilled parametric shapes for priors and noise › Unisensory data fits. ↔ manuscript_allplots.m, lines 1206–1303 · score 0.63 · exponential sensory noise, TwoGaussians, SingleGaussian, GaussianLaplace
[10] § Results › Distilled parametric shapes for priors and noise › Unisensory data fits. ↔ manuscript_allplots_old.m, lines 941–1038 · score 0.63 · exponential sensory noise, TwoGaussians, SingleGaussian, GaussianLaplace
[11] § Methods › Model recovery analysis ↔ manuscript_allplots.m, lines 373–496 · score 0.62 · ground truth generative, Model recovery, AIC, BIC
[12] § Results › Models ↔ analysis/nllfun_bc_parametric.m, lines 167–211 · score 0.61 · causal inference strategy, BC task, Model Selection, motor noise, posterior, Matching
[13] § Methods › Models and free parameters › Parametric model names. ↔ manuscript_allplots.m, lines 1206–1303 · score 0.55 · Gaussian Laplace prior, Exponential sensory noise, stimulus
[14] § Results › Distilled parametric shapes for priors and noise › All-tasks fits. ↔ manuscript_allplots_old.m, lines 607–696 · score 0.55 · lifted semiparametric fits, parametric model fit, BC, Exp, GaussianLaplace
[15] § Results › Lifting the semiparametric fits to all tasks ↔ manuscript_allplots_old.m, lines 607–696 · score 0.55 · best lifted semiparametric, model fits, AIC, BIC, causal
[16] § Results › Distilled parametric shapes for priors and noise › All-tasks fits. ↔ manuscript_allplots.m, lines 866–961 · score 0.54 · lifted SemiParametric fits, lifted SemiParametric model, model parameters, AIC, BIC, Exp
[17] § Methods › Models and free parameters › Parametric model fits on only unisensory data. ↔ create_composite_allfits_figures.ipynb, lines 133–193 · score 0.51 · unisensory parametric, Exp GaussianLaplace, sensory noise, Const, fitted, models
[18] § Results › Lifting the semiparametric fits to all tasks ↔ analysis/manuscript_bimodalavfits_visualization_resc.m, lines 176–313 · score 0.50 · stimulus location disparity, visual reliabilities, bias, ribbons, weight, stratified

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

MATLAB · 1,561 lines · 69 KB · no license · 7 matches

clear all; close all;
cd('C:\Users\liu_s\Audiovisual-causal-inference')
fig_maxwidth_inches = 7.5;
fig_maxheight_inches = 8.75;
set(0,'units','inches');
Inch_SS = get(0,'screensize');
set(0,'units','pixels');
figsize = get(0, 'ScreenSize');
Res = figsize(3)./Inch_SS(3);
set(groot,'DefaultAxesFontName','Arial')
%figsize_RespDistr = [0,0,figsize(4)*4/3, figsize(4)];
figsize_RespDistr = [0,0,fig_maxwidth_inches, fig_maxheight_inches] .* Res;
figformat = "svg";
figpath = "plots\"; %"newplots\"
fontsize=9; %9
png_dpi = 500;
plot_lapse = true;
lapse_type = "Uniform";
model_path = "modelfits\";
data_path = "data\";
analysis_path = "analysis\";
addpath(analysis_path,data_path,model_path,"utils\");
%% s_V, s_A generative distributions
linewidth = 1;
figure('Position', [0,0,figsize_RespDistr(3)./5,figsize_RespDistr(3)./8]);
tiledlayout(1,1,'TileSpacing','none', 'Padding','none'); set(gca,'TickDir','out'); hold on;
sV_vals = [-25:1:-20,-20:1:20,20:1:25]; p_sV_vals = [zeros(size(-25:1:-20)),repmat(1/40,1,length(-20:1:20)),zeros(size(20:1:25))];
plot(sV_vals, p_sV_vals, "k-", 'LineWidth',linewidth); area(sV_vals, p_sV_vals,'FaceColor','k', 'FaceAlpha',0.2);
xlim([-25,25]); ylim([0,0.03]); xticks(-20:20:20); yticks([0,1/40]); yticklabels(["0","1/40"])
xlabel("{\its}_V",'FontSize',fontsize+1); ylabel("p({\its}_V)",'FontSize',fontsize+1);
set(gca,'FontSize',fontsize)
exportgraphics(gcf,figpath+'sV_gen'+'.pdf',"ContentType","vector");
figure('Position', [0,0,figsize_RespDistr(3)./5,figsize_RespDistr(3)./8]);
tiledlayout(1,1,'TileSpacing','none', 'Padding','tight'); set(gca,'TickDir','out'); hold on;
sA_vals = [-15:5:15]; p_sA_vals = repmat(1/length(sA_vals), 1,length(sA_vals));
h=stem(sA_vals, p_sA_vals,"k-", 'LineWidth',linewidth); set(h, 'Marker', 'none')
xlim([-20,20]); ylim([0,0.16]); xticks(-15:15:15); xtickangle(0); yticks([0,1/length(sA_vals)]); yticklabels(["0","1/7"])
xlabel("{\its}_A",'FontSize',fontsize); ylabel("p({\its}_A)",'FontSize',fontsize);
set(gca,'FontSize',fontsize)
exportgraphics(gcf,figpath+'sA_gen'+'.pdf',"ContentType","vector");
%% UAV data visualized, without model prediction ribbons.
prior = "NaN";
noise = "NaN";
aud_rescale = "NaN";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type);
exportgraphics(gcf,figpath+'UAV_dataonly'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'UAV_dataonly'+'.pdf',"ContentType","vector");
%% UJoint parametric model response distribution visualization
prior = "SingleGaussian";
noise = "constant";
aud_rescale = "1";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type, false, true);
exportgraphics(gcf,figpath+'Const-SingleGaussian_rescaleaud1'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Const-SingleGaussian_rescaleaud1'+'.pdf',"ContentType","vector");
%%
prior = "SingleGaussian";
noise = "constant";
aud_rescale = "free";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type);
exportgraphics(gcf,figpath+'Const-SingleGaussian'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Const-SingleGaussian'+'.pdf',"ContentType","vector");
%%
prior = "GaussianLaplaceBothFixedZero";
noise = "exp";
aud_rescale = "free";
%manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, four_by_three_figsize);
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type);
exportgraphics(gcf,figpath+'Exp-GaussianLaplace'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Exp-GaussianLaplace'+'.pdf',"ContentType","vector");
%% Individual-level
prior = "GaussianLaplaceBothFixedZero";
noise = "exp";
aud_rescale = "free";
save_name = "Exp-GaussianLaplace";
plot_individual = true;
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type, plot_individual);
figure(1)
saveas(gca, figpath+save_name+'_Individualmean.fig')
exportgraphics(gcf,figpath+save_name+'_Individualmean'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+save_name + '_Individualmean'+'.pdf',"ContentType","vector");
figure(2)
saveas(gca, figpath+save_name+'_IndividualSD.fig')
exportgraphics(gcf,figpath+save_name+'_IndividualSD'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+save_name + '_IndividualSD'+'.pdf',"ContentType","vector");
%% Exemplary subject
subjidx=7;
fitted_on_all_data = false;
allindvsubjplots_to_onesubjplot(save_name, subjidx, fitted_on_all_data, 10, [0 0 figsize_RespDistr(3) figsize_RespDistr(3)*0.5], figpath)
exportgraphics(gcf,figpath+'Exp-GaussianLaplace_Individual_example'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Exp-GaussianLaplace_Individual_example'+'.pdf',"ContentType","vector");
%% No Lapse
prior = "GaussianLaplaceBothFixedZero";
noise = "exp";
aud_rescale = "free";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, false, lapse_type);
exportgraphics(gcf,figpath+'Exp-GaussianLaplace_nolapse'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Exp-GaussianLaplace_nolapse'+'.pdf',"ContentType","vector");
%% Gaussian
prior = "GaussianLaplaceBothFixedZero";
noise = "exp";
aud_rescale = "free";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, true, "Gaussian");
exportgraphics(gcf,figpath+'Exp-GaussianLaplace_Gaussianlapse'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Exp-GaussianLaplace_Gaussianlapse'+'.pdf',"ContentType","vector");
% Compute model comparison between uniform lapse model.
UnimodalData_ModelComparison_FinalTables_uniformgaussianlapse = unimodaldata_modelcomparison_visualize_uniformgaussianlapse(model_path, data_path);
%%
prior = "SingleGaussian";
noise = "exp";
aud_rescale = "free";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type);
exportgraphics(gcf,figpath+'Exp-SingleGaussian'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Exp-SingleGaussian'+'.pdf',"ContentType","vector");
prior = "GaussianLaplaceBothFixedZero";
noise = "constant";
aud_rescale = "free";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type);
exportgraphics(gcf,figpath+'Const-GaussianLaplace'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Const-GaussianLaplace'+'.pdf',"ContentType","vector");
prior = "TwoGaussiansBothFixedZero";
noise = "exp";
aud_rescale = "free";
manuscript_ujoint_respdistrvisualization(prior, noise, aud_rescale, fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type);
exportgraphics(gcf,figpath+'Exp-TwoGaussians'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Exp-TwoGaussians'+'.pdf',"ContentType","vector");
%% Unimodal semiparam model
% Response distribution visualization
% manuscript_ujoint_respdistrvisualization_semiparam(fontsize, four_by_three_figsize);
manuscript_ujoint_respdistrvisualization_semiparam(fontsize, figsize_RespDistr, model_path, plot_lapse, lapse_type);
exportgraphics(gcf,figpath+'Semiparam_FittedRespDistr'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Semiparam_FittedRespDistr'+'.pdf',"ContentType","vector");
%%
% sigma(s), p(s) visualization
semiparam_sigmafun_prior_visualization(fontsize+1, figsize_RespDistr, model_path);
exportgraphics(gcf,figpath+'Semiparam_FittedParams'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'Semiparam_FittedParams'+'.pdf',"ContentType","vector");
%% Unimodal data ModelComparison
priors = ["","GaussianLaplaceBothFixedZero","GaussianLaplaceBothFixedZero","GaussianLaplaceBothFixedZero", "SingleGaussian", "GaussianLaplaceBothFixedZero", "TwoGaussiansBothFixedZero","SingleGaussian","SingleGaussian","SingleGaussian"];
noises = ["","exp", "exp", "exp", "exp", "constant", "exp", "constant","constant","constant"];
rescales = ["","free", "4over3", "1", "free","free","free","free","4over3","1"];
model_types = ["semiparametric","exp-GaussianLaplace", "exp-GaussianLaplace\_4/3","exp-GaussianLaplace\_1", "exp-SingleGaussian", "const-GaussianLaplace","exp-TwoGaussians","const-SingleGaussian","const-SingleGaussian\_4/3","const-SingleGaussian\_1"];
num_params = [40,14,13,13,12,10,14,8,7,7];
%% Vanila 3 models on unimodal data only
figure('Position', [0 0 5.2*Res figsize_RespDistr(4)*0.4]);
set(gcf, 'Color', 'w')
UnimodalData_ModelComparison_FinalTables_Vanilla = unimodaldata_modelcomparison_visualize(priors((end-2):end), noises((end-2):end), rescales((end-2):end), model_types((end-2):end), num_params((end-2):end), true, fontsize+1, model_path, data_path, false);
exportgraphics(gcf,figpath+'UJoint_ModelSelection_vanilla'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'UJoint_ModelSelection_vanilla'+'.pdf',"ContentType","vector");
save(analysis_path+'unimodaldata_modelcomparison_finaltables_vanilla','UnimodalData_ModelComparison_FinalTables_Vanilla');
% All models on unimodal data
keep_modelidx = [2,5,6,7,8,1];
figure('Position', [0 0 5.2*Res figsize_RespDistr(3)*0.6*0.5]);
set(gcf, 'Color', 'w')
UnimodalData_ModelComparison_FinalTables = unimodaldata_modelcomparison_visualize(priors(keep_modelidx), noises(keep_modelidx), rescales(keep_modelidx), model_types(keep_modelidx), num_params(keep_modelidx), true, fontsize+1, model_path, data_path, true);
exportgraphics(gcf,figpath+'UJoint_ModelSelection_BIC'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'UJoint_ModelSelection_BIC'+'.pdf',"ContentType","vector");
keep_modelidx = [2,3,4,5,6,7,8,9,10,1];
figure('Position', [0 0 5.2*Res figsize_RespDistr(3)*0.6*0.7]);
set(gcf, 'Color', 'w')
UnimodalData_ModelComparison_FinalTables = unimodaldata_modelcomparison_visualize(priors(keep_modelidx), noises(keep_modelidx), rescales(keep_modelidx), model_types(keep_modelidx), num_params(keep_modelidx), true, fontsize+1, model_path, data_path, true);
exportgraphics(gcf,figpath+'UJoint_ModelSelection_BIC_full'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'UJoint_ModelSelection_BIC_full'+'.pdf',"ContentType","vector");
figure('Position', [0 0 figsize_RespDistr(3) figsize_RespDistr(3)]);
set(gcf, 'Color', 'w')
UnimodalData_ModelComparison_FinalTables = unimodaldata_modelcomparison_visualize(priors, noises, rescales, model_types, num_params, true, fontsize+1, model_path, data_path, false);
exportgraphics(gcf,figpath+'UJoint_ModelSelection'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'UJoint_ModelSelection'+'.pdf',"ContentType","vector");
save(analysis_path+'unimodaldata_modelcomparison_finaltables','UnimodalData_ModelComparison_FinalTables');
%% AllData LiftedSemiparam Response distribution visualization
causal_inf_strategy = "ProbMatching";
save_name = "PM";
manuscript_allfits_respdistrvisual_semiparaminsp_maintext(causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type);
manuscript_allfits_respdistrvisualization_semiparaminsp_resc(causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type);
%%
causal_inf_strategy = "ModelSelection";
save_name = "MS";
manuscript_allfits_respdistrvisualization_semiparaminsp_resc(causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type);
causal_inf_strategy = "ModelAveraging";
save_name = "MA";
manuscript_allfits_respdistrvisualization_semiparaminsp_resc(causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type);
%% AllData parametric model response distributions
prior_type = "GaussianLaplaceBothFixedZero";
hetero_type = "exp";
causal_inf_strategy = "ProbMatching";
save_name = "exp-GaussianLaplace-PM";
manuscript_allfits_respdistrvisualization_resc_maintext(prior_type, hetero_type, causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type)
manuscript_allfits_respdistrvisualization_resc(prior_type, hetero_type, causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type)
%% Individual-level plots for the above model
prior_type = "GaussianLaplaceBothFixedZero";
hetero_type = "exp";
causal_inf_strategy = "ProbMatching";
save_name = "exp-GaussianLaplace-PM";
manuscript_allfits_respdistrvisualization_resc(prior_type, hetero_type, causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type, true)
%% Exemplary subject
subjidx=7;
fitted_on_all_data = true;
allindvsubjplots_to_onesubjplot(save_name,subjidx, fitted_on_all_data, fontsize, figsize_RespDistr, figpath)
exportgraphics(gcf,figpath+'Exp-GaussianLaplace-PM_Individual_example'+'.pdf',"ContentType","vector");
exportgraphics(gcf,figpath+'Exp-GaussianLaplace-PM_Individual_example'+'.png','Resolution',png_dpi);
%%
prior_type = "GaussianLaplaceBothFixedZero";
hetero_type = "exp";
causal_inf_strategy = "ModelAveraging";
save_name = "exp-GaussianLaplace-MA";
manuscript_allfits_respdistrvisualization_resc(prior_type, hetero_type, causal_inf_strategy, fontsize, figsize_RespDistr, figpath, save_name, png_dpi, model_path, plot_lapse, lapse_type)
%% AllData ModelComparison
causal_inf_strategies = ["ModelSelection","ModelAveraging","ProbMatching"];
param_model_names = ["exp-GaussianLaplace","exp-SingleGaussian","const-GaussianLaplace","exp-TwoGaussians","const-SingleGaussian","paramBest", "LiftedSemiparam"];
% Only BIC for PM models
figure('Position', [0 0 5.2*Res figsize_RespDistr(3)*0.3]);
hold on;
set(gcf, 'Color', 'w')
alldata_modelcomparison_visualize([causal_inf_strategies(3)], param_model_names, true, fontsize, model_path, data_path, true, true, true);
alldata_modelcomparison_visualize([causal_inf_strategies(3)], param_model_names, true, fontsize, model_path, data_path, true, false, false);
exportgraphics(gcf,figpath+'SemiparamIndv_ModelSelection_BIC'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'SemiparamIndv_ModelSelection_BIC'+'.pdf',"ContentType","vector");
%% Only BIC for all models
figure('Position', [0 0 5.2*Res figsize_RespDistr(3)*0.6]);
set(gcf, 'Color', 'w')
alldata_modelcomparison_visualize(causal_inf_strategies, param_model_names, true, fontsize, model_path, data_path, true, true, true);
alldata_modelcomparison_visualize(causal_inf_strategies, param_model_names, true, fontsize, model_path, data_path, true, false, false);
exportgraphics(gcf,figpath+'SemiparamIndv_ModelSelection_BIC_full'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'SemiparamIndv_ModelSelection_BIC_full'+'.pdf',"ContentType","vector");
% NLL, AIC, BIC for all models
figure('Position', [0 0 figsize_RespDistr(3) figsize_RespDistr(3)]);
set(gcf, 'Color', 'w')
AllData_ModelComparison_FinalTables_liftedsemiparam_count = alldata_modelcomparison_visualize(causal_inf_strategies, param_model_names, true, fontsize, model_path, data_path, false, true, true);
AllData_ModelComparison_FinalTables = alldata_modelcomparison_visualize(causal_inf_strategies, param_model_names, true, fontsize, model_path, data_path, false, false, false);
writematrix(round(AllData_ModelComparison_FinalTables{1},2), figpath+"NLL.csv")
writematrix(round(AllData_ModelComparison_FinalTables{2},2), figpath+"AIC.csv")
writematrix(round(AllData_ModelComparison_FinalTables{3},2), figpath+"BIC.csv")
exportgraphics(gcf,figpath+'SemiparamIndv_ModelSelection'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'SemiparamIndv_ModelSelection'+'.pdf',"ContentType","vector");
save(analysis_path+'alldata_modelcomparison_finaltables','AllData_ModelComparison_FinalTables');
%% sigma(s) and p(s) examples
figure('Position', [0 0 figsize_RespDistr(3) figsize_RespDistr(3)]);
set(gcf, 'Color', 'w')
sigmafun_prior_examples(fontsize);
exportgraphics(gcf,figpath+'SensoryNoisePriorParamFamilies'+'.png','Resolution',png_dpi);
exportgraphics(gcf,figpath+'SensoryNoisePriorParamFamilies'+'.pdf',"ContentType","vector");
%% Parameter recovery
load(model_path+"fittedparams_UJoint_exp-GaussianLaplaceBothFixedZero_rescalefree_lapseUniform.mat")
theta_fitted_orig = theta_fitted;
F_vals_orig = F_vals;
load(model_path+"fittedparams_UJoint_exp-GaussianLaplaceBothFixedZero_rescalefree_lapseUniform__exp-GaussianLaplaceBothFixedZero_rescalefree_lapseUniform")
num_params_model = length(theta_fitted(1,:));
param_names = ["\sigma_{0,V}","k_{1,V}","k_{2,V}","\alpha_{med}","\alpha_{low}","\sigma_s","\lambda","b","w","\sigma_{motor}","\sigma_{0,A}","k_{1,A}","k_{2,A}", "\rho_A"];
figure('Position',[50,100,700,500]);
t=tiledlayout(3,ceil(num_params_model/3),'Padding', 'tight', 'TileSpacing', 'tight');
for param =1:num_params_model
nexttile(t); hold on;
theta_aug = [theta_fitted_orig(:,param); theta_fitted(:,param)];
scatter(theta_fitted_orig(:,param), theta_fitted(:,param), "k.")
plot([min(0,min(theta_aug)), max(theta_aug)], [min(0,min(theta_aug)), max(theta_aug)], "g-")
title("$"+param_names(param)+"$",'interpreter','latex','fontsize',12)
end
xlabel(t,'Ground-truth parameter value','FontSize',12)
ylabel(t,'Recovered parameter value','FontSize',12)
exportgraphics(gcf,figpath+"UAV_Exp-GaussianLaplace_ParamRecovery"+".pdf","ContentType","vector");
%% Check model recovery
prior_types = ["SingleGaussian","GaussianLaplaceBothFixedZero","SingleGaussian","SingleGaussian","GaussianLaplaceBothFixedZero","GaussianLaplaceBothFixedZero"]; % "SingleGaussian", "GaussianLaplaceBothFixedZero", or "TwoGaussiansBothFixedZero"
hetero_types = ["constant","exp","constant","exp","constant","exp"]; % "constant" or "exp";
lapse_types = repmat("Uniform",1,length(prior_types)); % "Uniform" or "Gaussian";
rescale_auds = ["free","free","1","free","free","1"]; % "1", "4/3", or "free";
prior_type_modelrecovdatas = prior_types;
hetero_type_modelrecovdatas = hetero_types; % "constant" or "exp";
lapse_type_modelrecovdatas = lapse_types; % "Uniform" or "Gaussian";
rescale_aud_modelrecovdatas = rescale_auds; % "1", "4/3", or "free";
% Model parameters
num_params = [8,14,7,12,10,13];
num_datasets = length(prior_type_modelrecovdatas);
num_models = length(prior_types);
num_subjects = 15;
load(data_path+"data_stratified_UV.mat");
load(data_path+"data_stratified_UA.mat");
data_UV = data_stratified_to_data(data_stratified_UV, false, true); % last argument is is_visual.
data_UA = data_stratified_to_data(data_stratified_UA, false, false);
n_data = zeros(1,15);
for subjidx=1:num_subjects
n_data(subjidx) = length(data_UA{subjidx}) +length(data_UV{subjidx});
end
% AIC, BIC
NLLs = zeros(num_datasets, num_models, num_subjects);
AICs = zeros(num_datasets, num_models, num_subjects);
BICs = zeros(num_datasets, num_models, num_subjects);
% load nonparam indv UJoint fits
for model_idx =1:6
prior_type = prior_types(model_idx);
hetero_type = hetero_types(model_idx);
lapse_type = lapse_types(model_idx);
rescale_aud = rescale_auds(model_idx);
for dataset_idx =1:6
prior_type_modelrecovdata = prior_type_modelrecovdatas(dataset_idx); % "SingleGaussian", "GaussianLaplaceBothFixedZero", or "TwoGaussiansBothFixedZero"
hetero_type_modelrecovdata = hetero_type_modelrecovdatas(dataset_idx); % "constant" or "exp";
lapse_type_modelrecovdata = lapse_type_modelrecovdatas(dataset_idx); % "Uniform" or "Gaussian";
rescale_aud_modelrecovdata = rescale_aud_modelrecovdatas(dataset_idx); % "1", "4/3", or "free";
filename = 'fittedparams_UJoint_'+hetero_type+"-"+prior_type+"_rescale"+rescale_aud+"_lapse"+lapse_type;
datafilename = hetero_type_modelrecovdata+"-"+prior_type_modelrecovdata+"_rescale"+rescale_aud_modelrecovdata+"_lapse"+lapse_type_modelrecovdata;
filename_final = filename + "__" + datafilename;
for model=1:(num_models)
load(model_path + filename_final+".mat")
[min_val, min_idx] = min(F_vals,[],2);
NLLs(dataset_idx, model_idx,:) = min_val';
AICs(dataset_idx, model_idx,:) = 2.*min_val' + 2.* num_params(model_idx);
BICs(dataset_idx, model_idx,:) = 2.*min_val' + num_params(model_idx).*log(n_data);
end
end
end
NLL_sum = sum(NLLs,3);
AIC_sum = sum(AICs,3);
BIC_sum = sum(BICs,3);
% figure; hold on;
NLL_sumvalues_diff = NLL_sum - diag(NLL_sum);
[~,min_NLL_model] = min(NLL_sum,[],2);
AIC_sumvalues_diff = AIC_sum - diag(AIC_sum);
[~,min_AIC_model] = min(AIC_sum,[],2);
BIC_sumvalues_diff = BIC_sum - diag(BIC_sum);
[~,min_BIC_model] = min(BIC_sum,[],2);
%% Manuscript polished figure
num_colors = 64;
half_colors = num_colors/2;
green_to_white = [linspace(0,1,half_colors)', linspace(0.5,1,half_colors)', linspace(0,1,half_colors)']; % dark blue→white
white_to_red = [linspace(1,1,half_colors)', linspace(1,0,half_colors)', linspace(1,0.5,half_colors)']; % white→dark red
diverging_cmap = flip([green_to_white; white_to_red]);
% Data stack (3 x 6 x 6)
matrices = zeros(3,6,6);
matrices(1,:,:) = NLL_sumvalues_diff;
matrices(2,:,:) = AIC_sumvalues_diff;
matrices(3,:,:) = BIC_sumvalues_diff;
% Panel letters
panel_letters = {'(a)','(b)','(c)'};
% Original names (as currently aligned with matrices rows/cols = 1..6)
model_names = {'Const-SingleGaussian', ...
'Exp-GaussianLaplace', ...
'Const-SingleGaussian-1', ...
'Exp-SingleGaussian', ...
'Const-GaussianLaplace', ...
'Exp-GaussianLaplace-1'};
nModels = numel(model_names);
% --- NEW desired order for both axes ---
new_order_names = { ...
'Exp-GaussianLaplace', ...
'Exp-GaussianLaplace-1', ...
'Exp-SingleGaussian', ...
'Const-GaussianLaplace', ...
'Const-SingleGaussian', ...
'Const-SingleGaussian-1'};
% Map names -> indices in the current matrices
new_idx = cellfun(@(nm) find(strcmp(model_names, nm), 1, 'first'), new_order_names);
% Precompute signed-log to unify color limits across all tiles (order doesn't affect max)
A_signedlog_all = zeros(size(matrices));
for i = 1:3
A0 = squeeze(matrices(i,:,:));
A_signedlog_all(i,:,:) = sign(A0) .* log10(1 + abs(A0));
end
clim_abs = max(abs(A_signedlog_all(:)));
% Figure & layout
figure('Position',[50,100,1200,400]);
t = tiledlayout(1,3,'TileSpacing','compact','Padding','tight'); % annotations use figure coords
ax = gobjects(1,3);
for subplot_idx = 1:3
ax(subplot_idx) = nexttile; hold on;
% --- Apply the row/col permutation by name ---
A0 = squeeze(matrices(subplot_idx,:,:));
A = A0(new_idx, new_idx); % reorder rows and columns
A_signedlog = sign(A) .* log10(1 + abs(A));
% Heatmap
imagesc(A_signedlog);
set(gca,'YDir','reverse');
axis image tight
% Unified color scale (no per-axes colorbar)
caxis([-clim_abs, clim_abs]);
% Tick setup:
set(gca, 'XTick', 1:nModels, 'XTickLabel', '', 'TickLabelInterpreter','none');
xtickangle(45);
% Y: show labels only on the first subplot to reduce clutter; others keep ticks but hide labels
set(gca, 'YTick', 1:nModels, 'YTickLabel', new_order_names, 'TickLabelInterpreter','none');
if subplot_idx > 1
set(gca, 'YTickLabel', []);
end
% Axis labels per subplot (no shared x-axis)
xlabel('Fitted model');
if subplot_idx == 1
ylabel('Ground-truth generative model');
end
% Overlay original (untransformed) values (matching the permuted grid)
[cols, rows] = deal(size(A,2), size(A,1));
[X, Y] = meshgrid(1:cols, 1:rows);
textStrings = compose('%.2f', A(:));
high_contrast = abs(A_signedlog(:)) > 0.6*clim_abs;
textColors = repmat([0 0 0], numel(textStrings), 1);
textColors(high_contrast,:) = repmat([1 1 1], sum(high_contrast), 1);
for k = 1:numel(textStrings)
text(X(k), Y(k), textStrings{k}, ...
'HorizontalAlignment','center', ...
'VerticalAlignment','middle', ...
'Color', textColors(k,:), ...
'FontSize', 9, 'FontWeight','bold');
end
set(gca,'TickDir','out','Box','on');
box off;
end
% One unified colorbar on the very right (compatible with older MATLAB)
colormap(gcf, diverging_cmap);
cb = colorbar(ax(end));
cb.Layout.Tile = 'east';
cb.Label.String = 'signed log_{10}(1 + |\Delta|) · sign(\Delta)';
% --- Panel letters ---
for i = 1:3
pos = ax(i).Position;
x = pos(1); y = pos(2) + pos(4);
dx = -0.03; dy = 0.05;
if i == 1
x = x - 0.01;
end
annotation('textbox', [x+dx, y+dy, 0.01, 0.01], ...
'String', panel_letters{i}, ...
'FontSize', 11, 'FontWeight','bold', ...
'LineStyle', 'none', 'HorizontalAlignment','left', 'VerticalAlignment','top');
end
exportgraphics(gcf,figpath+"UAV_ModelRecovery"+".pdf","ContentType","vector");
%% BELOW: Helper Functions
function [UnimodalData_ModelComparison_FinalTables] = unimodaldata_modelcomparison_visualize(priors, noises, rescales, model_types, num_params, is_plot, fontsize, model_path, data_path, plot_BIC_only)
num_models = length(priors);
num_subjects = 15;
load(data_path+"data_stratified_UV.mat");
load(data_path+"data_stratified_UA.mat");
data_UV = data_stratified_to_data(data_stratified_UV, false, true); % last argument is is_visual.
data_UA = data_stratified_to_data(data_stratified_UA, false, false);
n_data = zeros(1,15);
for subjidx=1:num_subjects
n_data(subjidx) = length(data_UA{subjidx}) +length(data_UV{subjidx});
end
% AIC, BIC
NLLs = zeros(num_models, num_subjects);
AICs = zeros(num_models, num_subjects);
BICs = zeros(num_models, num_subjects);
% load nonparam indv UJoint fits
for model=1:(num_models)
if(model_types(model)=="semiparametric")
load(model_path + "fittedparams_UJoint_Semiparam_rescalefree_lapseUniform.mat")
[min_val, min_idx] = min(F_vals,[],1);
[num_inits,num_params_nonparam]=size(theta_fitted);
NLLs(model,:) = min_val;
AICs(model,:) = 2.*min_val + 2.* num_params(model);
BICs(model,:) = 2.*min_val + num_params(model).*log(n_data);
else
load(model_path + "fittedparams_UJoint_"+noises(model)+"-"+priors(model)+"_rescale"+rescales(model)+"_lapseUniform.mat")
[min_val, min_idx] = min(F_vals,[],2);
NLLs(model,:) = min_val';
AICs(model,:) = 2.*min_val' + 2.* num_params(model);
BICs(model,:) = 2.*min_val' + num_params(model).*log(n_data);
end
end
NLL_sum = sum(NLLs,2);
AIC_sum = sum(AICs,2);
BIC_sum = sum(BICs,2);
NLL_sumvalues_diff = NLL_sum - min(NLL_sum);
[~,min_NLL_model] = min(NLL_sum);
AIC_sumvalues_diff = AIC_sum - min(AIC_sum);
[~,min_AIC_model] = min(AIC_sum);
BIC_sumvalues_diff = BIC_sum - min(BIC_sum);
[~,min_BIC_model] = min(BIC_sum);
% AIC/BIC bootstrapping
num_bootstrap_samps = 100000;
NLL_sum_bootstraps = zeros(num_models,num_bootstrap_samps);
AIC_sum_bootstraps = zeros(num_models,num_bootstrap_samps);
BIC_sum_bootstraps = zeros(num_models,num_bootstrap_samps);
rng("default")
rng(0)
for samp = 1:num_bootstrap_samps
sampled_subj = datasample(1:num_subjects,num_subjects);
for model=1:(num_models)
NLL_sum_bootstraps(model, samp) = sum(NLLs(model, sampled_subj));
AIC_sum_bootstraps(model, samp) = sum(AICs(model, sampled_subj));
BIC_sum_bootstraps(model, samp) = sum(BICs(model, sampled_subj));
end
end
NLL_bootstraps_diff = NLL_sum_bootstraps - NLL_sum_bootstraps(min_NLL_model,:);
AIC_bootstraps_diff = AIC_sum_bootstraps - AIC_sum_bootstraps(min_AIC_model,:);
BIC_bootstraps_diff = BIC_sum_bootstraps - BIC_sum_bootstraps(min_BIC_model,:);
NLL_bootstraps_errorbars = prctile(NLL_bootstraps_diff,[2.5, 97.5], 2);
AIC_bootstraps_errorbars = prctile(AIC_bootstraps_diff,[2.5, 97.5], 2);
BIC_bootstraps_errorbars = prctile(BIC_bootstraps_diff,[2.5, 97,5], 2);
% Plot
bootstraps_errorbars_allstats = {NLL_bootstraps_errorbars, AIC_bootstraps_errorbars, BIC_bootstraps_errorbars};
allstats = {NLL_sumvalues_diff, AIC_sumvalues_diff, BIC_sumvalues_diff};
stat_names = ["\DeltaNLL", "\DeltaAIC", "\DeltaBIC"];
if(is_plot)
if(plot_BIC_only)
statistics = 3;
else
statistics = 1:length(stat_names)
end
tiledlayout(length(statistics),1, 'TileSpacing', 'tight','Padding', 'none')
for statistic=statistics
nexttile
set(gca,'TickDir','out');
bootstraps_errorbars = bootstraps_errorbars_allstats{statistic};
mean_stat = allstats{statistic};
hold on
%for model=[1,2,3]
cats = model_types;
cats = insertBefore(cats,"_no","\");
C = categorical(cats);
C = reordercats(C,cellstr(C)');
barh(0:(length(C)-1),mean_stat,'FaceColor','k', 'FaceAlpha',0.2)
errorbar(mean_stat,0:(length(C)-1),mean_stat-squeeze(bootstraps_errorbars(:, 1)),squeeze(bootstraps_errorbars(:, 2))-mean_stat,'horizontal', 'k.')
yticks(0:(length(C)-1))
C_capitalized = C;
for c=1:length(C)
cat_char = char(string(C(c)));
C_capitalized(c) = convertCharsToStrings([upper(cat_char(1)), cat_char(2:end)]);
end
yticklabels(C_capitalized)
ylim([0-0.5, length(C)-0.5])
set(gca,'YDir','reverse')
%end
if(statistic~=length(stat_names))
xticklabels([]);
set(gca,'xtick',[])
end
set(gca,'FontSize',9)
xlabel(stat_names(statistic))
%set(gca,'xticklabel',["diff","max","ent"].')
end
end
% Create a .mat file with the delta NLL, AIC, and BIC tables
UnimodalData_ModelComparison_FinalTables = cell(1,3);
for stat=1:3
final_table = zeros(length(allstats{1}),3);
final_table(:,1) = bootstraps_errorbars_allstats{stat}(:,1); % 2.5% percentile
final_table(:,2) = allstats{stat}; % Sum
final_table(:,3) = bootstraps_errorbars_allstats{stat}(:,2); % 97.5% percentile
UnimodalData_ModelComparison_FinalTables{stat} = final_table;
end
end
%%
function [UnimodalData_ModelComparison_FinalTables] = unimodaldata_modelcomparison_visualize_uniformgaussianlapse(model_path, data_path)
num_models = 2;
num_subjects=15;
NLLs = [];
load(model_path + "fittedparams_UJoint_exp-GaussianLaplaceBothFixedZero_rescalefree_lapseUniform.mat")
NLLs = [NLLs; min(F_vals,[],2)'];
load(model_path + "fittedparams_UJoint_exp-GaussianLaplaceBothFixedZero_rescalefree_lapseGaussian.mat")
NLLs = [NLLs; min(F_vals,[],2)'];
delta_NLLs = NLLs - NLLs(1,:);
NLL_allmodels_sumdiff = sum(delta_NLLs,2);
num_model_params = [14;15];
AICs = 2.*(NLLs + num_model_params);
delta_AICs = AICs - AICs(1,:);
AIC_allmodels_sumdiff = sum(delta_AICs,2);
load(data_path+"data_stratified_UV.mat");
load(data_path+"data_stratified_UA.mat");
data_UV = data_stratified_to_data(data_stratified_UV, false, true); % last argument is is_visual.
data_UA = data_stratified_to_data(data_stratified_UA, false, false);
n_data = zeros(1,15);
for subjidx=1:15
n_data(subjidx) = length(data_UA{subjidx}) +length(data_UV{subjidx});
end
BICs = 2.*NLLs + num_model_params.*log(n_data);
delta_BICs = BICs - BICs(1,:);
BIC_allmodels_sumdiff = sum(delta_BICs,2);
% AIC/BIC bootstrapping
num_bootstrap_samps = 100000;
NLL_sum_bootstraps = zeros(num_models,num_bootstrap_samps);
AIC_sum_bootstraps = zeros(num_models,num_bootstrap_samps);
BIC_sum_bootstraps = zeros(num_models,num_bootstrap_samps);
rng('default')
rng(0)
for samp = 1:num_bootstrap_samps
sampled_subj = datasample(1:num_subjects,num_subjects);
for model=1:num_models
NLL_sum_bootstraps(model, samp) = sum(NLLs(model, sampled_subj));
AIC_sum_bootstraps(model, samp) = sum(AICs(model, sampled_subj));
BIC_sum_bootstraps(model, samp) = sum(BICs(model, sampled_subj));
end
end
NLL_bootstraps_diff = NLL_sum_bootstraps - NLL_sum_bootstraps(1,:);
AIC_bootstraps_diff = AIC_sum_bootstraps - AIC_sum_bootstraps(1,:);
BIC_bootstraps_diff = BIC_sum_bootstraps - BIC_sum_bootstraps(1,:);
NLL_bootstraps_errorbars = prctile(NLL_bootstraps_diff,[2.5, 97.5], 2);
AIC_bootstraps_errorbars = prctile(AIC_bootstraps_diff,[2.5, 97.5], 2);
BIC_bootstraps_errorbars = prctile(BIC_bootstraps_diff,[2.5, 97,5], 2);
% Plot
bootstraps_errorbars_allstats = {NLL_bootstraps_errorbars, AIC_bootstraps_errorbars, BIC_bootstraps_errorbars};
allstats = {NLL_allmodels_sumdiff, AIC_allmodels_sumdiff, BIC_allmodels_sumdiff};
stat_names = ["\DeltaNLL", "\DeltaAIC", "\DeltaBIC"];
if(false)
tiledlayout(length(stat_names),1, 'TileSpacing', 'tight','Padding', 'none')
for statistic=1:length(stat_names)
nexttile
set(gca,'TickDir','out');
bootstraps_errorbars = bootstraps_errorbars_allstats{statistic};
mean_stat = allstats{statistic};
hold on
bar(1:(num_models),mean_stat,'FaceColor','k', 'FaceAlpha',0.2)
errorbar(1:(num_models)',mean_stat,mean_stat-squeeze(bootstraps_errorbars(:, 1)),squeeze(bootstraps_errorbars(:, 2))-mean_stat, 'k.')
if(statistic~=length(stat_names))
xticklabels([])
set(gca,'xtick',[])
else
xticks(1:(3*num_models));
xticklabels(allmodel_xticklabels)
xtickangle(20)
end
set(gca,'FontSize',9)
ylabel(stat_names(statistic), 'FontSize',9)
xlim([0, num_models+0.7])
end
end
% Create a .mat file with the delta NLL, AIC, and BIC tables
UnimodalData_ModelComparison_FinalTables = cell(1,3);
for stat=1:3
final_table = zeros(num_models,3);
final_table(:,1) = bootstraps_errorbars_allstats{stat}(:,1); % 2.5% percentile
final_table(:,2) = allstats{stat}; % Sum
final_table(:,3) = bootstraps_errorbars_allstats{stat}(:,2); % 97.5% percentile
UnimodalData_ModelComparison_FinalTables{stat} = final_table;
end
end
%%
function [AllData_ModelComparison_FinalTables] = alldata_modelcomparison_visualize(causal_inf_strategies, param_model_names, is_plot, fontsize, model_path, data_path, plot_BIC_only, LiftedSemiparamParamAll, new_figure)
num_strategies = length(causal_inf_strategies);
num_models = length(param_model_names);
num_subjects = 15;
if(LiftedSemiparamParamAll)
barcolor = [1,1,1];
else
barcolor = [204,204,204]./255;
end
out_structs = cell(1,num_strategies);
NLL_allmodels = zeros(num_models*num_strategies, 15);
AIC_allmodels = zeros(num_models*num_strategies, 15);
BIC_allmodels = zeros(num_models*num_strategies, 15);
for causal_inf_strategy_idx=1:num_strategies
causal_inf_strategy_idx
out_struct = alldata_modelcomparison_visualize_helper(causal_inf_strategies(causal_inf_strategy_idx), param_model_names, model_path, data_path, LiftedSemiparamParamAll);
out_structs{causal_inf_strategy_idx} = out_struct;
NLL_allmodels(causal_inf_strategy_idx:num_strategies:end,:) = out_struct.NLLs;
AIC_allmodels(causal_inf_strategy_idx:num_strategies:end,:) = out_struct.AICs;
BIC_allmodels(causal_inf_strategy_idx:num_strategies:end,:) = out_struct.BICs;
end
% deltaNLL and deltaAIC across all 18 models
baseline_modelidx = find(causal_inf_strategies=="ProbMatching"); % Use the first model's PM version as baseline.
NLL_allmodels_sum = sum(NLL_allmodels,2);
NLL_allmodels_sumdiff = NLL_allmodels_sum - NLL_allmodels_sum(baseline_modelidx);
NLL_min_model = baseline_modelidx;
AIC_allmodels_sum = sum(AIC_allmodels,2);
AIC_allmodels_sumdiff = AIC_allmodels_sum - AIC_allmodels_sum(baseline_modelidx);
AIC_min_model = baseline_modelidx;
BIC_allmodels_sum = sum(BIC_allmodels,2);
BIC_allmodels_sumdiff = BIC_allmodels_sum - BIC_allmodels_sum(baseline_modelidx);
BIC_min_model = baseline_modelidx;
%causal_inf_strategies_abbrev = strrep(causal_inf_strategies,["ModelSelection", "ModelAveraging","ProbMatching"], "-"+["MS","MA","PM"]);
causal_inf_strategies_abbrev = [];
for strategy=1:num_strategies
if(causal_inf_strategies(strategy)=="ModelSelection")
causal_inf_strategies_abbrev = [causal_inf_strategies_abbrev,"-MS"];
elseif(causal_inf_strategies(strategy)=="ModelAveraging")
causal_inf_strategies_abbrev = [causal_inf_strategies_abbrev,"-MA"];
else(causal_inf_strategies(strategy)=="ProbMatching")
causal_inf_strategies_abbrev = [causal_inf_strategies_abbrev,"-PM"];
end
end
allmodel_xticklabels = param_model_names' + causal_inf_strategies_abbrev;
allmodel_xticklabels = allmodel_xticklabels';
allmodel_xticklabels = allmodel_xticklabels(:);
for lab_idx = 1:length(allmodel_xticklabels)
lab = char(allmodel_xticklabels(lab_idx));
allmodel_xticklabels(lab_idx) = convertCharsToStrings([upper(lab(1)) lab(2:end)]);
end
% AIC/BIC bootstrapping
num_bootstrap_samps = 100000;
NLL_sum_bootstraps = zeros(num_strategies*num_models,num_bootstrap_samps);
AIC_sum_bootstraps = zeros(num_strategies*num_models,num_bootstrap_samps);
BIC_sum_bootstraps = zeros(num_strategies*num_models,num_bootstrap_samps);
rng('default')
rng(0)
for samp = 1:num_bootstrap_samps
sampled_subj = datasample(1:num_subjects,num_subjects);
for model=1:(num_strategies*num_models)
NLL_sum_bootstraps(model, samp) = sum(NLL_allmodels(model, sampled_subj));
AIC_sum_bootstraps(model, samp) = sum(AIC_allmodels(model, sampled_subj));
BIC_sum_bootstraps(model, samp) = sum(BIC_allmodels(model, sampled_subj));
end
end
NLL_bootstraps_diff = NLL_sum_bootstraps - NLL_sum_bootstraps(NLL_min_model,:);
AIC_bootstraps_diff = AIC_sum_bootstraps - AIC_sum_bootstraps(AIC_min_model,:);
BIC_bootstraps_diff = BIC_sum_bootstraps - BIC_sum_bootstraps(BIC_min_model,:);
NLL_bootstraps_errorbars = prctile(NLL_bootstraps_diff,[2.5, 97.5], 2);
AIC_bootstraps_errorbars = prctile(AIC_bootstraps_diff,[2.5, 97.5], 2);
BIC_bootstraps_errorbars = prctile(BIC_bootstraps_diff,[2.5, 97,5], 2);
% Plot
bootstraps_errorbars_allstats = {NLL_bootstraps_errorbars, AIC_bootstraps_errorbars, BIC_bootstraps_errorbars};
allstats = {NLL_allmodels_sumdiff, AIC_allmodels_sumdiff, BIC_allmodels_sumdiff};
stat_names = ["\DeltaNLL", "\DeltaAIC", "\DeltaBIC"];
if(is_plot)
if(plot_BIC_only)
statistics = 3;
else
statistics = 1:length(stat_names);
end
if(new_figure)
tiledlayout(length(statistics),1, 'TileSpacing', 'tight','Padding', 'none')
end
stat_idx = 1;
for statistic=statistics
if(new_figure || length(statistics)>1)
nexttile(stat_idx)
end
set(gca,'TickDir','out');
bootstraps_errorbars = bootstraps_errorbars_allstats{statistic};
mean_stat = allstats{statistic};
hold on
barh(1:(num_strategies*num_models),mean_stat,'FaceColor',barcolor, 'FaceAlpha',1)
errorbar(mean_stat,(1:(num_strategies*num_models))', mean_stat-squeeze(bootstraps_errorbars(:, 1)),squeeze(bootstraps_errorbars(:, 2))-mean_stat,'horizontal', 'k.')
yticks(1:(num_strategies*num_models));
yticklabels(allmodel_xticklabels)
ylim([0.2, num_strategies*num_models+0.7])
set(gca, 'YDir','reverse')
set(gca,'FontSize',9)
xlabel(stat_names(statistic), 'FontSize', 9)
stat_idx=stat_idx+1;
end
end
% Create a .mat file with the delta NLL, AIC, and BIC tables
AllData_ModelComparison_FinalTables = cell(1,3);
for stat=1:3
final_table = zeros(num_models*num_strategies,3);
final_table(:,1) = bootstraps_errorbars_allstats{stat}(:,1); % 2.5% percentile
final_table(:,2) = allstats{stat}; % Sum
final_table(:,3) = bootstraps_errorbars_allstats{stat}(:,2); % 97.5% percentile
AllData_ModelComparison_FinalTables{stat} = final_table;
end
end
%%
function [out_struct_allmodels] = alldata_modelcomparison_visualize_helper(causal_inf_strategy, param_model_names, model_path, data_path, LiftedSemiparamParamAll)
if(nargin==0)
causal_inf_strategy = "ModelSelection";
is_plot=false;
elseif(nargin==1)
is_plot=false;
end
num_subjects = 15;
% Get number of observations for BIC
load(data_path+"BAV_data.mat");
load(data_path+"BC_data.mat");
load(data_path+"data_stratified_UV.mat");
load(data_path+"data_stratified_UA.mat");
data_UV = data_stratified_to_data(data_stratified_UV, false, true); % last argument is is_visual.
data_UA = data_stratified_to_data(data_stratified_UA, false, false);
UAV_data = cell(1,num_subjects);
for i=1:num_subjects
data_UA{i}(:,3) = 4;
UAV_data{i} = [data_UV{i}; data_UA{i}];
end
n_data = zeros(1,15);
n_data_bytasktype = zeros(4,15);
for subjidx=1:num_subjects
n_data(subjidx) = length(BAV_data{subjidx}) + length(BC_data{subjidx}) + length(UAV_data{subjidx});
n_data_bytasktype(1,subjidx) = length(data_UV{subjidx});
n_data_bytasktype(2,subjidx) = length(data_UA{subjidx});
n_data_bytasktype(3,subjidx) = length(BC_data{subjidx});
n_data_bytasktype(4,subjidx) = length(BAV_data{subjidx});
end
% save("NumTrials_allsubjects",'n_data_bytasktype','n_data')
colors = brewermap(10,"Accent");
colors = colors([1,2,3,5:10],:);
noises = ["exp","exp","constant","exp","constant"];
priors = ["GaussianLaplaceBothFixedZero", "SingleGaussian", "GaussianLaplaceBothFixedZero", "TwoGaussiansBothFixedZero","SingleGaussian"];
helper_model_order = strrep(noises+"-"+priors,"constant","const");
helper_model_order = strrep(helper_model_order,"BothFixedZero","");
helper_model_order = ["paramBest","LiftedSemiparam",helper_model_order]; % Order used by this helper function to read files.
num_models = length(priors);
NLLs_param = zeros(num_models,15);
num_parametric_model_params = zeros(num_models,1);
idx=0;
for i=1:num_models
idx = idx+1;
filename = "fittedparams_All_UBresc_"+noises(i)+"-"+priors(i)+"-"+causal_inf_strategy+"_rescalefree_lapseUniform.mat";
load(model_path + filename)
NLLs_param(idx,:) = min(F_vals,[],2);
num_parametric_model_params(idx) = length(theta_fitted(1,:));
if(noises(i)=="constant") % Remove place filler zeros in theta_fitted for k_vis, k_aud.
num_parametric_model_params(idx) = num_parametric_model_params(idx)-2;
end
end
[NLL_param_min, NLL_param_min_idx] = min(NLLs_param, [],1);
% LiftedSemiparam model for this causal_inf_strategy
filename_basis = "fittedparams_All_UBresc_SemiparamInspired_"+causal_inf_strategy+"_rescalefree_lapseUniform.mat";
load(model_path + filename_basis)
num_LiftedSemiparam_params = 9; %12+12+11
if(LiftedSemiparamParamAll) % Count in the sigma_v(s), sigma_a(s), p(s) pivot points from the semiparametric fits
num_LiftedSemiparam_params = num_LiftedSemiparam_params + 12+12+11;
end
[F_min_val, F_min_idx] = min(F_vals,[],1);
% Order: paramBest,LiftedSemiParam, ParamModels.
NLLs = [NLL_param_min;F_min_val;NLLs_param];
% AIC
AIC_LiftedSemiparam_best = 2.*(num_LiftedSemiparam_params+F_min_val);
AIC_param_best = 2.*(num_parametric_model_params(NLL_param_min_idx)'+NLL_param_min);
AICs = [AIC_param_best; AIC_LiftedSemiparam_best];
for model=1:num_models
AICs = [AICs; 2.*(num_parametric_model_params(repmat(model,1,num_subjects))'+NLLs_param(model,:))];
end
% BIC
%load("NumTrials_allsubjects")
BIC_LiftedSemiparam_best = 2.*F_min_val + num_LiftedSemiparam_params.*log(n_data);
BIC_param_best = 2.*NLL_param_min + num_parametric_model_params(NLL_param_min_idx)'.*log(n_data);
BICs = [BIC_param_best; BIC_LiftedSemiparam_best];
for model=1:num_models
BICs = [BICs; 2.*NLLs_param(model,:) + num_parametric_model_params(repmat(model,1,num_subjects))'.*log(n_data)];
end
[~,prompted_model_order] = ismember(param_model_names,helper_model_order);
out_struct_allmodels.NLLs = NLLs(prompted_model_order,:);
out_struct_allmodels.AICs = AICs(prompted_model_order,:);
out_struct_allmodels.BICs = BICs(prompted_model_order,:);
out_struct_allmodels.n_data = n_data;
end
%% Plot nonparam fitted sigma(s) and p(s) shapes
function [] = semiparam_sigmafun_prior_visualization(fontsize, figspec, model_path)
num_subjects = 15;
num_params = 40;
num_iters = 81*15;
num_inits = 81;
colors = brewermap(12,"Set3");
colors2 = brewermap(8,"Set2");
colors(2,:) = colors2(2,:);
colors = [colors; colors2([4,7,8],:)];
%colors = colors .* 0.9;
filename_basis = "fittedparams_UJoint_Semiparam_rescalefree_lapseUniform.mat";
load(model_path + filename_basis);
theta_fitted_cmaes = theta_fitted;
F_vals_cmaes = F_vals;
ModelComponents.SPivot = [0,0.1,0.3,1,2,4,6,8,10,15,20,45];
% Plot function shapes
s_pivot = ModelComponents.SPivot;
s_pivot_full = [-fliplr(s_pivot(2:end)), s_pivot];
s_fine = linspace(0,15,2^7);
s_fine_full = linspace(0,45,2^9);
num_pivots = length(s_pivot);
figure('Position', figspec);
set(gcf, 'Color', 'w')
T = tiledlayout(2,10,'TileSpacing','compact', "Padding","none");
linewidth = 1;
for fun_idx=1:3
t = nexttile(T,[1,5]);
set(t,'TickDir','out');
hold(t,'on')
pl = get(t, 'Position');
switch fun_idx
case 1
h = axes('Parent', gcf, 'Position', [pl(1)+pl(3)*.61 pl(2)+pl(4)*0.72 pl(3)*0.35 pl(3)*0.35.*3/4]);
case 2
h = axes('Parent', gcf, 'Position', [pl(1)+pl(3)*.65 pl(2)+.07 pl(3)*0.33 pl(3)*0.35.*3/4]);
case 3
h = axes('Parent', gcf, 'Position', [pl(1)+pl(3)*.61 pl(2)+pl(4)*0.75 pl(3)*0.35 pl(3)*0.35.*3/4]);
end
%box(h,'on');
hold(h,'on')
set(h,'TickDir','out');
for subj=1:num_subjects
%[fun_idx, subj]
theta=squeeze(theta_fitted_cmaes(subj,:));
color = colors(subj,:);
switch fun_idx
case 1
sigma_fun_vis_rel_high_pivots = cumsum(theta(1:num_pivots));
sigma_fun_vis_rel_high_pivots = [fliplr(sigma_fun_vis_rel_high_pivots(2:end)), sigma_fun_vis_rel_high_pivots];
sigma_fun_vis = @(s) min([repmat(45,length(s),1)' ; interp1(s_pivot_full, exp(sigma_fun_vis_rel_high_pivots), s, 'pchip')], [], 1);
p=plot(t,s_fine, sigma_fun_vis(s_fine),'-','Color', color, 'LineWidth',linewidth);
p.Color(4)=0.9;
scatter1 = scatter(t,s_pivot_full, exp(sigma_fun_vis_rel_high_pivots),'o','MarkerFaceColor',color,'MarkerEdgeColor',color);
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = linewidth.*5;
ylabel(t,"$\sigma_{\mathrm{V}}(s)$", 'interpreter','latex', 'FontSize', fontsize)
xlabel(t,"Visual stimulus location (\circ)", 'FontSize', fontsize)
xlim(t,[0,15])
ylim(t,[0,6.1])
xl = xlim(t); yl = ylim(t);
yticks(t,0:6)
t.XAxis.FontSize = 9;
t.YAxis.FontSize = 9;
p=plot(h, s_fine_full, sigma_fun_vis(s_fine_full),'-', 'Color',color);
p.Color(4)=0.9;
scatter1 = scatter(h, s_pivot_full, min(45,exp(sigma_fun_vis_rel_high_pivots)),'o','MarkerFaceColor',color,'MarkerEdgeColor',color);
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = 1;
rectangle('Position',[xl(1) yl(1) xl(2)-xl(1) yl(2)-yl(1)])
%ylabel(h,"$\sigma_{\mathrm{V}}(s)$", 'interpreter','latex', 'FontSize', fontsize)
xlim(h,[0,45])
ylim(h,[0,45])
xticks(h,0:15:45)
yticks(h,0:15:45)
h.XAxis.FontSize = 9;
h.YAxis.FontSize = 9;
case 2
sigma_fun_aud_pivots = cumsum(theta((num_pivots+1):(2*num_pivots)));
sigma_fun_aud_pivots = [fliplr(sigma_fun_aud_pivots(2:end)), sigma_fun_aud_pivots];
sigma_fun_aud = @(s) min([repmat(45,length(s),1)' ; interp1(s_pivot_full, exp(sigma_fun_aud_pivots), s, 'pchip')], [], 1);
p=plot(t,s_fine, sigma_fun_aud(s_fine),'-', 'Color',color, 'LineWidth', linewidth);
p.Color(4)=0.9;
scatter1 = scatter(t,s_pivot_full, exp(sigma_fun_aud_pivots),'o','MarkerFaceColor',color,'MarkerEdgeColor',color);
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = linewidth.*5;
ylabel(t,"$\sigma_{\mathrm{A}}(s)$", 'interpreter','latex', 'FontSize', fontsize)
xlabel(t,"Auditory stimulus location (\circ)", 'FontSize', fontsize)
xlim(t,[0,15])
ylim(t,[0,6.1])
xl = xlim(t); yl = ylim(t);
yticks(t,0:6)
t.XAxis.FontSize = 9;
t.YAxis.FontSize = 9;
p=plot(h, s_fine_full, sigma_fun_aud(s_fine_full),'-', 'Color',color);
p.Color(4)=0.9;
scatter1 = scatter(h, s_pivot_full, min(45,exp(sigma_fun_aud_pivots)),'o','MarkerFaceColor',color,'MarkerEdgeColor',color);
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = 1;
rectangle('Position',[xl(1) yl(1) xl(2)-xl(1) yl(2)-yl(1)])
%ylabel(h,"$\sigma_{\mathrm{A}}(s)$", 'interpreter','latex', 'FontSize', fontsize)
xlim(h,[0,45])
ylim(h,[0,45])
xticks(h,0:15:45)
yticks(h,0:15:45)
h.XAxis.FontSize = 9;
h.YAxis.FontSize = 9;
case 3
prior_pivots = cumsum([1,theta((2*num_pivots+1):(3*num_pivots-1))]);
prior_pivots = [fliplr(prior_pivots(2:end)), prior_pivots];
prior = @(s) exp(interp1(s_pivot_full, prior_pivots, s, 'pchip'));
s_fine_full_width = s_fine_full(2) - s_fine_full(1);
normalization_constant = 1./(qtrapz(prior(s_fine_full).*s_fine_full_width));
p=plot(t,s_fine_full, (prior(s_fine_full).*normalization_constant),'-', 'Color',color, 'LineWidth', linewidth);
p.Color(4)=0.9;
%scatter(s_pivot_full, prior_pivots,"o", "MarkerFaceAlpha",0.1)
scatter1 = scatter(t,s_pivot_full, exp(prior_pivots).*normalization_constant,'o','MarkerFaceColor',color,'MarkerEdgeColor',color);
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = linewidth.*5;
ylabel(t,"$p(s)$", 'interpreter','latex', 'FontSize', fontsize)
xlabel(t,"Visual/Auditory stimulus location (\circ)", 'FontSize', fontsize)
xlim(t,[0,3])
ylim(t,[0,5])
xl = xlim(t); yl = ylim(t);
xticks(t,0:1:3);
yticks(t, 0:1:5);
t.XAxis.FontSize = 9;
t.YAxis.FontSize = 9;
p=plot(h, s_fine_full, log(prior(s_fine_full).*normalization_constant),'-', 'Color',color);
p.Color(4)=0.9;
scatter1 = scatter(h, s_pivot_full, log(exp(prior_pivots).*normalization_constant),'o','MarkerFaceColor',color,'MarkerEdgeColor',color);
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = 1;
%rectangle('Position',[xl(1) yl(1) xl(2)-xl(1) yl(2)-yl(1)])
xlim(h,[0,45])
ylim(h,[-20,3])
xticks(h,0:15:45)
yticks(h, -20:10:0);
h.XAxis.FontSize = 9;
h.YAxis.FontSize = 9;
ylabel(h,"$\log p(s)$", 'interpreter','latex', 'FontSize', fontsize)
end
end
ax.XAxis.FontSize = fontsize;
ax.YAxis.FontSize = fontsize;
end
x_labels_pos = (3*num_pivots):num_params;
x_labels = "$"+{"\alpha_\mathrm{med}", "\alpha_\mathrm{low}", "\lambda","\sigma_\mathrm{motor}","\rho_\mathrm{A}"}+"$";
t=nexttile(T,[1,3]);
boxplot(theta_fitted_cmaes(:,x_labels_pos([1,2,5])), 'Color','k')
hold on
idx=0
for param=[1,2,5]
idx=idx+1
for subj=1:num_subjects
scatter1 = scatter(repmat(idx,1), theta_fitted_cmaes(subj,x_labels_pos(param)),'o','MarkerFaceColor','k','MarkerEdgeColor','none');
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = 5;
end
end
xticks(1:3);
xaxisproperties=get(gca, 'XAxis');
xaxisproperties.TickLabelInterpreter = 'latex';
xticklabels(x_labels([1,2,5]));
%xtickangle(45);
ax=gca;
ax.XAxis.FontSize = fontsize;
xlim([1-0.5,idx+0.5])
ylim([0, Inf])
yticks(0:1:5)
t.XAxis.FontSize = 9;
t.YAxis.FontSize = 9;
box off
set(gca,'TickDir','out');
t=nexttile(T);
boxplot(theta_fitted_cmaes(:,x_labels_pos(3)), 'Color','k')
hold on
for subj=1:num_subjects
scatter1 = scatter(1, theta_fitted_cmaes(subj,x_labels_pos(3)),'o','MarkerFaceColor','k','MarkerEdgeColor','none');
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = 5;
end
% scatter1 = scatter(repmat(1,num_subjects,1), theta_fitted_cmaes(:,x_labels_pos(3)),'o','MarkerFaceColor','k','MarkerEdgeColor','k');
% scatter1.SizeData = 5;
xticks([1])
xtickangle(0);
yticks(0:0.005:0.015)
xaxisproperties=get(gca, 'XAxis');
xaxisproperties.TickLabelInterpreter = 'latex';
xticklabels(x_labels(3));
ax=gca;
ax.XAxis.FontSize = fontsize;
xlim([0.5,1.5])
ylim([0, 0.015])
t.XAxis.FontSize = 9;
t.YAxis.FontSize = 9;
set(gca,'TickDir','out');
box off
t=nexttile(T);
boxplot(theta_fitted_cmaes(:,x_labels_pos(4)), 'Color','k')
hold on
for subj=1:num_subjects
scatter1 = scatter(1, theta_fitted_cmaes(subj,x_labels_pos(4)),'o','MarkerFaceColor','k','MarkerEdgeColor','none');
%scatter1.MarkerFaceAlpha = .2; scatter1.MarkerEdgeAlpha = .2;
scatter1.SizeData = 5;
end
% scatter1 = scatter(repmat(1,num_subjects,1), theta_fitted_cmaes(:,x_labels_pos(4)),'o','MarkerFaceColor','k','MarkerEdgeColor','k');
% scatter1.SizeData = 5;
xticks([1])
yticks(0:0.1:0.5);
xtickangle(0);
xaxisproperties=get(gca, 'XAxis');
xaxisproperties.TickLabelInterpreter = 'latex';
xticklabels(x_labels(4));
ax=gca;
ax.XAxis.FontSize = fontsize;
xlim([0.5,1.5])
ylim([0, 0.5])
t.XAxis.FontSize = 9;
t.YAxis.FontSize = 9;
set(gca,'TickDir','out');
box off
%set(gca,'fontsize', fontsize)
end
%%
function [] = sigmafun_prior_examples(fontsize)
s_grid = -45:0.1:45;
sigma0_vals = [0.5,1,2,3];
colors = brewermap(12,"Set1");
hetero_type = "exp";
sigma_fun_constant = heterotype_to_sigmafun("constant");
sigma_fun_exp = heterotype_to_sigmafun("exp");
subplot(2,3,1)
set(gca,'TickDir','out');
hold on
for i=1:length(sigma0_vals)
plot(s_grid, repmat(sigma0_vals(i), length(s_grid),1), "-", 'Color', colors(i,:));
end
lg = legend("$\sigma_0="+sigma0_vals+"$", 'Interpreter', 'latex', 'FontSize', fontsize);
set(lg,'Box','off')
ylim([0,6])
yticks(0:1:6)
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
ylabel("$\sigma(s)$", 'Interpreter', 'latex', 'FontSize', fontsize)
title("Constant sensory noise", 'FontSize', fontsize+1)
sigma0_vals = [0.5,1,1,1,3];
k1_vals = [1,1,2,1,1];
k2_vals = [0.1,0.1,0.1,0.5 0.5];
subplot(2,3,2)
set(gca,'TickDir','out');
hold on
for i=1:length(sigma0_vals)
plot(s_grid, sigma_fun_exp(s_grid,sigma0_vals(i), [k1_vals(i),k2_vals(i)]), "-", 'Color', colors(i,:));
end
lg = legend("$\sigma_0="+sigma0_vals+", k_1="+k1_vals+", k_2="+k2_vals+"$", 'Interpreter', 'latex', 'FontSize', fontsize);
set(lg,'Box','off')
lg.Position(1:2) = [0.63,0.8];
ylim([0,6])
yticks(0:1:6)
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
title("Exponential sensory noise", 'FontSize', fontsize+1)
% Priors
sigma_s_vals = [3, 5, 8, 10];
subplot(2,3,4)
set(gca,'TickDir','out');
hold on
for i=1:length(sigma_s_vals)
plot(s_grid, normpdf(s_grid, 0, sigma_s_vals(i)), "-", 'Color', colors(i,:))
end
lg = legend("$\sigma_s="+sigma_s_vals+"$", 'Interpreter', 'latex', 'FontSize', fontsize);
set(lg,'Box','off')
lg.Position(1) = 0.25;
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
ylabel("$p(s)$", 'Interpreter', 'latex', 'FontSize', fontsize)
title("SingleGaussian prior", 'FontSize', fontsize+1)
ylim([0,0.15])
yticks(0:0.05:0.15)
sigma_s_vals = [8, 8, 8, 15].*2;
b_vals = [1,2,1,2].*2;
w_vals = [0.3, 0.3, 0.5, 0.5];
subplot(2,3,6)
set(gca,'TickDir','out');
hold on
for i=1:length(sigma_s_vals)
plot(s_grid, (1-w_vals(i)).*normpdf(s_grid, 0, sigma_s_vals(i)) + w_vals(i).*1./(2.*b_vals(i)).*exp(-abs(s_grid)./b_vals(i)), "-", 'Color', colors(i,:))
end
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
lg = legend("$\sigma_s="+sigma_s_vals+", b="+b_vals+", w="+w_vals+"$", 'Interpreter', 'latex', 'FontSize', fontsize);
set(lg,'Box','off')
lg.Position(1:2) = [0.69,0.48];
title("GaussianLaplace prior", 'FontSize', fontsize+1)
ylim([0,0.15])
yticks(0:0.05:0.15)
sigma_s_vals = [3,5,8,10];
sigma_s2_vals = [8, 8, 8, 15].*2 - sigma_s_vals;
w_vals = [0.3, 0.3, 0.7, 0.7];
subplot(2,3,5)
set(gca,'TickDir','out');
hold on
for i=1:length(sigma_s_vals)
plot(s_grid, (1-w_vals(i)).*normpdf(s_grid, 0, sigma_s_vals(i)) + w_vals(i).*normpdf(s_grid, 0, sigma_s_vals(i)+sigma_s2_vals(i)), "-", 'Color', colors(i,:))
end
ylim([0,0.15])
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
lg = legend("$\sigma_s="+sigma_s_vals+", \sigma_{\Delta}="+sigma_s2_vals+", w="+w_vals+"$", 'Interpreter', 'latex', 'FontSize', fontsize);
set(lg,'Box','off')
lg.Position(1:2) = [0.42, 0.348];
title("TwoGaussians prior", 'FontSize', fontsize+1)
ylim([0,0.15])
yticks(0:0.05:0.15)
%set(gca,'fontsize', fontsize)
end
%%
function [] = allindvsubjplots_to_onesubjplot(save_name, subjidx, fitted_on_all_data, fontsize, figspecs, figpath)
% This function assumes that the individual-level plots have been saved
% as .fig files.
close all;
if(~fitted_on_all_data)
F1 = openfig(figpath + save_name + "_Individualmean.fig");
t1 = nexttile(subjidx);
ax1=gca;
F2 = openfig(figpath + save_name + "_IndividualSD.fig");
t2 = nexttile(subjidx);
ax2=gca;
figure('Position', figspecs);
set(gcf, 'Color', 'w')
T=tiledlayout(1,2,'Padding', 'tight', 'TileSpacing', 'tight');
t1 = nexttile(1);
set(gca,'TickDir','out');
hold on
plot([-20,20],[0,0],"k--",'HandleVisibility','off');
fig1 = get(ax1,'children');
copyobj(fig1, t1);
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
ylabel("Bias (\circ)", 'FontSize', fontsize)
ylim([-20,20])
xticks(-20:10:20)
set(gca,"FontSize",9)
t2 = nexttile(2);
fig2 = get(ax2,'children');
set(gca,'TickDir','out');
copyobj(fig2, t2);
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
ylabel("SD of location response (\circ)", 'FontSize', fontsize)
h = findall(gca, 'LineStyle', '-');
for i=1:4
h(i).HandleVisibility="off";
end
lg = legend("Visual (high reliability)","Visual (med. reliability)", "Visual (low reliability)", "Auditory");
set(lg,'Box','off')
lg.FontSize = max(9,fontsize-1);
lg.Location="northeast";
lg.ItemTokenSize(1) = 10;
xticks(-20:10:20)
set(gca,"FontSize",9)
else
F1 = openfig(figpath + save_name + "-UAV_Individualmean.fig");
t1 = nexttile(subjidx);
ax1=gca;
F2 = openfig(figpath + save_name + "-UAV_IndividualSD.fig");
t2 = nexttile(subjidx);
ax2=gca;
F3 = openfig(figpath + save_name + "-BC_Individual.fig");
ax3_center = F3.Children.Children((end-subjidx+1)).Children(2);
ax3_periphery = F3.Children.Children((end-subjidx+1)).Children(1);
F4 = openfig(figpath + save_name + "-BV_Individual.fig");
ax4_right = F4.Children.Children((end-subjidx+1)).Children(1);
ax4_center = F4.Children.Children((end-subjidx+1)).Children(2);
ax4_left = F4.Children.Children((end-subjidx+1)).Children(3);
F5 = openfig(figpath + save_name + "-BA_Individual.fig");
ax5_right = F5.Children.Children((end-subjidx+1)).Children(1);
ax5_center = F5.Children.Children((end-subjidx+1)).Children(2);
ax5_left = F5.Children.Children((end-subjidx+1)).Children(3);
%% Move to new plot
figure('Position', figspecs);
set(gcf, 'Color', 'w')
T=tiledlayout(2,12,'Padding', 'tight', 'TileSpacing', 'tight');
t12=tiledlayout(T,1,2, 'Padding','none','TileSpacing','tight');
t12.Layout.Tile = 1;
t12.Layout.TileSpan = [1 6];
t1 = nexttile(t12);
hold on
set(gca,'TickDir','out');
plot([-20,20],[0,0],"k--",'HandleVisibility','off');
fig1 = get(ax1,'children');
copyobj(fig1, t1);
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
ylabel("Bias of location response (\circ)", 'FontSize', fontsize)
ylim([-20,20])
ttl = title('(a)', "Fontsize", 10);
ttl.Units = 'Normalize';
ttl.Position(1) = -0.3; % use negative values (ie, -0.1) to move further left
ttl.HorizontalAlignment = 'left';
xticks(-20:10:20)
set(gca,"FontSize",9)
t2 = nexttile(t12);
set(gca,'TickDir','out');
fig2 = get(ax2,'children');
copyobj(fig2, t2);
xlabel("Stimulus location (\circ)", 'FontSize', fontsize)
ylabel("SD of location response (\circ)", 'FontSize', fontsize)
ylim([0,9])
h = findall(gca, 'LineStyle', '-');
for i=1:4
h(i).HandleVisibility="off";
end
lg = legend("Visual (high rel.)","Visual (med. rel.)", "Visual (low rel.)", "Auditory");
set(lg,'Box','off')
lg.FontSize = 9;
lg.ItemTokenSize(1) = 10;
xticks(-20:10:20)
set(gca,"FontSize",9)
% BC
%t3 = nexttile([1,2]);
t3=tiledlayout(T,1,2, 'Padding','none','TileSpacing','compact');
t3.Layout.Tile = 7;
t3.Layout.TileSpan = [1 6];
xlabel(t3,"Stimulus location disparity, {\its}_A– {\its}_V (\circ)", 'FontSize',fontsize)
ylabel(t3,{"{\rm \fontsize{9} {Proportion responding "+ '"'+'same'+ '"'+"}}"}, 'FontSize',fontsize);
BC_strat_names = ["Center", "Periphery"];
for strats=1:2
tt = nexttile(t3);
set(gca,'TickDir','out');
hold on;
if(strats==1)
fig31 = get(ax3_center,'children');
copyobj(fig31, tt);
ttl = title('(b)', "Fontsize", 10);
ttl.Units = 'Normalize';
ttl.Position(1) = -0.3; % use negative values (ie, -0.1) to move further left
ttl.HorizontalAlignment = 'left';
subtitle(tt,"Center",'Fontsize', fontsize, 'FontWeight','bold')
yticks(0:0.2:1)
else
fig32 = get(ax3_periphery,'children');
copyobj(fig32, tt);
subtitle(tt,"Periphery",'Fontsize', fontsize, 'FontWeight','bold')
yticks([])
end
h = findall(gca, 'LineStyle', '-');
for i=1:3
h(i).HandleVisibility="off";
end
xlim([-30,30])
ylim([0,1])
xticks(-30:15:30)
xtickangle(0)
set(gca,"FontSize",9)
lg = legend({"High vis. rel.","Med. vis. rel.","Low vis. rel."});
set(lg,'Box','off')
lg.FontSize = 9;
lg.Location="south";
lg.ItemTokenSize(1) = 10;
end
BAV_strat_names = ["Left","Center","Right"];
t4=tiledlayout(T,1,3, 'Padding','none','TileSpacing','compact');
t4.Layout.Tile = 13;
t4.Layout.TileSpan = [1 6];
xlabel(t4, "Stimulus location disparity, {\its}_A– {\its}_V (\circ)", 'FontSize',fontsize)
ylabel(t4,"{\rm \fontsize{10} {Visual bias (\circ)}}");
for strats=1:3
tt = nexttile(t4);
set(gca,'TickDir','out');
hold on;
if(strats==1)
fig4 = get(ax4_left,'children');
copyobj(fig4, tt);
ttl = title('(c)', "Fontsize", 10);
ttl.Units = 'Normalize';
ttl.Position(1) = -0.4; % use negative values (ie, -0.1) to move further left
ttl.HorizontalAlignment = 'left';
subtitle(tt,"Left",'Fontsize', fontsize, 'FontWeight','bold')
elseif(strats==2)
yticks([])
fig4 = get(ax4_center,'children');
copyobj(fig4, tt);
subtitle(tt,"Center",'Fontsize', fontsize, 'FontWeight','bold')
else
yticks([])
fig4 = get(ax4_right,'children');
copyobj(fig4, tt);
subtitle(tt,"Right",'Fontsize', fontsize, 'FontWeight','bold')
end
h = findall(gca, 'LineStyle', '-');
for i=1:3
h(i).HandleVisibility="off";
end
xlim([-35,35])
ylim([-15,15])
xticks(-30:15:30)
xtickangle(0)
set(gca,"FontSize",9)
lg = legend({"High vis. rel.","Med. vis. rel.","Low vis. rel."});
set(lg,'Box','off')
lg.FontSize = 9;
lg.Location="north";
lg.ItemTokenSize(1) = 10;
end
t5=tiledlayout(T,1,3, 'Padding','none','TileSpacing','compact');
t5.Layout.Tile = 19;
t5.Layout.TileSpan = [1 6];
xlabel(t5, "Stimulus location disparity, {\its}_A– {\its}_V (\circ)", 'FontSize',fontsize)
ylabel(t5,"{\rm \fontsize{10} {Auditory bias (\circ)}}");
for strats=1:3
tt = nexttile(t5);
set(gca,'TickDir','out');
hold on;
if(strats==1)
fig5 = get(ax5_left,'children');
copyobj(fig5, tt);
ttl = title('(d)', "Fontsize", 10);
ttl.Units = 'Normalize';
ttl.Position(1) = -0.4; % use negative values (ie, -0.1) to move further left
ttl.HorizontalAlignment = 'left';
subtitle(tt,"Left",'Fontsize', fontsize, 'FontWeight','bold')
elseif(strats==2)
yticks([])
fig5 = get(ax5_center,'children');
copyobj(fig5, tt);
subtitle(tt,"Center",'Fontsize', fontsize, 'FontWeight','bold')
else
yticks([])
fig5 = get(ax5_right,'children');
copyobj(fig5, tt);
subtitle(tt,"Right",'Fontsize', fontsize, 'FontWeight','bold')
end
h = findall(gca, 'LineStyle', '-');
for i=1:3
h(i).HandleVisibility="off";
end
set(gca,"FontSize",9)
xlim([-35,35])
ylim([-15,15])
xticks(-30:15:30)
xtickangle(0)
lg = legend({"High vis. rel.","Med. vis. rel.","Low vis. rel."});
set(lg,'Box','off')
lg.FontSize = 9;
lg.Location="north";
lg.ItemTokenSize(1) = 10;
switch strats
case 2
lg.Position(1) = 0.46;
lg.Position(2) = 0.7725;
case 3
lg.Position(1) = 0.74;
lg.Position(2) = 0.7725;
end
end
end
end

manuscript_allplots.m at commit b065a19, no license · at the source

Overview

Authors: Shuze Liu¹, Trevor Holland², Wei Ji Ma^2,3, Luigi Acerbi⁴

ORCID iDs: Shuze Liu, Luigi Acerbi

Institutions: Harvard University (United States); Baylor College of Medicine (United States); New York University (United States); University of Helsinki (Finland)

Journal: PLoS computational biology, volume 22, issue 5, article e1014251

Dates: received 13 June 2025; accepted 20 April 2026; published online 8 May 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1371/journal.pcbi.1014251 · PMID 42102227 · PMCID PMC13155690 · OpenAlex W4411170078

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.31234/osf.io/3cwm8_v1

Status: code verified

Categories: human (organism), cognitive (subfield)

Methods: Smoothing, state filtering, decompositions

MeSH: Brain*, Models, Neurological*, Perception*, Bayes Theorem, Causality, Humans, Neural Pathways, Sensation (* major topic)

Journal subjects: Biology and Life Sciences, Neuroscience, Cognitive Science, Cognitive Psychology, Perception, Sensory Perception, Psychology, Social Sciences, Vision, Physical Sciences, Mathematics, Statistics, Statistical Noise, Gaussian Noise, Computer and Information Sciences, Data Management, Data Visualization, Organisms, Eukaryota, Plants, Flowering Plants, Vanilla, Research and Analysis Methods, Separation Processes, Distillation, Sensory Cues, Animals, Vertebrates, Amniotes, Birds, Zoology

Topic: Multisensory perception and integration (Experimental and Cognitive Psychology, Psychology), according to OpenAlex

Funding: Research Council of Finland (356498, 358980)

Citations: not cited yet (Europe PMC); 111 references in the paper

Abstract

The perception of the external world relies on integrating information from multiple sensory modalities. To do this effectively, the brain must determine whether sensory signals come from a common source and, if so, combine them to reduce perceptual uncertainty. While Bayesian observer models have been successful in accounting for multisensory causal inference decisions by humans, they typically rely on simplifying assumptions that may not reflect the true complexity of human perception. In this study, we challenge two assumptions common in Bayesian multisensory perception models: homoskedastic (constant across space) sensory noise and Gaussian priors. We collected an auditory-visual perceptual dataset featuring both unisensory and bisensory tasks, where participants must either provide stimulus location estimates or same-different source judgments. Subsequently, we developed a flexible semiparametric approach that allowed us to infer the sensory noise and prior shapes from participants’ data, and subsequently ‘distill’ them into new model classes through visual inspection of the semiparametrically fitted function shapes. We find that human multisensory perception is best described by an eccentricity-dependent sensory noise that plateaus in the periphery and a prior distribution with a narrow central peak and smoother tails. We also found evidence for auditory range recalibration and increased sensory noise in multisensory conditions, suggesting complex interactions between sensory modalities. These findings deviate substantially from traditional modeling assumptions and highlight the value of data-driven rather than theory-driven modeling assumptions. Overall, our study demonstrates the value of systematically exploring model assumptions in multisensory research and provides a new set of modeling tools for perceptual causal inference.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 18 matches between paragraphs and lines of code.

LSZ2001/Audiovisual-causal-inference

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: b065a1984a7a2bbc9b6df5e06c9caddc927961ea, 20 January 2026

Languages: MATLAB (42), Jupyter (2)

Size: 288 files, 44 scripts

Software Heritage: not archived

Found in: “Data Availability”

Holds: README, 1 notebook

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Statistics and Machine Learning Toolbox (12 files)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

45 files

.ipynb_checkpoints/create_composite_allfits_figures-checkpoint.ipynb, Jupyter, 300 lines
analysis/complete_thetaua_for_ujointfits.m, MATLAB, 30 lines
analysis/data_stratified_to_data.m, MATLAB, 38 lines
analysis/fakedata_ujoint_respdistrvisualization.m, MATLAB, 221 lines
analysis/fit_alldatamodel_parametric_resc.m, MATLAB, 168 lines
analysis/fit_alldatamodel_semiparaminsp_resc.m, MATLAB, 220 lines
analysis/fit_ujointmodel_parametric.m, MATLAB, 155 lines
analysis/fit_ujointmodel_parametric_modelrecov.m, MATLAB, 151 lines
analysis/fit_ujointmodel_semiparam.m, MATLAB, 248 lines
analysis/heterotype_to_sigmafun.m, MATLAB, 16 lines
analysis/manuscript_allfits_respdistrvisual_semiparaminsp_maintext.m, MATLAB, 99 lines
analysis/manuscript_allfits_respdistrvisualization_resc.m, MATLAB, 209 lines
analysis/manuscript_allfits_respdistrvisualization_resc_maintext.m, MATLAB, 144 lines
analysis/manuscript_allfits_respdistrvisualization_semiparaminsp_resc.m, MATLAB, 102 lines
analysis/manuscript_bimodalavfits_visualization_resc.m, MATLAB, 339 lines, 1 match
analysis/manuscript_bimodalavfits_visualization_resc_maintext.m, MATLAB, 244 lines
analysis/manuscript_bimodalcfits_visualization_resc.m, MATLAB, 229 lines
analysis/manuscript_bimodalcfits_visualization_resc_maintext.m, MATLAB, 144 lines
analysis/manuscript_ujoint_respdistrvisualization.m, MATLAB, 182 lines
analysis/manuscript_ujoint_respdistrvisualization_semiparam.m, MATLAB, 98 lines
analysis/manuscript_unimodalfits_visualization.m, MATLAB, 504 lines
analysis/merge_ujoint_badsbounds.m, MATLAB, 33 lines
analysis/nllfun_bav_parametric.m, MATLAB, 541 lines
analysis/nllfun_bav_ubresc_semiparaminsp.m, MATLAB, 420 lines
analysis/nllfun_bc_parametric.m, MATLAB, 367 lines, 1 match
analysis/nllfun_bc_ubresc_semiparaminsp.m, MATLAB, 224 lines
analysis/nllfun_uav_parametric.m, MATLAB, 278 lines
analysis/nllfun_uav_semiparam.m, MATLAB, 265 lines
analysis/nllfun_uav_semiparaminsp.m, MATLAB, 211 lines
analysis/sigmafun_badsbounds_comprehensive.m, MATLAB, 108 lines
analysis/ujointmodel_parametric_modelrecov_datagen.m, MATLAB, 141 lines
create_composite_allfits_figures.ipynb, Jupyter, 300 lines, 1 match
data/parse_data.m, MATLAB, 61 lines
fast_fit_visualize.m, MATLAB, 77 lines
fit_models.m, MATLAB, 192 lines
fit_models_modelrecov_gpc.m, MATLAB, 217 lines
manuscript_allplots.m, MATLAB, 1,561 lines, 7 matches
manuscript_allplots_old.m, MATLAB, 1,296 lines, 6 matches
utils/brewermap.m, MATLAB, 511 lines
utils/cmaes.m, MATLAB, 3,106 lines, 1 match
utils/cmaes_modded.m, MATLAB, 3,070 lines, 1 match
utils/patchline.m, MATLAB, 120 lines
utils/qtrapz.m, MATLAB, 76 lines
utils/trandn.m, MATLAB, 104 lines
README.md, Text, 118 lines

At the source: github.com/LSZ2001/Audiovisual-causal-inference

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
44 scripts, each with its path and the digest of its content;
18 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

All data and analysis code are available at https://github.com/LSZ2001/Audiovisual-causal-inference.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 8 MeSH terms, 1 funder, 96 references.

Cite

This paper

Liu, S., Holland, T., Ma, W. J., & Acerbi, L. (2026). Distilling noise characteristics and prior expectations in multisensory causal inference. PLoS computational biology, 22(5), e1014251. https://doi.org/10.1371/journal.pcbi.1014251

BibTeX

@article{liu2026distilling, author = {Liu, Shuze and Holland, Trevor and Ma, Wei Ji and Acerbi, Luigi}, title = {{Distilling noise characteristics and prior expectations in multisensory causal inference}}, journal = {PLoS computational biology}, year = {2026}, month = may, volume = {22}, number = {5}, pages = {e1014251}, publisher = {PLOS}, issn = {1553-734X}, doi = {10.1371/journal.pcbi.1014251}, url = {https://doi.org/10.1371/journal.pcbi.1014251}, pmid = {42102227}, pmcid = {PMC13155690} }

RIS

TY - JOUR AU - Liu, Shuze AU - Holland, Trevor AU - Ma, Wei Ji AU - Acerbi, Luigi TI - Distilling noise characteristics and prior expectations in multisensory causal inference T2 - PLoS computational biology J2 - PLoS Comput Biol PY - 2026 DA - 2026/05/08 VL - 22 IS - 5 SP - e1014251 SN - 1553-734X PB - PLOS DO - 10.1371/journal.pcbi.1014251 UR - https://doi.org/10.1371/journal.pcbi.1014251 LA - en ER -

CSL-JSON

{ "id": "10.1371/journal.pcbi.1014251", "type": "article-journal", "title": "Distilling noise characteristics and prior expectations in multisensory causal inference", "container-title": "PLoS computational biology", "author": [ { "family": "Liu", "given": "Shuze" }, { "family": "Holland", "given": "Trevor" }, { "family": "Ma", "given": "Wei Ji" }, { "family": "Acerbi", "given": "Luigi" } ], "container-title-short": "PLoS Comput Biol", "volume": "22", "issue": "5", "page": "e1014251", "DOI": "10.1371/journal.pcbi.1014251", "PMID": "42102227", "PMCID": "PMC13155690", "ISSN": "1553-734X", "publisher": "PLOS", "URL": "https://doi.org/10.1371/journal.pcbi.1014251", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 8 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.