Skip to content

Commit f0d920d

Browse files
paulstaporpaulstapor
committed
restructuring confidence interval computation
1 parent 45ef6fd commit f0d920d

File tree

3 files changed

+191
-176
lines changed

3 files changed

+191
-176
lines changed

‎getConfidenceIntervals.m‎

Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,179 @@
1+
function pStruct = getConfidenceIntervals(pStruct, alpha, type, varargin)
2+
% getConfidenceIntervals() calculates the confidence intervals for the
3+
% model parameters or properties. This is done by four approaches:
4+
% The values of CI.local_PL and CI.PL are determined by the point on which
5+
% a threshold according to the confidence level alpha (calculated by a
6+
% chi2-distribution) is reached. local_PL computes this point by a local
7+
% approximation around the MAP estimate using the Hessian matrix, PL uses
8+
% the profile likelihoods instead.
9+
% The value of CI.local_B is computed by using the cummulative distribution
10+
% function of a local approximation of the profile based on the Hessian
11+
% matrix at the MAP estimate.
12+
% The value of CI.S is calculated using samples for the model parameters
13+
% and the according percentiles based on the confidence levels alpha.
14+
%
15+
% USAGE:
16+
% * pStruct = getConfidenceIntervals(pStruct, alpha)
17+
%
18+
% Parameters:
19+
% pStruct: parameter or properties struct
20+
% alpha: vector with desired confidence levels for the intervals
21+
% varargin:
22+
% options: A PestoOptions instance
23+
%
24+
% Return values:
25+
% pStruct: updated parameter or properties struct
26+
%
27+
% Generated fields of parameters:
28+
% CI: Information about confidence levels
29+
% * local_PL: Threshold based approach, uses a local approximation by
30+
% the Hessian matrix at the MAP estimate
31+
% (requires parameters.MS, e.g. from getMultiStarts)
32+
% * PL: Threshold based approach, uses profile likelihoods
33+
% (requires parameters.P, e.g. from getParameterProfiles)
34+
% * local_B: Mass based approach, uses a local approximation by
35+
% the Hessian matrix at the MAP estimate
36+
% (requires parameters.MS, e.g. from getMultiStarts)
37+
% * S: Bayesian approach, uses percentiles based on samples
38+
% (requires parameters.S, e.g. from getParameterSamples)
39+
40+
41+
%% Checking and assigning inputs
42+
if length(varargin) >= 1
43+
options = handleOptionArgument(varargin{1});
44+
else
45+
options = PestoOptions();
46+
end
47+
48+
% Maximum posterior index
49+
iMAP = options.MAP_index;
50+
if (isempty(iMAP))
51+
iMAP = 1;
52+
end
53+
54+
% set names for fields in pStruct
55+
if strcmp(type, 'par')
56+
p_index = 'parameter_index';
57+
else
58+
p_index = 'property_index';
59+
end
60+
61+
% parameter index
62+
if isempty(options.(p_index))
63+
options.(p_index) = 1 : pStruct.number;
64+
end
65+
66+
% Initialization
67+
pStruct.CI.alpha_levels = alpha;
68+
69+
% Loop: alpha levels
70+
for iConfLevel = 1:length(alpha)
71+
% Loop: Parameters
72+
for iP = options.(p_index)
73+
if isfield(pStruct,'MS')
74+
pStruct = getCIfromOptimization(pStruct, alpha(iConfLevel), type, iMAP, iP, iConfLevel, options);
75+
end
76+
77+
% Confidence intervals computed using profile likelihood
78+
if isfield(pStruct,'P')
79+
pStruct = getCIfromProfiles(pStruct, alpha(iConfLevel), type, iMAP, iP, iConfLevel);
80+
end
81+
82+
% Confidence intervals computed using sample
83+
if isfield(pStruct,'S')
84+
pStruct.CI.S(iP,:,iConfLevel) = prctile(pStruct.S.(type)(iP,:,1),50 + 100*[-alpha(iConfLevel)/2, alpha(iConfLevel)/2]);
85+
end
86+
end
87+
end
88+
89+
%% Output
90+
switch options.mode
91+
case 'visual'
92+
plotConfidenceIntervals(pStruct, alpha, [], options);
93+
disp('-> Calculation of confidence intervals for parameters FINISHED.');
94+
case 'text'
95+
disp('-> Calculation of confidence intervals for parameters FINISHED.');
96+
case 'silent' % no output
97+
end
98+
99+
end
100+
101+
102+
function pStruct = getCIfromOptimization(pStruct, alpha, type, iMAP, iP, iConfLevel, options)
103+
104+
if strcmp(type, 'par')
105+
% Inversion of Hessian
106+
if isempty(options.fixedParameters)
107+
Sigma = pinv(pStruct.MS.hessian(:,:,iMAP));
108+
else
109+
Sigma = nan(pStruct.number);
110+
ind = setdiff(1:pStruct.number,options.fixedParameters);
111+
Sigma(ind,ind) = pinv(pStruct.MS.hessian(ind,ind,iMAP));
112+
end
113+
else
114+
Sigma = pStruct.MS.prop_Sigma(:,:,iMAP);
115+
end
116+
117+
% Confidence intervals computed using local approximation and a
118+
% threshold (-> similar to PL-based confidence intervals)
119+
pStruct.CI.local_PL(iP,1,iConfLevel) = pStruct.MS.(type)(iP,iMAP) - sqrt(icdf('chi2',alpha,1)*Sigma(iP,iP));
120+
pStruct.CI.local_PL(iP,2,iConfLevel) = pStruct.MS.(type)(iP,iMAP) + sqrt(icdf('chi2',alpha,1)*Sigma(iP,iP));
121+
122+
% Confidence intervals computed using local approximation and the
123+
% probability mass (-> similar to Bayesian confidence intervals)
124+
pStruct.CI.local_B(iP,1,iConfLevel) = icdf('norm', (1-alpha)/2,pStruct.MS.(type)(iP,iMAP),sqrt(Sigma(iP,iP)));
125+
pStruct.CI.local_B(iP,2,iConfLevel) = icdf('norm',1-(1-alpha)/2,pStruct.MS.(type)(iP,iMAP),sqrt(Sigma(iP,iP)));
126+
127+
end
128+
129+
130+
function pStruct = getCIfromProfiles(pStruct, alpha, type, iMAP, iP, iConfLevel)
131+
132+
if ~isempty(pStruct.P(iP).(type))
133+
% left bound
134+
if strcmp(type, 'par')
135+
ind = find(pStruct.P(iP).(type)(iP,:) <= pStruct.MS.(type)(iP,iMAP));
136+
j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',alpha,1)/2),1,'last');
137+
if ~isempty(j)
138+
pStruct.CI.PL(iP,1,iConfLevel) = interp1(pStruct.P(iP).R(ind([j,j+1])),...
139+
pStruct.P(iP).(type)(iP,ind([j,j+1])),exp(-icdf('chi2',alpha,1)/2));
140+
else
141+
pStruct.CI.PL(iP,1,iConfLevel) = -inf;
142+
end
143+
else
144+
ind = find(pStruct.P(iP).(type) <= pStruct.MS.(type)(iP,1));
145+
j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',alpha,1)/2),1,'last');
146+
if ~isempty(j)
147+
pStruct.CI.PL(iP,1,iConfLevel) = interp1(pStruct.P(iP).R(ind([j,j+1])),...
148+
pStruct.P(iP).(type)(ind([j,j+1])),exp(-icdf('chi2',alpha,1)/2));
149+
else
150+
pStruct.CI.PL(iP,1,iConfLevel) = -inf;
151+
end
152+
end
153+
154+
% right bound
155+
if strcmp(type, 'par')
156+
ind = find(pStruct.P(iP).(type)(iP,:) >= pStruct.MS.(type)(iP,iMAP));
157+
j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',alpha,1)/2),1,'first');
158+
if ~isempty(j)
159+
pStruct.CI.PL(iP,2,iConfLevel) = interp1(pStruct.P(iP).R(ind([j-1,j])),...
160+
pStruct.P(iP).(type)(iP,ind([j-1,j])),exp(-icdf('chi2',alpha,1)/2));
161+
else
162+
pStruct.CI.PL(iP,2,iConfLevel) = inf;
163+
end
164+
else
165+
ind = find(pStruct.P(iP).(type) >= pStruct.MS.(type)(iP,1));
166+
j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',alpha,1)/2),1,'first');
167+
if ~isempty(j)
168+
pStruct.CI.PL(iP,2,iConfLevel) = interp1(pStruct.P(iP).R(ind([j-1,j])),...
169+
pStruct.P(iP).(type)(ind([j-1,j])),exp(-icdf('chi2',alpha,1)/2));
170+
else
171+
pStruct.CI.PL(iP,2,iConfLevel) = inf;
172+
end
173+
end
174+
175+
else
176+
pStruct.CI.PL(iP,[1,2],iConfLevel) = nan(1,2);
177+
end
178+
179+
end

‎getParameterConfidenceIntervals.m‎

Lines changed: 6 additions & 104 deletions
Original file line numberDiff line numberDiff line change
@@ -36,113 +36,15 @@
3636
% (requires parameters.MS, e.g. from getMultiStarts)
3737
% * S: Bayesian approach, uses percentiles based on samples
3838
% (requires parameters.S, e.g. from getParameterSamples)
39-
%
40-
% History:
41-
% * 2013/11/29 Jan Hasenauer
42-
% * 2016/12/01 Paul Stapor
43-
44-
%% Checking and assigning inputs
45-
if length(varargin) >= 1
46-
options = handleOptionArgument(varargin{1});
47-
else
48-
options = PestoOptions();
49-
end
50-
51-
% Maximum posterior index
52-
iMAP = options.MAP_index;
53-
if (isempty(iMAP))
54-
iMAP = 1;
55-
end
56-
57-
% parameter index
58-
if isempty(options.parameter_index)
59-
options.parameter_index = 1 : parameters.number;
60-
end
6139

62-
% Initialization
63-
parameters.CI.alpha_levels = alpha;
6440

65-
% Set default values for parameter profile
66-
if isfield(parameters,'P')
67-
for iPar = 1:parameters.number
68-
if iPar > length(parameters.P)
69-
parameters.P(iPar).par = [];
70-
parameters.P(iPar).logPost = [];
71-
parameters.P(iPar).R = [];
72-
parameters.P(iPar).t_cpu = [];
73-
parameters.P(iPar).optSteps = [];
74-
parameters.P(iPar).intSteps = [];
75-
parameters.P(iPar).reOptSteps = [];
76-
end
41+
%% Checking and assigning inputs
42+
if length(varargin) >= 1
43+
options = handleOptionArgument(varargin{1});
44+
else
45+
options = PestoOptions();
7746
end
78-
end
79-
80-
% Loop: alpha levels
81-
for k = 1:length(alpha)
82-
% Loop: Parameters
83-
for iPar = options.parameter_index
84-
if isfield(parameters,'MS')
85-
% Inversion of Hessian
86-
if isempty(options.fixedParameters)
87-
Sigma = pinv(parameters.MS.hessian(:,:,iMAP));
88-
else
89-
Sigma = nan(parameters.number);
90-
ind = setdiff(1:parameters.number,options.fixedParameters);
91-
Sigma(ind,ind) = pinv(parameters.MS.hessian(ind,ind,iMAP));
92-
end
93-
94-
% Confidence intervals computed using local approximation and a
95-
% threshold (-> similar to PL-based confidence intervals)
96-
parameters.CI.local_PL(iPar,1,k) = parameters.MS.par(iPar,iMAP) - sqrt(icdf('chi2',alpha(k),1)*Sigma(iPar,iPar));
97-
parameters.CI.local_PL(iPar,2,k) = parameters.MS.par(iPar,iMAP) + sqrt(icdf('chi2',alpha(k),1)*Sigma(iPar,iPar));
98-
99-
% Confidence intervals computed using local approximation and the
100-
% probability mass (-> similar to Bayesian confidence intervals)
101-
parameters.CI.local_B(iPar,1,k) = icdf('norm', (1-alpha(k))/2,parameters.MS.par(iPar,iMAP),sqrt(Sigma(iPar,iPar)));
102-
parameters.CI.local_B(iPar,2,k) = icdf('norm',1-(1-alpha(k))/2,parameters.MS.par(iPar,iMAP),sqrt(Sigma(iPar,iPar)));
103-
end
104-
105-
% Confidence intervals computed using profile likelihood
106-
if isfield(parameters,'P')
107-
if ~isempty(parameters.P(iPar).par)
108-
% left bound
109-
ind = find(parameters.P(iPar).par(iPar,:) <= parameters.MS.par(iPar,iMAP));
110-
j = find(parameters.P(iPar).R(ind) <= exp(-icdf('chi2',alpha(k),1)/2),1,'last');
111-
if ~isempty(j)
112-
parameters.CI.PL(iPar,1,k) = interp1(parameters.P(iPar).R(ind([j,j+1])),...
113-
parameters.P(iPar).par(iPar,ind([j,j+1])),exp(-icdf('chi2',alpha(k),1)/2));
114-
else
115-
parameters.CI.PL(iPar,1,k) = -inf;
116-
end
117-
% right bound
118-
ind = find(parameters.P(iPar).par(iPar,:) >= parameters.MS.par(iPar,iMAP));
119-
j = find(parameters.P(iPar).R(ind) <= exp(-icdf('chi2',alpha(k),1)/2),1,'first');
120-
if ~isempty(j)
121-
parameters.CI.PL(iPar,2,k) = interp1(parameters.P(iPar).R(ind([j-1,j])),...
122-
parameters.P(iPar).par(iPar,ind([j-1,j])),exp(-icdf('chi2',alpha(k),1)/2));
123-
else
124-
parameters.CI.PL(iPar,2,k) = inf;
125-
end
126-
else
127-
parameters.CI.PL(iPar,[1,2],k) = nan(1,2);
128-
end
129-
end
130-
131-
% Confidence intervals computed using sample
132-
if isfield(parameters,'S')
133-
parameters.CI.S(iPar,:,k) = prctile(parameters.S.par(iPar,:,1),50 + 100*[-alpha(k)/2, alpha(k)/2]);
134-
end
135-
end
136-
end
13747

138-
%% Output
139-
switch options.mode
140-
case 'visual'
141-
plotConfidenceIntervals(parameters, alpha, [], options);
142-
disp('-> Calculation of confidence intervals for parameters FINISHED.');
143-
case 'text'
144-
disp('-> Calculation of confidence intervals for parameters FINISHED.');
145-
case 'silent' % no output
146-
end
48+
parameters = getConfidenceIntervals(parameters, alpha, 'par', options);
14749

14850
end

‎getPropertyConfidenceIntervals.m‎

Lines changed: 6 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -32,81 +32,15 @@
3232
% * local_B: Mass based approach, uses a local approximation by
3333
% the Hessian matrix at the MAP estimate
3434
% (requires parameters.MS, e.g. from getMultiStarts)
35-
%
36-
% History:
37-
% * 2013/11/29 Jan Hasenauer
38-
% * 2016/12/01 Paul Stapor
39-
40-
%% Checking and assigning inputs
41-
if length(varargin) >= 1
42-
options = handleOptionArgument(varargin{1});
43-
else
44-
options = PestoOptions();
45-
end
4635

47-
% Initialization
48-
properties.CI.alpha_levels = alpha;
49-
50-
if isempty(options.property_index)
51-
options.property_index = 1 : properties.number;
52-
end
5336

54-
% Loop: alpha levels
55-
for k = 1:length(alpha)
56-
% Loop: properties
57-
for iProp = options.property_index
58-
% Hessian
59-
Sigma = properties.MS.prop_Sigma(:,:,1);
60-
61-
% Confidence intervals computed using local approximation and a
62-
% threshold (-> similar to PL-based confidence intervals)
63-
properties.CI.local_PL(iProp,1,k) = properties.MS.prop(iProp) - sqrt(icdf('chi2',alpha(k),1)*Sigma(iProp,iProp));
64-
properties.CI.local_PL(iProp,2,k) = properties.MS.prop(iProp) + sqrt(icdf('chi2',alpha(k),1)*Sigma(iProp,iProp));
65-
66-
% Confidence intervals computed using local approximation and the
67-
% probability mass (-> similar to Bayesian confidence intervals)
68-
properties.CI.local_B(iProp,1,k) = icdf('norm', (1-alpha(k))/2,properties.MS.prop(iProp),sqrt(Sigma(iProp,iProp)));
69-
properties.CI.local_B(iProp,2,k) = icdf('norm',1-(1-alpha(k))/2,properties.MS.prop(iProp),sqrt(Sigma(iProp,iProp)));
70-
71-
% Confidence intervals computed using profile likelihood
72-
if isfield(properties,'P')
73-
if ~isempty(properties.P(iProp).prop)
74-
% left bound
75-
ind = find(properties.P(iProp).prop <= properties.MS.prop(iProp,1));
76-
j = find(properties.P(iProp).R(ind) <= exp(-icdf('chi2',alpha(k),1)/2),1,'last');
77-
if ~isempty(j)
78-
properties.CI.PL(iProp,1,k) = interp1(properties.P(iProp).R(ind([j,j+1])),...
79-
properties.P(iProp).prop(ind([j,j+1])),exp(-icdf('chi2',alpha(k),1)/2));
80-
else
81-
properties.CI.PL(iProp,1,k) = -inf;
82-
end
83-
% right bound
84-
ind = find(properties.P(iProp).prop >= properties.MS.prop(iProp,1));
85-
j = find(properties.P(iProp).R(ind) <= exp(-icdf('chi2',alpha(k),1)/2),1,'first');
86-
if ~isempty(j)
87-
properties.CI.PL(iProp,2,k) = interp1(properties.P(iProp).R(ind([j-1,j])),...
88-
properties.P(iProp).prop(ind([j-1,j])),exp(-icdf('chi2',alpha(k),1)/2));
89-
else
90-
properties.CI.PL(iProp,2,k) = inf;
91-
end
92-
end
93-
end
94-
95-
% Confidence intervals computed using sample
96-
if isfield(properties,'S')
97-
properties.CI.S(iProp,:,k) = prctile(properties.S.prop(iProp,:),50 + 100*[-alpha(k)/2, alpha(k)/2]);
98-
end
37+
%% Checking and assigning inputs
38+
if length(varargin) >= 1
39+
options = handleOptionArgument(varargin{1});
40+
else
41+
options = PestoOptions();
9942
end
100-
end
10143

102-
%% Output
103-
switch options.mode
104-
case 'visual'
105-
plotConfidenceIntervals(properties, alpha, [], options);
106-
disp('-> Calculation of confidence intervals for properties FINISHED.');
107-
case 'text'
108-
disp('-> Calculation of confidence intervals for properties FINISHED.');
109-
case 'silent' % no output
110-
end
44+
properties = getConfidenceIntervals(properties, alpha, 'prop', options);
11145

11246
end

0 commit comments

Comments
 (0)