diff --git a/@PestoOptions/PestoOptions.m b/@PestoOptions/PestoOptions.m index 82b4d93..ef77e0a 100644 --- a/@PestoOptions/PestoOptions.m +++ b/@PestoOptions/PestoOptions.m @@ -172,7 +172,7 @@ % How should profiles be computed (if no more precise options are % set like profile_optim_index or profile_integ_index)? % Possibilities: {'optimization' (default), 'integration', 'mixed'} - profile_method = 'optimization'; + profile_method = 'default'; diff --git a/doc/PESTO-doc.pdf b/doc/PESTO-doc.pdf index 0f96c2a..59c90e3 100644 Binary files a/doc/PESTO-doc.pdf and b/doc/PESTO-doc.pdf differ diff --git a/examples/enzymatic_catalysis/getInitialConcentrations.m b/examples/enzymatic_catalysis/getInitialConcentrations.m deleted file mode 100644 index 16c2d1d..0000000 --- a/examples/enzymatic_catalysis/getInitialConcentrations.m +++ /dev/null @@ -1,6 +0,0 @@ -function con0 = getInitialConcentrations() - - con0 = nan(4, 1); - con0(:, 1) = [ 0.4505790; 0.4136356; 0.8565062; 1.9411520]; - - end \ No newline at end of file diff --git a/examples/enzymatic_catalysis/mainEnzymaticCatalysis.m b/examples/enzymatic_catalysis/mainEnzymaticCatalysis.m index 483e188..78e6d5d 100644 --- a/examples/enzymatic_catalysis/mainEnzymaticCatalysis.m +++ b/examples/enzymatic_catalysis/mainEnzymaticCatalysis.m @@ -9,10 +9,12 @@ % Demonstrates furthermore: % * how to do sampling without multi-start local optimization beforehand % * how to use the MEIGO toolbox for optimization (commented) -% * how to compute profile likelihoods via ODE integration +% * how to compute profile likelihoods via the integration method % * how to use plotting functions outside the get... routines % * the reliability of sampling and profiling in the case of % non-identifiabilites +% * how to use diagnosis tool (e.g. plotMCMCDiagnosis and +% plotMultiStartDiagnosis) % % This example provides a model for the reaction of a species X_1 to a % species X_4, which is catalyzed by an enzyme X_2. @@ -29,6 +31,11 @@ % optimization based on these measurements, demonstrating the use of % getMultiStarts(). % +% Parameter sampling is done first without prior information from +% optimization, then with information from optimization. +% +% Parameter optimization is done using multi-start local optimization. +% % The Profile likelihoods are calculated by integrating an ODE following % the profile path using getParameterProfiles with the option % optionsPesto.profile_method = 'integration'. @@ -99,8 +106,8 @@ % PT (with only 1 chain -> AM) specific options: optionsSampling.samplingAlgorithm = 'PT'; -optionsSampling.PT.nTemps = 5; -optionsSampling.PT.exponentT = 4; +optionsSampling.PT.nTemps = 6; +optionsSampling.PT.exponentT = 6; optionsSampling.PT.regFactor = 1e-8; optionsSampling.PT.temperatureAdaptionScheme = 'Lacki15'; %'Vousden16'; @@ -113,6 +120,8 @@ % Run the sampling parameters = getParameterSamples(parameters, objectiveFunction, optionsSampling); +% Use a diagnosis tool to see, how plotting worked (see burn-in etc.) +plotMCMCdiagnosis(parameters, 'parameters'); %% Calculate Confidence Intervals % Confidence Intervals for the Parameters are inferred from the local @@ -151,6 +160,9 @@ optionsPesto.n_starts = 5; parameters = getMultiStarts(parameters, objectiveFunction, optionsPesto); +% Use a diagnosis tool to see, how optimization worked +plotMultiStartDiagnosis(parameters); + %% Calculate Confidence Intervals % Confidence Intervals for the Parameters are inferred from the local diff --git a/getParameterConfidenceIntervals.m b/getParameterConfidenceIntervals.m index 2de6f95..866f172 100644 --- a/getParameterConfidenceIntervals.m +++ b/getParameterConfidenceIntervals.m @@ -108,7 +108,7 @@ % Confidence intervals computed using sample if isfield(parameters,'S') - parameters.CI.S(i,:,k) = prctile(parameters.S.par(i,:),50 + 100*[-alpha(k)/2, alpha(k)/2]); + parameters.CI.S(i,:,k) = prctile(parameters.S.par(i,:,1),50 + 100*[-alpha(k)/2, alpha(k)/2]); end end end diff --git a/getParameterProfiles.m b/getParameterProfiles.m index e330d81..ad77992 100644 --- a/getParameterProfiles.m +++ b/getParameterProfiles.m @@ -90,18 +90,33 @@ end % Process, which profiles should be computed in which manner +if strcmp(options.profile_method, 'default') + if (isempty(options.profile_optim_index) && isempty(options.profile_integ_index)) + options.profile_method = 'optimization'; + elseif (~isempty(options.profile_optim_index) && isempty(options.profile_integ_index)) + options.profile_method = 'optimization'; + elseif (isempty(options.profile_optim_index) && ~isempty(options.profile_integ_index)) + options.profile_method = 'integration'; + elseif (~isempty(options.profile_optim_index) && ~isempty(options.profile_integ_index)) + options.profile_method = 'mixed'; + end +end if isempty(options.parameter_index) switch options.profile_method case 'optimization' - options.profile_optim_index = 1:parameters.number; + if isempty(options.profile_optim_index) + options.profile_optim_index = 1:parameters.number; + end if ~isempty(options.profile_integ_index) - options.profile_optim_index(options.profile_integ_index) = []; + error('Some profiles seem to be computed twice. Please redefine consistent options!'); end case 'integration' - options.profile_integ_index = 1:parameters.number; + if isempty(options.profile_integ_index) + options.profile_integ_index = 1:parameters.number; + end if ~isempty(options.profile_optim_index) - options.profile_integ_index(options.profile_optim_index) = []; + error('Some profiles seem to be computed twice. Please redefine consistent options!'); end case 'mixed' @@ -113,6 +128,9 @@ if length(unique([options.profile_optim_index options.profile_integ_index])) < length([options.profile_optim_index options.profile_integ_index]) error('Some profiles seem to be computed twice. Please redefine consistent options!'); end + + otherwise + error('Unknown profile computationg method'); end options.parameter_index = sort(unique([options.profile_optim_index options.profile_integ_index])); else @@ -120,20 +138,24 @@ case 'optimization' options.profile_optim_index = options.parameter_index; if ~isempty(options.profile_integ_index) - options.profile_optim_index(options.profile_integ_index) = []; + error('Some profiles seem to be computed twice. Please redefine consistent options!'); end case 'integration' options.profile_integ_index = options.parameter_index; if ~isempty(options.profile_optim_index) - options.profile_integ_index(options.profile_optim_index) = []; + error('Some profiles seem to be computed twice. Please redefine consistent options!'); end case 'mixed' - if (length(unique([options.profile_optim_index, options.profile_integ_index])) < length(unique([options.parameter_index, options.profile_optim_index, options.profile_integ_index]))) + if (length(unique([options.profile_optim_index, options.profile_integ_index])) ~= length([options.profile_optim_index, options.profile_integ_index])) ... + || (length([options.profile_optim_index, options.profile_integ_index]) ~= length(options.parameter_index)) error('Inconsistent settings for indices in profile calculation.'); end options.parameter_index = sort(unique([options.profile_optim_index options.profile_integ_index])); + + otherwise + error('Unknown profile computationg method'); end end diff --git a/plotMCMCdiagnosis.m b/plotMCMCdiagnosis.m index ad8fb71..481cfdc 100644 --- a/plotMCMCdiagnosis.m +++ b/plotMCMCdiagnosis.m @@ -35,7 +35,7 @@ if nargin >= 2 && ~isempty(varargin{1}) type = varargin{1}; if ~max(strcmp({'parameters','log-posterior'},type)) - error('The ''type'' of plot is unknown. ''type'' can only be ''parameter'' or ''log-posterior''.') + error('The ''type'' of plot is unknown. ''type'' can only be ''parameters'' or ''log-posterior''.') end end diff --git a/private/performPT.m b/private/performPT.m index 893cca9..4298c3d 100644 --- a/private/performPT.m +++ b/private/performPT.m @@ -77,7 +77,7 @@ res.temperatures = nan(nIter, nTemps); beta = linspace(1,1/nTemps,nTemps).^exponentT; - if strcmp(temperatureAdaptionScheme,'Vousden16') + if strcmp(temperatureAdaptionScheme,'Vousden16') && nTemps > 1 beta(end) = 0; end T = ones(1,nTemps); diff --git a/tests/ConversionReactionTest.m b/tests/ConversionReactionTest.m index 1fb5da6..41e99f3 100644 --- a/tests/ConversionReactionTest.m +++ b/tests/ConversionReactionTest.m @@ -115,6 +115,7 @@ function testSamplingPT(testCase) optionsSampling = PestoSamplingOptions(); optionsSampling.nIterations = 100; + optionsSampling.mode = 'silent'; optionsSampling.samplingAlgorithm = 'PT'; optionsSampling.PT.nTemps = 3;