diff --git a/ddToolbox/CODA/CODA.m b/ddToolbox/CODA/CODA.m index 6d3af37e..cb938b7f 100644 --- a/ddToolbox/CODA/CODA.m +++ b/ddToolbox/CODA/CODA.m @@ -173,12 +173,6 @@ function densityplot(obj, targetAxisHandle, samples) if ~isempty(targetAxisHandle) subplot(targetAxisHandle) end - - -% % using my plot tools package -% mcmc.UnivariateDistribution(samples',... -% 'plotStyle','hist',... -% 'plotHDI',false); univariateObject = Stochastic('name_here'); univariateObject.addSamples(samples); @@ -208,7 +202,7 @@ function plot_bivariate_distribution(obj, targetAxisHandle, x_var_name, y_var_na 'axisSquare', true); end - function plotUnivariateSummaries(obj, variables, plotOptions, modelFilename, idNames) + function plotUnivariateSummaries(obj, variables, plotOptions, idNames) % create a figure with multiple subplots (boxplots) for each input variable privided % create subplots @@ -238,14 +232,6 @@ function plotUnivariateSummaries(obj, variables, plotOptions, modelFilename, idN % screen_size = get(0,'ScreenSize'); % fig_width = min(screen_size(3), 100+numel(participantNames)*20); % set(gcf,'Position',[100 200 fig_width 1000]) - - % Export - if plotOptions.shouldExportPlots - myExport(plotOptions.savePath,... - 'UnivariateSummary',... - 'suffix', modelFilename,... - 'formats', plotOptions.exportFormats) - end end diff --git a/ddToolbox/CODA/computeStats.m b/ddToolbox/CODA/computeStats.m index c78cd0ff..dc759d7c 100644 --- a/ddToolbox/CODA/computeStats.m +++ b/ddToolbox/CODA/computeStats.m @@ -7,105 +7,121 @@ disp('CODA: Calculating statistics') assert(isstruct(all_samples)) -variable_names = fieldnames(all_samples); +variable_list = asRowVector(fieldnames(all_samples)); stats = struct('Rhat',[], 'mean', [], 'median', [], 'std', [],... 'ci_low' , [] , 'ci_high' , [],... 'hdi_low', [] , 'hdi_high' , []); -for v=1:length(variable_names) - var_name = variable_names{v}; - var_samples = all_samples.(var_name); - - sz = size(var_samples); - Nchains = sz(1); - Nsamples = sz(2); - dims = ndims(var_samples); - - % Calculate stats - switch dims - case{2} % scalar - stats.mode.(var_name) = calcMode( var_samples(:) ); - stats.median.(var_name) = median( var_samples(:) ); - stats.mean.(var_name) = mean( var_samples(:) ); - stats.std.(var_name) = std( var_samples(:) ); - [stats.ci_low.(var_name), stats.ci_high.(var_name)] = calcCI(var_samples(:)); - [stats.hdi_low.(var_name), stats.hdi_high.(var_name)] = calcHDI(var_samples(:)); - - case{3} % vector - for n=1:sz(3) - stats.mode.(var_name)(n) = calcMode( vec(var_samples(:,:,n)) ); - stats.median.(var_name)(n) = median( vec(var_samples(:,:,n)) ); - stats.mean.(var_name)(n) = mean( vec(var_samples(:,:,n)) ); - stats.std.(var_name)(n) = std( vec(var_samples(:,:,n)) ); - [stats.ci_low.(var_name)(n),... - stats.ci_high.(var_name)(n)] = calcCI( vec(var_samples(:,:,n)) ); - [stats.hdi_low.(var_name)(n),... - stats.hdi_high.(var_name)(n)] = calcHDI( vec(var_samples(:,:,n)) ); - end - case{4} % 2D matrix - for a=1:sz(3) - for b=1:sz(4) - stats.mode.(var_name)(a,b) = calcMode( vec(var_samples(:,:,a,b)) ); - stats.median.(var_name)(a,b) = median( vec(var_samples(:,:,a,b)) ); - stats.mean.(var_name)(a,b) = mean( vec(var_samples(:,:,a,b)) ); - stats.std.(var_name)(a,b) = std( vec(var_samples(:,:,a,b)) ); - [stats.ci_low.(var_name)(a,b),... - stats.ci_high.(var_name)(a,b)] = calcCI( vec(var_samples(:,:,a,b)) ); - [stats.hdi_low.(var_name)(a,b),... - stats.hdi_high.(var_name(a,b))] = calcHDI( vec(var_samples(:,:,a,b)) ); - end - end +for variable_name = variable_list + + variable_samples = all_samples.(variable_name{:}); + + switch ndims(variable_samples) + case{2} + stats = calcStatsScalar(stats, variable_samples, variable_name{:}); + case{3} + stats = calcStatsVector(stats, variable_samples, variable_name{:}); + case{4} + stats = calcStatsMatrix(stats, variable_samples, variable_name{:}); otherwise - warning('calculation of stats not supported for >2D variables. You could implement it and send a pull request.') - stats.mode.(var_name) = []; - stats.median.(var_name) = []; - stats.mean.(var_name) = []; - stats.std.(var_name) = []; - stats.ci_low.(var_name) = []; - stats.ci_high.(var_name) = []; - stats.hdi_low.(var_name) = []; - stats.hdi_high.(var_name) = []; + stats = calcStatsTensor3(stats, variable_samples, variable_name{:}); end - - %% "estimated potential scale reduction" statistics due to Gelman and Rubin. - Rhat = calcRhat(); + + Rhat = calcRhat(variable_samples); if ~isnan(Rhat) - stats.Rhat.(var_name) = squeeze(Rhat); + stats.Rhat.(variable_name{:}) = squeeze(Rhat); end + +end end - function Rhat = calcRhat() - st_mean_per_chain = mean(var_samples, 2); - st_mean_overall = mean(st_mean_per_chain, 1); +function stats = calcStatsScalar(stats, var_samples, var_name) +stats.mode.(var_name) = calcMode( var_samples(:) ); +stats.median.(var_name) = median( var_samples(:) ); +stats.mean.(var_name) = mean( var_samples(:) ); +stats.std.(var_name) = std( var_samples(:) ); +[stats.ci_low.(var_name), stats.ci_high.(var_name)] = calcCI(var_samples(:)); +[stats.hdi_low.(var_name), stats.hdi_high.(var_name)] = calcHDI(var_samples(:)); +end - if Nchains > 1 - B = (Nsamples/Nchains-1) * ... - sum((st_mean_per_chain - repmat(st_mean_overall, [Nchains,1])).^2); - varPerChain = var(var_samples, 0, 2); - W = (1/Nchains) * sum(varPerChain); - vhat = ((Nsamples-1)/Nsamples) * W + (1/Nsamples) * B; - Rhat = sqrt(vhat./(W+eps)); - else - Rhat = nan; - end - end +function stats = calcStatsVector(stats, var_samples, var_name) +sz = size(var_samples); +for n=1:sz(3) + stats.mode.(var_name)(n) = calcMode( vec(var_samples(:,:,n)) ); + stats.median.(var_name)(n) = median( vec(var_samples(:,:,n)) ); + stats.mean.(var_name)(n) = mean( vec(var_samples(:,:,n)) ); + stats.std.(var_name)(n) = std( vec(var_samples(:,:,n)) ); + [stats.ci_low.(var_name)(n),... + stats.ci_high.(var_name)(n)] = calcCI( vec(var_samples(:,:,n)) ); + [stats.hdi_low.(var_name)(n),... + stats.hdi_high.(var_name)(n)] = calcHDI( vec(var_samples(:,:,n)) ); +end +end - function [low, high] = calcCI(reshaped_samples) - % get the 95% interval of the posterior - ci_samples_overall = prctile( reshaped_samples , [ 2.5 97.5 ] , 1 ); - ci_samples_overall_low = ci_samples_overall( 1,: ); - ci_samples_overall_high = ci_samples_overall( 2,: ); - low = squeeze(ci_samples_overall_low); - high = squeeze(ci_samples_overall_high); +function stats = calcStatsMatrix(stats, var_samples, var_name) +sz = size(var_samples); +for a=1:sz(3) + for b=1:sz(4) + stats.mode.(var_name)(a,b) = calcMode( vec(var_samples(:,:,a,b)) ); + stats.median.(var_name)(a,b) = median( vec(var_samples(:,:,a,b)) ); + stats.mean.(var_name)(a,b) = mean( vec(var_samples(:,:,a,b)) ); + stats.std.(var_name)(a,b) = std( vec(var_samples(:,:,a,b)) ); + [stats.ci_low.(var_name)(a,b),... + stats.ci_high.(var_name)(a,b)] = calcCI( vec(var_samples(:,:,a,b)) ); + [stats.hdi_low.(var_name)(a,b),... + stats.hdi_high.(var_name(a,b))] = calcHDI( vec(var_samples(:,:,a,b)) ); end +end +end - function [low, high] = calcHDI(reshaped_samples) - % get the 95% highest density intervals of the posterior - [hdi] = HDIofSamples(reshaped_samples, 0.95); - low = squeeze(hdi(1)); - high = squeeze(hdi(2)); - end +function stats = calcStatsTensor3(stats, var_samples, var_name) +warning('calculation of stats not supported for >2D matricies. You could implement it and send a pull request.') +stats.mode.(var_name) = []; +stats.median.(var_name) = []; +stats.mean.(var_name) = []; +stats.std.(var_name) = []; +stats.ci_low.(var_name) = []; +stats.ci_high.(var_name) = []; +stats.hdi_low.(var_name) = []; +stats.hdi_high.(var_name) = []; +end + +function Rhat = calcRhat(var_samples) +% "estimated potential scale reduction" statistics due to Gelman and Rubin +sz = size(var_samples); +Nchains = sz(1); +Nsamples = sz(2); + +st_mean_per_chain = mean(var_samples, 2); +st_mean_overall = mean(st_mean_per_chain, 1); + +if Nchains > 1 + B = (Nsamples/Nchains-1) * ... + sum((st_mean_per_chain - repmat(st_mean_overall, [Nchains,1])).^2); + varPerChain = var(var_samples, 0, 2); + W = (1/Nchains) * sum(varPerChain); + vhat = ((Nsamples-1)/Nsamples) * W + (1/Nsamples) * B; + Rhat = sqrt(vhat./(W+eps)); +else + Rhat = nan; +end +end + +function [low, high] = calcCI(reshaped_samples) +% get the 95% interval of the posterior +ci_samples_overall = prctile( reshaped_samples , [ 2.5 97.5 ] , 1 ); +ci_samples_overall_low = ci_samples_overall( 1,: ); +ci_samples_overall_high = ci_samples_overall( 2,: ); +low = squeeze(ci_samples_overall_low); +high = squeeze(ci_samples_overall_high); +end + +function [low, high] = calcHDI(reshaped_samples) +% get the 95% highest density intervals of the posterior +[hdi] = HDIofSamples(reshaped_samples, 0.95); +low = squeeze(hdi(1)); +high = squeeze(hdi(2)); end diff --git a/ddToolbox/Data/Data.m b/ddToolbox/Data/Data.m index b8c8729f..3cc80c46 100644 --- a/ddToolbox/Data/Data.m +++ b/ddToolbox/Data/Data.m @@ -14,7 +14,7 @@ properties (Dependent) totalTrials - groupTable % table of A§, DA, B, DB, R, ID, PA, PB + groupTable % table of A, DA, B, DB, R, ID, PA, PB end % NOTE TO SELF: These public methods need to be seen as interfaces to @@ -234,9 +234,13 @@ function exportGroupDataFileToDisk(obj) end end - function output = getEverythingAboutAnExperiment(obj, ind) - % return a structure of everything about the data file 'ind' - + function [samples] = getGroupLevelSamples(obj, fieldsToGet) + if ~obj.isUnobservedPartipantPresent() + error('Looks like we don''t have group level estimates.') + else + index = obj.data.getIndexOfUnobservedParticipant(); + samples = obj.coda.getSamplesAtIndex_asStruct(index, fieldsToGet); + end end function participantIndexList = getParticipantIndexList(obj) diff --git a/ddToolbox/Data/DataImporter.m b/ddToolbox/Data/DataImporter.m index 098348aa..fec7e50e 100644 --- a/ddToolbox/Data/DataImporter.m +++ b/ddToolbox/Data/DataImporter.m @@ -18,13 +18,13 @@ obj.path = path; obj.fnames = fnames; obj.nFiles = numel(fnames); - obj = obj.import(); - + % do importing + obj.dataArray = obj.import(); disp('The following data files were imported:') disp(fnames') end - function obj = import(obj) + function dataArray = import(obj) for n=1:obj.nFiles %% do the actual file import experimentTable = readtable(... @@ -47,7 +47,6 @@ dataArray(n) = DataFile(experimentTable); end - obj.dataArray = dataArray; end function dataArray = getData(obj) diff --git a/ddToolbox/DeterministicFunction/DiscountFunction.m b/ddToolbox/DeterministicFunction/DiscountFunction.m index 91a3bfe5..2f3f24f0 100644 --- a/ddToolbox/DeterministicFunction/DiscountFunction.m +++ b/ddToolbox/DeterministicFunction/DiscountFunction.m @@ -16,7 +16,7 @@ function plot(obj, pointEstimateType, dataPlotType, timeUnits) timeUnitFunction = str2func(timeUnits); N_SAMPLES_FROM_POSTERIOR = 100; - delays = obj.determineDelayValues(); + delays = obj.getDelayValues(); if verLessThan('matlab','9.1') % backward compatability delaysDuration = delays; else @@ -59,39 +59,36 @@ function plot(obj, pointEstimateType, dataPlotType, timeUnits) drawnow end - - function delayValues = determineDelayValues(obj) - % TODO: remove this stupid special-case handling of group-level - % participant with no data - try - maxDelayRange = max( obj.data.getDelayRange() )*1.2; - catch - % default (happens when there is no data, ie group level - % observer). - maxDelayRange = 365; - end - delayValues = linspace(0, maxDelayRange, 1000); - end - - function nansPresent = anyNaNsPresent(obj) - nansPresent = false; - for field = fields(obj.theta)' - if any(isnan(obj.theta.(field{:}).samples)) - nansPresent = true; - warning('NaN''s detected in theta') - break - end - end - end - + end - methods (Abstract) - - - - end - - - + + methods (Access = private) + + function delayValues = getDelayValues(obj) + % TODO: remove this stupid special-case handling of group-level + % participant with no data + try + maxDelayRange = max( obj.data.getDelayRange() )*1.2; + catch + % default (happens when there is no data, ie group level + % observer). + maxDelayRange = 365; + end + delayValues = linspace(0, maxDelayRange, 1000); + end + + function nansPresent = anyNaNsPresent(obj) + nansPresent = false; + for field = fields(obj.theta)' + if any(isnan(obj.theta.(field{:}).samples)) + nansPresent = true; + warning('NaN''s detected in theta') + break + end + end + end + + end + end diff --git a/ddToolbox/DeterministicFunction/PsychometricFunction.m b/ddToolbox/DeterministicFunction/PsychometricFunction.m index 8233de1d..dd60590e 100644 --- a/ddToolbox/DeterministicFunction/PsychometricFunction.m +++ b/ddToolbox/DeterministicFunction/PsychometricFunction.m @@ -6,6 +6,7 @@ function obj = PsychometricFunction(varargin) obj = obj@DeterministicFunction(varargin{:}); + % TODO: this violates dependency injection, so we may want to pass these Stochastic objects in obj.theta.alpha = Stochastic('alpha'); obj.theta.epsilon = Stochastic('epsilon'); diff --git a/ddToolbox/PosteriorPrediction.m b/ddToolbox/PosteriorPrediction.m index 511bbd80..453ee013 100644 --- a/ddToolbox/PosteriorPrediction.m +++ b/ddToolbox/PosteriorPrediction.m @@ -2,7 +2,7 @@ %PosteriorPrediction Summary of this class goes here % Detailed explanation goes here - properties + properties (SetAccess = private, GetAccess = private) postPred end @@ -35,19 +35,32 @@ end function plot(obj, plotOptions, modelFilename) - % loop over all experiments/people, prodicing plot figures + % loop over all experiments/people, producing plot figures N = length(obj.postPred); for n = 1:N obj.posterior_prediction_figure(n, plotOptions, modelFilename) + + % Exporting + prefix_string = obj.postPred(n).IDname{:}; + obj.exportFigure(plotOptions, prefix_string, modelFilename) end - end + end + + % PUBLIC GETTERS -------------------------------------------------- + + function score = getScores(obj) + score = [obj.postPred(:).score]'; + end + + function pp = getPercentPredictedDistribution(obj) + pp = {obj.postPred(:).percentPredictedDistribution}; + end end methods (Access = private) function posterior_prediction_figure(obj, n, plotOptions, modelFilename) - pp = obj.postPred(n); % Sort figure figure(1), colormap(gray), clf latex_fig(16, 9, 6) @@ -56,27 +69,20 @@ function posterior_prediction_figure(obj, n, plotOptions, modelFilename) subplot(h(1)), obj.pp_plotTrials(n) subplot(h(2)), obj.pp_plotGOFdistribution(n, plotOptions) subplot(h(3)), obj.pp_plotPercentPredictedDistribution(n, plotOptions) - % Export figure drawnow - if plotOptions.shouldExportPlots - myExport(plotOptions.savePath, 'PosteriorPredictive',... - 'prefix', pp.IDname{:},... - 'suffix', modelFilename,... - 'formats', plotOptions.exportFormats) - end end function pp_plotGOFdistribution(obj, n, plotOptions) - pp = obj.postPred(n); - uni = mcmc.UnivariateDistribution(pp.GOF_distribtion(:),... + uni = mcmc.UnivariateDistribution(... + obj.postPred(n).GOF_distribtion(:),... 'xLabel', 'goodness of fit score',... 'plotStyle','hist',... 'pointEstimateType', plotOptions.pointEstimateType); end function pp_plotPercentPredictedDistribution(obj, n, plotOptions) - pp = obj.postPred(n); - uni = mcmc.UnivariateDistribution(pp.percentPredictedDistribution(:),... + uni = mcmc.UnivariateDistribution(... + obj.postPred(n).percentPredictedDistribution(:),... 'xLabel', '$\%$ proportion responses accounted for',... 'plotStyle','hist',... 'pointEstimateType', plotOptions.pointEstimateType); @@ -86,22 +92,20 @@ function pp_plotPercentPredictedDistribution(obj, n, plotOptions) end function pp_plotTrials(obj, n) - pp = obj.postPred(n); % plot predicted probability of choosing delayed - bar(pp.responses_predicted, 'BarWidth',1) + bar(obj.postPred(n).responses_predicted, 'BarWidth',1) % formatting box off axis tight % plot response data hold on - plot([1:numel(pp.responses_actual)], pp.responses_actual, '+') + plot([1:numel(obj.postPred(n).responses_actual)], obj.postPred(n).responses_actual, '+') % formatting xlabel('trial') ylabel('response') legend('prediction','response', 'Location','East') end - end methods (Access = private, Static) @@ -118,7 +122,7 @@ function pp_plotTrials(obj, n) end function percentResponsesPredicted = calcPercentResponsesCorrectlyPredicted(responses_predictedMCMC, responses_actual) - %% Calculate % responses predicted by the model + % Calculate % responses predicted by the model totalSamples = size(responses_predictedMCMC,2); nQuestions = numel(responses_actual); modelPrediction = zeros(size(responses_predictedMCMC)); @@ -132,11 +136,21 @@ function pp_plotTrials(obj, n) % Calculate log posterior odds of data under the model and a % control model where prob of responding is 0.5. responses_control_model = ones(size(responses_predicted)).*0.5; - score = calcLogOdds(... calcDataLikelihood(responses_actual, responses_predicted'),... calcDataLikelihood(responses_actual, responses_control_model')); end + + function exportFigure(plotOptions, prefix_string, modelFilename) + % TODO: Exporting is not the responsibility of PosteriorPrediction class, so we need to extract this up to Model subclasses. They call it as: obj.postPred.plot(obj.plotOptions, obj.modelFilename) + if plotOptions.shouldExportPlots + myExport(plotOptions.savePath,... + 'PosteriorPredictive',... + 'prefix', prefix_string,... + 'suffix', modelFilename,... + 'formats', plotOptions.exportFormats) + end + end end diff --git a/ddToolbox/ResultsExporter.m b/ddToolbox/ResultsExporter.m index b970a265..79e0e6d2 100644 --- a/ddToolbox/ResultsExporter.m +++ b/ddToolbox/ResultsExporter.m @@ -85,7 +85,7 @@ function export(obj, savePath, pointEstimateType) end function postPredTable = makePostPredTable(obj) - postPredTable = table([obj.postPred(:).score]',... + postPredTable = table(obj.postPred.getScores(),... obj.calc_percent_predicted_point_estimate(),... obj.any_percent_predicted_warnings(),... 'RowNames', obj.data.getIDnames('experiments'),... @@ -106,7 +106,7 @@ function export(obj, savePath, pointEstimateType) % estimate type that the user specified pointEstFunc = str2func(obj.plotOptions.pointEstimateType); percentPredicted = cellfun(pointEstFunc,... - {obj.postPred.percentPredictedDistribution})'; + obj.postPred.getPercentPredictedDistribution())'; end function pp_warning = any_percent_predicted_warnings(obj) @@ -115,7 +115,7 @@ function export(obj, savePath, pointEstimateType) warningFunc = @(x) x(1) < ppLowerThreshold; warnOnHDI = @(x) warningFunc( hdiFunc(x) ); pp_warning = cellfun( warnOnHDI,... - {obj.postPred.percentPredictedDistribution})'; + obj.postPred.getPercentPredictedDistribution())'; end end end diff --git a/ddToolbox/models/Model.m b/ddToolbox/models/Model.m index 10416a8a..2ac7d3d4 100644 --- a/ddToolbox/models/Model.m +++ b/ddToolbox/models/Model.m @@ -12,16 +12,12 @@ properties (SetAccess = protected, GetAccess = protected) dfClass % function handle to DiscountFunction class samplerType - discountFuncType postPred - parameterEstimateTable - experimentFigPlotFuncs mcmcParams % structure of user-supplied params observedData % User supplied preferences modelFilename % string (ie modelFilename.jags, or modelFilename.stan) varList - plotFuncs % structure of function handles plotOptions timeUnits % string whose name must be a function to create a Duration. end @@ -31,11 +27,6 @@ plot() initialiseChainValues() experimentMultiPanelFigure() - plot_discount_functions_in_grid() - - % TODO: sort these out, #172... probably want to rejig so it's done differently - %getExperimentData - %extractDiscountFunctionSamples end methods (Access = public) @@ -128,7 +119,7 @@ function export(obj) obj.data.getIDnames('all')) exporter = ResultsExporter(obj.coda, obj.data,... - obj.postPred.postPred,... + obj.postPred,... obj.varList,... obj.plotOptions); exporter.printToScreen(); @@ -141,6 +132,10 @@ function export(obj) function obj = plotMCMCchains(obj,vars) obj.coda.plotMCMCchains(vars); end + + function [samples] = getGroupLevelSamples(obj, fieldsToGet) + [samples] = obj.data.getGroupLevelSamples(fieldsToGet); + end end @@ -152,37 +147,24 @@ function export(obj) nChains = obj.mcmcParams.nchains; end - function [samples] = getGroupLevelSamples(obj, fieldsToGet) - if ~obj.data.isUnobservedPartipantPresent() - % exit if we don't have any group level inference - error('Looks like we don''t have group level estimates.') - else - index = obj.data.getIndexOfUnobservedParticipant(); - samples = obj.coda.getSamplesAtIndex_asStruct(index, fieldsToGet); - end - end + function [predicted_subjective_values] = get_inferred_present_subjective_values(obj) %% calculate point estimates % get point estimates of present subjective values. These will % be vectors. Each value corresponds to one trial in the % overall dataset - VA_point_estimate = obj.coda.getStats(obj.plotOptions.pointEstimateType, 'VA'); - VB_point_estimate = obj.coda.getStats(obj.plotOptions.pointEstimateType, 'VB'); - assert(isvector(VA_point_estimate)) - assert(isvector(VB_point_estimate)) - + + %% return point estimates of present subjective values... all_data_table = obj.data.groupTable; - all_data_table.VA = VA_point_estimate; - all_data_table.VB = VB_point_estimate; + % 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; - %% Return full posterior distributions of present subjective values - % TODO + %% TODO: Return full posterior distributions of present subjective values % predicted_subjective_values.A_full_posterior = % predicted_subjective_values.B_full_posterior = - - %% return point estimates of present subjectiv values... - predicted_subjective_values.point_estimates = all_data_table; end end @@ -211,7 +193,7 @@ function export(obj) observedData.nRealExperimentFiles = obj.data.getNRealExperimentFiles(); observedData.totalTrials = height(all_data); % protected method which can be over-ridden by model sub-classes - observedData = obj.addititional_model_specific_ObservedData(observedData); + observedData = obj.additional_model_specific_ObservedData(observedData); end function obj = calcDerivedMeasures(obj) @@ -238,7 +220,8 @@ function plotAllExperimentFigures(obj) drawnow if obj.plotOptions.shouldExportPlots - myExport(obj.plotOptions.savePath, 'expt',... + myExport(obj.plotOptions.savePath,... + 'expt',... 'prefix', names{experimentIndex},... 'suffix', obj.modelFilename,... 'formats', obj.plotOptions.exportFormats); @@ -247,13 +230,38 @@ function plotAllExperimentFigures(obj) close(fh); end end + + function plot_discount_functions_in_grid(obj) + latex_fig(12, 11,11) + clf, 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) + obj.plot_discount_function(subplot_handles(n), n) + title(names{n}, 'FontSize',10) + set(gca,'FontSize',10) + end + drawnow + end end methods (Static, Access = protected) - function observedData = addititional_model_specific_ObservedData(observedData) + 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? diff --git a/ddToolbox/models/nonparametric_models/DF_NonParametric.m b/ddToolbox/models/nonparametric_models/DF_NonParametric.m index 73ee35ec..d6a791d0 100644 --- a/ddToolbox/models/nonparametric_models/DF_NonParametric.m +++ b/ddToolbox/models/nonparametric_models/DF_NonParametric.m @@ -33,7 +33,9 @@ % We will estimate AUC AUC_samples = obj.calcAUC(); + % create AUC as a Stochastic object + % TODO: this violates dependency injection, so we may want to pass these Stochastic objects in obj.AUC = Stochastic('AUC'); obj.AUC.addSamples(AUC_samples); end diff --git a/ddToolbox/models/nonparametric_models/DF_SLICE_PsychometricFunction.m b/ddToolbox/models/nonparametric_models/DF_SLICE_PsychometricFunction.m index 0fdbacb1..a51e879d 100644 --- a/ddToolbox/models/nonparametric_models/DF_SLICE_PsychometricFunction.m +++ b/ddToolbox/models/nonparametric_models/DF_SLICE_PsychometricFunction.m @@ -9,6 +9,7 @@ obj = obj@DeterministicFunction(); % create Stochastic objects + % TODO: this violates dependency injection, so we may want to pass these Stochastic objects in obj.theta.indifference = Stochastic('indifference'); obj.theta.alpha = Stochastic('alpha'); obj.theta.epsilon = Stochastic('epsilon'); diff --git a/ddToolbox/models/nonparametric_models/NonParametric.m b/ddToolbox/models/nonparametric_models/NonParametric.m index e1239f96..881439bf 100644 --- a/ddToolbox/models/nonparametric_models/NonParametric.m +++ b/ddToolbox/models/nonparametric_models/NonParametric.m @@ -3,7 +3,6 @@ properties (Access = private) AUC_DATA - getDiscountRate % function handle end methods (Access = public) @@ -14,21 +13,11 @@ % Create variables obj.varList.participantLevel = {'Rstar'}; obj.varList.monitored = {'Rstar', 'alpha', 'epsilon', 'Rpostpred', 'P'}; + + obj.varList.discountFunctionParams(1).name = 'Rstar'; + obj.varList.discountFunctionParams(1).label = 'Rstar'; end - end - - - methods (Access = protected) - - function obj = calcDerivedMeasures(obj) - end - - end - - - - methods (Access = public) function plot(obj, varargin) % overriding from Model base class close all @@ -40,6 +29,13 @@ function plot(obj, varargin) % overriding from Model base class p.parse(varargin{:}); obj.plot_discount_functions_in_grid(); + % Export + if obj.plotOptions.shouldExportPlots + myExport(obj.plotOptions.savePath,... + 'grid_discount_functions',... + 'suffix', obj.modelFilename,... + 'formats', obj.plotOptions.exportFormats); + end % EXPERIMENT PLOT ================================================== obj.psychometric_plots(); @@ -64,24 +60,24 @@ function experimentMultiPanelFigure(obj, ind) latex_fig(12, 14, 3) h = layout([1 2 3 3]); - %% Set up psychometric function - psycho = PsychometricFunction('samples', obj.coda.getSamplesAtIndex_asStruct(ind,{'alpha','epsilon'})); + opts.pointEstimateType = obj.plotOptions.pointEstimateType; + opts.timeUnits = obj.timeUnits; - %% plot bivariate distribution of alpha, epsilon - subplot(h(1)) - samples = obj.coda.getSamplesAtIndex_asStruct(ind,{'alpha','epsilon'}); - mcmc.BivariateDistribution(... - samples.epsilon(:),... - samples.alpha(:),... - 'xLabel','error rate, $\epsilon$',... - 'ylabel','comparison accuity, $\alpha$',... - 'pointEstimateType',obj.plotOptions.pointEstimateType,... - 'plotStyle', 'hist',... - 'axisSquare', true); - + responseErrorVariables = {obj.varList.responseErrorParams.name}; + + %% Set up psychometric function + psycho = PsychometricFunction('samples', obj.coda.getSamplesAtIndex_asStruct(ind,responseErrorVariables)); + + %% PLOT: density plot of (alpha, epsilon) + obj.coda.plot_bivariate_distribution(h(1),... + responseErrorVariables(1),... + responseErrorVariables(2),... + ind,... + opts) + %% Set up discount function personInfo = obj.getExperimentData(ind); - discountFunction = DF_NonParametric('delays',personInfo.delays,... + discountFunction = obj.dfClass('delays',personInfo.delays,... 'theta', personInfo.dfSamples); % inject a DataFile object into the discount function discountFunction.data = obj.data.getExperimentObject(ind); @@ -92,60 +88,27 @@ function experimentMultiPanelFigure(obj, ind) xlim([0 2]) %% plot discount function - subplot(h(3)) - discountFunction.plot(); + obj.plot_discount_function(h(3), ind) end - - - function plot_discount_functions_in_grid(obj) - latex_fig(12, 11,11) - - % 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'); - - clf, drawnow - - % create grid layout - N = numel(names); - subplot_handles = create_subplots(N, 'square'); - - % Iterate over files, plotting - disp('Plotting...') - - for n = 1:numel(names) - subplot(subplot_handles(n)) - % ~~~~~~~~~~~~~~~~~~ - plot_df(n) - % ~~~~~~~~~~~~~~~~~~ - end - drawnow - - if obj.plotOptions.shouldExportPlots - myExport(obj.plotOptions.savePath, 'grid_discount_functions',... - 'suffix', obj.modelFilename,... - 'formats', obj.plotOptions.exportFormats); - end - - function plot_df(ind) - % Set up discount function - personInfo = obj.getExperimentData(ind); - discountFunction = DF_NonParametric('delays',personInfo.delays,... - 'theta', personInfo.dfSamples); - discountFunction.data = obj.data.getExperimentObject(ind); - - discountFunction.plot(); - end - - end - - - - + + + % TODO: work to be able to move this method up to Model base class from both Parametric and NonParamtric + function plot_discount_function(obj, subplot_handle, ind) + discountFunctionVariables = {obj.varList.discountFunctionParams.name}; + subplot(subplot_handle) + + personInfo = obj.getExperimentData(ind); % TODO: do we really need this? + + + % TODO: DF_NonParametric should have same interface as DF_Hyperbolic1 etc + discountFunction = obj.dfClass('delays',personInfo.delays,... + 'theta', personInfo.dfSamples); + discountFunction.data = obj.data.getExperimentObject(ind); + + discountFunction.plot(); + % TODO #166 avoid having to parse these args in here + end function personStruct = getExperimentData(obj, p) @@ -214,8 +177,13 @@ function plot_df(ind) end + + methods (Access = protected) + function obj = calcDerivedMeasures(obj) + end + function psychometric_plots(obj) % TODO: plot data on these figures @@ -240,7 +208,6 @@ function psychometric_plots(obj) psycho.plot(); %% plot response data TODO: move this to Data ~~~~~~~~~ hold on - %pTable = obj.data.getRawDataTableForParticipant(ind); AoverB = personStruct.data.A ./ personStruct.data.B; R = personStruct.data.R; % grab just for this delay @@ -253,7 +220,8 @@ function psychometric_plots(obj) end drawnow if obj.plotOptions.shouldExportPlots - myExport(obj.plotOptions.savePath, 'expt_psychometric',... + myExport(obj.plotOptions.savePath,... + 'expt_psychometric',... 'prefix', names{ind},... 'suffix', obj.modelFilename,... 'formats', obj.plotOptions.exportFormats ); @@ -266,7 +234,7 @@ function psychometric_plots(obj) methods (Static, Access = protected) - function observedData = addititional_model_specific_ObservedData(observedData) + function observedData = additional_model_specific_ObservedData(observedData) observedData.uniqueDelays = sort(unique(observedData.DB))'; observedData.delayLookUp = calcDelayLookup(); diff --git a/ddToolbox/models/parametric_models/Parametric.m b/ddToolbox/models/parametric_models/Parametric.m index 06ba009d..7561e915 100644 --- a/ddToolbox/models/parametric_models/Parametric.m +++ b/ddToolbox/models/parametric_models/Parametric.m @@ -1,9 +1,5 @@ classdef (Abstract) Parametric < Model - - properties (Access = private) - end - methods (Access = public) function obj = Parametric(data, varargin) @@ -21,18 +17,42 @@ function plot(obj, varargin) variables = obj.varList.participantLevel; obj.coda.plotUnivariateSummaries(variables,... obj.plotOptions,... - obj.modelFilename,... obj.data.getParticipantNames()); + + % Export + if obj.plotOptions.shouldExportPlots + myExport(obj.plotOptions.savePath,... + 'UnivariateSummary',... + 'suffix', obj.modelFilename,... + 'formats', obj.plotOptions.exportFormats) + end + % summary figure of core discounting parameters clusterPlot(... obj.coda,... obj.data,... [1 0 0],... - obj.modelFilename,... obj.plotOptions,... - obj.varList.discountFunctionParams) - + obj.varList.discountFunctionParams) + + % Export + if obj.plotOptions.shouldExportPlots + myExport(obj.plotOptions.savePath,... + 'summary_plot',... + 'suffix', obj.modelFilename,... + 'formats', obj.plotOptions.exportFormats) + end + + obj.plot_discount_functions_in_grid(); + % Export + if obj.plotOptions.shouldExportPlots + myExport(obj.plotOptions.savePath,... + 'grid_discount_functions',... + 'suffix', obj.modelFilename,... + 'formats', obj.plotOptions.exportFormats); + end + %% Plots, one per data file =================================== obj.plotAllExperimentFigures(); @@ -56,70 +76,85 @@ function experimentMultiPanelFigure(obj, ind) discountFunctionVariables = {obj.varList.discountFunctionParams.name}; responseErrorVariables = {obj.varList.responseErrorParams.name}; - %% PLOT: density plot of (alpha, epsilon) - obj.coda.plot_bivariate_distribution(h(1),... - responseErrorVariables(1),... - responseErrorVariables(2),... - ind,... - opts) - - %% Plot the psychometric function ---------------------------------- - subplot(h(2)) - psycho = PsychometricFunction('samples', obj.coda.getSamplesAtIndex_asStruct(ind, responseErrorVariables)); - psycho.plot(obj.plotOptions.pointEstimateType) - - %% Plot the discount function parameters --------------------------- - switch numel(discountFunctionVariables) - case{1} - obj.coda.plot_univariate_distribution(h(3),... - discountFunctionVariables(1),... - ind,... - opts) - case{2} - obj.coda.plot_bivariate_distribution(h(3),... - discountFunctionVariables(1),... - discountFunctionVariables(2),... - ind,... - opts) - otherwise - error('Currently only set up to plot univariate or bivariate distributions, ie discount functions 1 or 2 params.') - end - - %% Plot the discount function parameters ---------------------- - subplot(h(4)) - discountFunction = obj.dfClass(... - 'samples', obj.coda.getSamplesAtIndex_asStruct(ind, discountFunctionVariables),... - 'data', obj.data.getExperimentObject(ind)); - discountFunction.plot(obj.plotOptions.pointEstimateType,... - obj.plotOptions.dataPlotType,... - obj.timeUnits); - % TODO #166 avoid having to parse these args in here + plot_density_alpha_epsilon(h(1)) + plot_psychometric_function(h(2)) + plot_discount_function_parameters(h(3)) + obj.plot_discount_function(h(4), ind) + + function plot_density_alpha_epsilon(subplot_handle) + obj.coda.plot_bivariate_distribution(subplot_handle,... + responseErrorVariables(1),... + responseErrorVariables(2),... + ind,... + opts) + end + + function plot_psychometric_function(subplot_handle) + subplot(subplot_handle) + psycho = PsychometricFunction('samples', obj.coda.getSamplesAtIndex_asStruct(ind, responseErrorVariables)); + psycho.plot(obj.plotOptions.pointEstimateType) + end + + function plot_discount_function_parameters(subplot_handle) + switch numel(discountFunctionVariables) + case{1} + obj.coda.plot_univariate_distribution(subplot_handle,... + discountFunctionVariables(1),... + ind,... + opts) + case{2} + obj.coda.plot_bivariate_distribution(subplot_handle,... + discountFunctionVariables(1),... + discountFunctionVariables(2),... + ind,... + opts) + otherwise + error('Currently only set up to plot univariate or bivariate distributions, ie discount functions 1 or 2 params.') + end + end end - - function plot_discount_functions_in_grid(obj) - error('Implement this! See #170') + + + % TODO: work to be able to move this method up to Model base class from both Parametric and NonParamtric + function plot_discount_function(obj, subplot_handle, ind) + discountFunctionVariables = {obj.varList.discountFunctionParams.name}; + + subplot(subplot_handle) + discountFunction = obj.dfClass(... + 'samples', obj.coda.getSamplesAtIndex_asStruct(ind, discountFunctionVariables),... + 'data', obj.data.getExperimentObject(ind)); + discountFunction.plot(obj.plotOptions.pointEstimateType,... + obj.plotOptions.dataPlotType,... + obj.timeUnits); + % TODO #166 avoid having to parse these args in here end end methods (Access = private) + % TODO: this should be moved to CODA function plotAllTriPlots(obj, plotOptions, modelFilename) + + pVariableNames = obj.varList.participantLevel; + for p = 1:obj.data.getNExperimentFiles figure(87), clf - pVariableNames = obj.varList.participantLevel; posteriorSamples = obj.coda.getSamplesAtIndex_asMatrix(p, pVariableNames); mcmc.TriPlotSamples(posteriorSamples,... pVariableNames,... 'pointEstimateType', plotOptions.pointEstimateType); + id_prefix = obj.data.getIDnames(p); + + % Export if plotOptions.shouldExportPlots - id = obj.data.getIDnames(p); - myExport(plotOptions.savePath, 'triplot',... - 'prefix', id{:},... + myExport(plotOptions.savePath,... + 'triplot',... + 'prefix', id_prefix{:},... 'suffix', modelFilename,... 'formats', plotOptions.exportFormats); end diff --git a/ddToolbox/models/parametric_models/bens_new_model/ExponentialPower.m b/ddToolbox/models/parametric_models/bens_new_model/ExponentialPower.m index 50a97c46..29f1e6f5 100644 --- a/ddToolbox/models/parametric_models/bens_new_model/ExponentialPower.m +++ b/ddToolbox/models/parametric_models/bens_new_model/ExponentialPower.m @@ -1,9 +1,5 @@ classdef (Abstract) ExponentialPower < Parametric - properties (Access = private) - getDiscountRate % function handle - end - methods (Access = public) function obj = ExponentialPower(data, varargin) diff --git a/ddToolbox/models/parametric_models/exponential_models/DF_Exponential1.m b/ddToolbox/models/parametric_models/exponential_models/DF_Exponential1.m index 0544493a..5560d88a 100644 --- a/ddToolbox/models/parametric_models/exponential_models/DF_Exponential1.m +++ b/ddToolbox/models/parametric_models/exponential_models/DF_Exponential1.m @@ -10,6 +10,7 @@ function obj = DF_Exponential1(varargin) obj = obj@DiscountFunction(varargin{:}); + % TODO: this violates dependency injection, so we may want to pass these Stochastic objects in obj.theta.k = Stochastic('k'); obj = obj.parse_for_samples_and_data(varargin{:}); diff --git a/ddToolbox/models/parametric_models/exponential_models/Exponential1.m b/ddToolbox/models/parametric_models/exponential_models/Exponential1.m index 1252f5ec..6758b7dd 100644 --- a/ddToolbox/models/parametric_models/exponential_models/Exponential1.m +++ b/ddToolbox/models/parametric_models/exponential_models/Exponential1.m @@ -1,9 +1,5 @@ classdef (Abstract) Exponential1 < Parametric - properties (Access = private) - getDiscountRate % function handle - end - methods (Access = public) function obj = Exponential1(data, varargin) diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/DF_HyperbolicMagnitudeEffect.m b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/DF_HyperbolicMagnitudeEffect.m index 4e4677f4..715e159c 100644 --- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/DF_HyperbolicMagnitudeEffect.m +++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/DF_HyperbolicMagnitudeEffect.m @@ -12,6 +12,7 @@ obj = obj@DiscountFunction(varargin{:}); obj.theta = []; % clear anything from superclass + % TODO: this violates dependency injection, so we may want to pass these Stochastic objects in obj.theta.m = Stochastic('m'); obj.theta.c = Stochastic('c'); diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m index 0525ede2..0383ca06 100644 --- a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m +++ b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/Hyperbolic1MagEffect.m @@ -1,10 +1,6 @@ classdef (Abstract) Hyperbolic1MagEffect < Parametric %Hyperbolic1MagEffect Hyperbolic1MagEffect is a subclass of Model for examining the 1-parameter hyperbolic discounting function. - properties (Access = private) - getDiscountRate % function handle - end - methods (Access = public) function obj = Hyperbolic1MagEffect(data, varargin) @@ -109,6 +105,27 @@ function experimentMultiPanelFigure(obj, ind) end + + + % OVERRIDING : TODO: remove the need for this by making all + % discountFunction.plot methods have the same interface / input + % arguments. + function plot_discount_function(obj, subplot_handle, ind) + discountFunctionVariables = {obj.varList.discountFunctionParams.name}; + + subplot(subplot_handle) + discountFunction = obj.dfClass(... + 'samples', obj.coda.getSamplesAtIndex_asStruct(ind, discountFunctionVariables),... + 'data', obj.data.getExperimentObject(ind)); + + discountFunction.plot(obj.plotOptions.pointEstimateType,... + obj.plotOptions.dataPlotType,... + obj.timeUnits,... + obj.data.getMaxRewardValue(ind),... + obj.data.getMaxDelayValue(ind)); + % TODO #166 avoid having to parse these args in here + end + % MIDDLE-MAN METHOD function logk = getLogDiscountRate(obj, reward, index, varargin) diff --git a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME index 8ae9f87d..42f7eed8 100755 Binary files a/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME and b/ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchicalME differ diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/DF_Hyperbolic1.m b/ddToolbox/models/parametric_models/hyperbolic_models/DF_Hyperbolic1.m index 67fd64b9..7d425b32 100644 --- a/ddToolbox/models/parametric_models/hyperbolic_models/DF_Hyperbolic1.m +++ b/ddToolbox/models/parametric_models/hyperbolic_models/DF_Hyperbolic1.m @@ -10,6 +10,7 @@ function obj = DF_Hyperbolic1(varargin) obj = obj@DiscountFunction(varargin{:}); + % TODO: this violates dependency injection, so we may want to pass these Stochastic objects in obj.theta.logk = Stochastic('logk'); obj = obj.parse_for_samples_and_data(varargin{:}); diff --git a/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m b/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m index c2c08f25..d95ac78b 100644 --- a/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m +++ b/ddToolbox/models/parametric_models/hyperbolic_models/Hyperbolic1.m @@ -1,10 +1,6 @@ classdef (Abstract) Hyperbolic1 < Parametric %Hyperbolic1 Hyperbolic1 is a subclass of Model for examining the 1-parameter hyperbolic discounting function. - properties (Access = private) - getDiscountRate % function handle - end - methods (Access = public) function obj = Hyperbolic1(data, varargin) diff --git a/ddToolbox/utils-plot/clusterPlot.m b/ddToolbox/utils-plot/clusterPlot.m index 2c540cd7..da4e387d 100644 --- a/ddToolbox/utils-plot/clusterPlot.m +++ b/ddToolbox/utils-plot/clusterPlot.m @@ -1,4 +1,4 @@ -function clusterPlot(mcmcContainer, data, col, modelType, plotOptions, vars) +function clusterPlot(mcmcContainer, data, col, plotOptions, vars) % deal with cluster plotting of either univariate or bivariate % distributions @@ -6,9 +6,9 @@ function clusterPlot(mcmcContainer, data, col, modelType, plotOptions, vars) varNames = {vars.name}; if numel(varNames) == 1 - plot1Dclusters(mcmcContainer, data, col, modelType, plotOptions, vars); + plot1Dclusters(mcmcContainer, data, col, plotOptions, vars); elseif numel(varNames) == 2 - plot2Dclusters(mcmcContainer, data, col, modelType, plotOptions, vars); + plot2Dclusters(mcmcContainer, data, col, plotOptions, vars); else error('can only deal with plotting univariate or bivariate distributions') end diff --git a/ddToolbox/utils-plot/plot1Dclusters.m b/ddToolbox/utils-plot/plot1Dclusters.m index 63dfd7d3..c64ad847 100644 --- a/ddToolbox/utils-plot/plot1Dclusters.m +++ b/ddToolbox/utils-plot/plot1Dclusters.m @@ -1,4 +1,4 @@ -function plot1Dclusters(mcmcContainer, data, col, modelType, plotOptions, varInfo) +function plot1Dclusters(mcmcContainer, data, col, plotOptions, varInfo) % TODO: plotExpclusters.m and plotLOGKclusters are pretty much identical @@ -55,22 +55,16 @@ function plot1Dclusters(mcmcContainer, data, col, modelType, plotOptions, varInf % 'YAxisLocation','origin') drawnow -if plotOptions.shouldExportPlots - myExport(plotOptions.savePath, 'summary_plot',... - 'suffix', modelType,... - 'formats', {'png'}) -end - - function plotOpts = definePlotOptions4Participant(col) - plotOpts = {'FaceAlpha', 0.2,... + function properties = definePlotOptions4Participant(col) + properties = {'FaceAlpha', 0.2,... 'FaceColor', col,... 'LineStyle', 'none'}; end - function plotOpts = definePlotOptions4Group(col) - plotOpts = {'FaceColor', 'none',... + function properties = definePlotOptions4Group(col) + properties = {'FaceColor', 'none',... 'EdgeColor', col/2,... 'LineWidth', 3}; end -end \ No newline at end of file +end diff --git a/ddToolbox/utils-plot/plot2Dclusters.m b/ddToolbox/utils-plot/plot2Dclusters.m index 0a233480..9555c900 100644 --- a/ddToolbox/utils-plot/plot2Dclusters.m +++ b/ddToolbox/utils-plot/plot2Dclusters.m @@ -1,4 +1,4 @@ -function plot2Dclusters(mcmcContainer, data, col, modelType, plotOptions, varInfo) +function plot2Dclusters(mcmcContainer, data, col, plotOptions, varInfo) % TODO: % remove input: @@ -21,18 +21,35 @@ function plot2Dclusters(mcmcContainer, data, col, modelType, plotOptions, varInf x(:,p) = tempx; y(:,p) = tempy; end + + %plot(x(:,p),y(:,p),'.'), pause end %% plot all actual participants -mcBivariateParticipants = mcmc.BivariateDistribution(... - x(:,[1:data.getNRealExperimentFiles()]),... - y(:,[1:data.getNRealExperimentFiles()]),... - 'xLabel', varInfo(1).label,... - 'yLabel', varInfo(2).label,... - 'plotStyle','contour',... - 'probMass',probMass,... - 'pointEstimateType','mode',... - 'patchProperties',definePlotOptions4Participant(col)); +% TODO: there is a plotting error when calling this code... +% Think it's a plotting bug in mcmc.BivariateDistribution + +% mcBivariateParticipants = mcmc.BivariateDistribution(... +% x(:,[1:data.getNRealExperimentFiles()]),... +% y(:,[1:data.getNRealExperimentFiles()]),... +% 'xLabel', varInfo(1).label,... +% 'yLabel', varInfo(2).label,... +% 'plotStyle','contour',... +% 'probMass',probMass,... +% 'pointEstimateType','mode',... +% 'patchProperties',definePlotOptions4Participant(col)); +for p = 1:size(x,2) + mcBivariateParticipants = mcmc.BivariateDistribution(... + x(:,p),... + y(:,p),... + 'xLabel', varInfo(1).label,... + 'yLabel', varInfo(2).label,... + 'plotStyle','contour',... + 'probMass',probMass,... + 'pointEstimateType','mode',... + 'patchProperties',definePlotOptions4Participant(col)); +end + % TODO: enable this functionality in BivariateDistribution % % plot numbers @@ -73,21 +90,16 @@ function plot2Dclusters(mcmcContainer, data, col, modelType, plotOptions, varInf 'YAxisLocation','origin') drawnow -if plotOptions.shouldExportPlots - myExport(plotOptions.savePath, 'summary_plot',... - 'suffix', modelType,... - 'formats', {'png'}) -end - function plotOpts = definePlotOptions4Participant(col) - plotOpts = {'FaceAlpha', 0.1,... + function properties = definePlotOptions4Participant(col) + properties = {'FaceAlpha', 0.1,... 'FaceColor', col,... 'LineStyle', 'none'}; end - function plotOpts = definePlotOptions4Group(col) - plotOpts = {'FaceColor', 'none',... + function properties = definePlotOptions4Group(col) + properties = {'FaceColor', 'none',... 'EdgeColor', col,... 'LineWidth', 2}; end -end \ No newline at end of file +end diff --git a/ddToolbox/utils/asRowVector.m b/ddToolbox/utils/asRowVector.m new file mode 100644 index 00000000..df5de204 --- /dev/null +++ b/ddToolbox/utils/asRowVector.m @@ -0,0 +1,6 @@ +function output = asRowVector(input) +% coerce into a row vector +if iscolumn(input) + output = input'; +end +end