diff --git a/ddToolbox/CODA/CODA.m b/ddToolbox/CODA/CODA.m index 2ab299ef..91d65208 100644 --- a/ddToolbox/CODA/CODA.m +++ b/ddToolbox/CODA/CODA.m @@ -2,23 +2,22 @@ %CODA This class is useful for storing, ploting, and getting %information about MCMC samples. % Detailed explanation goes here - - % TODO: make this general, and create a seperate, project-specific subclass with additional methods if I need to do that - + + % TODO: make this general, and create a seperate, project-specific subclass with additional methods if I need to do that + properties (GetAccess = private, SetAccess = private) samples % structure. Fields correspond to variables. stats % structure. Fields correspond to stats, subfield correspond to variables end - + properties (GetAccess = public, SetAccess = private) variableNames % cell array of variables end - + methods (Access = public) - - % TODO: be able to create just from samples. - % TODO: add a makeFromJAGS alternative constructor, just like we have one for Stan - + + % TODO: be able to create just from samples. + function obj = CODA(samples, stats) % constructor assert(isstruct(samples)) assert(isstruct(stats)) @@ -27,10 +26,10 @@ obj.variableNames = fieldnames(samples); % TODO: Check presence of my mcmc-utils code as the plotting relies upon it. end - + function paramEstimateTable = buildParameterEstimateTable(obj,... variablesRequested, rowNames, savePath, pointEstimateType, varargin) - + p = inputParser; p.FunctionName = mfilename; p.addRequired('variablesRequested',@iscellstr); @@ -41,14 +40,14 @@ p.addParameter('includeCI',false, @islogical); p.parse(variablesRequested, rowNames, savePath, pointEstimateType, varargin{:}); - + % TODO: act on includeCI preference. Ie get, or do not get CI's. - + colHeaderNames = createColumnHeaders(... p.Results.variablesRequested,... p.Results.includeCI,... p.Results.pointEstimateType); - + % TODO: FIX THIS FAFF TO DEAL WITH POSSIBLE VECTOR/MATRIX VARIABLES errorFlag = false; tableEntries = NaN(numel(rowNames), numel(colHeaderNames)); @@ -64,7 +63,7 @@ tableEntries([1:numel(vals)],n) = vals; end end - + if ~errorFlag paramEstimateTable = array2table(tableEntries,... 'VariableNames',colHeaderNames,... @@ -74,7 +73,7 @@ % return an empty table paramEstimateTable= table(); end - + function colHeaderNames = createColumnHeaders(varNames, getCI, pointEstimateType) colHeaderNames = {}; for k = 1:numel(varNames) @@ -86,31 +85,31 @@ end end end - - + + % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %% Plotting methods - % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - + % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + function trellisplots(obj, variablesToPlot) - % plot mcmc chains (left column) and density plots (right column) + % plot mcmc chains (left column) and density plots (right column) assert(iscellstr(variablesToPlot)) - + for n = 1:numel(variablesToPlot) % TODO: REMOVE THIS LOOP BY A MAP ? - + % sort figure out figure latex_fig(16, 12,10) - + % get info mcmcsamples = obj.getSamples(variablesToPlot(n)); mcmcsamples = mcmcsamples.(variablesToPlot{n}); [chains,Nsamples,rows] = size(mcmcsamples); rhat_all = obj.getStats('Rhat', variablesToPlot{n}); - - % TODO #166 Unclear what this code is doing - % TODO #166 Exact this into a method. Can we use layout() or create_subplots() ? + + % TODO #166 Unclear what this code is doing + % TODO #166 Exact this into a method. Can we use layout() or create_subplots() ? % create geometry for row = 1:rows % traceplot @@ -119,50 +118,50 @@ function trellisplots(obj, variablesToPlot) % density plot axis_handle(row, 2) = subplot(rows,6,row*6); end - + % call the plot functions for row = 1:rows obj.traceplot(axis_handle(row, 1),... mcmcsamples(:,:,row),... variablesToPlot{n},... rhat_all(row)); - + obj.densityplot(axis_handle(row, 2),... mcmcsamples(:,:,row)) end - + %x tick only on bottom plot set(axis_handle([1:end-1], 1),'XTick',[]) axis_handle(rows, 1).XLabel.String = 'MCMC sample'; - + % link y-axes of all traceplots linkaxes(axis_handle(:,1),'xy') - + % link x-axis of all density plots linkaxes(axis_handle(:,2),'x') end end - - + + function traceplot(obj, targetAxisHandle, samples, paramString, rhat) - % Plot mcmc chains - + % Plot mcmc chains + % TODO: make targetAxisHandle an optional input - + assert(ischar(paramString)) assert(isscalar(rhat)) assert(ishandle(targetAxisHandle)) assert(size(samples,3)==1) % plot - subplot(targetAxisHandle) + subplot(targetAxisHandle) h = plot(samples', 'LineWidth',0.5); % format - box off + box off ylabel(sprintf('$$ %s $$', paramString), 'Interpreter','latex') if ~isempty(rhat), addRhatStringToFigure(targetAxisHandle, rhat), end end - - + + function densityplot(obj, targetAxisHandle, samples) % TODO: MAKE targetAxisHandle OPTIONAL @@ -173,10 +172,10 @@ function densityplot(obj, targetAxisHandle, samples) if ~isempty(targetAxisHandle) subplot(targetAxisHandle) end - - univariateObject = Stochastic('name_here'); - univariateObject.addSamples(samples); - univariateObject.plot; + + univariateObject = Stochastic('name_here'); + univariateObject.addSamples(samples); + univariateObject.plot; end @@ -201,22 +200,22 @@ function plot_bivariate_distribution(obj, targetAxisHandle, x_var_name, y_var_na 'plotStyle', 'hist',... 'axisSquare', true); end - - function plotUnivariateSummaries(obj, variables, plotOptions, idNames) - % create a figure with multiple subplots (boxplots) for each input variable privided - - % create subplots - N = numel(variables); - figure(8) - subplot_handles = create_subplots(N, 'col'); - - % call plotting function for each variable (subplot) - for n = 1:numel(variables) - subplot( subplot_handles(n) ) + + function plotUnivariateSummaries(obj, variables, plotOptions, idNames) + % create a figure with multiple subplots (boxplots) for each input variable privided + + % create subplots + N = numel(variables); + figure(8) + subplot_handles = create_subplots(N, 'col'); + + % call plotting function for each variable (subplot) + for n = 1:numel(variables) + subplot( subplot_handles(n) ) hold off - - pointEstVals = obj.getStats(plotOptions.pointEstimateType, variables{n}); - hdi =[obj.getStats('hdi_low',variables{n}),... + + pointEstVals = obj.getStats(plotOptions.pointEstimateType, variables{n}); + hdi =[obj.getStats('hdi_low',variables{n}),... obj.getStats('hdi_high',variables{n})]; % clunky manual checking: sometimes we will not have a @@ -224,25 +223,25 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) % off the idNames idNamesCropped = idNames(1:numel(pointEstVals)); - plotErrorBars(idNamesCropped,... - pointEstVals,... - hdi,... - variables{n}); - end - - % % Scale width of figure - % screen_size = get(0,'ScreenSize'); - % fig_width = min(screen_size(3), 100+numel(participantNames)*20); - % set(gcf,'Position',[100 200 fig_width 1000]) - end - - - % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + plotErrorBars(idNamesCropped,... + pointEstVals,... + hdi,... + variables{n}); + end + + % % Scale width of figure + % screen_size = get(0,'ScreenSize'); + % fig_width = min(screen_size(3), 100+numel(participantNames)*20); + % set(gcf,'Position',[100 200 fig_width 1000]) + end + + + % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %% Get methods - % TODO #103 rethink all these get methods. - % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - + % TODO #103 rethink all these get methods. + % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + function data = grabParamEstimates(obj, varNames, getCI, pointEstimateType) assert(islogical(getCI)) data=[]; @@ -254,7 +253,7 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) end end end - + function [samples] = getSamplesAtIndex_asStruct(obj, index, fieldsToGet) assert(iscellstr(fieldsToGet),'arguments needs to be a cell array of strings') assert(isnumeric(index), 'argument needs to be numeric') @@ -263,7 +262,7 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) % 1. mcmc chain number % 2. mcmc sample number % 3. index of variable, meaning depends upon context of the model - + [flatSamples] = obj.flattenChains(obj.samples, fieldsToGet); for n = 1:numel(fieldsToGet) try @@ -302,30 +301,30 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) function [samplesMatrix] = getSamplesAtIndex_asMatrix(obj, index, fieldsToGet) samples = getSamplesAtIndex_asStruct(obj, index, fieldsToGet); -% % Checks -% % We need to check and remove fields for which there are no samples (either empty, or NaN) -% fields = fieldnames(samples); -% for n = 1:numel(fields) -% if any(isnan(samples.(fields{n}))) || isempty(samples.(fields{n})) -% samples = rmfield(samples, fields{n}); -% end -% end - + % % Checks + % % We need to check and remove fields for which there are no samples (either empty, or NaN) + % fields = fieldnames(samples); + % for n = 1:numel(fields) + % if any(isnan(samples.(fields{n}))) || isempty(samples.(fields{n})) + % samples = rmfield(samples, fields{n}); + % end + % end + [samplesMatrix] = struct2Matrix(samples); end - + function [samples] = getSamples(obj, fieldsToGet) % Doesn't flatten across chains assert(iscellstr(fieldsToGet)) fieldsToGet = ismember(fieldnames(obj.samples), fieldsToGet); samples = filterFields(fieldsToGet, obj.samples); end - + function [samplesMatrix] = getSamplesAsMatrix(obj, fieldsToGet) % TODO: this makes assumptions, which are not true. Add checks or robustify. samplesMatrix = struct2Matrix( obj.flattenChains(obj.samples, fieldsToGet) ); end - + function [columnVector] = getStats(obj, field, variable) % TODO: check requested field exists in stats try @@ -338,7 +337,7 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) columnVector =[]; end end - + function pointEstimates = getPointEstimates(obj, n, variableNames) assert(iscellstr(variableNames)) for var = each(variableNames) @@ -346,7 +345,7 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) pointEstimates.(var) = temp(n); end end - + function [predicted] = getParticipantPredictedResponses(obj, ind) % ind is a binary valued vector indicating the trials % corresponding to a particular participant @@ -356,33 +355,24 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) % Calculate predicted response probability predicted = sum(participantRpostpredSamples,1) ./ size(participantRpostpredSamples,1); end - + function [PChooseDelayed] = getPChooseDelayed(obj, pInd) PChooseDelayed = obj.samples.P(:,:,pInd); PChooseDelayed = collapseFirstTwoColumnsOfMatrix(PChooseDelayed)'; end - - end - - %% Alternate constructors - methods (Static) - function obj = buildFromStanFit(stanFitObject) - % TODO: add an assert about the type of object being passed in - samples = stanFitObject.extract('collapseChains', false, 'permuted', false); - stats = computeStats(samples); - obj = CODA(samples, stats); - end end - - + + + + %% PRIVATE METHODS ==================================================== % Not to be covered by tests, unless it is useful during development. % But we do not need tests to constrain the way how these % implementation details work. - + methods(Static, Access = private) - + function [new_samples] = flattenChains(samples, fieldsToGet) assert(isstruct(samples)) assert(iscellstr(fieldsToGet)) @@ -390,7 +380,7 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames) samples = filterFields(fieldsToGet, samples); new_samples = structfun(@collapseFirstTwoColumnsOfMatrix, samples, 'UniformOutput', false); end - + end - + end diff --git a/ddToolbox/Data/DataImporter.m b/ddToolbox/Data/DataImporter.m index c7b97c65..fac261bd 100644 --- a/ddToolbox/Data/DataImporter.m +++ b/ddToolbox/Data/DataImporter.m @@ -92,8 +92,9 @@ newR(t,1) = 1; end end + experimentTable.R = newR; end - experimentTable.R = newR; + %% Add ID column experimentTable = appendTableColOfVals(experimentTable, n); diff --git a/ddToolbox/dependencies.txt b/ddToolbox/dependencies.txt index b09f69a8..a8ff84c8 100644 --- a/ddToolbox/dependencies.txt +++ b/ddToolbox/dependencies.txt @@ -1,5 +1,3 @@ https://github.com/drbenvincent/mcmc-utils-matlab https://github.com/altmany/export_fig https://github.com/drbenvincent/matjags -https://github.com/brian-lau/MatlabProcessManager -https://github.com/drbenvincent/MatlabStan diff --git a/ddToolbox/models/Model.m b/ddToolbox/models/Model.m index bd2b30e9..f3152ea3 100644 --- a/ddToolbox/models/Model.m +++ b/ddToolbox/models/Model.m @@ -1,13 +1,13 @@ classdef (Abstract) Model %Model Base class model. - + % Allow acces to these via Model, but we still only get access to these % class's public interface. properties (Hidden = true, SetAccess = protected, GetAccess = public) coda % handle to coda object data % handle to Data class end - + %% Private properties properties (SetAccess = protected, GetAccess = protected) dfClass % function handle to DiscountFunction class @@ -16,13 +16,13 @@ mcmcParams % structure of user-supplied params observedData % User supplied preferences - modelFilename % string (ie modelFilename.jags, or modelFilename.stan) + modelFilename % string (ie modelFilename.jags) varList plotOptions timeUnits % string whose name must be a function to create a Duration. auc % array of Stochastic objects end - + % methods that subclasses must implement methods (Abstract, Access = public) plot() @@ -32,10 +32,10 @@ methods (Abstract, Access = protected) initialiseChainValues() end - - + + methods (Access = public) - + function obj = Model(data, varargin) % Input parsing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ p = inputParser; @@ -45,7 +45,7 @@ % Required p.addRequired('data', @(x) isa(x,'Data')); % Optional inference related parameters - p.addParameter('samplerType', 'jags', @(x) any(strcmp(x,{'jags','stan'}))); + p.addParameter('samplerType', 'jags', @(x) any(strcmp(x,{'jags'}))); p.addParameter('mcmcParams', struct, @isstruct) % Define the time units. This must correspond to Duration % creation function, such as hours, days, etc. See `help @@ -63,29 +63,29 @@ obj.mcmcParams = obj.parse_mcmcparams(obj.mcmcParams); obj.plotOptions = obj.parse_plot_options(varargin{:}); % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + obj.varList.responseErrorParams(1).name = 'alpha'; obj.varList.responseErrorParams(1).label = 'comparison accuity, $\alpha$'; - + obj.varList.responseErrorParams(2).name = 'epsilon'; obj.varList.responseErrorParams(2).label = 'error rate, $\epsilon$'; end - - + + function export(obj) %model.export Exports summary information about the model % model.EXPORT() exports the following: % - extensive MCMC convergence information % - an experiment-level table of point estimates and model % checking quantities - + % TODO: This should be a method of CODA % TODO: better function name. What does it do? Calculate or % export, or both? convergenceSummary(obj.coda.getStats('Rhat',[]),... obj.plotOptions.savePath,... obj.data.getIDnames('all')) - + exporter = ResultsExporter(obj.coda,... obj.data,... obj.getAUCasTable,... @@ -96,14 +96,14 @@ function export(obj) exporter.export(obj.plotOptions.savePath, obj.plotOptions.pointEstimateType); % TODO ^^^^ avoid this duplicate use of pointEstimateType end - + end - + methods (Hidden = true) - + function obj = conductInference(obj) % Not intended to be called directly by user - + % pre-sampling preparation obj.observedData = obj.constructObservedDataForMCMC(); path_of_model_file = makeProbModelsPath(obj.modelFilename, obj.samplerType); @@ -119,11 +119,11 @@ function export(obj) obj = obj.postSamplingActivities(); end end - + %% GETTERS - + methods - + function [predicted_subjective_values] = getInferredPresentSubjectiveValues(obj) % info = model.getInferredPresentSubjectiveValues Returns information % on the dataset along with inferred present subjective values of @@ -136,38 +136,38 @@ function export(obj) % % NOTE: We do have access to full distributions of inferred present % subjective values, but currently we return the point estimates. - + %% return point estimates of present subjective values... all_data_table = obj.data.groupTable; % add new columns for present subjective value (VA, VB) all_data_table.VA = obj.coda.getStats(obj.plotOptions.pointEstimateType, 'VA'); all_data_table.VB = obj.coda.getStats(obj.plotOptions.pointEstimateType, 'VB'); predicted_subjective_values.point_estimates = all_data_table; - + %% TODO: Return full posterior distributions of present subjective values % predicted_subjective_values.A_full_posterior = % predicted_subjective_values.B_full_posterior = end - + function [auc] = getAUCasStochasticArray(obj) %getAUCasStochasticArray() returns an array of Stochastic objects each describing a posterior distribution of AUC. - + auc = obj.auc; end - + function [T] = getAUCasTable(obj) %getAUCasTable() returns a Table with AUC information. - + stochasticArray = obj.auc; - + T = table; - + % catch cases where we don't have any AUC information (eg with % 2 dimensional discount functions) and return an empty table if isempty(cell2mat({stochasticArray(:).samples})) return end - + %% grab point estimate pointEstimateType = obj.plotOptions.pointEstimateType; T.(['auc_' pointEstimateType]) = [stochasticArray(:).(pointEstimateType)]'; @@ -181,51 +181,51 @@ function export(obj) % possible error with reordering for example. T.Properties.RowNames = obj.data.getIDnames('experiments'); end - + end - + methods (Hidden = true) - + function discountFunctionVariables = getDiscountFunctionVariables(obj) discountFunctionVariables = {obj.varList.discountFunctionParams.name}; end - + function responseErrorVariables = getResponseErrorVariables(obj) responseErrorVariables = {obj.varList.responseErrorParams.name}; end - + end - - + + %% Protected methods - + methods (Access = protected) - + function nChains = get_nChains(obj) nChains = obj.mcmcParams.nchains; end - + function [samples] = getGroupLevelSamples(obj, fieldsToGet) [samples] = obj.data.getGroupLevelSamples(fieldsToGet); end - + function obj = postSamplingActivities(obj) %% Post-sampling activities (for model sub-classes) ----------- % If a model has additional measures that need to be calculated % from the MCMC samples, then we can do by overriding this % method in the model sub-classes obj = obj.calcDerivedMeasures(); - + %% Post-sampling activities (common to all models) ------------ % posterior prediction calculation obj.postPred = PosteriorPrediction(obj.coda,... obj.data,... obj.observedData,... obj.plotOptions.pointEstimateType); - + % calc and export convergence summary and parameter estimates obj.export(); - + % plot or not if ~strcmp(obj.plotOptions.shouldPlot,'no') % TODO: Allow public calls of obj.plot to specify options. @@ -233,10 +233,10 @@ function export(obj) % object construction obj.plot() end - + obj.displayPublicMethods() end - + function observedData = constructObservedDataForMCMC(obj) % This function can be overridden by model subclasses, however % we still expect them to call this model baseclass method to @@ -249,43 +249,43 @@ function export(obj) % protected method which can be over-ridden by model sub-classes observedData = obj.additional_model_specific_ObservedData(observedData); end - + function obj = calcDerivedMeasures(obj) - + %% Calculate AUC MAX_DELAY = 365; obj.auc = obj.calcAreaUnderCurveForAll(MAX_DELAY); - + end - + function auc = calcAreaUnderCurveForAll(obj, MAX_DELAY) % Calculate Area Under Curve. % Returns an array of Stochastic objects - + % TODO: store experiment filenames along with the auc data to % ensure we don't mix up ordering - + N = obj.data.getNRealExperimentFiles(); - + for ind = 1:N % loop over files discountFunctionVariables = {obj.varList.discountFunctionParams.name}; samples = obj.coda.getSamplesAtIndex_asStochastic(ind, discountFunctionVariables); discountFunction = obj.dfClass('samples', samples); - + % grab a distribution over AUC (a Stochastic object) uniqueDelays = obj.data.getUniqueDelaysForThisParticant(ind); auc(ind) = discountFunction.calcAUC(MAX_DELAY, uniqueDelays); end - + % TODO: This is a both solution. We don't calculate AUC for % group level, so if it is present, then return an empty % Stochastic object if obj.data.getNRealExperimentFiles() ~= obj.data.getNExperimentFiles() auc = [auc Stochastic('AUC')]; end - + end - + function displayPublicMethods(obj) display('For more info (assuming your model is named ''model'') then you can type') display(' >> help model.methodName') @@ -294,32 +294,32 @@ function displayPublicMethods(obj) methodsAvailable = methods(obj); disp(methodsAvailable) end - + function obj = addUnobservedParticipant(obj, str) % TODO: Check we need this obj.data = obj.data.add_unobserved_participant(str); % add name (eg 'GROUP') end - + end - - + + methods (Static, Access = protected) - + function observedData = additional_model_specific_ObservedData(observedData) % KEEP THIS HERE. IT IS OVER-RIDDEN IN SOME MODEL SUB-CLASSES - + % TODO: can we move this to NonParamtric abstract class? end - + function samplerFunction = selectSampler(samplerType) switch samplerType case{'jags'} samplerFunction = @sampleWithMatjags; - case{'stan'} - samplerFunction = @sampleWithMatlabStan; + otherwise + error('samplerType must be ''jags'' '); end end - + function mcmcparams = parse_mcmcparams(mcmcParams) defaultMCMCParams.doparallel = 1; defaultMCMCParams.nburnin = 1000; @@ -331,20 +331,20 @@ function displayPublicMethods(obj) end mcmcparams = kwargify(defaultMCMCParams, mcmcParams); end - + end - - - - + + + + % ========================================================================== % ========================================================================== % PLOTTING % ========================================================================== % ========================================================================== - + methods (Access = public) - + function plotDiscountFunction(obj, ind, varargin) %model.PLOTDISCOUNTFUNCTION(N) plots a discount % function where H is a handle to a subplot, and IND is the nth @@ -353,33 +353,33 @@ function plotDiscountFunction(obj, ind, varargin) % Optional arguments as key/value pairs % 'axisHandle' - handle to axes % 'figureHandle' - handle to figure - + parseFigureAndAxisRequested(varargin{:}); - + p = inputParser; p.FunctionName = mfilename; p.KeepUnmatched = true; p.addParameter('plot_mode', 'full', @isstr); p.parse(varargin{:}); - + discountFunctionVariables = {obj.varList.discountFunctionParams.name}; - + samples = obj.coda.getSamplesAtIndex_asStochastic(ind, discountFunctionVariables); - + discountFunction = obj.dfClass(... 'samples', samples,... 'data', obj.data.getExperimentObject(ind)); - + plotOptions.pointEstimateType = obj.plotOptions.pointEstimateType; plotOptions.dataPlotType = obj.plotOptions.dataPlotType; plotOptions.timeUnits = obj.timeUnits; plotOptions.plotMode = p.Results.plot_mode; plotOptions.maxRewardValue = obj.data.getMaxRewardValue(ind); plotOptions.maxDelayValue = obj.data.getMaxDelayValue(ind); - + discountFunction.plot(plotOptions); end - + function plotDiscountFunctionGrid(obj, varargin) %plotDiscountFunctionGrid Plots a montage of discount functions % model.PLOTDISCOUNTFUNCTIONGRID() plots discount functions for @@ -388,24 +388,24 @@ function plotDiscountFunctionGrid(obj, varargin) % Optional arguments as key/value pairs % 'axisHandle' - handle to axes % 'figureHandle' - handle to figure - + [figureHandle, ~] = parseFigureAndAxisRequested(varargin{:}); - + figure(figureHandle), clf latex_fig(12, 11,11) drawnow - + % TODO: extract the grid formatting stuff to be able to call % any plot function we want % USE: apply_plot_function_to_subplot_handle.m ?? - + %fh = figure('Name', names{experimentIndex}); names = obj.data.getIDnames('all'); - + % create grid layout N = numel(names); subplot_handles = create_subplots(N, 'square'); - + % Iterate over files, plotting disp('Plotting...') for n = 1:numel(names) @@ -415,8 +415,8 @@ function plotDiscountFunctionGrid(obj, varargin) end drawnow end - - + + function plotDiscountFunctionsOverlaid(obj, varargin) %plotDiscountFunctionsOverlaid Plots all discount functions in one figure % model.PLOTDISCOUNTFUNCTIONSOVERLAID() plots discount functions for @@ -425,16 +425,16 @@ function plotDiscountFunctionsOverlaid(obj, varargin) % Optional arguments as key/value pairs % 'axisHandle' - handle to axes % 'figureHandle' - handle to figure - + latex_fig(12, 8,6) clf drawnow - + [figureHandle, axisHandle] = parseFigureAndAxisRequested(varargin{:}); - + % don't want the group level estimate, so not asking for 'all' names = obj.data.getIDnames('experiments'); - + % Iterate over files, plotting disp('Plotting...') for n = 1:numel(names) @@ -450,8 +450,8 @@ function plotDiscountFunctionsOverlaid(obj, varargin) % set(gca,'FontSize',10) drawnow end - - + + function plotPosteriorAUC(obj, ind, varargin) %model.plotPosteriorAUC(N, varargin) plots the posterior distribution % of AUC values where N is the Nth experiment. @@ -459,17 +459,17 @@ function plotPosteriorAUC(obj, ind, varargin) % Optional arguments as key/value pairs % 'axisHandle' - handle to axes % 'figureHandle' - handle to figure - + parseFigureAndAxisRequested(varargin{:}) - + obj.auc(ind).plot(); - + % set min of x-axis equal to 0 a = get(gca, 'XLim'); set(gca, 'XLim', [0, a(2)]); end - - + + function plotUnivarateSummary(obj, varargin) %plotUnivarateSummary(varargin) Produces a univariate summary plot % Creates a forest plot of univariate summary stats @@ -477,12 +477,12 @@ function plotUnivarateSummary(obj, varargin) % Optional arguments as key/value pairs % 'axisHandle' - handle to axes % 'figureHandle' - handle to figure - + p = inputParser; p.FunctionName = mfilename; p.addParameter('variablesToPlot', 'all', @(x) isstr(x)|iscellstr(x)); p.parse(varargin{:}); - + if iscellstr(p.Results.variablesToPlot) % user is providing a cell array of strings, which we will assume % corresponds to specific variable names. @@ -496,16 +496,16 @@ function plotUnivarateSummary(obj, varargin) error('unrecognised input provided for ''variablesToPlot''.') end end - + obj.coda.plotUnivariateSummaries(variables,... obj.plotOptions,... obj.data.getParticipantNames()); - + plot_savename = 'UnivariateSummary'; obj.pleaseExportFigure(plot_savename) - + end - + function plotPosteriorClusterPlot(obj, varargin) %plotPosteriorClusterPlot(H, varargin) Plots posterior distributions % for all experiments. @@ -513,9 +513,9 @@ function plotPosteriorClusterPlot(obj, varargin) % Optional arguments as key/value pairs % 'axisHandle' - handle to axes % 'figureHandle' - handle to figure - + parseFigureAndAxisRequested(varargin{:}) - + % TODO: put clusterPlot function code here rather than have it as a separate function? clusterPlot(... obj.coda,... @@ -524,11 +524,11 @@ function plotPosteriorClusterPlot(obj, varargin) obj.plotOptions,... obj.varList.discountFunctionParams) end - + end - + methods (Access = protected) - + function pleaseExportFigure(obj, savename) if obj.plotOptions.shouldExportPlots myExport(obj.plotOptions.savePath,... @@ -537,17 +537,17 @@ function pleaseExportFigure(obj, savename) 'formats', obj.plotOptions.exportFormats) end end - + function plotAllExperimentFigures(obj) % this is a wrapper function to loop over all data files, producing multi-panel figures. This is implemented by the plotExperimentOverviewFigure method, which may be overridden by subclasses if need be. names = obj.data.getIDnames('all'); - + for experimentIndex = 1:numel(names) fh = figure('Name', names{experimentIndex}); - + obj.plotExperimentOverviewFigure(experimentIndex); drawnow - + if obj.plotOptions.shouldExportPlots myExport(obj.plotOptions.savePath,... 'expt',... @@ -555,12 +555,12 @@ function plotAllExperimentFigures(obj) 'suffix', obj.modelFilename,... 'formats', obj.plotOptions.exportFormats); end - + close(fh); end end - - + + function plotPosteriorErrorParams(obj, ind, varargin) %plotPosteriorErrorParams(H, N) Plots the posterior distributions over % response error related parameters, for experiment N. The plot is @@ -569,69 +569,69 @@ function plotPosteriorErrorParams(obj, ind, varargin) % Optional arguments as key/value pairs % 'axisHandle' - handle to axes % 'figureHandle' - handle to figure - + [figureHandle, axisHandle] = parseFigureAndAxisRequested(varargin{:}); - + responseErrorVariables = obj.getResponseErrorVariables(); - + % TODO: remove duplication of "opts" in mulitple places, but also should perhaps be a single coherent structure in the first place. opts.pointEstimateType = obj.plotOptions.pointEstimateType; opts.timeUnits = obj.timeUnits; opts.dataPlotType = obj.plotOptions.dataPlotType; - + obj.coda.plot_bivariate_distribution(axisHandle,... responseErrorVariables(1),... responseErrorVariables(2),... ind,... opts) end - + end - + methods (Static, Access = protected) - + function plotOptions = parse_plot_options(varargin) p = inputParser; p.StructExpand = false; p.KeepUnmatched = true; p.FunctionName = mfilename; - + p.addParameter('exportFormats', {'png'}, @iscellstr); p.addParameter('savePath',tempname, @isstr); p.addParameter('shouldPlot', 'no', @(x) any(strcmp(x,{'yes','no'}))); p.addParameter('shouldExportPlots', true, @islogical); p.addParameter('pointEstimateType','mode',... @(x) any(strcmp(x,{'mean','median','mode'}))); - + p.parse(varargin{:}); - + plotOptions = p.Results; end - + end - + methods (Hidden = true) - + function disp(obj) - + disp('Discounting model object') linebreak fprintf('Model class: %s\n', class(obj)) obj.dispModelInfo - + linebreak disp(obj.data) - + linebreak disp('MCMC options used in parameter estimation:') disp(obj.mcmcParams) - + linebreak disp('Model parameters (Discounting related):') disp(obj.getDiscountFunctionVariables()) disp('Model parameters (Response error related):') disp(obj.getResponseErrorVariables()) - + linebreak disp('Point estimates of parameters:') exporter = ResultsExporter(obj.coda,... @@ -641,10 +641,10 @@ function disp(obj) obj.varList,... obj.plotOptions); exporter.printToScreen(); - + linebreak obj.displayPublicMethods() end end - + end diff --git a/ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.stan b/ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.stan deleted file mode 100644 index 074cdab7..00000000 --- a/ddToolbox/models/parametric_models/exponential_models/hierarchicalExp1.stan +++ /dev/null @@ -1,74 +0,0 @@ -// RANDOM FACTORS: k[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: k[p], epsilon[p], alpha[p] - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector df_exponential1(vector reward, vector k, vector delay){ - return reward .* exp(-k .* delay); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - real k_mu; - real k_sigma; - vector[nRealExperimentFiles+1] k; // +1 for unobserved participant - - real alpha_mu; - real alpha_sigma; - vector[nRealExperimentFiles+1] alpha; // +1 for unobserved participant - - real omega; - real kappa; - vector[nRealExperimentFiles+1] epsilon; // +1 for unobserved participants -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - VA = df_exponential1(A, k[ID], DA); - VB = df_exponential1(B, k[ID], DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - k_mu ~ normal(0.01, 0.5); // TODO: pick this in a more meaningul manner - k_sigma ~ inv_gamma(0.01,0.01); // TODO: pick this in a more meaningul manner - k ~ normal(k_mu, k_sigma); - - alpha_mu ~ uniform(0,100); - alpha_sigma ~ inv_gamma(0.01,0.01); - alpha ~ normal(alpha_mu, alpha_sigma); - - omega ~ beta(1.1, 10.9); // mode for lapse rate - kappa ~ gamma(0.1,0.1); // concentration parameter - epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); - - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/exponential_models/mixedExp1.stan b/ddToolbox/models/parametric_models/exponential_models/mixedExp1.stan deleted file mode 100644 index 0a88752f..00000000 --- a/ddToolbox/models/parametric_models/exponential_models/mixedExp1.stan +++ /dev/null @@ -1,73 +0,0 @@ -// RANDOM FACTORS: k[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: epsilon[p], alpha[p] - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector df_exponential1(vector reward, vector k, vector delay){ - return reward .* exp(-k .* delay); - } - -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - real alpha_mu; - real alpha_sigma; - vector[nRealExperimentFiles+1] alpha; - - real omega; - real kappa; - vector[nRealExperimentFiles+1] epsilon; - - vector[nRealExperimentFiles] k; // No hierarchical, so no +1 -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - VA = df_exponential1(A, k[ID], DA); - VB = df_exponential1(B, k[ID], DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - alpha_mu ~ uniform(0,100); - alpha_sigma ~ inv_gamma(0.01,0.01); - alpha ~ lognormal(alpha_mu, alpha_sigma); // positive values for alpha - - omega ~ beta(1.1, 10.9); // mode for lapse rate - kappa ~ gamma(0.1,0.1); // concentration parameter - epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); - - // no hierarchical inference for k - k ~ normal(0.01, 0.5^2); - - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/exponential_models/separateExp1.stan b/ddToolbox/models/parametric_models/exponential_models/separateExp1.stan deleted file mode 100644 index 5ebf4d6b..00000000 --- a/ddToolbox/models/parametric_models/exponential_models/separateExp1.stan +++ /dev/null @@ -1,59 +0,0 @@ -// RANDOM FACTORS: k[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: none - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector df_exponential1(vector reward, vector k, vector delay){ - return reward .* exp(-k .* delay); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - vector[nRealExperimentFiles] k; - vector[nRealExperimentFiles] alpha; - vector[nRealExperimentFiles] epsilon; -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - VA = df_exponential1(A, k[ID], DA); - VB = df_exponential1(B, k[ID], DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - // no hierarchical inference for k, alpha, epsilon - k ~ normal(0.01, 0.5^2); - alpha ~ exponential(0.01); - epsilon ~ beta(1.1, 10.9); - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK ? - int Rpostpred[totalTrials]; - - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower b/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower deleted file mode 100755 index 6eec848b..00000000 Binary files a/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower and /dev/null differ diff --git a/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.stan b/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.stan deleted file mode 100644 index dd70a107..00000000 --- a/ddToolbox/models/parametric_models/exponential_power_models/mixedExpPower.stan +++ /dev/null @@ -1,78 +0,0 @@ -// RANDOM FACTORS: k[p], tau[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: epsilon[p], alpha[p] - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector df_exp_power(vector reward, vector delay, vector k, vector tau){ - vector[rows(reward)] V; - for (t in 1:rows(reward)){ - V[t] = reward[t] *( exp( -k[t] * (delay[t]^tau[t]) ) ); - } - return V; - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - // Discounting parameters - vector[nRealExperimentFiles] k; - vector[nRealExperimentFiles] tau; - - // Psychometric function parameters - real alpha_mu; - real alpha_sigma; - vector[nRealExperimentFiles+1] alpha; - - real omega; - real kappa; - vector[nRealExperimentFiles+1] epsilon; -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - VA = df_exp_power(A, DA, k[ID], tau[ID]); - VB = df_exp_power(B, DB, k[ID], tau[ID]); - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - alpha_mu ~ uniform(0,100); - alpha_sigma ~ inv_gamma(0.01,0.01); - alpha ~ normal(alpha_mu, alpha_sigma); - - omega ~ beta(1.1, 10.9); // mode for lapse rate - kappa ~ gamma(0.1,0.1); // concentration parameter - epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); - - // no hierarchical inference for k or tau - k ~ normal(0, 2); - tau ~ exponential(0.01); - - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower b/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower deleted file mode 100755 index 16181ee3..00000000 Binary files a/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower and /dev/null differ diff --git a/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.stan b/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.stan deleted file mode 100644 index bd61413d..00000000 --- a/ddToolbox/models/parametric_models/exponential_power_models/separateExpPower.stan +++ /dev/null @@ -1,60 +0,0 @@ -// RANDOM FACTORS: k[p], tau[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: none - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - real df_exp_power(real reward, real delay, real k, real tau){ - return reward *( exp( -k * (delay^tau) ) ); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - vector[nRealExperimentFiles] k; - vector[nRealExperimentFiles] tau; - vector[nRealExperimentFiles] alpha; - vector[nRealExperimentFiles] epsilon; -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - for (t in 1:totalTrials){ - VA[t] = df_exp_power(A[t], DA[t], k[ID[t]], tau[ID[t]]); - VB[t] = df_exp_power(B[t], DB[t], k[ID[t]], tau[ID[t]]); - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - // no hierarchical inference for k, tau, alpha, epsilon - k ~ normal(0, 0.01); # sigma = 0.1 - tau ~ exponential(0.01); - - alpha ~ exponential(1); - epsilon ~ beta(1+1, 1+1000); - - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK ? - int Rpostpred[totalTrials]; - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME deleted file mode 100755 index e8d30257..00000000 Binary files a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME and /dev/null differ diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.stan b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.stan deleted file mode 100644 index f775d755..00000000 --- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME.stan +++ /dev/null @@ -1 +0,0 @@ -// RANDOM FACTORS: m[p], c[p], epsilon[p], alpha[p] // HYPER-PRIORS ON: none functions { real psychometric_function(real alpha, real epsilon, real VA, real VB){ // returns probability of choosing B (delayed reward) return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); } vector magnitude_effect(vector m, vector c, vector reward){ return m .* log(reward) + c; // we assume reward is positive } vector df_hyperbolic1(vector reward, vector logk, vector delay){ return reward ./ (1+(exp(logk).*delay)); } } data { int totalTrials; int nRealExperimentFiles; vector[totalTrials] A; vector[totalTrials] B; vector[totalTrials] DA; vector[totalTrials] DB; int R[totalTrials]; int ID[totalTrials]; } transformed data { vector[totalTrials] Aabs; vector[totalTrials] Babs; for (t in 1:totalTrials){ Aabs[t] = fabs(A[t]); Babs[t] = fabs(B[t]); } } parameters { real m_mu; real m_sigma; vector[nRealExperimentFiles+1] m; real c_mu; real c_sigma; vector[nRealExperimentFiles+1] c; // real alpha_mu; // real alpha_sigma; real alpha_mode; real alpha_std; vector[nRealExperimentFiles+1] alpha; // real omega; // real kappa; real epsilon_mean; real epsilon_sample_size; vector[nRealExperimentFiles+1] epsilon; } transformed parameters { real epsilon_alpha; real epsilon_beta; real alpha_rate; real alpha_shape; vector[totalTrials] logkA; vector[totalTrials] logkB; vector[totalTrials] VA; vector[totalTrials] VB; vector[totalTrials] P; logkA = magnitude_effect(m[ID], c[ID], Aabs); logkB = magnitude_effect(m[ID], c[ID], Babs); VA = df_hyperbolic1(A, logkA, DA); VB = df_hyperbolic1(B, logkB, DB); for (t in 1:totalTrials){ P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); } // reparameterisation for epsilon epsilon_alpha = epsilon_mean * epsilon_sample_size; epsilon_beta = (1-epsilon_mean) * epsilon_sample_size; // reparameterisation for alpha alpha_rate = (alpha_mode + sqrt(alpha_mode^2 + 4* alpha_std^2))/(2*alpha_std); alpha_shape = 1 + alpha_mode * alpha_rate; } model { m_mu ~ normal(-0.243, 0.27); m_sigma ~ normal( 0.072, 0.25); // convert to inv_gamma(0.01,0.01) ?????????? m ~ normal(m_mu, m_sigma); c_mu ~ normal(0, 10); c_sigma ~ inv_gamma(0.01,0.01); c ~ normal(c_mu, c_sigma); // alpha_mu ~ uniform(0, 100); // alpha_sigma ~ uniform(0, 100); alpha_mode ~ exponential(0.1); alpha_std ~ inv_gamma(0.1, 0.1); alpha ~ gamma(alpha_shape, alpha_rate); // omega ~ beta(1.1, 10.9); // mode for lapse rate // kappa ~ gamma(0.01, 0.01); // concentration parameter //epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); epsilon_mean ~ exponential(1); epsilon_sample_size ~ normal(200,10); epsilon ~ beta(epsilon_alpha , epsilon_beta); R ~ bernoulli(P); } generated quantities { // see page 76 of manual // NO VECTORIZATION IN THIS BLOCK int Rpostpred[totalTrials]; for (t in 1:totalTrials){ Rpostpred[t] = bernoulli_rng(P[t]); } } \ No newline at end of file diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.stan b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.stan deleted file mode 100644 index 496c3851..00000000 --- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalMEupdated.stan +++ /dev/null @@ -1,99 +0,0 @@ -// RANDOM FACTORS: m[p], c[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: p[p], c[p], epsilon[p], alpha[p] - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector magnitude_effect(vector m, vector c, vector reward){ - return m .* log(reward) + c; // we assume reward is positive - } - - vector df_hyperbolic1(vector reward, vector logk, vector delay){ - return reward ./ (1+(exp(logk).*delay)); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -transformed data { - vector[totalTrials] Aabs; - vector[totalTrials] Babs; - for (t in 1:totalTrials){ - Aabs[t] = fabs(A[t]); - Babs[t] = fabs(B[t]); - } -} - -parameters { - real m_mu; - real m_sigma; - vector[nRealExperimentFiles+1] m; - - real c_mu; - real c_sigma; - vector[nRealExperimentFiles+1] c; - - real alpha_mu; - real alpha_sigma; - vector[nRealExperimentFiles+1] alpha; - - real omega; - real kappa; - vector[nRealExperimentFiles+1] epsilon; -} - -transformed parameters { - vector[totalTrials] logkA; - vector[totalTrials] logkB; - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - logkA = magnitude_effect(m[ID], c[ID], Aabs); - logkB = magnitude_effect(m[ID], c[ID], Babs); - VA = df_hyperbolic1(A, logkA, DA); - VB = df_hyperbolic1(B, logkB, DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - m_mu ~ normal(-0.243, 0.27); - m_sigma ~ normal(0,5); - m ~ normal(m_mu, m_sigma); - - c_mu ~ normal(0, 100); - c_sigma ~ normal(0,10); - c ~ normal(c_mu, c_sigma); - - alpha_mu ~ exponential(0.01); - alpha_sigma ~ normal(0,5); - alpha ~ normal(alpha_mu, alpha_sigma); - - omega ~ beta(1+1, 1+100); // mode for lapse rate - kappa ~ gamma(0.1, 0.1); // concentration parameter - epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); - - R ~ bernoulli(P); -} - -generated quantities { // see page 76 of manual // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME deleted file mode 100755 index 88ab3f13..00000000 Binary files a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME and /dev/null differ diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.stan b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.stan deleted file mode 100644 index ff1f57aa..00000000 --- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/mixedME.stan +++ /dev/null @@ -1,89 +0,0 @@ -// RANDOM FACTORS: m[p], c[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: epsilon[p], alpha[p] -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector magnitude_effect(vector m, vector c, vector reward){ - return m .* log(reward) + c; // we assume reward is positive - } - - vector df_hyperbolic1(vector reward, vector logk, vector delay){ - return reward ./ (1+(exp(logk).*delay)); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -transformed data { - vector[totalTrials] Aabs; - vector[totalTrials] Babs; - for (t in 1:totalTrials){ - Aabs[t] = fabs(A[t]); - Babs[t] = fabs(B[t]); - } -} - -parameters { - vector[nRealExperimentFiles] m; // No hierarchical, so no +1 - vector[nRealExperimentFiles] c; // No hierarchical, so no +1 - - real alpha_mu; - real alpha_sigma; - vector[nRealExperimentFiles+1] alpha; - - real omega; - real kappa; - vector[nRealExperimentFiles+1] epsilon; -} - -transformed parameters { - vector[totalTrials] logkA; - vector[totalTrials] logkB; - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - logkA = magnitude_effect(m[ID], c[ID], Aabs); - logkB = magnitude_effect(m[ID], c[ID], Babs); - VA = df_hyperbolic1(A, logkA, DA); - VB = df_hyperbolic1(B, logkB, DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - m ~ normal(-0.243, 0.5); - c ~ normal(0, 10); - - alpha_mu ~ uniform(0, 100); - alpha_sigma ~ uniform(0, 100); - alpha ~ normal(alpha_mu, alpha_sigma); - - omega ~ beta(1.1, 10.9); // mode for lapse rate - kappa ~ gamma(0.01, 0.01); // concentration parameter - epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); - - R ~ bernoulli(P); -} - -generated quantities { // see page 76 of manual // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME deleted file mode 100755 index 09f381db..00000000 Binary files a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME and /dev/null differ diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.stan b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.stan deleted file mode 100644 index ef6a51ac..00000000 --- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/separateME.stan +++ /dev/null @@ -1,79 +0,0 @@ -// RANDOM FACTORS: m[p], c[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: none - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector magnitude_effect(vector m, vector c, vector reward){ - return m .* log(reward) + c; // we assume reward is positive - } - - vector df_hyperbolic1(vector reward, vector logk, vector delay){ - return reward ./ (1+(exp(logk).*delay)); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -transformed data { - vector[totalTrials] Aabs; - vector[totalTrials] Babs; - - for (t in 1:totalTrials){ - Aabs[t] = fabs(A[t]); - Babs[t] = fabs(B[t]); - } -} - -parameters { - vector[nRealExperimentFiles] m; - vector[nRealExperimentFiles] c; - vector[nRealExperimentFiles] alpha; - vector[nRealExperimentFiles] epsilon; -} - -transformed parameters { - vector[totalTrials] logkA; - vector[totalTrials] logkB; - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - logkA = magnitude_effect(m[ID], c[ID], Aabs); - logkB = magnitude_effect(m[ID], c[ID], Babs); - VA = df_hyperbolic1(A, logkA, DA); - VB = df_hyperbolic1(B, logkB, DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - m ~ normal(-0.243, 0.5); - c ~ normal(0, 10); - epsilon ~ beta(1+1, 1+10); - alpha ~ exponential(0.01); - - R ~ bernoulli(P); -} - -generated quantities { // see page 76 of manual // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK b/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK deleted file mode 100755 index 6a2f6c31..00000000 Binary files a/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK and /dev/null differ diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.stan b/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.stan deleted file mode 100644 index f67a5ee4..00000000 --- a/ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK.stan +++ /dev/null @@ -1,74 +0,0 @@ -// RANDOM FACTORS: logk[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: logk[p], epsilon[p], alpha[p] - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector df_hyperbolic1(vector reward, vector logk, vector delay){ - return reward ./ (1+(exp(logk).*delay)); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - real logk_mu; - real logk_sigma; - vector[nRealExperimentFiles+1] logk; // +1 for unobserved participant - - real alpha_mu; - real alpha_sigma; - vector[nRealExperimentFiles+1] alpha; // +1 for unobserved participant - - real omega; - real kappa; - vector[nRealExperimentFiles+1] epsilon; // +1 for unobserved participants -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - VA = df_hyperbolic1(A, logk[ID], DA); - VB = df_hyperbolic1(B, logk[ID], DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - logk_mu ~ normal(-3.9120,2.5); - logk_sigma ~ inv_gamma(0.01,0.01); - logk ~ normal(logk_mu, logk_sigma); - - alpha_mu ~ uniform(0,100); - alpha_sigma ~ inv_gamma(0.01,0.01); - alpha ~ normal(alpha_mu, alpha_sigma); - - omega ~ beta(1.1, 10.9); // mode for lapse rate - kappa ~ gamma(0.1,0.1); // concentration parameter - epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); - - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK b/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK deleted file mode 100755 index 76c04c66..00000000 Binary files a/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK and /dev/null differ diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.stan b/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.stan deleted file mode 100644 index 0e7b6757..00000000 --- a/ddToolbox/models/parametric_models/hyperbolic_models/mixedLogK.stan +++ /dev/null @@ -1,72 +0,0 @@ -// RANDOM FACTORS: logk[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: epsilon[p], alpha[p] - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector df_hyperbolic1(vector reward, vector logk, vector delay){ - return reward ./ (1+(exp(logk).*delay)); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - real alpha_mu; - real alpha_sigma; - vector[nRealExperimentFiles+1] alpha; - - real omega; - real kappa; - vector[nRealExperimentFiles+1] epsilon; - - vector[nRealExperimentFiles] logk; // No hierarchical, so no +1 -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - VA = df_hyperbolic1(A, logk[ID], DA); - VB = df_hyperbolic1(B, logk[ID], DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - alpha_mu ~ uniform(0,100); - alpha_sigma ~ inv_gamma(0.01,0.01); - alpha ~ normal(alpha_mu, alpha_sigma); - - omega ~ beta(1.1, 10.9); // mode for lapse rate - kappa ~ gamma(0.1,0.1); // concentration parameter - epsilon ~ beta(omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ); - - // no hierarchical inference for logk - logk ~ normal(log(1.0/50.0), 2.5); - - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK - int Rpostpred[totalTrials]; - - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK b/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK deleted file mode 100755 index bc4d5174..00000000 Binary files a/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK and /dev/null differ diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK.stan b/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK.stan deleted file mode 100644 index 97b0b3bc..00000000 --- a/ddToolbox/models/parametric_models/hyperbolic_models/separateLogK.stan +++ /dev/null @@ -1,59 +0,0 @@ -// RANDOM FACTORS: logk[p], epsilon[p], alpha[p] -// HYPER-PRIORS ON: none - -functions { - real psychometric_function(real alpha, real epsilon, real VA, real VB){ - // returns probability of choosing B (delayed reward) - return epsilon + (1-2*epsilon) * Phi( (VB-VA) / alpha); - } - - vector df_hyperbolic1(vector reward, vector logk, vector delay){ - return reward ./ (1+(exp(logk).*delay)); - } -} - -data { - int totalTrials; - int nRealExperimentFiles; - vector[totalTrials] A; - vector[totalTrials] B; - vector[totalTrials] DA; - vector[totalTrials] DB; - int R[totalTrials]; - int ID[totalTrials]; -} - -parameters { - vector[nRealExperimentFiles] logk; - vector[nRealExperimentFiles] alpha; - vector[nRealExperimentFiles] epsilon; -} - -transformed parameters { - vector[totalTrials] VA; - vector[totalTrials] VB; - vector[totalTrials] P; - - VA = df_hyperbolic1(A, logk[ID], DA); - VB = df_hyperbolic1(B, logk[ID], DB); - - for (t in 1:totalTrials){ - P[t] = psychometric_function(alpha[ID[t]], epsilon[ID[t]], VA[t], VB[t]); - } -} - -model { - // no hierarchical inference for logk, alpha, epsilon - logk ~ normal(log(1.0/50.0), 2.5); - alpha ~ exponential(0.01); - epsilon ~ beta(1.1, 10.9); - R ~ bernoulli(P); -} - -generated quantities { // NO VECTORIZATION IN THIS BLOCK ? - int Rpostpred[totalTrials]; - - for (t in 1:totalTrials){ - Rpostpred[t] = bernoulli_rng(P[t]); - } -} diff --git a/ddToolbox/sampleWithMatlabStan.m b/ddToolbox/sampleWithMatlabStan.m deleted file mode 100644 index 9c8c2a67..00000000 --- a/ddToolbox/sampleWithMatlabStan.m +++ /dev/null @@ -1,44 +0,0 @@ -function codaObject = sampleWithMatlabStan(... - modelFilename, observedData, mcmcparams, initialParameters, ~) - -assert(ischar(modelFilename)) -assert(isstruct(observedData)) -assert(isstruct(mcmcparams)) - -%% sampler-specific preparation ++++++++++++++++++++ -% NOTE: mstan.stan_home() function defined in the repo 'MatlabStan' which was -% downloaded to your default matlab directory. This is findable by typing: -% >> getenv('HOME') -stan_model = StanModel('file', modelFilename,... - 'stan_home', mstan.stan_home()); -display('COMPILING STAN MODEL...') -tic -stan_model.compile(); -toc -% ++++++++++++++++++++++++++++++++++++++++++++++++++ - -% % convert initialParameters. % TODO: make this general -% init.groupW = [initialParameters(:).groupW]'; -% init.groupALPHAmu = [initialParameters(:).groupALPHAmu]'; -% init.groupALPHAsigma = [initialParameters(:).groupALPHAsigma]'; - -%% Get our sampler to sample -display('SAMPLING STAN MODEL...') -tic -obj.stanFit = stan_model.sampling(... - 'data', observedData,... %'init', init,... - 'warmup', mcmcparams.nburnin,... % warmup = burn-in - 'iter', ceil( mcmcparams.nsamples / mcmcparams.nchains),... % iter = number of MCMC samples - 'chains', mcmcparams.nchains,... - 'verbose', true,... - 'stan_home', mstan.stan_home() ); -% block command window access until sampling finished -obj.stanFit.block(); -toc - -% Uncomment this line if you want auditory feedback -%speak('sampling complete') - -codaObject = CODA.buildFromStanFit(obj.stanFit); - -end \ No newline at end of file diff --git a/ddToolbox/utils/stanExitHandler.m b/ddToolbox/utils/stanExitHandler.m deleted file mode 100644 index 86468e54..00000000 --- a/ddToolbox/utils/stanExitHandler.m +++ /dev/null @@ -1,8 +0,0 @@ -function stanExitHandler(src,data) - fprintf('\n'); - beep; - fprintf('Listener notified!\n'); - fprintf('Stan finished. Chains exited with exitValue = \n'); - disp(src.exit_value) - fprintf('\n'); -end \ No newline at end of file diff --git a/demo/demo_compare_samplers.m b/demo/demo_compare_samplers.m deleted file mode 100644 index e779c411..00000000 --- a/demo/demo_compare_samplers.m +++ /dev/null @@ -1,75 +0,0 @@ -function [models] = demo_compare_samplers(modelType, varargin) -% [models] = demo_compare_samplers('ModelHierarchicalLogK') - -% SEE ISSUE #94 - - -%% Input parsing -defaultMcmcParams = struct('nsamples', 50000,... - 'nchains', 4,... - 'nburnin', 5000); - -p = inputParser; -p.StructExpand = false; -p.FunctionName = mfilename; -% Required -p.addRequired('modelType', @isstr); -p.addParameter('mcmcParams', defaultMcmcParams, @isstruct) -% parse inputs -p.parse(modelType, varargin{:}); -mcmcParams = p.Results.mcmcParams; - - -%% Setup -samplers = {'jags', 'stan'}; -modelFunction = str2func(modelType); - -% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -path_of_this_mfile = strrep(which(mfilename),[mfilename '.m'],''); -toolbox_path = fullfile(path_of_this_mfile,'..','ddToolbox'); -datapath = fullfile(path_of_this_mfile,'datasets','kirby'); -% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -% Run setup routine -addpath(toolbox_path) -ddAnalysisSetUp(); - - -%% Do parameter estimation for each of the sampler types -for sampler = samplers - models.(sampler{:}) = modelFunction(... - Data(datapath, 'files', allFilesInFolder(datapath, 'txt')),... - 'savePath', fullfile(pwd,'output','my_analysis'),... - 'pointEstimateType', 'median',... - 'sampler', sampler{:},... - 'shouldPlot', 'no',... - 'shouldExportPlots', false,... - 'mcmcParams', mcmcParams); -end - - -%% Compare group level parameter estimates (repeated-measures) -assert(numel(samplers)==2,'this is only implemented for comparing TWO samplers against each other') - -% get posteriors -varName = 'logk'; % TODO: Automatically do this for all variables -p=1; % TODO: Automatically do it for all data files -for sampler = samplers - posterior.(sampler{:}) = models.(sampler{:}).coda.getSamplesAtIndex_asStruct(p, {varName}); -end - -% calculate differences -varNames = fieldnames(posterior.(sampler{1})); -for varName = varNames - posteriorDifference.(varName{:}) = minus(... - posterior.(samplers{1}).(varName{:}),... - posterior.(samplers{2}).(varName{:}) ); -end - -% plot distribution of differences -figure -for varName = varNames - mcmc.UnivariateDistribution(posteriorDifference.(varName{:}),... - 'XLabel', [samplers{1} ' - ' samplers{2}],... - 'plotHDI', true) -end diff --git a/demo/experimental_test_nonparametric.m b/demo/experimental_test_nonparametric.m index c8609176..96042d32 100644 --- a/demo/experimental_test_nonparametric.m +++ b/demo/experimental_test_nonparametric.m @@ -25,9 +25,9 @@ 'savePath', 'ModelGaussianRandomWalkSimple',... 'pointEstimateType','median'); -% Do some Bayesian inference with JAGS or STAN -grwModel = grwModel.conductInference('jags',... % {'jags', 'stan'} +% Do some Bayesian inference with JAGS +grwModel = grwModel.conductInference('jags',... 'shouldPlot', 'no',... 'mcmcSamples', 10^3); % TODO: add mcmcparams over-ride -grwModel.plot() \ No newline at end of file +grwModel.plot() diff --git a/demo/run_me.m b/demo/run_me.m index 1435f4f7..002fda0d 100644 --- a/demo/run_me.m +++ b/demo/run_me.m @@ -59,8 +59,8 @@ % If you analysed your data with a model which accounts for the magnitude % effect, then you may want to work out what the discount rate, log(k), % might be for a given reward magnitude. You can do this by: -% >> logk = model.getLogDiscountRate(100,1) % <-------------TODO: MAKE IMPROVEMENTS -% >> logk.plot() % <-------------TODO: MAKE IMPROVEMENTS +% >> logk = model.getLogDiscountRate(100,1) % <-------------TODO: MAKE IMPROVEMENTS +% >> logk.plot() % <-------------TODO: MAKE IMPROVEMENTS % % You can get access to samples using code such as the following. They will % be returned into a structure: @@ -122,8 +122,7 @@ % running proper analyses. I have provided small numbers here just to % confirm the code is working without having to wait a long time. % - you can change the point estimate type to mean, median, or mode -% - the sampler can be 'jags' or 'stan', although stan models are not yet -% full complete +% - the sampler can be 'jags' % If we didn't ask for plots when we ran the model, then we do that diff --git a/tests/test_AllNonParametricModels.m b/tests/test_AllNonParametricModels.m index 65128a8e..77d015ba 100644 --- a/tests/test_AllNonParametricModels.m +++ b/tests/test_AllNonParametricModels.m @@ -4,36 +4,36 @@ % http://uk.mathworks.com/help/matlab/matlab_prog/create-basic-parameterized-test.html % So we don't have to laborously write down test methods for all the % model class types. - + properties data saveName = 'temp.mat' savePath = tempname(); % matlab auto-generated temp folders end - + properties (TestParameter) model = {'ModelSeparateNonParametric'} pointEstimateType = {'mean','median','mode'} - sampler = {'jags', 'stan'} % TODO: ADD STAN + sampler = {'jags'} chains = {2,3,4} end - - + + %% CLASS-LEVEL SETUP/TEARDOWN ----------------------------------------- - + methods (TestClassSetup) function setupData(testCase) % assuming this is running on my maching addpath('~/git-local/delay-discounting-analysis/ddToolbox') datapath = '~/git-local/delay-discounting-analysis/demo/datasets/non_parametric'; - + % only analyse 2 people, for speed of running tests filesToAnalyse = allFilesInFolder(datapath, 'txt'); filesToAnalyse = filesToAnalyse(1:2); testCase.data = Data(datapath, 'files', filesToAnalyse); end end - + methods(TestClassTeardown) function remove_temp_folder(testCase) rmdir(testCase.savePath,'s') @@ -45,17 +45,17 @@ function close_figs(testCase) close all end end - - + + %% THE ACTUAL TESTS --------------------------------------------------- - - + + methods (Test) - + function defaultSampler(testCase, model) % make model makeModelFunction = str2func(model); - + model = makeModelFunction(testCase.data,... 'savePath', testCase.savePath,... 'mcmcParams', struct('nsamples', get_numer_of_samples_for_tests(),... @@ -65,7 +65,7 @@ function defaultSampler(testCase, model) 'shouldExport',false); % TODO: DO AN ACTUAL TEST HERE !!!!!!!!!!!!!!!!!!!!!! end - + function nChains(testCase, model, chains) % make model makeModelFunction = str2func(model); @@ -75,13 +75,13 @@ function nChains(testCase, model, chains) 'nchains', chains,... 'nburnin', get_burnin_for_tests()),... 'shouldPlot','no'); - + % removed this because I've removed get_nChains. It was an % uncessary public method %testCase.verifyEqual(chains, model.get_nChains()) end - - + + function specifiedSampler(testCase, model, sampler) % make model makeModelFunction = str2func(model); @@ -94,11 +94,11 @@ function specifiedSampler(testCase, model, sampler) 'nburnin', get_burnin_for_tests())); % TODO: DO AN ACTUAL TEST HERE !!!!!!!!!!!!!!!!!!!!!! end - + function plotting(testCase, model) % make model makeModelFunction = str2func(model); - + model = makeModelFunction(testCase.data,... 'savePath', testCase.savePath,... 'mcmcParams', struct('nsamples', get_numer_of_samples_for_tests(),... @@ -108,7 +108,7 @@ function plotting(testCase, model) 'shouldExport',false); % TODO: DO AN ACTUAL TEST HERE !!!!!!!!!!!!!!!!!!!!!! end - + function model_disp_function(testCase, model, sampler) % make model makeModelFunction = str2func(model); @@ -122,7 +122,7 @@ function model_disp_function(testCase, model, sampler) % Can we call the disp function without error? disp(modelFitted) end - + end - + end diff --git a/tests/test_AllParametricModels.m b/tests/test_AllParametricModels.m index 1d110ea0..c4889a5e 100644 --- a/tests/test_AllParametricModels.m +++ b/tests/test_AllParametricModels.m @@ -11,7 +11,7 @@ end properties (TestParameter) - sampler = {'jags', 'stan'}; + sampler = {'jags'}; model = getAllParametricModelNames(); pointEstimateType = {'mean','median','mode'} chains = {2,3,4} @@ -24,7 +24,7 @@ function setup(testCase) addpath('~/git-local/delay-discounting-analysis/ddToolbox') datapath = '~/git-local/delay-discounting-analysis/demo/datasets/kirby'; - % only analyse 2 people, for speed of running tests + % only analyse 2 people, for speed of running tests filesToAnalyse = allFilesInFolder(datapath, 'txt'); filesToAnalyse = filesToAnalyse(1:2); testCase.data = Data(datapath, 'files', filesToAnalyse); @@ -97,8 +97,8 @@ function specifiedSampler(testCase, model, sampler) % TODO: DO AN ACTUAL TEST HERE !!!!!!!!!!!!!!!!!!!!!! end - - + + function getting_predicted_values(testCase, model, sampler) % make model makeModelFunction = str2func(model); @@ -109,20 +109,20 @@ function getting_predicted_values(testCase, model, sampler) 'nchains', 2,... 'nburnin', get_burnin_for_tests()),... 'shouldPlot','no'); - + % Get inferred present subjective values of [predicted_subjective_values] =... modelFitted.getInferredPresentSubjectiveValues(); testCase.assertTrue(isstruct(predicted_subjective_values)) - + % TODO: one test per method? % tests for point estimates being a table testCase.assertTrue(istable(predicted_subjective_values.point_estimates)) - + end - - + + function model_disp_function(testCase, model, sampler) % make model makeModelFunction = str2func(model); @@ -136,7 +136,7 @@ function model_disp_function(testCase, model, sampler) % Can we call the disp function without error? disp(modelFitted) end - + end diff --git a/tests/test_PlotAllModels.m b/tests/test_PlotAllModels.m index 540bd719..d641ea4e 100644 --- a/tests/test_PlotAllModels.m +++ b/tests/test_PlotAllModels.m @@ -1,48 +1,43 @@ classdef test_PlotAllModels < matlab.unittest.TestCase - + properties data savePath = tempname(); % matlab auto-generated temp folders end - + properties (TestParameter) model = getAllParametricModelNames; pointEstimateType = {'mean','median','mode'}; dataset = {'kirby', 'non_parametric'}; chains = {2,3} end - + %% CLASS-LEVEL SETUP/TEARDOWN ----------------------------------------- % methods (TestClassSetup) % function setup(testCase) % end % end - - + + methods(TestClassTeardown) - + function remove_temp_directory(testCase) rmdir(testCase.savePath,'s') end - - function delete_stan_related_outputs(testCase) - delete('*.data.R') - delete('*.csv') - end - + function close_figures(testCase) close all end - + end - - + + %% THE ACTUAL TESTS --------------------------------------------------- - + methods (Test) - + function plot(testCase, model, dataset) - + %% Set up data datapath = ['~/git-local/delay-discounting-analysis/demo/datasets/' dataset]; % assuming this is running on my maching @@ -51,7 +46,7 @@ function plot(testCase, model, dataset) filesToAnalyse = allFilesInFolder(datapath, 'txt'); filesToAnalyse = filesToAnalyse(1:2); testCase.data = Data(datapath, 'files', filesToAnalyse); - + %% Run model makeModelFunction = str2func(model); model = makeModelFunction(testCase.data,... @@ -62,11 +57,11 @@ function plot(testCase, model, dataset) 'nchains', 2,... 'nburnin', get_burnin_for_tests() )); %model.plot(); - + close all % TODO: DO AN ACTUAL TEST HERE !!!!!!!!!!!!!!!!!!!!!! end - + end - + end diff --git a/tests/test_public_plot_methods.m b/tests/test_public_plot_methods.m index 34f13fa6..b2eac911 100644 --- a/tests/test_public_plot_methods.m +++ b/tests/test_public_plot_methods.m @@ -1,25 +1,25 @@ classdef test_public_plot_methods < matlab.unittest.TestCase - + properties data model savePath = tempname(); % matlab auto-generated temp folders end - + properties (TestParameter) % % just test with a couple of selected models % %model = {'ModelHierarchicalME', 'ModelHierarchicalLogK', 'ModelSeparateNonParametric'}; % model = {'ModelHierarchicalME'}; % dataset = {'non_parametric'}; end - + %% CLASS-LEVEL SETUP/TEARDOWN ----------------------------------------- methods (TestClassSetup) function setup(testCase) - + model = 'ModelHierarchicalLogK'; dataset = 'non_parametric'; - + %% Set up data datapath = ['~/git-local/delay-discounting-analysis/demo/datasets/' dataset]; % assuming this is running on my maching @@ -28,7 +28,7 @@ function setup(testCase) filesToAnalyse = allFilesInFolder(datapath, 'txt'); filesToAnalyse = filesToAnalyse(1:2); testCase.data = Data(datapath, 'files', filesToAnalyse); - + %% Run model makeModelFunction = str2func(model); testCase.model = makeModelFunction(testCase.data,... @@ -40,61 +40,56 @@ function setup(testCase) 'nburnin', get_burnin_for_tests() )); end end - - + + methods(TestClassTeardown) - + function remove_temp_directory(testCase) rmdir(testCase.savePath,'s') end - - function delete_stan_related_outputs(testCase) - delete('*.data.R') - delete('*.csv') - end - + function close_figures(testCase) close all end - + end - - + + %% THE ACTUAL TESTS --------------------------------------------------- - + methods (Test) - + function plotDiscountFunction(testCase) expt_to_plot = 1; testCase.model.plotDiscountFunction(expt_to_plot); close all end - + function plotDiscountFunctionGrid(testCase) testCase.model.plotDiscountFunctionGrid(); close all end - + function plotDiscountFunctionsOverlaid(testCase) testCase.model.plotDiscountFunctionsOverlaid(); close all end - + function plotPosteriorAUC(testCase) expt_to_plot = 1; testCase.model.plotPosteriorAUC(expt_to_plot); close all end - + function plotPosteriorClusterPlot(testCase) testCase.model.plotPosteriorClusterPlot(); close all end - + function plotUnivarateSummary(testCase) testCase.model.plotUnivarateSummary(); close all end end - + end