Identify the Best Model
*************

Lesson
================

.. article-info::
    :avatar: UCLA_Suit.jpg
    :avatar-link: https://www.decisionneurosciencelab.com/elijahgalvan
    :author: Elijah Galvan
    :date: September 1, 2023
    :read-time: 15 min read
    :class-container: sd-p-2 sd-outline-muted sd-rounded-1

Goal During this Stage
---------------

To determine which model is best. 

.. Note::

    Now if you've been following along, you've probably realized that we've only talked so far about a single model - I'll briefly explain why. 
    We use Social Utility Models to encapsulate individual differences in value based preferences and exisiting research on this topic essentially shows that **people are usually different in most conceivable ways**. 
    Thus, we always try to design experiments around trying to see this individual differences in behavior and our "favored" model - the one we think will best fit the data - is the one which captures all of these differences.
    Since the favored model must capture all meaningful potential differences, it is therefore the most complex model that we will test.

    So why are we reiterating this now? 
    Well, the answer is that any *other types* of models you test on value-based decision-making behavior in non-forced, non-probabilistic, social choice tasks are almost always going to perform much worse. 
    These tasks are designed to elicit consistent, well-defined preferences with a very high signal to noise ratio: therefore, models capturing other aspects of the decision-making process other than these preferences are theoretically and pragmatically inferior in these scenarios. 
    All of this to say, we are going to test this favored model - which is the most complex model - against models which are simplifications of this model.
    Thus, we can use the same :bdg-secondary:`Construct Value` computations, removing one of these :bdg-secondary:`Construct Value` computations from the utility equation and a :bdg-success:`Free Parameter` per model.

How to Achieve this Goal
------------

.. dropdown:: Create :bdg-warning:`Utility` Equations for our Alternative Models

    Generally, you should first enumerate all of the reasonable potential alternative models. 
    This means that if you have a model with 3 different :bdg-secondary:`Construct Value` terms (as most models that we will talk about indeed do) then the alternative model set would include models with only 2 or 1 :bdg-secondary:`Construct Value` terms.
    Are any models not worth testing? 
    Then don't include them. 
    For every model we think is worth testing, let's talk about how you're going to go about testing them.

    Any models that we want to test with one :bdg-secondary:`Construct Value` term are really, really simple - :bdg-warning:`Utility` = :bdg-secondary:`Construct Value`. 
    Notice that there are no :bdg-success:`Free Parameters` in this Equation - since only one :bdg-secondary:`Construct Value` determines :bdg-warning:`Utility`, there is no way that it can capture individual differences. 

    The models that we want to test with two :bdg-secondary:`Construct Value` terms are also pretty simple - we want to capture how people prioritize one norm over another. 
    Thus, we express this as :bdg-warning:`Utility` = (:bdg-secondary:`Construct Value 1` x :bdg-success:`Free Parameter`) + (:bdg-secondary:`Construct Value 2` x (1 - :bdg-success:`Free Parameter`)).

.. dropdown:: Preallocating and Defining Functions 

    .. tab-set::

        .. tab-item:: Plain English

            We're going to need several items for this section:

            1. :bdg-warning:`Utility` Functions for each each Model
            2. Objective Functions for each :bdg-warning:`Utility` Equation with :bdg-success:`Free Parameters`
            3. New Optimizer Inputs for each :bdg-warning:`Utility` Equation with :bdg-success:`Free Parameters`
            4. Outputs for the Alternative Models
        
        .. tab-item:: R

            :: 

                utility_alt1 = function(construct1){
                    return(utility)
                }
                utility_alt2 = function(construct1, construct2, parameter1){
                    return(utility)
                }

                #alternative model with 1 construct has 0 free parameters: doesn't need an objective function

                obj_function_alt2 = function(param, decisions, method = "OLS") {
                    Parameter1 = param[1]
                    
                    predicted_utility = vector('numeric', length(trialList[,1]))
                    chosen = decisions + 1
                    for (k in 1:length(trialList[,1])){
                        IV = trialList[k, 1]
                        Constant = trialList[k, 2]
                        Choices = seq(0, (I * M), 1)
                        
                        Utility = vector('numeric', length(Choices))
                        for (n in 1:length(Choices)){
                        Utility[n] = utility_alt2(Parameter1, construct1(IV, Constant, Choices[n]))
                        }
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[chosen[k]]
                    }
                    if (method == "OLS"){
                        return(sum((predicted_utility - observed_utility)**2))
                    } else if (method == "MLE"){
                        return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
                    }
                }

                initial_param_alt2 = #something
                lower_bound_alt2 = #something
                upper_bound_alt2 = #something

                altSubjectData = data.frame()

        .. tab-item:: MatLab

                ::

                    function utility = utility_alt1(construct1)
                        % Your utility calculation for utility_alt1 goes here
                        utility = 0; % Placeholder value, replace with actual calculation
                    end

                    function utility = utility_alt2(construct1, construct2, parameter1)
                        % Your utility calculation for utility_alt2 goes here
                        utility = 0; % Placeholder value, replace with actual calculation
                    end

                    function obj_value = obj_function_alt2(param, decisions, method)
                        Parameter1 = param(1);
                        
                        trialList = []; % Define trialList here
                        
                        predicted_utility = zeros(size(trialList, 1), 1);
                        chosen = decisions + 1;
                        
                        for k = 1:size(trialList, 1)
                            IV = trialList(k, 1);
                            Constant = trialList(k, 2);
                            Choices = 0:(I * M); % Define I and M
                            
                            Utility = zeros(size(Choices));
                            
                            for n = 1:length(Choices)
                                Utility(n) = utility_alt2(Parameter1, construct1(IV, Constant, Choices(n)));
                            end
                            
                            predicted_utility(k) = max(Utility);
                            observed_utility(k) = Utility(chosen(k));
                        end
                        
                        if strcmp(method, 'OLS')
                            obj_value = sum((predicted_utility - observed_utility).^2);
                        elseif strcmp(method, 'MLE')
                            obj_value = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
                        end
                    end

                    % Define initial_param_alt2, lower_bound_alt2, and upper_bound_alt2 here

                    altSubjectData = table(); % Create an empty table for altSubjectData

        .. tab-item:: Python

                ::

                    from scipy.stats import norm

                    def utility_alt1(construct1):
                        # Your utility calculation for utility_alt1 goes here
                        utility = 0  # Placeholder value, replace with actual calculation
                        return utility

                    def utility_alt2(construct1, construct2, parameter1):
                        # Your utility calculation for utility_alt2 goes here
                        utility = 0  # Placeholder value, replace with actual calculation
                        return utility

                    def obj_function_alt2(param, decisions, method="OLS"):
                        Parameter1 = param[0]  # MATLAB-style indexing
                        
                        trialList = np.array([])  # Define trialList here
                        
                        predicted_utility = np.zeros(len(trialList))
                        chosen = decisions + 1
                        
                        for k in range(len(trialList)):
                            IV = trialList[k, 0]
                            Constant = trialList[k, 1]
                            Choices = np.arange(0, I * M + 1)  # Define I and M
                            
                            Utility = np.zeros(len(Choices))
                            
                            for n in range(len(Choices)):
                                Utility[n] = utility_alt2(Parameter1, construct1(IV, Constant, Choices[n]))
                            
                            predicted_utility[k] = np.max(Utility)
                            observed_utility[k] = Utility[chosen[k]]
                        
                        if method == "OLS":
                            obj_value = np.sum((predicted_utility - observed_utility)**2)
                        elif method == "MLE":
                            obj_value = -np.sum(np.log(norm.pdf(observed_utility, predicted_utility, sd)))
                        
                        return obj_value

                    # Define initial_param_alt2, lower_bound_alt2, and upper_bound_alt2 here

                    altSubjectData = pd.DataFrame()  # Create an empty DataFrame for altSubjectData


.. dropdown:: Recover :bdg-success:`Free Parameters` for Each Alternative Model, Per :bdg-success:`Subject`

    .. tab-set::

        .. tab-item:: Plain English

           For the alternative models with :bdg-success:`Free Parameters`, we'll need to recover these :bdg-success:`Free Parameters` in order to generate model predictions. 
           Let's do that quickly in the same way that we did for the other model, leaving a demand to subsequently determine model predictions.
           
           .. Note::
            
                Models with only one :bdg-success:`Free Parameter` may require a different optimzer than models with multiple :bdg-success:`Free Parameters`. 
                However, you won't have to change anything about how your objective function works or anything like that so don't worry!
        
        .. tab-item:: R

            :: 

                for (i in 1:length(included_subjects)){
                    datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
                    df = read.csv2(datafile)
                    reorder = df$trialsTask.thisIndex + 1

                    result_alt2 = optim(obj_function_alt2, par = initial_param_alt2, lower = lower_bound_alt2, upper = upper_bound_alt2, decisions = df$Decisions)

                    # Determine Predictions
                }

        .. tab-item:: MatLab

                    ::

                        for i = 1:length(included_subjects)
                            datafile = strcat(parentfolder, included_subjects{i}, restoffilepath); % produces a character vector 'parentfolder/included_subjects{i}**.filetype'
                            df = readtable(datafile);
                            reorder = df.trialsTask.thisIndex + 1;

                            result_alt2 = fmincon(@obj_function_alt2, initial_param_alt2, [], [], [], [], lower_bound_alt2, upper_bound_alt2, [], optimset('Display', 'off'), df.Decisions);

                            % Determine Predictions
                        end


        .. tab-item:: Python

                    ::

                        for i in range(len(included_subjects)):
                            datafile = parentfolder + included_subjects[i] + restoffilepath  # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
                            df = pd.read_csv(datafile, sep='\t')
                            reorder = df['trialsTask.thisIndex'] + 1

                            result_alt2 = optimize.minimize(obj_function_alt2, initial_param_alt2, bounds=list(zip(lower_bound_alt2, upper_bound_alt2)), args=(df['Decisions'],))

                            # Determine Predictions


.. dropdown:: Determine Predicted :bdg-danger:`Decisions` for Each Alternative Model, Per :bdg-success:`Subject`

    .. tab-set::

            .. tab-item:: Plain English

                Now, we are going to answer the Determine Predictions demand placed on us.
                We have found the :bdg-success:`Subject`'s :bdg-success:`Free Parameters` so we need to specifically know what it is that our model predicts that they will do.
                In the previous step, we could have cut a corner and gotten the predictions from the closest point we simulated data for. 
                In all likelihood, the model predictions would be indistinguishable from these, but for the sake of being punctual let's get these predictions in the same way we did with our favored model!

            .. tab-item:: R

                ::

                    for (i in 1:length(included_subjects)){
                        datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
                        df = read.csv2(datafile)
                        reorder = df$trialsTask.thisIndex + 1

                        result_alt2 = optim(obj_function_alt2, par = initial_param_alt2, lower = lower_bound_alt2, upper = upper_bound_alt2, decisions = df$Decisions)

                        #Just Added
                        
                        df$PredictionAlt1 = vector('numeric')
                        df$PredictionAlt2 = df$PredictionAlt1
                        for (k in 1:length(df$Decisions)){
                            UtilityAlt1 = vector('numeric', length(Choices))
                            UtilityAlt2 = vector('numeric', length(Choices))
                            for (n in 1:length(Choices)){
                                UtilityAlt1[n] = utility_alt1(construct1 = construct1(df$IV[k], df$Constant[k], Choices[n]))
                                UtilityAlt2[n] = utility_alt2(parameter1 = result_alt2$par[1],
                                                              construct1 = construct1(df$IV[k], df$Constant[k], Choices[n]),
                                                              construct2 = construct2(df$IV[k], df$Constant[k], Choices[n]))
                            }
                            correct_choice_alt1 = which(UtilityAlt1 == max(UtilityAlt1))
                            correct_choice_alt2 = which(UtilityAlt2 == max(UtilityAlt2))
                            if (length(correct_choice) > 1){
                                correct_choice = correct_choice[sample(correct_choice, 1)]
                            }
                            df$PredictionAlt1[k] = Choices[correct_choice_alt1]
                            df$PredictionAlt2[k] = Choices[correct_choice_alt2]
                        }

                        model_NLL_Alt1 = -2 * log(sum(dnorm(df$Decision, mean = df$PredictionAlt1)))
                        model_SS_Alt1 = sum((df$Decision - df$PredictionAlt1)**2)
                        model_NLL_Alt2 = -2 * log(sum(dnorm(df$Decision, mean = df$PredictionAlt2)))
                        model_SS_Alt2 = sum((df$Decision - df$PredictionAlt2)**2)

                        altSubjectData[i, 1:6] = c(included_subjects[i], result_alt2$par[1], model_NLL_Alt1, model_SS_Alt1, model_NLL_Alt2, model_SS_Alt2) 
                    }
                    colnames(altSubjectData) = c('SubjectID', 'Parameter1_Alt2', 'alt1_modelNLL', 'alt1_modelSS', 'alt2_modelNLL', 'alt2_modelSS')

            .. tab-item:: MatLab

                ::

                    for i = 1:length(included_subjects)
                        datafile = strcat(parentfolder, included_subjects{i}, restoffilepath);
                        df = readtable(datafile);
                        reorder = df.trialsTask.thisIndex + 1;

                        result_alt2 = fmincon(@(x) obj_function_alt2(x, df.Decisions), initial_param_alt2, [], [], [], [], lower_bound_alt2, upper_bound_alt2);

                        % Just Added
                        df.PredictionAlt1 = zeros(size(df.Decisions));
                        df.PredictionAlt2 = df.PredictionAlt1;
                        
                        for k = 1:length(df.Decisions)
                            UtilityAlt1 = zeros(size(Choices));
                            UtilityAlt2 = zeros(size(Choices));
                            
                            for n = 1:length(Choices)
                                UtilityAlt1(n) = utility_alt1(construct1(df.IV(k), df.Constant(k), Choices(n)));
                                UtilityAlt2(n) = utility_alt2(result_alt2(1), construct1(df.IV(k), df.Constant(k), Choices(n)), construct2(df.IV(k), df.Constant(k), Choices(n)));
                            end
                            
                            [~, correct_choice_alt1] = max(UtilityAlt1);
                            [~, correct_choice_alt2] = max(UtilityAlt2);
                            
                            if length(correct_choice_alt1) > 1
                                correct_choice_alt1 = correct_choice_alt1(randi(length(correct_choice_alt1)));
                            end
                            
                            df.PredictionAlt1(k) = Choices(correct_choice_alt1);
                            df.PredictionAlt2(k) = Choices(correct_choice_alt2);
                        end

                        model_NLL_Alt1 = -2 * sum(log(normpdf(df.Decisions, df.PredictionAlt1)));
                        model_SS_Alt1 = sum((df.Decisions - df.PredictionAlt1).^2);
                        model_NLL_Alt2 = -2 * sum(log(normpdf(df.Decisions, df.PredictionAlt2)));
                        model_SS_Alt2 = sum((df.Decisions - df.PredictionAlt2).^2);

                        altSubjectData(i, 1:6) = [included_subjects{i}, result_alt2(1), model_NLL_Alt1, model_SS_Alt1, model_NLL_Alt2, model_SS_Alt2];
                    end

                    altSubjectData.Properties.VariableNames = {'SubjectID', 'Parameter1_Alt2', 'alt1_modelNLL', 'alt1_modelSS', 'alt2_modelNLL', 'alt2_modelSS'};


            .. tab-item:: Python
                
                ::

                    for i in range(len(included_subjects)):
                        datafile = parentfolder + included_subjects[i] + restoffilepath
                        df = pd.read_csv(datafile, delimiter=',')
                        reorder = df['trialsTask.thisIndex'] + 1

                        result_alt2 = minimize(lambda x: obj_function_alt2(x, df['Decisions']), initial_param_alt2, bounds=list(zip(lower_bound_alt2, upper_bound_alt2)))

                        # Just Added
                        df['PredictionAlt1'] = np.zeros(len(df['Decisions']))
                        df['PredictionAlt2'] = df['PredictionAlt1']

                        for k in range(len(df['Decisions'])):
                            UtilityAlt1 = np.zeros(len(Choices))
                            UtilityAlt2 = np.zeros(len(Choices))

                            for n in range(len(Choices)):
                                UtilityAlt1[n] = utility_alt1(construct1(df['IV'][k], df['Constant'][k], Choices[n]))
                                UtilityAlt2[n] = utility_alt2(result_alt2.x[0], construct1(df['IV'][k], df['Constant'][k], Choices[n]), construct2(df['IV'][k], df['Constant'][k], Choices[n]))

                            correct_choice_alt1 = np.argmax(UtilityAlt1)
                            correct_choice_alt2 = np.argmax(UtilityAlt2)

                            if len(correct_choice_alt1) > 1:
                                correct_choice_alt1 = np.random.choice(correct_choice_alt1, 1)

                            df['PredictionAlt1'][k] = Choices[correct_choice_alt1]
                            df['PredictionAlt2'][k] = Choices[correct_choice_alt2]

                        model_NLL_Alt1 = -2 * np.sum(np.log(np.random.normal(df['Decisions'], df['PredictionAlt1'])))
                        model_SS_Alt1 = np.sum((df['Decisions'] - df['PredictionAlt1'])**2)
                        model_NLL_Alt2 = -2 * np.sum(np.log(np.random.normal(df['Decisions'], df['PredictionAlt2'])))
                        model_SS_Alt2 = np.sum((df['Decisions'] - df['PredictionAlt2'])**2)

                        altSubjectData[i, 1:6] = [included_subjects[i], result_alt2.x[0], model_NLL_Alt1, model_SS_Alt1, model_NLL_Alt2, model_SS_Alt2]

                    altSubjectData.columns = ['SubjectID', 'Parameter1_Alt2', 'alt1_modelNLL', 'alt1_modelSS', 'alt2_modelNLL', 'alt2_modelSS']


.. dropdown:: Compute Model Fit Index for Each :bdg-success:`Subject`, for Each Alternative Model

    .. tab-set::

        .. tab-item:: Plain English

           Now that we have the model error - either the sum of squared errors or the negative log likelihood of the real :bdg-success:`Subjects` :bdg-danger:`Decisions` versus the model's predicted :bdg-danger:`Decisions`.
        
        .. tab-item:: R

            :: 

                altSubjectData$modelAlt1AIC = N * log(altSubjectData$modelSS_Alt1/N) + 2*0
                altSubjectData$modelAlt2AIC = N * log(altSubjectData$modelSS_Alt2/N) + 2*1

        .. tab-item:: MatLab

            ::

                altSubjectData.modelAlt1AIC = N * log(altSubjectData.modelSS_Alt1/N) + 2*0;
                altSubjectData.modelAlt2AIC = N * log(altSubjectData.modelSS_Alt2/N) + 2*1;


        .. tab-item:: Python

            ::

                altSubjectData['modelAlt1AIC'] = N * log(altSubjectData['modelSS_Alt1']/N) + 2*0
                altSubjectData['modelAlt2AIC'] = N * log(altSubjectData['modelSS_Alt2']/N) + 2*1


.. dropdown:: Compare Model Performance

    .. tab-set::

        .. tab-item:: Plain English

            Now we simply want to identify which model is best. 
            Thus, we're going to create a vector with the MFI for each model averaged across all :bdg-success:`Subjects` and select the model with the lowest MFI. 
            This approach tells us which model provides the best average fit for :bdg-success:`Subjects` in our sample. 
            Importantly, if any :bdg-success:`Subjects` data are fully explained by the model (i.e. observed :bdg-danger:`Decisions` always equal :bdg-danger:`Decisions` predicted by the model) then these :bdg-success:`Subjects` must be excluded from your analysis since their MFIs are negative infinity. 

            Another approach would be to compute MFIs for the entire dataset - this approach does not require that you exclude :bdg-success:`Subjects` from analysis. 
            It is more appropriate to use the first approach if you are focused on individual differences (i.e. trying to characterize how people are different) rather than general trends in behavior (i.e. tring to characterize how various factors affect decision-making within a person).
        
        .. tab-item:: R

            :: 

                excluded = which(is.infinite(subjectData$modelAIC) | is.infinite(altSubjectData$modelAlt1AIC) | is.infinite(altSubjectData$modelAlt2AIC))
                averageAIC = c(mean(subjectData$modelAIC[-excluded]), mean(altSubjectData$modelAlt1AIC[-excluded]), mean(altSubjectData$modelAlt2AIC[-excluded]))
                fullAIC = length(trialData$SubjectID) * log(sum(subjectData$modelSS)/length(trialData$SubjectID)) + (2 * k * length(subjectData$SubjectID))
                fullAICAlt1 = length(trialData$SubjectID) * log(sum(altSubjectData$modelAlt1SS)/length(trialData$SubjectID)) + (2 * 0 * length(subjectData$SubjectID))
                fullAICAlt2 = length(trialData$SubjectID) * log(sum(altSubjectData$modelAlt2SS)/length(trialData$SubjectID)) + (2 * 0 * length(subjectData$SubjectID))

                bestModel = c("Favored Model", "Alternative Model 1", "Alternative Model 2")[which(averageAIC == min(averageAIC))] #best model based on average performance per subject
                bestModelFullDataset = c("Favored Model", "Alternative Model 1", "Alternative Model 2")[which(c(fullAIC, fullAICAlt1, fullAICAlt2) == min(c(fullAIC, fullAICAlt1, fullAICAlt2)))] #best model based on all data observations

        .. tab-item:: MatLab

            ::

                excluded = find(isinf(subjectData.modelAIC) | isinf(altSubjectData.modelAlt1AIC) | isinf(altSubjectData.modelAlt2AIC));
                averageAIC = [mean(subjectData.modelAIC(~excluded)), mean(altSubjectData.modelAlt1AIC(~excluded)), mean(altSubjectData.modelAlt2AIC(~excluded))];
                fullAIC = length(trialData.SubjectID) * log(sum(subjectData.modelSS)/length(trialData.SubjectID)) + (2 * k * length(subjectData.SubjectID));
                fullAICAlt1 = length(trialData.SubjectID) * log(sum(altSubjectData.modelAlt1SS)/length(trialData.SubjectID)) + (2 * 0 * length(subjectData.SubjectID));
                fullAICAlt2 = length(trialData.SubjectID) * log(sum(altSubjectData.modelAlt2SS)/length(trialData.SubjectID)) + (2 * 0 * length(subjectData.SubjectID));

                [~, bestModel] = min(averageAIC); % best model based on average performance per subject
                [~, bestModelFullDataset] = min([fullAIC, fullAICAlt1, fullAICAlt2]); % best model based on all data observations
                bestModel = {'Favored Model', 'Alternative Model 1', 'Alternative Model 2'}{bestModel};
                bestModelFullDataset = {'Favored Model', 'Alternative Model 1', 'Alternative Model 2'}{bestModelFullDataset};


        .. tab-item:: Python

            ::
                
                excluded = np.where(np.isinf(subjectData['modelAIC']) | np.isinf(altSubjectData['modelAlt1AIC']) | np.isinf(altSubjectData['modelAlt2AIC']))[0]
                averageAIC = [np.mean(subjectData['modelAIC'][~excluded]), np.mean(altSubjectData['modelAlt1AIC'][~excluded]), np.mean(altSubjectData['modelAlt2AIC'][~excluded])]
                fullAIC = len(trialData['SubjectID']) * np.log(np.sum(subjectData['modelSS']) / len(trialData['SubjectID'])) + (2 * k * len(subjectData['SubjectID']))
                fullAICAlt1 = len(trialData['SubjectID']) * np.log(np.sum(altSubjectData['modelAlt1SS']) / len(trialData['SubjectID'])) + (2 * 0 * len(subjectData['SubjectID']))
                fullAICAlt2 = len(trialData['SubjectID']) * np.log(np.sum(altSubjectData['modelAlt2SS']) / len(trialData['SubjectID'])) + (2 * 0 * len(subjectData['SubjectID']))

                bestModel = ['Favored Model', 'Alternative Model 1', 'Alternative Model 2'][np.argmin(averageAIC)]  # best model based on average performance per subject
                bestModelFullDataset = ['Favored Model', 'Alternative Model 1', 'Alternative Model 2'][np.argmin([fullAIC, fullAICAlt1, fullAICAlt2])]  # best model based on all data observations


Tutorials
==========

Tutorial 1 - van Baar, Chang, & Sanfey, 2019
----------------------

.. dropdown:: Create :bdg-warning:`Utility` Equations for our Alternative Models

    .. tab-set::

        .. tab-item:: R

            ::

                utility_greed = function(greed){
                    return(greed)
                }
                utility_guilt = function(theta, greed, guilt){
                    return(theta * greed - (1 - theta) * guilt)
                }
                utility_inequity = function(theta, greed, inequity){
                    return(theta * greed - (1 - theta) * inequity)
                }

        .. tab-item:: MatLab

            ::

                function value = utility_greed(greed)
                    value = greed;
                end

                function value = utility_guilt(theta, greed, guilt)
                    value = theta * greed - (1 - theta) * guilt;
                end

                function value = utility_inequity(theta, greed, inequity)
                    value = theta * greed - (1 - theta) * inequity;
                end

        .. tab-item:: Python

            ::

                def utility_greed(greed):
                    return greed

                def utility_guilt(theta, greed, guilt):
                    return theta * greed - (1 - theta) * guilt

                def utility_inequity(theta, greed, inequity):
                    return theta * greed - (1 - theta) * inequity

.. dropdown:: Preallocating and Defining Functions

    .. tab-set::

        .. tab-item:: R

            ::

                obj_function_guilt = function(params, df, method = "OLS") {
                    Theta = params[1]
                    
                    predicted_utility = vector('numeric', length(df[,1]))
                    observed_utility = vector('numeric', length(df[,1]))
                    chosen = as.numeric(df[,4]) + 1
                    for (k in 1:length(df[,1])){
                        I = df[k, 2]
                        M = df[k, 3]
                        B = 4
                        E = 10
                        if (I > 10) {Choices = seq(0, (I*M), round((I*M)/10))} else {Choices = seq(0, (I * M), 1)}
                        
                        Utility = vector('numeric', length(Choices))
                        for (n in 1:length(Choices)){
                        Utility[n] = utility_guilt(theta = Theta, 
                                                    greed = payout_maximization(I, M, Choices[n]), 
                                                    guilt = guilt(I, B, Choices[n], M))
                        }
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[chosen[k]]
                    }
                    if (method == "OLS"){
                        return(sum((predicted_utility - observed_utility)**2))
                    } else if (method == "MLE"){
                        return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
                    }
                } 
                obj_function_inequity = function(params, df, method = "OLS") {
                    Theta = params[1]
                    
                    predicted_utility = vector('numeric', length(df[,1]))
                    observed_utility = vector('numeric', length(df[,1]))
                    chosen = as.numeric(df[,4]) + 1
                    for (k in 1:length(df[,1])){
                        I = df[k, 2]
                        M = df[k, 3]
                        B = 4
                        E = 10
                        if (I > 10) {Choices = seq(0, (I*M), round((I*M)/10))} else {Choices = seq(0, (I * M), 1)}
                        
                        Utility = vector('numeric', length(Choices))
                        for (n in 1:length(Choices)){
                        Utility[n] = utility_inequity(theta = Theta, 
                                                        greed = payout_maximization(I, M, Choices[n]), 
                                                        inequity = inequity(I, M, Choices[n], E))
                        }
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[chosen[k]]
                    }
                    if (method == "OLS"){
                        return(sum((predicted_utility - observed_utility)**2))
                    } else if (method == "MLE"){
                        return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
                    }
                } 

                altSubjectData = data.frame()

        .. tab-item:: MatLab

            ::

                function result = obj_function_guilt(params, df, method)
                    Theta = params(1);

                    predicted_utility = zeros(length(df), 1);
                    observed_utility = zeros(length(df), 1);
                    chosen = df(:, 4) + 1;

                    for k = 1:length(df)
                        I = df(k, 2);
                        M = df(k, 3);
                        B = 4;
                        E = 10;

                        if I > 10
                            Choices = 0:round((I * M) / 10):(I * M);
                        else
                            Choices = 0:1:(I * M);
                        end

                        Utility = zeros(size(Choices));

                        for n = 1:length(Choices)
                            Utility(n) = utility_guilt(Theta, payout_maximization(I, M, Choices(n)), guilt(I, B, Choices(n), M));
                        end

                        predicted_utility(k) = max(Utility);
                        observed_utility(k) = Utility(chosen(k));
                    end

                    if strcmp(method, 'OLS')
                        result = sum((predicted_utility - observed_utility).^2);
                    elseif strcmp(method, 'MLE')
                        result = -1 * sum(log(normpdf(observed_utility, predicted_utility, sd)));
                    end
                end

                function result = obj_function_inequity(params, df, method)
                    Theta = params(1);

                    predicted_utility = zeros(length(df), 1);
                    observed_utility = zeros(length(df), 1);
                    chosen = df(:, 4) + 1;

                    for k = 1:length(df)
                        I = df(k, 2);
                        M = df(k, 3);
                        B = 4;
                        E = 10;

                        if I > 10
                            Choices = 0:round((I * M) / 10):(I * M);
                        else
                            Choices = 0:1:(I * M);
                        end

                        Utility = zeros(size(Choices));

                        for n = 1:length(Choices)
                            Utility(n) = utility_inequity(Theta, payout_maximization(I, M, Choices(n)), inequity(I, M, Choices(n), E));
                        end

                        predicted_utility(k) = max(Utility);
                        observed_utility(k) = Utility(chosen(k));
                    end

                    if strcmp(method, 'OLS')
                        result = sum((predicted_utility - observed_utility).^2);
                    elseif strcmp(method, 'MLE')
                        result = -1 * sum(log(normpdf(observed_utility, predicted_utility, sd)));
                    end
                end

                altSubjectData = table();

        .. tab-item:: Python

            ::

                def obj_function_guilt(params, df, method='OLS'):
                    Theta = params[0]

                    predicted_utility = np.zeros(len(df))
                    observed_utility = np.zeros(len(df))
                    chosen = df[:, 3].astype(int)

                    for k in range(len(df)):
                        I = df[k, 1]
                        M = df[k, 2]
                        B = 4
                        E = 10

                        if I > 10:
                            Choices = np.arange(0, (I * M) + 1, round((I * M) / 10))
                        else:
                            Choices = np.arange(0, (I * M) + 1, 1)

                        Utility = np.zeros(len(Choices))

                        for n in range(len(Choices)):
                            Utility[n] = utility_guilt(Theta, payout_maximization(I, M, Choices[n]), guilt(I, B, Choices[n], M))

                        predicted_utility[k] = np.max(Utility)
                        observed_utility[k] = Utility[chosen[k] - 1]

                    if method == 'OLS':
                        result = np.sum((predicted_utility - observed_utility)**2)
                    elif method == 'MLE':
                        result = -1 * np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))

                    return result


                def obj_function_inequity(params, df, method='OLS'):
                    Theta = params[0]

                    predicted_utility = np.zeros(len(df))
                    observed_utility = np.zeros(len(df))
                    chosen = df[:, 3].astype(int)

                    for k in range(len(df)):
                        I = df[k, 1]
                        M = df[k, 2]
                        B = 4
                        E = 10

                        if I > 10:
                            Choices = np.arange(0, (I * M) + 1, round((I * M) / 10))
                        else:
                            Choices = np.arange(0, (I * M) + 1, 1)

                        Utility = np.zeros(len(Choices))

                        for n in range(len(Choices)):
                            Utility[n] = utility_inequity(Theta, payout_maximization(I, M, Choices[n]), inequity(I, M, Choices[n], E))

                        predicted_utility[k] = np.max(Utility)
                        observed_utility[k] = Utility[chosen[k] - 1]

                    if method == 'OLS':
                        result = np.sum((predicted_utility - observed_utility)**2)
                    elif method == 'MLE':
                        result = -1 * np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))

                    return result

                altSubjectData = pd.DataFrame()

.. dropdown:: Recover :bdg-success:`Free Parameters` for Each Alternative Model, Per :bdg-success:`Subject`

    .. tab-set::

        .. tab-item:: R

            ::

                for (i in 1:length(included_subjects)){
                    df = trialData[which(included_subjects[i] == trialData$Subject), ]
  
                    result_guilt = optim(fn = obj_function_guilt, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                    result_inequity = optim(fn = obj_function_inequity, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')

                    #Determine Model Predictions using result_guilt$pars and result_inequity$pars

                }

        .. tab-item:: MatLab

            ::

                for i = 1:length(included_subjects)
                    subject = included_subjects(i);
                    df = trialData(trialData.Subject == subject, :);

                    % Optimize for guilt
                    result_guilt = fmincon(@(theta) obj_function_guilt(theta, df), 0.5, [], [], [], [], 0, 1);

                    % Optimize for inequity
                    result_inequity = fmincon(@(theta) obj_function_inequity(theta, df), 0.5, [], [], [], [], 0, 1);

                    % Determine Model Predictions using result_guilt.pars and result_inequity.pars
                    
                end

        .. tab-item:: Python

            ::

                for i = 1:length(included_subjects)
                    subject = included_subjects(i);
                    df = trialData(trialData.Subject == subject, :);

                    % Optimize for guilt
                    result_guilt = fmincon(@(theta) obj_function_guilt(theta, df), 0.5, [], [], [], [], 0, 1);

                    % Optimize for inequity
                    result_inequity = fmincon(@(theta) obj_function_inequity(theta, df), 0.5, [], [], [], [], 0, 1);

                end

.. dropdown:: Determine Predicted :bdg-danger:`Decisions` for Each Alternative Model, Per :bdg-success:`Subject`

    .. tab-set::

        .. tab-item:: R

            ::

                for (i in 1:length(included_subjects)){
                    df = trialData[which(included_subjects[i] == trialData$Subject), ]
  
                    result_guilt = optim(fn = obj_function_guilt, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                    result_inequity = optim(fn = obj_function_inequity, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')

                    #Just Added

                    df$PredictionGreed = vector('numeric', length(df$Subject))
                    df$PredictionGuilt = df$PredictionAlt1
                    df$PredictionInequity = df$PredictionAlt1
                    for (k in 1:length(df$Returned)){
                        I = df$Investment[k]
                        M = df$Multiplier[k]
                        R = df$Returned[k]
                        B = 4
                        E = 10
                        if (I > 10) {Choices = seq(0, (I*M), round((I*M)/10))} else {Choices = seq(0, (I * M), 1)}
                        UtilityGreed = vector('numeric', length(Choices))
                        UtilityGuilt = vector('numeric', length(Choices))
                        UtilityInequity = vector('numeric', length(Choices))
                        for (n in 1:length(Choices)){
                            UtilityGreed[n] = utility_greed(greed = payout_maximization(I, M, Choices[n]))
                            UtilityGuilt[n] = utility_guilt(theta = result_guilt$par[1],
                                                            greed = payout_maximization(I, M, Choices[n]),
                                                            guilt = guilt(I, B, Choices[n], M))
                            UtilityInequity[n] = utility_inequity(theta = result_guilt$par[1],
                                                                greed = payout_maximization(I, M, Choices[n]),
                                                                inequity = inequity(I, M, Choices[n], E))
                        }
                        correct_choice_greed = which(UtilityGreed == max(UtilityGreed))
                        correct_choice_guilt = which(UtilityGuilt == max(UtilityGuilt))
                        correct_choice_inequity = which(UtilityInequity == max(UtilityInequity))
                        if (length(correct_choice_greed) > 1){
                            correct_choice_greed = correct_choice_greed[sample(correct_choice_greed, 1)]
                        }
                        if (length(correct_choice_guilt) > 1){
                            correct_choice_guilt = correct_choice_guilt[sample(correct_choice_guilt, 1)]
                        }
                        if (length(correct_choice_inequity) > 1){
                            correct_choice_inequity = correct_choice_inequity[sample(correct_choice_inequity, 1)]
                        }
                        df$PredictionGreed[k] = Choices[correct_choice_greed]
                        df$PredictionGuilt[k] = Choices[correct_choice_guilt]
                        df$PredictionInequity[k] = Choices[correct_choice_inequity]
                    }
  
                    model_NLL_Greed = -2 * log(sum(dnorm(as.numeric(df$Returned), mean = df$PredictionGreed)))
                    model_SS_Greed = sum((as.numeric(df$Returned) - df$PredictionGreed)**2)
                    model_NLL_Guilt = -2 * log(sum(dnorm(as.numeric(df$Returned), mean = df$PredictionGuilt)))
                    model_SS_Guilt = sum((as.numeric(df$Returned) - df$PredictionGuilt)**2)
                    model_NLL_Inequity = -2 * log(sum(dnorm(as.numeric(df$Returned), mean = df$PredictionInequity)))
                    model_SS_Inequity = sum((as.numeric(df$Returned) - df$PredictionInequity)**2)
                    
                    altSubjectData[i, 1:9] = c(included_subjects[i], result_guilt$par[1], result_inequity$par[1], model_NLL_Greed, model_SS_Greed, model_NLL_Guilt, model_SS_Guilt, model_NLL_Inequity, model_SS_Inequity)
                }
                colnames(altSubjectData) = c('SubjectID', 'guilt_theta', 'inequity_theta', 'greed_modelNLL', 'greed_modelSS', 'guilt_modelNLL', 'guilt_modelSS', 'inequity_ModelNLL', 'inequity_ModelSS')


        .. tab-item:: MatLab

            ::

                % Loop over included_subjects
                for i = 1:length(included_subjects)
                    subject = included_subjects(i);
                    df = trialData(trialData.Subject == subject, :);

                    % Optimize for guilt
                    result_guilt = fmincon(@(theta) obj_function_guilt(theta, df), 0.5, [], [], [], [], 0, 1);

                    % Optimize for inequity
                    result_inequity = fmincon(@(theta) obj_function_inequity(theta, df), 0.5, [], [], [], [], 0, 1);

                    % Determine Predictions
                    df.PredictionGreed = df.PredictionAlt1;
                    df.PredictionGuilt = df.PredictionAlt1;
                    df.PredictionInequity = df.PredictionAlt1;

                    for k = 1:height(df)
                        I = df.Investment(k);
                        M = df.Multiplier(k);
                        R = df.Returned(k);
                        B = 4;
                        E = 10;

                        if I > 10
                            Choices = 0:round((I * M) / 10, 3):(I * M);
                        else
                            Choices = 0:1:(I * M);
                        end

                        UtilityGreed = arrayfun(@(choice) utility_greed(payout_maximization(I, M, choice)), Choices);
                        UtilityGuilt = arrayfun(@(choice) utility_guilt(result_guilt, payout_maximization(I, M, choice), guilt(I, B, choice, M)), Choices);
                        UtilityInequity = arrayfun(@(choice) utility_inequity(result_inequity, payout_maximization(I, M, choice), inequity(I, M, choice, E)), Choices);

                        [~, correct_choice_greed] = max(UtilityGreed);
                        [~, correct_choice_guilt] = max(UtilityGuilt);
                        [~, correct_choice_inequity] = max(UtilityInequity);

                        df.PredictionGreed(k) = Choices(correct_choice_greed);
                        df.PredictionGuilt(k) = Choices(correct_choice_guilt);
                        df.PredictionInequity(k) = Choices(correct_choice_inequity);
                    end

                    % Calculate model statistics
                    model_NLL_Greed = -2 * sum(log(normpdf(df.Returned, df.PredictionGreed, 1)));
                    model_SS_Greed = sum((df.Returned - df.PredictionGreed).^2);

                    model_NLL_Guilt = -2 * sum(log(normpdf(df.Returned, df.PredictionGuilt, 1)));
                    model_SS_Guilt = sum((df.Returned - df.PredictionGuilt).^2);

                    model_NLL_Inequity = -2 * sum(log(normpdf(df.Returned, df.PredictionInequity, 1)));
                    model_SS_Inequity = sum((df.Returned - df.PredictionInequity).^2);

                    altSubjectData(i, :) = [subject, result_guilt, result_inequity, model_NLL_Greed, model_SS_Greed, model_NLL_Guilt, model_SS_Guilt, model_NLL_Inequity, model_SS_Inequity];
                end

                % Set column names
                colnames_altSubjectData = {'SubjectID', 'guilt_theta', 'inequity_theta', 'greed_modelNLL', 'greed_modelSS', 'guilt_modelNLL', 'guilt_modelSS', 'inequity_ModelNLL', 'inequity_ModelSS'};
                altSubjectData = array2table(altSubjectData, 'VariableNames', colnames_altSubjectData);

        .. tab-item:: Python

            ::

                for i, subject in enumerate(included_subjects):
                    df = trialData[trialData['Subject'] == subject]

                    # Optimize for guilt
                    result_guilt = minimize(obj_function_guilt, x0=[0.5], args=(df,), bounds=[(0, 1)], method='L-BFGS-B')

                    # Optimize for inequity
                    result_inequity = minimize(obj_function_inequity, x0=[0.5], args=(df,), bounds=[(0, 1)], method='L-BFGS-B')

                    # Determine Predictions
                    df['PredictionGreed'] = df['PredictionAlt1']
                    df['PredictionGuilt'] = df['PredictionAlt1']
                    df['PredictionInequity'] = df['PredictionAlt1']

                    for k in range(len(df['Returned'])):
                        I, M, R = df.loc[k, ['Investment', 'Multiplier', 'Returned']]
                        B, E = 4, 10
                        Choices = np.arange(0, I * M + 0.001, round((I * M) / 10, 3)) if I > 10 else np.arange(0, I * M + 0.001, 1)
                        
                        UtilityGreed = np.array([utility_greed(payout_maximization(I, M, choice)) for choice in Choices])
                        UtilityGuilt = np.array([utility_guilt(result_guilt.x[0], payout_maximization(I, M, choice), guilt(I, B, choice, M)) for choice in Choices])
                        UtilityInequity = np.array([utility_inequity(result_inequity.x[0], payout_maximization(I, M, choice), inequity(I, M, choice, E)) for choice in Choices])

                        correct_choice_greed = np.argmax(UtilityGreed)
                        correct_choice_guilt = np.argmax(UtilityGuilt)
                        correct_choice_inequity = np.argmax(UtilityInequity)

                        df.at[k, 'PredictionGreed'] = Choices[correct_choice_greed]
                        df.at[k, 'PredictionGuilt'] = Choices[correct_choice_guilt]
                        df.at[k, 'PredictionInequity'] = Choices[correct_choice_inequity]

                    # Calculate model statistics
                    model_NLL_Greed = -2 * np.log(np.sum(norm.pdf(df['Returned'], loc=df['PredictionGreed'], scale=1)))
                    model_SS_Greed = np.sum((df['Returned'] - df['PredictionGreed'])**2)
                    
                    model_NLL_Guilt = -2 * np.log(np.sum(norm.pdf(df['Returned'], loc=df['PredictionGuilt'], scale=1)))
                    model_SS_Guilt = np.sum((df['Returned'] - df['PredictionGuilt'])**2)

                    model_NLL_Inequity = -2 * np.log(np.sum(norm.pdf(df['Returned'], loc=df['PredictionInequity'], scale=1)))
                    model_SS_Inequity = np.sum((df['Returned'] - df['PredictionInequity'])**2)

                    altSubjectData[i, :] = [subject, result_guilt.x[0], result_inequity.x[0], model_NLL_Greed, model_SS_Greed, model_NLL_Guilt, model_SS_Guilt, model_NLL_Inequity, model_SS_Inequity]

                # Set column names
                colnames_altSubjectData = ['SubjectID', 'guilt_theta', 'inequity_theta', 'greed_modelNLL', 'greed_modelSS', 'guilt_modelNLL', 'guilt_modelSS', 'inequity_ModelNLL', 'inequity_ModelSS']
                altSubjectData = pd.DataFrame(altSubjectData, columns=colnames_altSubjectData)

.. dropdown:: Compute Model Fit Index for Each :bdg-success:`Subject`, for Each Alternative Model

    .. tab-set::

        .. tab-item:: R

            ::

                N = length(altSubjectData)

                altSubjectData$modelAICGreed = N * log(altSubjectData$greed_modelSS/N) + 2*0
                altSubjectData$modelAICGuilt = N * log(altSubjectData$guilt_modelSS/N) + 2*1
                altSubjectData$modelAICInequity = N * log(altSubjectData$inequity_ModelSS/N) + 2*1

        .. tab-item:: MatLab

            ::

                N = length(altSubjectData);

                altSubjectData.modelAICGreed = N * log(altSubjectData.greed_modelSS / N) + 2 * 0;
                altSubjectData.modelAICGuilt = N * log(altSubjectData.guilt_modelSS / N) + 2 * 1;
                altSubjectData.modelAICInequity = N * log(altSubjectData.inequity_ModelSS / N) + 2 * 1;

        .. tab-item:: Python

            ::

                N = len(altSubjectData)

                altSubjectData['modelAICGreed'] = N * np.log(altSubjectData['greed_modelSS'] / N) + 2 * 0
                altSubjectData['modelAICGuilt'] = N * np.log(altSubjectData['guilt_modelSS'] / N) + 2 * 1
                altSubjectData['modelAICInequity'] = N * np.log(altSubjectData['inequity_ModelSS'] / N) + 2 * 1

.. dropdown:: Compare Model Performance

    .. tab-set::

        .. tab-item:: R

            ::

                #subjects excluded from analyses not included in dataset
                averageAIC = c(mean(subjectData$modelAIC), mean(altSubjectData$modelAICGreed), mean(altSubjectData$modelAICGuilt), mean(altSubjectData$modelAICInequity))
                fullAIC = length(trialData$Subject) * log(sum(subjectData$modelSS)/length(trialData$Subject)) + (2 * k * length(subjectData$Subject))
                fullAICGreed = length(trialData$Subject) * log(sum(altSubjectData$greed_modelSS)/length(trialData$Subject)) + (2 * 0 * length(subjectData$Subject))
                fullAICGuilt = length(trialData$Subject) * log(sum(altSubjectData$guilt_modelSS)/length(trialData$Subject)) + (2 * 1 * length(subjectData$Subject))
                fullAICInequity = length(trialData$Subject) * log(sum(altSubjectData$inequity_ModelSS)/length(trialData$Subject)) + (2 * 1 * length(subjectData$Subject))

                bestModel = c("Moral Strategies Model", "Greed Model", "Guilt Model", "Inequity Model")[which(averageAIC == min(averageAIC))] #best model based on average performance per subject
                bestModelFullDataset = c("Moral Strategies Model", "Greed Model", "Guilt Model", "Inequity Model")[which(c(fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity) == min(c(fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity)))] #best model based on all data observations

        .. tab-item:: MatLab

            ::

                averageAIC = [mean(subjectData.modelAIC), mean(altSubjectData.modelAICGreed), mean(altSubjectData.modelAICGuilt), mean(altSubjectData.modelAICInequity)];
                fullAIC = length(trialData.Subject) * log(sum(subjectData.modelSS)/length(trialData.Subject)) + (2 * k * length(subjectData.Subject));
                fullAICGreed = length(trialData.Subject) * log(sum(altSubjectData.greed_modelSS)/length(trialData.Subject)) + (2 * 0 * length(subjectData.Subject));
                fullAICGuilt = length(trialData.Subject) * log(sum(altSubjectData.guilt_modelSS)/length(trialData.Subject)) + (2 * 1 * length(subjectData.Subject));
                fullAICInequity = length(trialData.Subject) * log(sum(altSubjectData.inequity_ModelSS)/length(trialData.Subject)) + (2 * 1 * length(subjectData.Subject));

                [~, bestModelIndex] = min(averageAIC);
                bestModel = {'Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'}{bestModelIndex}; % best model based on average performance per subject

                [~, bestModelFullDatasetIndex] = min([fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity]);
                bestModelFullDataset = {'Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'}{bestModelFullDatasetIndex}; % best model based on all data observations

        .. tab-item:: Python

            ::

                averageAIC = [np.mean(subjectData['modelAIC']), np.mean(altSubjectData['modelAICGreed']), np.mean(altSubjectData['modelAICGuilt']), np.mean(altSubjectData['modelAICInequity'])]
                fullAIC = len(trialData['Subject']) * np.log(np.sum(subjectData['modelSS']) / len(trialData['Subject'])) + (2 * k * len(subjectData['Subject']))
                fullAICGreed = len(trialData['Subject']) * np.log(np.sum(altSubjectData['greed_modelSS']) / len(trialData['Subject'])) + (2 * 0 * len(subjectData['Subject']))
                fullAICGuilt = len(trialData['Subject']) * np.log(np.sum(altSubjectData['guilt_modelSS']) / len(trialData['Subject'])) + (2 * 1 * len(subjectData['Subject']))
                fullAICInequity = len(trialData['Subject']) * np.log(np.sum(altSubjectData['inequity_ModelSS']) / len(trialData['Subject'])) + (2 * 1 * len(subjectData['Subject']))

                bestModelIndex = np.argmin(averageAIC)
                bestModel = ['Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'][bestModelIndex]  # best model based on average performance per subject

                bestModelFullDatasetIndex = np.argmin([fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity])
                bestModelFullDataset = ['Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'][bestModelFullDatasetIndex]  # best model based on all data observations

Tutorial 2 - Galvan & Sanfey, 2024
-------------------

.. dropdown:: Create :bdg-warning:`Utility` Equations for our Alternative Models

    .. tab-set::

        .. tab-item:: R

            ::

                utilityP = function(Payout){
                    return(Payout)
                }
                utilityL = function(Equality){
                    return(Equality)
                }
                utilityT = function(Equity){
                    return(Equity)
                }
                utilityPL = function(Payout, Equality, theta){
                    return(theta * Payout + (1 - theta) * Equality)
                }
                utilityPT = function(Payout, Equity, theta){
                    return(theta * Payout + (1 - theta) * Equity)
                }
                utilityLT = function(Equity, Equality, phi){
                    return(phi * Equality + (1 - phi) * Equity)
                }

        .. tab-item:: MatLab

            ::

                function u = utilityP(Payout)
                    u = Payout;
                end

                function u = utilityL(Equality)
                    u = Equality;
                end

                function u = utilityT(Equity)
                    u = Equity;
                end

                function u = utilityPL(Payout, Equality, theta)
                    u = theta * Payout + (1 - theta) * Equality;
                end

                function u = utilityPT(Payout, Equity, theta)
                    u = theta * Payout + (1 - theta) * Equity;
                end

                function u = utilityLT(Equity, Equality, phi)
                    u = phi * Equality + (1 - phi) * Equity;
                end

        .. tab-item:: Python

            ::

                def utilityP(Payout):
                    return Payout

                def utilityL(Equality):
                    return Equality

                def utilityT(Equity):
                    return Equity

                def utilityPL(Payout, Equality, theta):
                    return theta * Payout + (1 - theta) * Equality

                def utilityPT(Payout, Equity, theta):
                    return theta * Payout + (1 - theta) * Equity

                def utilityLT(Equity, Equality, phi):
                    return phi * Equality + (1 - phi) * Equity

.. dropdown:: Preallocating and Defining Functions

    .. tab-set::

        .. tab-item:: R

            ::

                obj_functionPL = function(params, df, method = "OLS") {
                    Theta = params[1]

                    predicted_utility = vector('numeric', length(df[,1]))
                    observed_utility = vector('numeric', length(df[,1]))
                    choices = seq(0, 1, 0.1)

                    for (k in 1:length(df[,1])){
                        Utility = vector('numeric', length(choices))
                        for (n in 1:length(choices)){
                            Utility[n] = utilityPL(theta = Theta,
                                                   Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
                                                   Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
                        }
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[which(df[k, 11] == choices)]
                    }
                    if (method == "OLS"){
                        return(sum((predicted_utility - observed_utility)**2))
                    } else if (method == "MLE"){
                        return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
                    }
                }

                obj_functionPT = function(params, df, method = "OLS") {
                    Theta = params[1]

                    predicted_utility = vector('numeric', length(df[,1]))
                    observed_utility = vector('numeric', length(df[,1]))
                    choices = seq(0, 1, 0.1)

                    for (k in 1:length(df[,1])){
                        Utility = vector('numeric', length(choices))
                        for (n in 1:length(choices)){
                            Utility[n] = utilityPT(theta = Theta,
                                                   Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
                                                   Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
                        }
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[which(df[k, 11] == choices)]
                    }
                    if (method == "OLS"){
                        return(sum((predicted_utility - observed_utility)**2))
                    } else if (method == "MLE"){
                        return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
                    }
                }

                obj_functionLT = function(params, df, method = "OLS") {
                    Phi = params[1]

                    predicted_utility = vector('numeric', length(df[,1]))
                    observed_utility = vector('numeric', length(df[,1]))
                    choices = seq(0, 1, 0.1)

                    for (k in 1:length(df[,1])){
                        Utility = vector('numeric', length(choices))
                        for (n in 1:length(choices)){
                            Utility[n] = utilityLT(phi = Phi,
                                                   Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
                                                   Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
                        }
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[which(df[k, 11] == choices)]
                    }
                    if (method == "OLS"){
                        return(sum((predicted_utility - observed_utility)**2))
                    } else if (method == "MLE"){
                        return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
                    }
                }

                altSubjectData = data.frame()

        .. tab-item:: MatLab

            ::

                function objValue = obj_functionPL(params, df, method)
                    Theta = params(1);

                    predicted_utility = zeros(length(df(:, 1)), 1);
                    observed_utility = zeros(length(df(:, 1)), 1);
                    choices = 0:0.1:1;

                    for k = 1:length(df(:, 1))
                        Utility = zeros(length(choices), 1);
                        for n = 1:length(choices)
                            Utility(n) = utilityPL(Theta, ...
                                                equality(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)), ...
                                                payout(new_value(df(k, 1), choices(n)), df(k, 1), choices(n)));
                        end
                        predicted_utility(k) = max(Utility);
                        observed_utility(k) = Utility(df(k, 11) == choices);
                    end

                    if strcmp(method, 'OLS')
                        objValue = sum((predicted_utility - observed_utility).^2);
                    elseif strcmp(method, 'MLE')
                        objValue = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
                    end
                end

                function objValue = obj_functionPT(params, df, method)
                    Theta = params(1);

                    predicted_utility = zeros(length(df(:, 1)), 1);
                    observed_utility = zeros(length(df(:, 1)), 1);
                    choices = 0:0.1:1;

                    for k = 1:length(df(:, 1))
                        Utility = zeros(length(choices), 1);
                        for n = 1:length(choices)
                            Utility(n) = utilityPT(Theta, ...
                                                equity(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)), ...
                                                payout(new_value(df(k, 1), choices(n)), df(k, 1), choices(n)));
                        end
                        predicted_utility(k) = max(Utility);
                        observed_utility(k) = Utility(df(k, 11) == choices);
                    end

                    if strcmp(method, 'OLS')
                        objValue = sum((predicted_utility - observed_utility).^2);
                    elseif strcmp(method, 'MLE')
                        objValue = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
                    end
                end

                function objValue = obj_functionLT(params, df, method)
                    Phi = params(1);

                    predicted_utility = zeros(length(df(:, 1)), 1);
                    observed_utility = zeros(length(df(:, 1)), 1);
                    choices = 0:0.1:1;

                    for k = 1:length(df(:, 1))
                        Utility = zeros(length(choices), 1);
                        for n = 1:length(choices)
                            Utility(n) = utilityLT(Phi, ...
                                                equity(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)), ...
                                                equality(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)));
                        end
                        predicted_utility(k) = max(Utility);
                        observed_utility(k) = Utility(df(k, 11) == choices);
                    end

                    if strcmp(method, 'OLS')
                        objValue = sum((predicted_utility - observed_utility).^2);
                    elseif strcmp(method, 'MLE')
                        objValue = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
                    end
                end

                altSubjectData = table();

        .. tab-item:: Python

            ::

                def obj_functionPL(params, df, method="OLS"):
                    Theta = params[0]

                    predicted_utility = np.zeros(len(df), dtype=float)
                    observed_utility = np.zeros(len(df), dtype=float)
                    choices = np.arange(0, 1.1, 0.1)

                    for k in range(len(df)):
                        Utility = np.zeros(len(choices), dtype=float)
                        for n in range(len(choices)):
                            Utility[n] = utilityPL(Theta,
                                                equality(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]),
                                                payout(new_value(df.iloc[k, 0], choices[n]), df.iloc[k, 0], choices[n]))
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[df.iloc[k, 10] == choices]

                    if method == "OLS":
                        obj_value = np.sum((predicted_utility - observed_utility)**2)
                    elif method == "MLE":
                        obj_value = -np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))

                    return obj_value


                def obj_functionPT(params, df, method="OLS"):
                    Theta = params[0]

                    predicted_utility = np.zeros(len(df), dtype=float)
                    observed_utility = np.zeros(len(df), dtype=float)
                    choices = np.arange(0, 1.1, 0.1)

                    for k in range(len(df)):
                        Utility = np.zeros(len(choices), dtype=float)
                        for n in range(len(choices)):
                            Utility[n] = utilityPT(Theta,
                                                equity(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]),
                                                payout(new_value(df.iloc[k, 0], choices[n]), df.iloc[k, 0], choices[n]))
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[df.iloc[k, 10] == choices]

                    if method == "OLS":
                        obj_value = np.sum((predicted_utility - observed_utility)**2)
                    elif method == "MLE":
                        obj_value = -np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))

                    return obj_value


                def obj_functionLT(params, df, method="OLS"):
                    Phi = params[0]

                    predicted_utility = np.zeros(len(df), dtype=float)
                    observed_utility = np.zeros(len(df), dtype=float)
                    choices = np.arange(0, 1.1, 0.1)

                    for k in range(len(df)):
                        Utility = np.zeros(len(choices), dtype=float)
                        for n in range(len(choices)):
                            Utility[n] = utilityLT(Phi,
                                                equity(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]),
                                                equality(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]))
                        predicted_utility[k] = max(Utility)
                        observed_utility[k] = Utility[df.iloc[k, 10] == choices]

                    if method == "OLS":
                        obj_value = np.sum((predicted_utility - observed_utility)**2)
                    elif method == "MLE":
                        obj_value = -np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))

                    return obj_value


                altSubjectData = pd.DataFrame()

.. dropdown:: Recover :bdg-success:`Free Parameters` for Each Alternative Model, Per :bdg-success:`Subject` Per :bdg-primary:`Condition`

    .. tab-set::

        .. tab-item:: R

            ::

                for (i in 1:length(included_subjects)){
                    datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
                    fullData = read.csv2(datafile)

                    parametersPerCondition = vector('numeric', length(conditions) * 3)
                    SSPerCondition = vector('numeric', length(conditions) * 6)

                    for (j in 1:length(conditions)){
                        df = fullData[which(fullData$condition == conditions[j]), c(49, 40:48), 33] #49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
                        df$redistributionRate = df$redistributionRate/100 #converting to a decimal from a percent

                        resultPL = optim(obj_functionPL, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                        resultPT = optim(obj_functionPT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                        resultLT = optim(obj_functionLT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                        
                        parametersPerCondition[(((j - 1) * 3) + 1):(j * 3)] = c(resultPL$par[1], resultPT$par[1], resultLT$par[1])

                        #Determine Predictions
                    }
                }

        .. tab-item:: MatLab

            ::

                parametersPerCondition = zeros(length(conditions) * 3, 1);
                SSPerCondition = zeros(length(conditions) * 6, 1);

                for j = 1:length(conditions)
                    df = fullData(fullData.condition == conditions{j}, [49, 40:48, 33]); % 49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
                    df.redistributionRate = df.redistributionRate / 100; % converting to a decimal from a percent

                    resultPL = fmincon(@(params) obj_functionPL(params, df, 'OLS'), 0.5, [], [], [], [], 0, 1);
                    resultPT = fmincon(@(params) obj_functionPT(params, df, 'OLS'), 0.5, [], [], [], [], 0, 1);
                    resultLT = fmincon(@(params) obj_functionLT(params, df, 'OLS'), 0.5, [], [], [], [], 0, 1);

                    parametersPerCondition(((j - 1) * 3) + 1 : j * 3) = [resultPL; resultPT; resultLT];

                    % Determine Predictions
                end

        .. tab-item:: Python

            ::

                from scipy.optimize import minimize

                parametersPerCondition = np.zeros(len(conditions) * 3)
                SSPerCondition = np.zeros(len(conditions) * 6)

                for j in range(len(conditions)):
                    df = fullData[fullData['condition'] == conditions[j]][[47, * range(39, 48), 32]]  # 49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
                    df['redistributionRate'] = df['redistributionRate'] / 100  # converting to a decimal from a percent

                    resultPL = minimize(lambda params: obj_functionPL(params, df, 'OLS'), 0.5, bounds=[(0, 1)])
                    resultPT = minimize(lambda params: obj_functionPT(params, df, 'OLS'), 0.5, bounds=[(0, 1)])
                    resultLT = minimize(lambda params: obj_functionLT(params, df, 'OLS'), 0.5, bounds=[(0, 1)])

                    parametersPerCondition[((j - 1) * 3):(j * 3)] = [resultPL.x[0], resultPT.x[0], resultLT.x[0]]

                    # Determine Predictions

.. dropdown:: Determine Predicted :bdg-danger:`Decisions` for Each Alternative Model, Per :bdg-success:`Subject` Per :bdg-primary:`Condition`

    .. tab-set::

        .. tab-item:: R

            ::

                for (i in 1:length(included_subjects)){
                    datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
                    fullData = read.csv2(datafile)

                    parametersPerCondition = vector('numeric', length(conditions) * 3)
                    SSPerCondition = vector('numeric', length(conditions) * 6)

                    for (j in 1:length(conditions)){
                        df = fullData[which(fullData$condition == conditions[j]), c(49, 40:48), 33] #49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
                        df$redistributionRate = df$redistributionRate/100 #converting to a decimal from a percent

                        resultPL = optim(obj_functionPL, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                        resultPT = optim(obj_functionPT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                        resultLT = optim(obj_functionLT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
                        
                        parametersPerCondition[(((j - 1) * 3) + 1):(j * 3)] = c(resultPL$par[1], resultPT$par[1], resultLT$par[1])

                        #Just Added

                        df$PredictionP = vector('numeric')
                        df$PredictionL = vector('numeric')
                        df$PredictionT = vector('numeric')
                        df$PredictionPT = vector('numeric')
                        df$PredictionPL = vector('numeric')
                        df$PredictionLT = vector('numeric')
                        for (k in 1:length(df$Decisions)){
                            UtilityP = vector('numeric', length(Choices))
                            UtilityL = vector('numeric', length(Choices))
                            UtilityT = vector('numeric', length(Choices))
                            UtilityPT = vector('numeric', length(Choices))
                            UtilityPL = vector('numeric', length(Choices))
                            UtilityLT = vector('numeric', length(Choices))
                            for (n in 1:length(Choices)){
                                UtilityP[n] = utilityP(Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
                                UtilityL[n] = utilityL(Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
                                UtilityT[n] = utilityT(Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
                                UtilityPT[n] = utilityPT(theta = resultPT$par[1],
                                                         Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
                                                         Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
                                UtilityPL[n] = utilityPL(theta = resultPL$par[1],
                                                         Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
                                                         Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
                                UtilityLT[n] = utilityLT(phi = resultLT$par[1],
                                                         Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
                                                         Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
                            }
                            correct_choiceP = which(UtilityP == max(UtilityP))
                            correct_choiceL = which(UtilityL == max(UtilityL))
                            correct_choiceT = which(UtilityT == max(UtilityT))
                            correct_choicePT = which(UtilityPT == max(UtilityPT))
                            correct_choicePL = which(UtilityPL == max(UtilityPL))
                            correct_choiceLT = which(UtilityLT == max(UtilityLT))

                            df$predictedP[k] = new_value(df$myself[k], Choices[correct_choiceP[sample(length(correct_choiceP), 1)]])
                            df$predictedL[k] = new_value(df$myself[k], Choices[correct_choiceL[sample(length(correct_choiceL), 1)]])
                            df$predictedT[k] = new_value(df$myself[k], Choices[correct_choiceT[sample(length(correct_choiceT), 1)]])
                            df$predictedPL[k] = new_value(df$myself[k], Choices[correct_choicePT[sample(length(correct_choicePT), 1)]])
                            df$predictedPT[k] = new_value(df$myself[k], Choices[correct_choicePL[sample(length(correct_choiceTL), 1)]])
                            df$predictedLT[k] = new_value(df$myself[k], Choices[correct_choiceLT[sample(length(correct_choiceLT), 1)]])
                        }
                        df$Outcome = new_value(df$myself, df$redistributionRate)
                        SSPerCondition[(((j - 1) * 6) + 1):(j*6)] = c(sum((df$Outcome - df$predictedP)**2), 
                                                                      sum((df$Outcome - df$predictedL)**2),
                                                                      sum((df$Outcome - df$predictedT)**2),
                                                                      sum((df$Outcome - df$predictedPL)**2),
                                                                      sum((df$Outcome - df$predictedPT)**2),
                                                                      sum((df$Outcome - df$predictedLT)**2))
                    }
                    altSubjectData[i, 1:37] = c(included_subjects[i], parametersPerCondition, SSPerCondition)
                }
                colnames(altSubjectData) = c('SubjectID', 
                                             'thetaPLMerit', 'thetaPTMerit', 'phiLTMerit', 
                                             'thetaPLEntitlement', 'thetaPTEntitlement', 'phiLTEntitlement', 
                                             'thetaPLCorruption', 'thetaPTCorruption', 'phiLTCorruption', 
                                             'thetaPLLuck', 'thetaPTLuck', 'phiLTLuck',
                                             'SSPMerit', 'SSLMerit', 'SSTMerit', 'SSPLMerit', 'SSPTMerit', 'SSLTMerit',
                                             'SSPEntitlement', 'SSLEntitlement', 'SSTEntitlement', 'SSPLEntitlement', 'SSPTEntitlement', 'SSLTEntitlement',
                                             'SSPCorruption', 'SSLCorruption', 'SSTCorruption', 'SSPLCorruption', 'SSPTCorruption', 'SSLTCorruption',
                                             'SSPLuck', 'SSLLuck', 'SSTLuck', 'SSPLLuck', 'SSPLuck', 'SSLTLuck')

        .. tab-item:: MatLab

            ::

                for i = 1:length(included_subjects)
                    datafile = strcat(parentfolder, included_subjects{i}, restoffilepath);
                    fullData = readtable(datafile);

                    parametersPerCondition = zeros(length(conditions) * 3, 1);
                    SSPerCondition = zeros(length(conditions) * 6, 1);

                    for j = 1:length(conditions)
                        df = fullData(fullData.condition == conditions(j), [49, 40:48, 33]);
                        df.redistributionRate = df.redistributionRate / 100;

                        resultPL = fmincon(@(params) obj_functionPL(params, df), 0.5, [], [], [], [], 0, 1);
                        resultPT = fmincon(@(params) obj_functionPT(params, df), 0.5, [], [], [], [], 0, 1);
                        resultLT = fmincon(@(params) obj_functionLT(params, df), 0.5, [], [], [], [], 0, 1);

                        parametersPerCondition(((j - 1) * 3) + 1:(j * 3)) = [resultPL(1), resultPT(1), resultLT(1)];

                        df.PredictionP = zeros(height(df), 1);
                        df.PredictionL = zeros(height(df), 1);
                        df.PredictionT = zeros(height(df), 1);
                        df.PredictionPT = zeros(height(df), 1);
                        df.PredictionPL = zeros(height(df), 1);
                        df.PredictionLT = zeros(height(df), 1);

                        for k = 1:length(df.Decisions)
                            UtilityP = zeros(length(Choices), 1);
                            UtilityL = zeros(length(Choices), 1);
                            UtilityT = zeros(length(Choices), 1);
                            UtilityPT = zeros(length(Choices), 1);
                            UtilityPL = zeros(length(Choices), 1);
                            UtilityLT = zeros(length(Choices), 1);

                            for n = 1:length(Choices)
                                UtilityP(n) = utilityP(payout(new_value(df{k, 1}, Choices{n}), df{k, 1}, Choices{n}));
                                UtilityL(n) = utilityL(equality(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}));
                                UtilityT(n) = utilityT(equity(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}));
                                UtilityPT(n) = utilityPT(resultPT(1), equity(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}), ...
                                                        payout(new_value(df{k, 1}, Choices{n}), df{k, 1}, Choices{n}));
                                UtilityPL(n) = utilityPL(resultPL(1), equality(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}), ...
                                                        payout(new_value(df{k, 1}, Choices{n}), df{k, 1}, Choices{n}));
                                UtilityLT(n) = utilityLT(resultLT(1), equity(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}), ...
                                                        equality(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}));
                            end

                            [~, correct_choiceP] = max(UtilityP);
                            [~, correct_choiceL] = max(UtilityL);
                            [~, correct_choiceT] = max(UtilityT);
                            [~, correct_choicePT] = max(UtilityPT);
                            [~, correct_choicePL] = max(UtilityPL);
                            [~, correct_choiceLT] = max(UtilityLT);

                            df.PredictionP(k) = new_value(df.myself{k}, Choices{correct_choiceP(randi(length(correct_choiceP)))});
                            df.PredictionL(k) = new_value(df.myself{k}, Choices{correct_choiceL(randi(length(correct_choiceL)))});
                            df.PredictionT(k) = new_value(df.myself{k}, Choices{correct_choiceT(randi(length(correct_choiceT)))});
                            df.PredictionPL(k) = new_value(df.myself{k}, Choices{correct_choicePT(randi(length(correct_choicePT)))});
                            df.PredictionPT(k) = new_value(df.myself{k}, Choices{correct_choicePL(randi(length(correct_choicePL)))});
                            df.PredictionLT(k) = new_value(df.myself{k}, Choices{correct_choiceLT(randi(length(correct_choiceLT)))});
                        end

                        df.Outcome = new_value(df.myself, df.redistributionRate);
                        SSPerCondition(((j - 1) * 6) + 1:(j * 6)) = [sum((df.Outcome - df.PredictionP).^2), ...
                                                                    sum((df.Outcome - df.PredictionL).^2), ...
                                                                    sum((df.Outcome - df.PredictionT).^2), ...
                                                                    sum((df.Outcome - df.PredictionPL).^2), ...
                                                                    sum((df.Outcome - df.PredictionPT).^2), ...
                                                                    sum((df.Outcome - df.PredictionLT).^2)];
                    end

                    altSubjectData(i, 1:37) = [included_subjects{i}, parametersPerCondition', SSPerCondition'];
                end

                colnames_altSubjectData = {'SubjectID', 'thetaPLMerit', 'thetaPTMerit', 'phiLTMerit', 'thetaPLEntitlement', 'thetaPTEntitlement', 'phiLTEntitlement', ...
                                            'thetaPLCorruption', 'thetaPTCorruption', 'phiLTCorruption', 'thetaPLLuck', 'thetaPTLuck', 'phiLTLuck', ...
                                            'SSPMerit', 'SSLMerit', 'SSTMerit', 'SSPLMerit', 'SSPTMerit', 'SSLTMerit', ...
                                            'SSPEntitlement', 'SSLEntitlement', 'SSTEntitlement', 'SSPLEntitlement', 'SSPTEntitlement', 'SSLTEntitlement', ...
                                            'SSPCorruption', 'SSLCorruption', 'SSTCorruption', 'SSPLCorruption', 'SSPTCorruption', 'SSLTCorruption', ...
                                            'SSPLuck', 'SSLLuck', 'SSTLuck', 'SSPLLuck', 'SSPLuck', 'SSLTLuck'};

                colnames(altSubjectData) = colnames_altSubjectData;

        .. tab-item:: Python

            ::

                for i in range(len(included_subjects)):
                    datafile = parentfolder + included_subjects[i] + restoffilepath
                    fullData = pd.read_csv(datafile)

                    parametersPerCondition = np.zeros(len(conditions) * 3)
                    SSPerCondition = np.zeros(len(conditions) * 6)

                    for j in range(len(conditions)):
                        df = fullData[fullData['condition'] == conditions[j]].iloc[:, [48] + list(range(39, 48)) + [32]]
                        df['redistributionRate'] = df['redistributionRate'] / 100

                        resultPL = minimize(lambda params: obj_functionPL(params, df), 0.5, bounds=[(0, 1)])
                        resultPT = minimize(lambda params: obj_functionPT(params, df), 0.5, bounds=[(0, 1)])
                        resultLT = minimize(lambda params: obj_functionLT(params, df), 0.5, bounds=[(0, 1)])

                        parametersPerCondition[((j - 1) * 3):(j * 3)] = [resultPL.x[0], resultPT.x[0], resultLT.x[0]]

                        df['PredictionP'] = np.nan
                        df['PredictionL'] = np.nan
                        df['PredictionT'] = np.nan
                        df['PredictionPT'] = np.nan
                        df['PredictionPL'] = np.nan
                        df['PredictionLT'] = np.nan

                        for k in range(len(df['Decisions'])):
                            UtilityP = np.zeros(len(Choices))
                            UtilityL = np.zeros(len(Choices))
                            UtilityT = np.zeros(len(Choices))
                            UtilityPT = np.zeros(len(Choices))
                            UtilityPL = np.zeros(len(Choices))
                            UtilityLT = np.zeros(len(Choices))

                            for n in range(len(Choices)):
                                UtilityP[n] = utilityP(Payout=payout(new_value(df.iloc[k, 0], Choices[n]), df.iloc[k, 0], Choices[n]))
                                UtilityL[n] = utilityL(Equality=equality(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]))
                                UtilityT[n] = utilityT(Equity=equity(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]))
                                UtilityPT[n] = utilityPT(theta=resultPT.x[0], Equity=equity(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]),
                                                        Payout=payout(new_value(df.iloc[k, 0], Choices[n]), df.iloc[k, 0], Choices[n]))
                                UtilityPL[n] = utilityPL(theta=resultPL.x[0], Equality=equality(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]),
                                                        Payout=payout(new_value(df.iloc[k, 0], Choices[n]), df.iloc[k, 0], Choices[n]))
                                UtilityLT[n] = utilityLT(phi=resultLT.x[0], Equality=equality(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]),
                                                        Equity=equity(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]))

                            correct_choiceP = np.argmax(UtilityP)
                            correct_choiceL = np.argmax(UtilityL)
                            correct_choiceT = np.argmax(UtilityT)
                            correct_choicePT = np.argmax(UtilityPT)
                            correct_choicePL = np.argmax(UtilityPL)
                            correct_choiceLT = np.argmax(UtilityLT)

                            df.at[k, 'PredictionP'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceP])
                            df.at[k, 'PredictionL'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceL])
                            df.at[k, 'PredictionT'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceT])
                            df.at[k, 'PredictionPL'] = new_value(df.iloc[k, 'myself'], Choices[correct_choicePT])
                            df.at[k, 'PredictionPT'] = new_value(df.iloc[k, 'myself'], Choices[correct_choicePL])
                            df.at[k, 'PredictionLT'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceLT])

                        df['Outcome'] = new_value(df['myself'], df['redistributionRate'])
                        SSPerCondition[((j - 1) * 6):(j * 6)] = [np.sum((df['Outcome'] - df['PredictionP'])**2),
                                                                np.sum((df['Outcome'] - df['PredictionL'])**2),
                                                                np.sum((df['Outcome'] - df['PredictionT'])**2),
                                                                np.sum((df['Outcome'] - df['PredictionPL'])**2),
                                                                np.sum((df['Outcome'] - df['PredictionPT'])**2),
                                                                np.sum((df['Outcome'] - df['PredictionLT'])**2)]

                    altSubjectData[i, 0:36] = np.concatenate(([included_subjects[i]], parametersPerCondition, SSPerCondition))

                colnames = ['SubjectID', 'thetaPLMerit', 'thetaPTMerit', 'phiLTMerit', 'thetaPLEntitlement', 'thetaPTEntitlement', 'phiLTEntitlement',
                            'thetaPLCorruption', 'thetaPTCorruption', 'phiLTCorruption', 'thetaPLLuck', 'thetaPTLuck', 'phiLTLuck',
                            'SSPMerit', 'SSLMerit', 'SSTMerit', 'SSPLMerit', 'SSPTMerit', 'SSLTMerit',
                            'SSPEntitlement', 'SSLEntitlement', 'SSTEntitlement', 'SSPLEntitlement', 'SSPTEntitlement', 'SSLTEntitlement',
                            'SSPCorruption', 'SSLCorruption', 'SSTCorruption', 'SSPLCorruption', 'SSPTCorruption', 'SSLTCorruption',
                            'SSPLuck', 'SSLLuck', 'SSTLuck', 'SSPLLuck', 'SSPLuck', 'SSLTLuck']

                altSubjectData = pd.DataFrame(altSubjectData, columns=colnames)

.. dropdown:: Compute Model Fit Index for Each :bdg-success:`Subject`, for Each Alternative Model Per :bdg-primary:`Condition`

    .. tab-set::

        .. tab-item:: R

            ::

                N = length(df[, 1])
                k = 2

                altSubjectData$AICPMerit = N * log(altSubjectData$SSPMerit/N) + 2*k
                altSubjectData$AICLMerit = N * log(altSubjectData$SSLMerit/N) + 2*k
                altSubjectData$AICTMerit = N * log(altSubjectData$SSTMerit/N) + 2*k
                altSubjectData$AICPLMerit = N * log(altSubjectData$SSPLMerit/N) + 2*k
                altSubjectData$AICPTMerit = N * log(altSubjectData$SSPTMerit/N) + 2*k
                altSubjectData$AICLTMerit = N * log(altSubjectData$SSLTMerit/N) + 2*k

                altSubjectData$AICPEntitlement = N * log(altSubjectData$SSPEntitlement/N) + 2*k
                altSubjectData$AICLEntitlement = N * log(altSubjectData$SSLEntitlement/N) + 2*k
                altSubjectData$AICTEntitlement = N * log(altSubjectData$SSTEntitlement/N) + 2*k
                altSubjectData$AICPLEntitlement = N * log(altSubjectData$SSPLEntitlement/N) + 2*k
                altSubjectData$AICPTEntitlement = N * log(altSubjectData$SSPTEntitlement/N) + 2*k
                altSubjectData$AICLTEntitlement = N * log(altSubjectData$SSLTEntitlement/N) + 2*k
                
                altSubjectData$AICPCorruption = N * log(altSubjectData$SSPCorruption/N) + 2*k
                altSubjectData$AICLCorruption = N * log(altSubjectData$SSLCorruption/N) + 2*k
                altSubjectData$AICTCorruption = N * log(altSubjectData$SSTCorruption/N) + 2*k
                altSubjectData$AICPLCorruption = N * log(altSubjectData$SSPLCorruption/N) + 2*k
                altSubjectData$AICPTCorruption = N * log(altSubjectData$SSPTCorruption/N) + 2*k
                altSubjectData$AICLTCorruption = N * log(altSubjectData$SSLTCorruption/N) + 2*k
                
                altSubjectData$AICPLuck = N * log(altSubjectData$SSPLuck/N) + 2*k
                altSubjectData$AICLLuck = N * log(altSubjectData$SSLLuck/N) + 2*k
                altSubjectData$AICTLuck = N * log(altSubjectData$SSTLuck/N) + 2*k
                altSubjectData$AICPLLuck = N * log(altSubjectData$SSPLLuck/N) + 2*k
                altSubjectData$AICPTLuck = N * log(altSubjectData$SSPTLuck/N) + 2*k
                altSubjectData$AICLTLuck = N * log(altSubjectData$SSLTLuck/N) + 2*k

        .. tab-item:: MatLab

            ::

                N = length(df(:, 1));
                k = 2;

                altSubjectData.AICPMerit = N * log(altSubjectData.SSPMerit / N) + 2 * k;
                altSubjectData.AICLMerit = N * log(altSubjectData.SSLMerit / N) + 2 * k;
                altSubjectData.AICTMerit = N * log(altSubjectData.SSTMerit / N) + 2 * k;
                altSubjectData.AICPLMerit = N * log(altSubjectData.SSPLMerit / N) + 2 * k;
                altSubjectData.AICPTMerit = N * log(altSubjectData.SSPTMerit / N) + 2 * k;
                altSubjectData.AICLTMerit = N * log(altSubjectData.SSLTMerit / N) + 2 * k;

                altSubjectData.AICPEntitlement = N * log(altSubjectData.SSPEntitlement / N) + 2 * k;
                altSubjectData.AICLEntitlement = N * log(altSubjectData.SSLEntitlement / N) + 2 * k;
                altSubjectData.AICTEntitlement = N * log(altSubjectData.SSTEntitlement / N) + 2 * k;
                altSubjectData.AICPLEntitlement = N * log(altSubjectData.SSPLEntitlement / N) + 2 * k;
                altSubjectData.AICPTEntitlement = N * log(altSubjectData.SSPTEntitlement / N) + 2 * k;
                altSubjectData.AICLTEntitlement = N * log(altSubjectData.SSLTEntitlement / N) + 2 * k;

                altSubjectData.AICPCorruption = N * log(altSubjectData.SSPCorruption / N) + 2 * k;
                altSubjectData.AICLCorruption = N * log(altSubjectData.SSLCorruption / N) + 2 * k;
                altSubjectData.AICTCorruption = N * log(altSubjectData.SSTCorruption / N) + 2 * k;
                altSubjectData.AICPLCorruption = N * log(altSubjectData.SSPLCorruption / N) + 2 * k;
                altSubjectData.AICPTCorruption = N * log(altSubjectData.SSPTCorruption / N) + 2 * k;
                altSubjectData.AICLTCorruption = N * log(altSubjectData.SSLTCorruption / N) + 2 * k;

                altSubjectData.AICPLuck = N * log(altSubjectData.SSPLuck / N) + 2 * k;
                altSubjectData.AICLLuck = N * log(altSubjectData.SSLLuck / N) + 2 * k;
                altSubjectData.AICTLuck = N * log(altSubjectData.SSTLuck / N) + 2 * k;
                altSubjectData.AICPLLuck = N * log(altSubjectData.SSPLLuck / N) + 2 * k;
                altSubjectData.AICPTLuck = N * log(altSubjectData.SSPTLuck / N) + 2 * k;
                altSubjectData.AICLTLuck = N * log(altSubjectData.SSLTLuck / N) + 2 * k;

        .. tab-item:: Python

            ::

                N = len(df)
                k = 2

                altSubjectData['AICPMerit'] = N * np.log(altSubjectData['SSPMerit'] / N) + 2 * k
                altSubjectData['AICLMerit'] = N * np.log(altSubjectData['SSLMerit'] / N) + 2 * k
                altSubjectData['AICTMerit'] = N * np.log(altSubjectData['SSTMerit'] / N) + 2 * k
                altSubjectData['AICPLMerit'] = N * np.log(altSubjectData['SSPLMerit'] / N) + 2 * k
                altSubjectData['AICPTMerit'] = N * np.log(altSubjectData['SSPTMerit'] / N) + 2 * k
                altSubjectData['AICLTMerit'] = N * np.log(altSubjectData['SSLTMerit'] / N) + 2 * k

                altSubjectData['AICPEntitlement'] = N * np.log(altSubjectData['SSPEntitlement'] / N) + 2 * k
                altSubjectData['AICLEntitlement'] = N * np.log(altSubjectData['SSLEntitlement'] / N) + 2 * k
                altSubjectData['AICTEntitlement'] = N * np.log(altSubjectData['SSTEntitlement'] / N) + 2 * k
                altSubjectData['AICPLEntitlement'] = N * np.log(altSubjectData['SSPLEntitlement'] / N) + 2 * k
                altSubjectData['AICPTEntitlement'] = N * np.log(altSubjectData['SSPTEntitlement'] / N) + 2 * k
                altSubjectData['AICLTEntitlement'] = N * np.log(altSubjectData['SSLTEntitlement'] / N) + 2 * k

                altSubjectData['AICPCorruption'] = N * np.log(altSubjectData['SSPCorruption'] / N) + 2 * k
                altSubjectData['AICLCorruption'] = N * np.log(altSubjectData['SSLCorruption'] / N) + 2 * k
                altSubjectData['AICTCorruption'] = N * np.log(altSubjectData['SSTCorruption'] / N) + 2 * k
                altSubjectData['AICPLCorruption'] = N * np.log(altSubjectData['SSPLCorruption'] / N) + 2 * k
                altSubjectData['AICPTCorruption'] = N * np.log(altSubjectData['SSPTCorruption'] / N) + 2 * k
                altSubjectData['AICLTCorruption'] = N * np.log(altSubjectData['SSLTCorruption'] / N) + 2 * k

                altSubjectData['AICPLuck'] = N * np.log(altSubjectData['SSPLuck'] / N) + 2 * k
                altSubjectData['AICLLuck'] = N * np.log(altSubjectData['SSLLuck'] / N) + 2 * k
                altSubjectData['AICTLuck'] = N * np.log(altSubjectData['SSTLuck'] / N) + 2 * k
                altSubjectData['AICPLLuck'] = N * np.log(altSubjectData['SSPLLuck'] / N) + 2 * k
                altSubjectData['AICPTLuck'] = N * np.log(altSubjectData['SSPTLuck'] / N) + 2 * k
                altSubjectData['AICLTLuck'] = N * np.log(altSubjectData['SSLTLuck'] / N) + 2 * k

.. dropdown:: Compare Model Performance

    .. tab-set::

        .. tab-item:: R

            ::

                excludedM = which(is.infinite(subjectData$AICMerit)) #any perfectly predicted subjects by alternative models should also be perfectly predicted by the most complex model as well
                excludedE = which(is.infinite(subjectData$AICEntitlement))
                excludedC = which(is.infinite(subjectData$AICCorruption))
                excludedL = which(is.infinite(subjectData$AICLuck))

                models = c('Favored Model', 'Equality Only', 'Equity Only', 'Payout Only', 'Payout-Equality', 'Payout-Equity', 'Equality-Equity')

                bestModelMerit = models[which(c(mean(subjectData$AICMerit[-excluded]),
                                                mean(altSubjectData$AICPMerit[-excluded]),
                                                mean(altSubjectData$AICLMerit[-excluded]),
                                                mean(altSubjectData$AICTMerit[-excluded]),
                                                mean(altSubjectData$AICPLMerit[-excluded]),
                                                mean(altSubjectData$AICPTMerit[-excluded]),
                                                mean(altSubjectData$AICLTMerit[-excluded])) == max(c(mean(subjectData$AICMerit[-excluded]),
                                                                                                     mean(altSubjectData$AICPMerit[-excluded]),
                                                                                                     mean(altSubjectData$AICLMerit[-excluded]),
                                                                                                     mean(altSubjectData$AICTMerit[-excluded]),
                                                                                                     mean(altSubjectData$AICPLMerit[-excluded]),
                                                                                                     mean(altSubjectData$AICPTMerit[-excluded]),
                                                                                                     mean(altSubjectData$AICLTMerit[-excluded]))))]

                bestModelEntitlement = models[which(c(mean(subjectData$AICEntitlement[-excluded]),
                                                      mean(altSubjectData$AICPEntitlement[-excluded]),
                                                      mean(altSubjectData$AICLEntitlement[-excluded]),
                                                      mean(altSubjectData$AICTEntitlement[-excluded]),
                                                      mean(altSubjectData$AICPLEntitlement[-excluded]),
                                                      mean(altSubjectData$AICPTEntitlement[-excluded]),
                                                      mean(altSubjectData$AICLTEntitlement[-excluded])) == max(c(mean(subjectData$AICEntitlement[-excluded]),
                                                                                                                 mean(altSubjectData$AICPEntitlement[-excluded]),
                                                                                                                 mean(altSubjectData$AICLEntitlement[-excluded]),
                                                                                                                 mean(altSubjectData$AICTEntitlement[-excluded]),
                                                                                                                 mean(altSubjectData$AICPLEntitlement[-excluded]),
                                                                                                                 mean(altSubjectData$AICPTEntitlement[-excluded]),
                                                                                                                 mean(altSubjectData$AICLTEntitlement[-excluded]))))]

                bestModelCorruption = models[which(c(mean(subjectData$AICCorruption[-excluded]),
                                                     mean(altSubjectData$AICPCorruption[-excluded]),
                                                     mean(altSubjectData$AICLCorruption[-excluded]),
                                                     mean(altSubjectData$AICTCorruption[-excluded]),
                                                     mean(altSubjectData$AICPLCorruption[-excluded]),
                                                     mean(altSubjectData$AICPTCorruption[-excluded]),
                                                     mean(altSubjectData$AICLTCorruption[-excluded])) == max(c(mean(subjectData$AICCorruption[-excluded]),
                                                                                                               mean(altSubjectData$AICPCorruption[-excluded]),
                                                                                                               mean(altSubjectData$AICLCorruption[-excluded]),
                                                                                                               mean(altSubjectData$AICTCorruption[-excluded]),
                                                                                                               mean(altSubjectData$AICPLCorruption[-excluded]),
                                                                                                               mean(altSubjectData$AICPTCorruption[-excluded]),
                                                                                                               mean(altSubjectData$AICLTCorruption[-excluded]))))]

                bestModelLuck = models[which(c(mean(subjectData$AICLuck[-excluded]),
                                               mean(altSubjectData$AICPLuck[-excluded]),
                                               mean(altSubjectData$AICLLuck[-excluded]),
                                               mean(altSubjectData$AICTLuck[-excluded]),
                                               mean(altSubjectData$AICPLLuck[-excluded]),
                                               mean(altSubjectData$AICPTLuck[-excluded]),
                                               mean(altSubjectData$AICLTLuck[-excluded])) == max(c(mean(subjectData$AICLuck[-excluded]),
                                                                                                   mean(altSubjectData$AICPLuck[-excluded]),
                                                                                                   mean(altSubjectData$AICLLuck[-excluded]),
                                                                                                   mean(altSubjectData$AICTLuck[-excluded]),
                                                                                                   mean(altSubjectData$AICPLLuck[-excluded]),
                                                                                                   mean(altSubjectData$AICPTLuck[-excluded]),
                                                                                                   mean(altSubjectData$AICLTLuck[-excluded]))))]

        .. tab-item:: MatLab

            ::

                excludedM = find(isinf(subjectData.AICMerit)); % any perfectly predicted subjects by alternative models should also be perfectly predicted by the most complex model as well
                excludedE = find(isinf(subjectData.AICEntitlement));
                excludedC = find(isinf(subjectData.AICCorruption));
                excludedL = find(isinf(subjectData.AICLuck));

                models = {'Favored Model', 'Equality Only', 'Equity Only', 'Payout Only', 'Payout-Equality', 'Payout-Equity', 'Equality-Equity'};

                meanAICMerit = mean(subjectData.AICMerit(~excludedM));
                meanAICPEMerit = mean(altSubjectData.AICPMerit(~excludedM));
                meanAICLEMerit = mean(altSubjectData.AICLMerit(~excludedM));
                meanAICTEMerit = mean(altSubjectData.AICTMerit(~excludedM));
                meanAICPLEMerit = mean(altSubjectData.AICPLMerit(~excludedM));
                meanAICPTEMerit = mean(altSubjectData.AICPTMerit(~excludedM));
                meanAICLTEMerit = mean(altSubjectData.AICLTMerit(~excludedM));

                meanAICs = [meanAICMerit, meanAICPEMerit, meanAICLEMerit, meanAICTEMerit, meanAICPLEMerit, meanAICPTEMerit, meanAICLTEMerit];

                [~, bestModelIndex] = max(meanAICs);
                bestModelMerit = models{bestModelIndex};meanAICEntitlement = mean(subjectData.AICEntitlement(~excludedE));
                meanAICPEEntitlement = mean(altSubjectData.AICPEntitlement(~excludedE));
                meanAICLEEntitlement = mean(altSubjectData.AICLEntitlement(~excludedE));
                meanAICTEEntitlement = mean(altSubjectData.AICTEntitlement(~excludedE));
                meanAICPLEEntitlement = mean(altSubjectData.AICPLEntitlement(~excludedE));
                meanAICPTEEntitlement = mean(altSubjectData.AICPTEntitlement(~excludedE));
                meanAICLTEEntitlement = mean(altSubjectData.AICLTEntitlement(~excludedE));

                meanAICsEntitlement = [meanAICEntitlement, meanAICPEEntitlement, meanAICLEEntitlement, meanAICTEEntitlement, meanAICPLEEntitlement, meanAICPTEEntitlement, meanAICLTEEntitlement];

                [~, bestModelIndexEntitlement] = max(meanAICsEntitlement);
                bestModelEntitlement = models{bestModelIndexEntitlement};

                meanAICCorruption = mean(subjectData.AICCorruption(~excludedC));
                meanAICPCorruption = mean(altSubjectData.AICPCorruption(~excludedC));
                meanAICLCorruption = mean(altSubjectData.AICLCorruption(~excludedC));
                meanAICTCorruption = mean(altSubjectData.AICTCorruption(~excludedC));
                meanAICPLCorruption = mean(altSubjectData.AICPLCorruption(~excludedC));
                meanAICPTCorruption = mean(altSubjectData.AICPTCorruption(~excludedC));
                meanAICLTCorruption = mean(altSubjectData.AICLTCorruption(~excludedC));

                meanAICsCorruption = [meanAICCorruption, meanAICPCorruption, meanAICLCorruption, meanAICTCorruption, meanAICPLCorruption, meanAICPTCorruption, meanAICLTCorruption];

                [~, bestModelIndexCorruption] = max(meanAICsCorruption);
                bestModelCorruption = models{bestModelIndexCorruption};

                meanAICLuck = mean(subjectData.AICLuck(~excludedL));
                meanAICPLuck = mean(altSubjectData.AICPLuck(~excludedL));
                meanAICLLuck = mean(altSubjectData.AICLLuck(~excludedL));
                meanAICTLuck = mean(altSubjectData.AICTLuck(~excludedL));
                meanAICPLLuck = mean(altSubjectData.AICPLLuck(~excludedL));
                meanAICPTLuck = mean(altSubjectData.AICPTLuck(~excludedL));
                meanAICLTLuck = mean(altSubjectData.AICLTLuck(~excludedL));

                meanAICsLuck = [meanAICLuck, meanAICPLuck, meanAICLLuck, meanAICTLuck, meanAICPLLuck, meanAICPTLuck, meanAICLTLuck];

                [~, bestModelIndexLuck] = max(meanAICsLuck);
                bestModelLuck = models{bestModelIndexLuck};

        .. tab-item:: Python

            ::

                # Find indices of subjects with infinite AIC values
                excludedM = np.where(np.isinf(subjectData['AICMerit']))[0]
                excludedE = np.where(np.isinf(subjectData['AICEntitlement']))[0]
                excludedC = np.where(np.isinf(subjectData['AICCorruption']))[0]
                excludedL = np.where(np.isinf(subjectData['AICLuck']))[0]

                models = ['Favored Model', 'Equality Only', 'Equity Only', 'Payout Only', 'Payout-Equality', 'Payout-Equity', 'Equality-Equity']

                # Calculate means for each model and find the best model for each category
                bestModelMerit = models[np.argmax([np.mean(subjectData['AICMerit'][-excludedM]),
                                                   np.mean(altSubjectData['AICPMerit'][-excludedM]),
                                                   np.mean(altSubjectData['AICLMerit'][-excludedM]),
                                                   np.mean(altSubjectData['AICTMerit'][-excludedM]),
                                                   np.mean(altSubjectData['AICPLMerit'][-excludedM]),
                                                   np.mean(altSubjectData['AICPTMerit'][-excludedM]),
                                                   np.mean(altSubjectData['AICLTMerit'][-excludedM])])]

                bestModelEntitlement = models[np.argmax([np.mean(subjectData['AICEntitlement'][-excludedE]),
                                                         np.mean(altSubjectData['AICPEntitlement'][-excludedE]),
                                                         np.mean(altSubjectData['AICLEntitlement'][-excludedE]),
                                                         np.mean(altSubjectData['AICTEntitlement'][-excludedE]),
                                                         np.mean(altSubjectData['AICPLEntitlement'][-excludedE]),
                                                         np.mean(altSubjectData['AICPTEntitlement'][-excludedE]),
                                                         np.mean(altSubjectData['AICLTEntitlement'][-excludedE])])]

                bestModelCorruption = models[np.argmax([np.mean(subjectData['AICCorruption'][-excludedC]),
                                                        np.mean(altSubjectData['AICPCorruption'][-excludedC]),
                                                        np.mean(altSubjectData['AICLCorruption'][-excludedC]),
                                                        np.mean(altSubjectData['AICTCorruption'][-excludedC]),
                                                        np.mean(altSubjectData['AICPLCorruption'][-excludedC]),
                                                        np.mean(altSubjectData['AICPTCorruption'][-excludedC]),
                                                        np.mean(altSubjectData['AICLTCorruption'][-excludedC])])]

                bestModelLuck = models[np.argmax([np.mean(subjectData['AICLuck'][-excludedL]),
                                                  np.mean(altSubjectData['AICPLuck'][-excludedL]),
                                                  np.mean(altSubjectData['AICLLuck'][-excludedL]),
                                                  np.mean(altSubjectData['AICTLuck'][-excludedL]),
                                                  np.mean(altSubjectData['AICPLLuck'][-excludedL]),
                                                  np.mean(altSubjectData['AICPTLuck'][-excludedL]),
                                                  np.mean(altSubjectData['AICLTLuck'][-excludedL])])]



Tutorial 3 - Crockett et al., 2014
-------------------

.. dropdown:: Preallocating and Defining Functions

    .. tab-set::

        .. tab-item:: R

            ::

                #adjusting
                optimize = function(of, using, df){
                    tryCatch({
                        fmincon(of, x0 = initial_params[using], lb = lower_bounds[using], ub = upper_bounds[using], df = df, optimMethod = "MLE")
                    }, error = function(e){
                        fmincon(of, x0 = initial_params[using], lb = lower_bounds[using], ub = upper_bounds[using], df = df, optimMethod = "OLS")
                    })
                }

                obj_functionNoGamma = function(params, df, optimMethod = "MLE") { 
                    #we want kappa and beta per condition, so we'll assign those based on the condition
                    #we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
                    Epsilon = params[5]
                    Gamma = 0 #0 means no bias, everything else is the same
                    
                    ProbHarm = vector('numeric', length(df[,1]))
                    choices = as.numeric(df[, 11])
                    for (k in 1:length(df[,1])){
                        isLeft = as.logical(df[k, 8])
                        isSelf = as.logical(df[k, 10])
                        if (isSelf){Kappa = params[1]; Beta = params[3]} else {Kappa = params[2]; Beta = params[4]}
                        moneyMore = as.numeric(df[k, 4])
                        moneyLess = as.numeric(df[k, 2])
                        shocksMore = as.numeric(df[k, 7])
                        shocksLess = as.numeric(df[k, 5])
                        
                        Utility = utility(Harm = harm(shocksMore, shocksLess),
                                        Payout = payout(moneyMore, moneyLess),
                                        kappa = Kappa)
                        ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
                    }
                    if (optimMethod == "OLS"){
                        return(sum((choices - ProbHarm)**2))
                    } else if (optimMethod == "MLE"){
                        return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
                    }
                }
                obj_functionNoEpsilon = function(params, df, optimMethod = "MLE") { 
                    #we want kappa and beta per condition, so we'll assign those based on the condition
                    #we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
                    Epsilon = 0 #0 means no noise
                    Gamma = 0 #and 0 means no bias, everything else is the same
                    
                    ProbHarm = vector('numeric', length(df[,1]))
                    choices = as.numeric(df[, 11])
                    for (k in 1:length(df[,1])){
                        isLeft = as.logical(df[k, 8])
                        isSelf = as.logical(df[k, 10])
                        if (isSelf){Kappa = params[1]; Beta = params[3]} else {Kappa = params[2]; Beta = params[4]}
                        moneyMore = as.numeric(df[k, 4])
                        moneyLess = as.numeric(df[k, 2])
                        shocksMore = as.numeric(df[k, 7])
                        shocksLess = as.numeric(df[k, 5])
                        
                        Utility = utility(Harm = harm(shocksMore, shocksLess),
                                        Payout = payout(moneyMore, moneyLess),
                                        kappa = Kappa)
                        ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
                    }
                    if (optimMethod == "OLS"){
                        return(sum((choices - ProbHarm)**2))
                    } else if (optimMethod == "MLE"){
                        return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
                    }
                }
                obj_functionNoConditions = function(params, df, optimMethod = "MLE") { 
                    #we want kappa and beta per condition, so we'll assign those based on the condition
                    #we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
                    Kappa = params[1] #moved up here because not changing depending on trial
                    Beta = params[2] #same as above
                    Epsilon = params[3] #change from 5 since 2 fewer parameters
                    Gamma = params[4] #change from 6 for same reason
                    
                    ProbHarm = vector('numeric', length(df[,1]))
                    choices = as.numeric(df[, 11])
                    for (k in 1:length(df[,1])){
                        isLeft = as.logical(df[k, 8])
                        isSelf = as.logical(df[k, 10])
                        moneyMore = as.numeric(df[k, 4])
                        moneyLess = as.numeric(df[k, 2])
                        shocksMore = as.numeric(df[k, 7])
                        shocksLess = as.numeric(df[k, 5])
                        
                        Utility = utility(Harm = harm(shocksMore, shocksLess),
                                        Payout = payout(moneyMore, moneyLess),
                                        kappa = Kappa)
                        ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
                    }
                    if (optimMethod == "OLS"){
                        return(sum((choices - ProbHarm)**2))
                    } else if (optimMethod == "MLE"){
                        return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
                    }
                }
                obj_functionNoKappa = function(params, df, optimMethod = "MLE") { 
                    #we want kappa and beta per condition, so we'll assign those based on the condition
                    #we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
                    Kappa = 0 #nonfactor
                    Beta = 0 #nonfactor
                    Epsilon = 0.5 #the only factor is noise and bias
                    Gamma = params[1] #bias can vary
                    
                    ProbHarm = vector('numeric', length(df[,1]))
                    choices = as.numeric(df[, 11])
                    for (k in 1:length(df[,1])){
                        isLeft = as.logical(df[k, 8])
                        isSelf = as.logical(df[k, 10])
                        moneyMore = as.numeric(df[k, 4])
                        moneyLess = as.numeric(df[k, 2])
                        shocksMore = as.numeric(df[k, 7])
                        shocksLess = as.numeric(df[k, 5])
                        
                        Utility = utility(Harm = harm(shocksMore, shocksLess),
                                        Payout = payout(moneyMore, moneyLess),
                                        kappa = Kappa)
                        ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
                    }
                    if (optimMethod == "OLS"){
                        return(sum((choices - ProbHarm)**2))
                    } else if (optimMethod == "MLE"){
                        return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
                    }
                }
                NoGamma = seq(1, 5)
                NoEpsilon = seq(1, 4)
                NoConditions = c(1, 3, 5, 6)
                NoKappa = 6

                altSubjectData = data.frame()

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python

            ::
.. dropdown:: Recover :bdg-success:`Free Parameters` and Determine Predicted :bdg-danger:`Decisions` 

    .. tab-set::

        .. tab-item:: R

            ::

                for (i in 1:length(included_subjects)){
                    df = grab_data(included_subjects[i])
                    if (is.character(df)){next}
                    
                    resultNoGamma = optimize(obj_functionNoGamma, NoGamma, df)
                    resultNoEpsilon = optimize(obj_functionNoEpsilon, NoEpsilon, df)
                    resultNoConditions = optimize(obj_functionNoConditions, NoConditions, df)
                    resultNoKappa = optim(obj_functionNoKappa,par = initial_params[NoKappa], lower = lower_bounds[NoKappa], upper = upper_bounds[NoKappa], df = df, method = 'L-BFGS-B')
                    
                    dfNoGamma = generate_predictions(df, data.frame(par = c(resultNoGamma$par[1:5], 0)))
                    dfNoEpsilon = generate_predictions(df, data.frame(par = c(resultNoEpsilon$par[1:4], 0, 0)))
                    dfNoConditions = generate_predictions(df, data.frame(par = c(resultNoConditions$par[c(1, 1, 2, 2, 3, 4)])))
                    dfNoKappa = generate_predictions(df, data.frame(par = c(0, 0, 0, 0, 0.5, resultNoKappa$par[1])))
                    
                    modelNG_SS = sum((dfNoGamma$choseHarm - dfNoGamma$ProbHarm)**2)
                    modelNG_NLL = -2*sum(dfNoGamma$choseHarm * log(dfNoGamma$ProbHarm) + (1 - dfNoGamma$choseHarm) * log(1 - dfNoGamma$ProbHarm))
                    modelNE_SS = sum((dfNoEpsilon$choseHarm - dfNoEpsilon$ProbHarm)**2)
                    modelNE_NLL = -2*sum(dfNoEpsilon$choseHarm * log(dfNoEpsilon$ProbHarm) + (1 - dfNoEpsilon$choseHarm) * log(1 - dfNoEpsilon$ProbHarm))
                    modelNC_SS = sum((dfNoConditions$choseHarm - dfNoConditions$ProbHarm)**2)
                    modelNC_NLL = -2*sum(dfNoConditions$choseHarm * log(dfNoConditions$ProbHarm) + (1 - dfNoConditions$choseHarm) * log(1 - dfNoConditions$ProbHarm))
                    modelNK_SS = sum((dfNoKappa$choseHarm - dfNoKappa$ProbHarm)**2)
                    modelNK_NLL = -2*sum(dfNoKappa$choseHarm * log(dfNoKappa$ProbHarm) + (1 - dfNoKappa$choseHarm) * log(1 - dfNoKappa$ProbHarm))
                    
                    altSubjectData[i, 1:23] = c(included_subjects[i], 
                                                resultNoGamma$par[1], resultNoGamma$par[2], resultNoGamma$par[3], resultNoGamma$par[4], resultNoGamma$par[5], 
                                                resultNoEpsilon$par[1], resultNoEpsilon$par[2], resultNoEpsilon$par[3], resultNoEpsilon$par[4], 
                                                resultNoConditions$par[1], resultNoConditions$par[2], resultNoConditions$par[3], resultNoConditions$par[4], 
                                                resultNoKappa$par[1],
                                                modelNG_SS, modelNG_NLL, modelNE_SS, modelNE_NLL, modelNC_SS, modelNC_NLL, modelNK_SS, modelNK_NLL)
                    
                    
                    df = cbind(data.frame(SubjectID = rep(included_subjects[i], length(df$trialNumber))), df)
                    df$ProbHarmNG = dfNoGamma$ProbHarm
                    df$ProbHarmNE = dfNoEpsilon$ProbHarm
                    df$ProbHarmNC = dfNoConditions$ProbHarm
                    df$ProbHarmNK = dfNoKappa$ProbHarm
                    if (i == 1) {altTrialData = df} else {altTrialData = rbind(altTrialData, df)}
                }

                colnames(altSubjectData) = c('SubjectID', 'KappaSelfM2', 'KappaOtherM2', 'BetaSelfM2', 'BetaOtherM2', 'EpsilonM2',
                                            'KappaSelfM3', 'KappaOtherM3', 'BetaSelfM3', 'BetaOtherM3', 
                                            'KappaM4', 'BetaM4', 'EpsilonM4', 'GammaM4', 
                                            'GammaM5', 
                                            'SSM2', 'DevianceM2', 'SSM3', 'DevianceM3', 'SSM4', 'DevianceM4', 'SSM5', 'DevianceM5')

                for (i in 2:ncol(altSubjectData)){altSubjectData[,i] = as.numeric(altSubjectData[,i])}

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python

            ::
.. dropdown:: Compute Model Fit Index for Each :bdg-success:`Subject`, for Each Alternative Model

    .. tab-set::

        .. tab-item:: R

            ::

                altSubjectData$BICM2 = as.numeric(altSubjectData$DevianceM2) + log(152) * 5
                altSubjectData$BICM3 = as.numeric(altSubjectData$DevianceM3) + log(152) * 4
                altSubjectData$BICM4 = as.numeric(altSubjectData$DevianceM4) + log(152) * 4
                altSubjectData$BICM5 = as.numeric(altSubjectData$DevianceM5) + log(152) * 1

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python

            ::
.. dropdown:: Compare Model Performance

    .. tab-set::

        .. tab-item:: R

            ::

                modelBIC = c(sum(subjectData$BIC[-20]), sum(altSubjectData$BICM2[-20]), sum(altSubjectData$BICM3[-20]), sum(altSubjectData$BICM4[-20]), sum(altSubjectData$BICM5[-20]))
                which(modelBIC == min(modelBIC))

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python

            ::

Tutorial 4 - Li et al., 2022
-------------------

.. dropdown:: Preallocating and Defining Functions

    .. tab-set::

        .. tab-item:: R

            ::

                #adjusting

                of_alphaOnly = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(1, 4:6)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_deltaOnly = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(2, 4:6)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_rhoOnly = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(3:6)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_ad = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(1:2, 4:6)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_ar = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(1, 3:6)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_dr = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(2:6)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_noEpsilon = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(1:4)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_noGamma = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[c(1:5)] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }
                of_GammaOnly = function(params, df, optimMethod = "MLE"){
                    params_new = rep(0, 6); params_new[5] = 0.5; params_new[6] = params
                    Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
                    if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2)) 
                    } else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
                }

                ofs = list(of_alphaOnly, of_deltaOnly, of_rhoOnly, of_ad, of_ar, of_dr, of_noEpsilon, of_noGamma, of_GammaOnly)
                idxs = list(c(1,4:6), c(2,4:6), c(3:6), c(1:2, 4:6), c(1, 3:6), c(2:6), c(1:4), c(1:5), c(6))

                altSubjectData = data.frame()
                altTrialData = trialData[,-9]
                altTrialData$alphaOnly_Prob1 = 0
                altTrialData$deltaOnly_Prob1 = 0
                altTrialData$rhoOnly_Prob1 = 0
                altTrialData$ad_Prob1 = 0
                altTrialData$ar_Prob1 = 0
                altTrialData$dr_Prob1 = 0
                altTrialData$noEpsilon_Prob1 = 0
                altTrialData$noGamma_Prob1 = 0
                altTrialData$gammaOnly_Prob1 = 0

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python

            ::
.. dropdown:: Recover :bdg-success:`Free Parameters` and Determine Predicted :bdg-danger:`Decisions`

    .. tab-set::

        .. tab-item:: R

            ::

                for (i in 1:length(included_subjects)){
                df = grab_data(included_subjects[i])
                outputs = vector('numeric', length = sum(length(c(c(1,4:6), c(2,4:6), c(3:6), c(1:2, 4:6), c(1, 3:6), c(2:6), c(1:4), c(1:5), c(6)))) + 2*9)
                j = 0
                for (k in 1:length(idxs)){
                        idx = idxs[k][[1]]; initials = initial_params[idx]; uppers = upper_bounds[idx]; lowers = lower_bounds[idx]; of = ofs[k][[1]]
                        if (length(idx) > 1){
                            result = optimize(of, initials, lowers, uppers, df)
                        } else {
                            result = optim(par = initials, fn = of, method = "L-BFGS-B", lower = lowers, upper = uppers, df = df)
                        }
                        pars = rep(0, times = 6)
                        pars[idx] = result$par
                        df$Prob1 = 0; df$Prob1 = generatePredictions(pars, df)
                        
                        model_SS = sum((df$Chose1 - df$Prob1)**2)
                        model_NLL = -2*sum(df$Chose1 * log(df$Prob1) + (1 - df$Chose1) * log(1 - df$Prob1))
                        outputs[(j+1):(j+2+length(result$par))] = c(result$par, model_SS, model_NLL)
                        j = j+2+length(result$par)
                        altTrialData[which(altTrialData$SubjectID == included_subjects[i]), (8+k)] = df$Prob1
                    }
                
                altSubjectData[i, 1:56] = c(included_subjects[i], outputs)
                }

                colnames(altSubjectData) = c('SubjectID', 
                                            'Alpha_M1', 'Beta_M1', 'Epsilon_M1', 'Gamma_M1', 'SS_M1', 'Deviance_M1',
                                            'Delta_M2', 'Beta_M2', 'Epsilon_M2', 'Gamma_M2', 'SS_M2', 'Deviance_M2', 
                                            'Rho_M3', 'Beta_M3', 'Epsilon_M3', 'Gamma_M3', 'SS_M3', 'Deviance_M3', 
                                            'Alpha_M4', 'Delta_M4', 'Beta_M4', 'Epsilon_M4', 'Gamma_M4', 'SS_M4', 'Deviance_M4',
                                            'Alpha_M5', 'Rho_M5', 'Beta_M5', 'Epsilon_M5', 'Gamma_M5', 'SS_M5', 'Deviance_M5',
                                            'Delta_M6', 'Rho_M6', 'Beta_M6', 'Epsilon_M6', 'Gamma_M6', 'SS_M6', 'Deviance_M6',
                                            'Alpha_M7', 'Delta_M7', 'Rho_M7', 'Beta_M7', 'SS_M7', 'Deviance_M7',
                                            'Alpha_M8', 'Delta_M8', 'Rho_M8', 'Beta_M8', 'Epsilon_M8', 'SS_M8', 'Deviance_M8',
                                            'Gamma_M9', 'SS_M9', 'Deviance_M9')

                for (i in 2:ncol(altSubjectData)){altSubjectData[,i] = as.numeric(altSubjectData[,i])}

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python

            ::
.. dropdown:: Compute Model Fit Index for Each :bdg-success:`Subject`, for Each Alternative Model

    .. tab-set::

        .. tab-item:: R

            ::

                altSubjectData$BIC_M1 = (altSubjectData$Deviance_M1) + log(65) * 4
                altSubjectData$BIC_M2 = (altSubjectData$Deviance_M2) + log(65) * 4
                altSubjectData$BIC_M3 = (altSubjectData$Deviance_M3) + log(65) * 4
                altSubjectData$BIC_M4 = (altSubjectData$Deviance_M4) + log(65) * 5
                altSubjectData$BIC_M5 = (altSubjectData$Deviance_M5) + log(65) * 5
                altSubjectData$BIC_M6 = (altSubjectData$Deviance_M6) + log(65) * 5
                altSubjectData$BIC_M7 = (altSubjectData$Deviance_M7) + log(65) * 4
                altSubjectData$BIC_M8 = (altSubjectData$Deviance_M8) + log(65) * 5
                altSubjectData$BIC_M9 = (altSubjectData$Deviance_M9) + log(65) * 1

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python

            ::
.. dropdown:: Compare Model Performance

    .. tab-set::

        .. tab-item:: R

            ::

                modelBIC = c(sum(subjectData$BIC), sum(altSubjectData$BIC_M1), sum(altSubjectData$BIC_M2), sum(altSubjectData$BIC_M3), sum(altSubjectData$BIC_M4), sum(altSubjectData$BIC_M5), sum(altSubjectData$BIC_M6), sum(altSubjectData$BIC_M7), sum(altSubjectData$BIC_M8), sum(altSubjectData$BIC_M9))
                which(modelBIC == min(modelBIC))

        .. tab-item:: MatLab

            ::

        .. tab-item:: Python