To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit c17f87de authored by dandreea's avatar dandreea
Browse files

Initial commit

parent ba8880d8
function get_multiple_conditions(iSubjectArray, doPlotFigures,iResponseModel)
% computes HGF for given subjects and creates parametric modulators for
% concatenated design matrix, plus base regressors for event onsets
%
if nargin < 1
iSubjectArray = setdiff([3:47], [9 14 25 32 33 34 37]);
% 6,7 = noisy; 9
end
if nargin < 2
doPlotFigures = 1;
end
if nargin < 3
iResponseModel = 1;
end
typeDesign = 'ModelBased';
errorSubjects = {};
errorIds = {};
for iSubj = iSubjectArray
%% Load Model and inputs
iD = iSubj;
% try % continuation with new subjects, if error
idDesign = 13;
paths = get_paths_wagad(iSubj,1,idDesign);
if ismac
doFitModel = false;
else
doFitModel = false;
end
addpath(paths.code.model);
input_u = load(fullfile(paths.code.model, 'final_inputs_advice_reward.txt'));% input structure: is this the input structure?
y = [];
%% Load Onsets
% construct output matrix from behavioral log files:
% outputmatrix=[onsets1 onsets2 onsets3 choice onsets_resp RS' inmatrix(:,17)];
outputmatrix = [];
for iRun = 1:2
% try whether run 1 and 2 (male adviser) exist
fileBehav = fullfile(paths.behav, ...
sprintf('%sperblock_IOIO_run%d.mat', paths.idSubjBehav, iRun));
if ~exist(fileBehav)
% we use run 5+6 (female adviser)
fileBehav = fullfile(paths.behav, ...
sprintf('%sperblock_IOIO_run%d.mat', paths.idSubjBehav, iRun+4));
end
load(fileBehav);
trigger = SOC.param(2).scanstart;
fileTrigger = fullfile(paths.behav, sprintf('scanner_trigger_%d.txt', iRun));
save(fileTrigger,'trigger','-ascii','-tabs');
% later runs are offset by duration of previous runs for
% concatentation
offsetRunSeconds = 0 + ...
sum(paths.scanInfo.TR(1:iRun-1).*paths.scanInfo.nVols(1:iRun-1));
outputmatrixSession{iRun} = apply_trigger(fileTrigger, ...
SOC.Session(2).exp_data, offsetRunSeconds);
choice = outputmatrixSession{iRun}(:,4);
wager = outputmatrixSession{iRun}(:,7);
y = [y; choice wager];
outputmatrix = [outputmatrix; outputmatrixSession{iRun}];
end
save(paths.fnBehavMatrix,'outputmatrix','-mat');
%% Run Inversion
for iRsp= 1 % 1:numel(paths.fileResponseModels)
if doFitModel
est=fitModel(y,input_u,'hgf_binary3l_reward_social_config',...
paths.fileResponseModels{iRsp});
save(paths.fnFittedModel{iRsp}, 'est');
else
load(paths.fnFittedModel{iResponseModel},'est','-mat'); % Select the winning model only
hgf_plotTraj_reward_social(est);
end
%% Create Parametric Modulators / Define Conditions
if strcmp(typeDesign,'ModelBased')==1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First: Arbitration (Precision Ratio),
% Advice Prediction (Mu1hat_A), Reward Prediction (Mu1hat_R)
% Time-Locked to Prediction Phase
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pmod(1,1).name = {'Arbitration','Social Weighting',...
'Card Weighting','Arbitration1','Precision_Advice'};
advice_card_space = input_u(:,3);
x_a = est.traj.muhat_a(:,2);
x_r = est.traj.muhat_r(:,2);
transformed_x_r = x_r.^advice_card_space.*(-1.*x_r).^(1-advice_card_space);
sa2hat_a = est.traj.sahat_a(:,2);
sa2hat_r = est.traj.sahat_r(:,2);
px = 1./sa2hat_a;
pc = 1./sa2hat_r;
ze1 = est.p_obs.ze1;
wx = ze1.*px./(ze1.*px + pc);
wc = pc./(ze1.*px + pc);
arbitr = px./(px + pc);
% 1st level precision
px1 = 1./est.traj.sahat_a(:,1);
pc1 = 1./est.traj.sahat_r(:,1);
wx1 = ze1.*px1./(ze1.*px1 + pc1);
Social_weighting = [0; wx.*x_a];
Card_weighting = [0; wc.*transformed_x_r];
Arbitration = [0.5; arbitr];
pmod(1,1).param = {[Arbitration],[Social_weighting],[Card_weighting],[0.5;wx1],[1; px1]}; % Precision (Model-based wager)
pmod(1,1).poly={[1],[1],[1],[1],[1]};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second: Wager (Belief Precision), Belief
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pmod(1,2).name = {'BeliefPrecision','Belief','Alpha','Wager_Amount'}; % Belief Precision
mu2b = wx.*x_a + wc.*transformed_x_r;
b = tapas_sgm(mu2b,1);
pib = 1./(b.*(1-b));
alpha = sgm((pib-4),1);
predict_wager = (2.*alpha -1).*10;
pmod(1,2).param = {[4; pib],[0.5; b],[0.5; alpha],[y(:,2)]}; % Precision (Model-based wager)
pmod(1,2).poly={[1],[1],[1],[1]};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Third: Social and Reward PEs
% Social and Reward Volatility PEs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pmod(1,3).name = {'Epsilon2 Adv','Epsilon2 Cue','Epsilon3_Adv','Epsilon3_Cue'}; % PEs
Epsilon2.Advice = est.traj.sa_a(:,2).*est.traj.da_a(:,1);
Epsilon2.Reward = abs(est.traj.sa_r(:,2).*est.traj.da_r(:,1));
Epsilon3.Advice = est.traj.sa_a(:,3).*est.traj.da_a(:,2);
Epsilon3.Reward = est.traj.sa_r(:,3).*est.traj.da_r(:,2);
pmod(1,3).param = {[0; Epsilon2.Advice],[0; Epsilon2.Reward],[0; Epsilon3.Advice],[0; Epsilon3.Reward]};
pmod(1,3).poly={[1], [1],[1], [1]};
%% Plot
probabilityReward = [ones(25,1)'.*0.9, ones(15,1)'.*0.60, ones(30,1)'.*0.5, ones(25,1)'.*0.4, ...
ones(25,1)'.*0.9,ones(15,1)'.*0.9, ones(25,1)'.*0.1];
probabilityAdvice = [ones(25,1)'.*0.9, ones(15,1)'.*0.90, ones(30,1)'.*0.6, ones(25,1)'.*0.1, ...
ones(25,1)'.*0.6, ones(15,1)'.*0.1, ones(25,1)'.*0.5];
PlotFigureA = 1;
PlotFigureB = 0;
if doPlotFigures
if PlotFigureA
figure;
% Subplots
subplot(5,1,1);
plot(probabilityReward,'b', 'LineWidth', 4);
hold on
plot(probabilityAdvice,'r', 'LineWidth', 4);
ylabel('Probabilities');
subplot(5,1,2);
plot(pmod(1,1).param{1}, 'm', 'LineWidth', 4);
ylabel(pmod(1,1).name{1});
subplot(5,1,3);
plot(pmod(1,1).param{2} , 'r', 'LineWidth', 4);
ylabel(pmod(1,1).name{2});
subplot(5,1,4);
plot(pmod(1,3).param{1}, 'c', 'LineWidth', 4);
ylabel(pmod(1,3).name{1});
subplot(5,1,5);
plot(pmod(1,3).param{2}, 'b', 'LineWidth', 4);
ylabel(pmod(1,3).name{2});
hold on;
xlabel('Trial number');
subplot(5,1,1);
hold on;
title([sprintf('cscore = %d', SOC.cscore), ' with \zeta=', ...
num2str(est.p_obs.ze1), ...
' for subject ', paths.idSubjBehav], ...
'FontWeight', 'bold');
elseif PlotFigureB
figure;
% Subplots
subplot(4,1,1);
plot(pmod(1,2).param{1}, 'm', 'LineWidth', 4);
hold on;
plot(ones(170,1).*0.25,'k','LineWidth', 1,'LineStyle','-.');
subplot(4,1,2);
plot(pmod(1,2).param{2} , 'r', 'LineWidth', 4);
hold on;
plot(ones(170,1).*0.25,'k','LineWidth', 1,'LineStyle','-.');
subplot(4,1,3);
plot(pmod(1,2).param{3}, 'c', 'LineWidth', 4);
hold on;
plot(ones(170,1).*1,'k','LineWidth', 1,'LineStyle','-.');
subplot(4,1,4);
plot(pmod(1,2).param{4}, 'b', 'LineWidth', 4);
hold on;
plot(predict_wager , 'r', 'LineWidth', 4);
plot(ones(170,1).*0.5,'k','LineWidth', 1,'LineStyle','-.');
xlabel('Trial number');
subplot(4,1,1);
hold on;
title([sprintf('cscore = %d', SOC.cscore), ' with \zeta=', ...
num2str(est.p_obs.ze1), ...
' for subject ', paths.idSubjBehav], ...
'FontWeight', 'bold');
else
figure;
% Subplots
subplot(4,1,1);
plot(pmod(1,3).param{1}, 'm', 'LineWidth', 4);
ylabel(pmod(1,3).name{1});
subplot(4,1,2);
plot(pmod(1,3).param{2} , 'r', 'LineWidth', 4);
ylabel(pmod(1,3).name{2});
subplot(4,1,3);
plot(pmod(1,3).param{3}, 'c', 'LineWidth', 4);
ylabel(pmod(1,3).name{3});
subplot(4,1,4);
plot(pmod(1,3).param{4}, 'b', 'LineWidth', 4);
ylabel(pmod(1,3).name{1});
hold on;
xlabel('Trial number');
subplot(4,1,1);
hold on;
title([sprintf('cscore = %d', SOC.cscore), ' with \zeta=', ...
num2str(est.p_obs.ze1), ...
' for subject ', paths.idSubjBehav], ...
'FontWeight', 'bold');
end
end
onsets{1} = outputmatrix(:,1);
onsets{2} = outputmatrix(:,2);
onsets{3} = outputmatrix(:,3);
% Switch off orthogonalization for each condition separately
orth{1} = 0;
orth{2} = 0;
orth{3} = 0;
names={'Advice','Wager','Outcome'};
durations{1} = 2;
durations{2} = 0;
durations{3} = 0;
save(paths.fnMultipleConditions, 'onsets', 'names', 'durations', 'names', 'pmod', 'orth', '-mat');
elseif strcmp(typeDesign,'ModelFree')==1
AdviceCodingUnstable=[zeros(25,1)' zeros(15,1)' ones(30,1)' zeros(25,1)' ones(25,1)' zeros(15,1)' ones(25,1)'];
RewardCodingUnstable=[zeros(25,1)' ones(15,1)' ones(30,1)' ones(25,1)' zeros(25,1)' zeros(15,1)' zeros(25,1)'];
AdviceCodingStable = [ones(25,1)' ones(15,1)' zeros(30,1)' ones(25,1)' zeros(25,1)' ones(15,1)' zeros(25,1)'];
RewardCodingStable = [ones(25,1)' zeros(15,1)' zeros(30,1)' zeros(25,1)' ones(25,1)' ones(15,1)' ones(25,1)'];
names={'RewardStable','RewardUnstable','AdviceStable','AdviceUnstable'};
% onsets{1} = outputmatrix(RewardCodingStable==1,1);
% onsets{2} = outputmatrix(RewardCodingUnstable==1,1);
% onsets{3} = outputmatrix(AdviceCodingStable==1,1);
% onsets{4} = outputmatrix(AdviceCodingUnstable==1,1);
% durations{1} = 0;
% durations{2} = 0;
% durations{3} = 0;
% durations{4} = 0;
% save(paths.fnModelFreePredictionConditions, 'onsets', 'names', 'durations', '-mat');
%
% clear onsets;
% clear durations;
onsets{1} = outputmatrix(RewardCodingStable==1,2);
onsets{2} = outputmatrix(RewardCodingUnstable==1,2);
onsets{3} = outputmatrix(AdviceCodingStable==1,2);
onsets{4} = outputmatrix(AdviceCodingUnstable==1,2);
durations{1} = 0;
durations{2} = 0;
durations{3} = 0;
durations{4} = 0;
save(paths.fnModelFreeWagerConditions, 'onsets', 'names', 'durations', '-mat');
end
end
% catch err
% errorSubjects{end+1,1}.id = iD;
% errorSubjects{end}.error = err;
% errorIds{end+1} = iD;
% end
end
\ No newline at end of file
function loop_analyze_behaviour_local(options)
% Loops over subjects in WAGADAP study with a local loop and executes all analysis steps %%
%
% USAGE
% loop_analyze_subject_local(options);
%
% IN
% options (subject-independent) analysis pipeline options,
% retrieve via options wagadap_options
%
for idCell = options.subjectIDs
id = char(idCell);
wagadap_analyze_subject(id,options);
end
function [behaviour_variables] = wagad_extract_behaviour(details, options,fileBehav)
% Get responses, performance accuracy, cumulative score
% Behaviour statistics collected here
% iValid is the row of trials where a response was made; 1 = response made;
% 0 = miss
% logical array
sub = details.dirSubject;
input_u = load(fullfile(options.task.inputs));% input structure
if ~exist(fileBehav, 'file')
warning('Behavioral logfile for subject %s not found', sub)
else
subjectData = load(fileBehav);
[outputmatrixSession,iValid,Probe] = wagadap_extract_from_rawdata(subjectData.SOC.Session(2).exp_data);
choice = outputmatrixSession(:,1);
wager = outputmatrixSession(:,2);
RT = outputmatrixSession(:,3);
y = [choice wager];
probeSelection = Probe;
y(1,:) = []; % Remove the first trial since it is not cued by the card
iValid(1,:) = []; % Remove the first trial since it is not cued by the card
RT(1,:) = [];
probeCategories = [probeSelection(1,1) probeSelection(14,1), ...
probeSelection(49,1) probeSelection(73,1),...
probeSelection(110,1)];
probes = [];
for iProbe = 1:numel(probeCategories)
switch probeCategories(iProbe)
case 0
probeValue = 0.5;
case 1
probeValue = 0.8;
case 2
probeValue = 0.2;
case 3
probeValue = 0.5;
end
probeValue(iProbe) = probeValue;
probes = [probes; probeValue(iProbe)];
end
end
iValid = logical(iValid);
totalTrials = size(y(iValid,1),1);
congruenceBehaviourAdvice = double(y(iValid,1)==input_u(iValid,1));
overall_perf_acc = sum(congruenceBehaviourAdvice)./totalTrials; % percentage of correct trials, disregarding misses
overall_wager = y(end,2);
take_adv_overall = sum(y(iValid,1))./totalTrials;
temp1 = (congruenceBehaviourAdvice).*2;
cScore = sum((temp1+(ones(size(y(iValid,1),1),1).*-1)));
AccuracyStableCard = sum(congruenceBehaviourAdvice.*options.task.design.stableCardPhase(iValid))./sum(options.task.design.stableCardPhase(iValid));
AccuracyVolatileCard = sum(congruenceBehaviourAdvice.*options.task.design.volatileCardPhase(iValid))./sum(options.task.design.volatileCardPhase(iValid));
AccuracyStableAdvice = sum(congruenceBehaviourAdvice.*options.task.design.stableAdvicePhase(iValid))./sum(options.task.design.stableAdvicePhase(iValid));
AccuracyVolatileAdvice= sum(congruenceBehaviourAdvice.*options.task.design.volatileAdvicePhase(iValid))./sum(options.task.design.volatileAdvicePhase(iValid));
AdviceStableCard = sum(y(iValid,1).*options.task.design.stableCardPhase(iValid))./sum(options.task.design.stableCardPhase(iValid));
AdviceVolatileCard = sum(y(iValid,1).*options.task.design.volatileCardPhase(iValid))./sum(options.task.design.volatileCardPhase(iValid));
AdviceStableAdvice = sum(y(iValid,1).*options.task.design.stableAdvicePhase(iValid))./sum(options.task.design.stableAdvicePhase(iValid));
AdviceVolatileAdvice= sum(y(iValid,1).*options.task.design.volatileAdvicePhase(iValid))./sum(options.task.design.volatileAdvicePhase(iValid));
WagerStableCard = sum(y(iValid,2).*options.task.design.stableCardPhase(iValid));
WagerVolatileCard = sum(y(iValid,2).*options.task.design.volatileCardPhase(iValid));
WagerStableAdvice = sum(y(iValid,2).*options.task.design.stableAdvicePhase(iValid));
WagerVolatileAdvice= sum(y(iValid,2).*options.task.design.volatileAdvicePhase(iValid));
RTStableCard = mean(RT(iValid,1).*options.task.design.stableCardPhase(iValid));
RTVolatileCard = mean(RT(iValid,1).*options.task.design.volatileCardPhase(iValid));
RTStableAdvice = mean(RT(iValid,1).*options.task.design.stableAdvicePhase(iValid));
RTVolatileAdvice= mean(RT(iValid,1).*options.task.design.volatileAdvicePhase(iValid));
behaviour_variables = [];
behaviour_variables.overall_perf_acc = overall_perf_acc;
behaviour_variables.overall_wager = overall_wager;
behaviour_variables.cScore = cScore;
behaviour_variables.take_adv_overall = take_adv_overall;
behaviour_variables.AccuracyStableCard = AccuracyStableCard;
behaviour_variables.AccuracyVolatileCard = AccuracyVolatileCard;
behaviour_variables.AccuracyStableAdvice = AccuracyStableAdvice;
behaviour_variables.AccuracyVolatileAdvice = AccuracyVolatileAdvice;
behaviour_variables.AdviceStableCard = AdviceStableCard;
behaviour_variables.AdviceVolatileCard = AdviceVolatileCard;
behaviour_variables.AdviceStableAdvice = AdviceStableAdvice;
behaviour_variables.AdviceVolatileAdvice = AdviceVolatileAdvice;
behaviour_variables.WagerStableCard = WagerStableCard;
behaviour_variables.WagerVolatileCard = WagerVolatileCard;
behaviour_variables.WagerStableAdvice = WagerStableAdvice;
behaviour_variables.WagerVolatileAdvice = WagerVolatileAdvice;
behaviour_variables.RTStableCard = RTStableCard;
behaviour_variables.RTVolatileCard = RTVolatileCard;
behaviour_variables.RTStableAdvice = RTStableAdvice;
behaviour_variables.RTVolatileAdvice = RTVolatileAdvice;
behaviour_variables.probes = probes;
end
\ No newline at end of file
function wagadap_analyze_subject(id, options)
doinvertSubject = ismember('inversion', options.pipe.executeStepsPerSubject);
doplotTrajectory = ismember('plotting', options.pipe.executeStepsPerSubject);
dobehavVariables = ismember('behaviour', options.pipe.executeStepsPerSubject);
% Inverts the model for a given subject
if doinvertSubject
wagadap_invert_subject(id, options);
end
% Plots trajectories of a given subject
if doplotTrajectory
wagadap_plot_subject(id, options);
end
% Plots trajectories of a given subject
if dobehavVariables
wagadap_get_behaviour_data(id, options);
end
end
function [outputmatrix,isValidTrial,probe_selection]=wagadap_extract_from_rawdata(InputFile)
%read inputfile1 into a matrix
inmatrix = InputFile;
isNotValid = inmatrix(:,8)<0;
RT = inmatrix(:,10)./1000;
RT(find(isNotValid))=NaN;
isValidTrial = inmatrix(:,8)>=0;
adviceBlue = mod(inmatrix(:,4),2);
resp = inmatrix(:,8);
respBlue = mod(resp,2); % blue = 1, green = 2
choice_congr = (adviceBlue == respBlue);
choice = double(choice_congr);
probe_selection = inmatrix(:,6);
wager_amount = inmatrix(:,17);
outputmatrix=[choice wager_amount RT];
function wagadap_get_behaviour_data(id,options)
% extracts behaviour variables, computes HGF for given subjects for
% concatenated design matrix, plus base regressors for event onsets
%
details = wagadap_subjects(id, options);
fileBehav = details.behav.fileRawBehav;
[behaviour_variables] = wagad_extract_behaviour(details,options,fileBehav);
save(details.behav.fileBehav,'behaviour_variables','-mat');
function wagadap_invert_subject(id, options)
% IN
% id subject id string, only number (e.g. '0002')
% options general analysis options%
% options = wagadap_options
details = wagadap_subjects(id, options);
fileBehav = details.behav.fileRawBehav;
input_u = load(fullfile(options.task.inputs));% input structure
perceptualModels = options.model.perceptualModels;
responseModels = options.model.responseModels';
iCombPercResp = zeros(6,2);
iCombPercResp(1:3,1) = 1;
iCombPercResp(4:6,1) = 2;
iCombPercResp(1,2) = 1;
iCombPercResp(2,2) = 2;
iCombPercResp(3,2) = 3;
iCombPercResp(4,2) = 1;
iCombPercResp(5,2) = 2;
iCombPercResp(6,2) = 3;
nModels = size(iCombPercResp,1);
for iModel = 1:nModels
subjectData = load(fileBehav);
[outputmatrixSession,~,~] = wagadap_extract_from_rawdata(subjectData.SOC.Session(2).exp_data);
choice = outputmatrixSession(:,1);
wager = outputmatrixSession(:,2);
y_responses = [choice wager];
if isempty(input_u) % cue and cue advice are not subject-dependent, have to check inputs
error('Behavioral data for subject %s not found', sub);
end
est_wagadap=tapas_fitModel(y_responses,input_u,perceptualModels{iCombPercResp(iModel,1)},...
responseModels{iCombPercResp(iModel,2)});
save(fullfile(details.behav.pathResults,[perceptualModels{iCombPercResp(iModel,1)}, ...
responseModels{iCombPercResp(iModel,2)},'.mat']), 'est_wagadap','-mat');
end
end
function wagadap_plot_subject(id, options)
details = wagadap_subjects(id, options);
perceptualModels = options.model.perceptualModels;
responseModels = options.model.responseModels';
iCombPercResp = zeros(6,2);
iCombPercResp(1:3,1) = 1;
iCombPercResp(4:6,1) = 2;
iCombPercResp(1,2) = 1;
iCombPercResp(2,2) = 2;
iCombPercResp(3,2) = 3;
iCombPercResp(4,2) = 1;
iCombPercResp(5,2) = 2;
iCombPercResp(6,2) = 3;
nModels = size(iCombPercResp,1);
for iModel=1:nModels
load(fullfile(details.behav.pathResults,[perceptualModels{iCombPercResp(iModel,1)}, ...
responseModels{iCombPercResp(iModel,2)},'.mat']), 'est_wagadap','-mat');
hgf_plotTraj_reward_social(est_wagadap);
tapas_fit_plotCorr(est_wagadap);
end
end
\ No newline at end of file
1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 0 1 0 0 0 0 0 1 1 1 1 1 0 1 0 1 1 1 0 1 0 1 1 1 0 1 0 0 1 0 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 0 1 0 1 1 0 0 0 1 1 0 0 1 0 0 1 0 0 1 0 1 1 1 1 0 1 1 0 0 0 0 1 1 0 0 1 1 0 1 0 1 0 0 1 1 0 0 0 0 1 1 1 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 0 0 0 0 1 1 0 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 0 0 1 1 0 0 1 0 0 1 0 0 1 0 0 1 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 1 1 0 1 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 1 1 1 1 1 1 1 1 0 1 0 0 1 0 0 1 0 1 1 1 1 1 0 1 0 1 1 1 1 0 1 0 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 0 1 0 0 1 0 0 1 1 1 0 0 0 0 1 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 0 0 1 0 0 1 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 0 1 1 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 1 1
\ No newline at end of file
function [predicted_wager] = calculate_predicted_wager(est,paths)
%% Inferred states
mu1hat_a = est.traj.muhat_a(:,1);
mu1hat_r = est.traj.muhat_r(:,1);
mu2hat_a = est.traj.muhat_a(:,2);
mu2hat_r = est.traj.muhat_r(:,2);
sa2hat_r = est.traj.sahat_r(:,2);
sa2hat_a = est.traj.sahat_a(:,2);
mu3hat_r = est.traj.muhat_r(:,3);
mu3hat_a = est.traj.muhat_a(:,3);
ze = est.p_obs.ze;