|
1 | | -function pStruct = getConfidenceIntervals(pStruct, alpha, type, varargin) |
| 1 | +function pStruct = getConfidenceIntervals(pStruct, confLevels, type, varargin) |
2 | 2 | % getConfidenceIntervals() calculates the confidence intervals for the |
3 | 3 | % model parameters or properties. This is done by four approaches: |
4 | 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 |
| 5 | + % a threshold according to the confidence level (calculated by a |
6 | 6 | % chi2-distribution) is reached. local_PL computes this point by a local |
7 | 7 | % approximation around the MAP estimate using the Hessian matrix, PL uses |
8 | 8 | % the profile likelihoods instead. |
9 | 9 | % The value of CI.local_B is computed by using the cummulative distribution |
10 | 10 | % function of a local approximation of the profile based on the Hessian |
11 | 11 | % matrix at the MAP estimate. |
12 | 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. |
| 13 | + % and the according percentiles based on the confidence levels. |
14 | 14 | % |
15 | 15 | % USAGE: |
16 | | - % * pStruct = getConfidenceIntervals(pStruct, alpha) |
| 16 | + % * pStruct = getConfidenceIntervals(pStruct, confLevels) |
17 | 17 | % |
18 | 18 | % Parameters: |
19 | 19 | % pStruct: parameter or properties struct |
20 | | - % alpha: vector with desired confidence levels for the intervals |
| 20 | + % confLevels: vector with desired confidence levels for the intervals |
21 | 21 | % varargin: |
22 | 22 | % options: A PestoOptions instance |
23 | 23 | % |
|
64 | 64 | end |
65 | 65 |
|
66 | 66 | % Initialization |
67 | | - pStruct.CI.alpha_levels = alpha; |
| 67 | + pStruct.CI.confLevels = confLevels; |
68 | 68 |
|
69 | | - % Loop: alpha levels |
70 | | - for iConfLevel = 1:length(alpha) |
| 69 | + % Loop: confidence levels |
| 70 | + for iConfLevel = 1:length(confLevels) |
71 | 71 | % Loop: Parameters |
72 | 72 | for iP = options.(p_index) |
73 | 73 | if isfield(pStruct,'MS') |
74 | | - pStruct = getCIfromOptimization(pStruct, alpha(iConfLevel), type, iMAP, iP, iConfLevel, options); |
| 74 | + pStruct = getCIfromOptimization(pStruct, confLevels(iConfLevel), type, iMAP, iP, iConfLevel, options); |
75 | 75 | end |
76 | 76 |
|
77 | 77 | % Confidence intervals computed using profile likelihood |
78 | 78 | if isfield(pStruct,'P') |
79 | | - pStruct = getCIfromProfiles(pStruct, alpha(iConfLevel), type, iMAP, iP, iConfLevel); |
| 79 | + pStruct = getCIfromProfiles(pStruct, confLevels(iConfLevel), type, iMAP, iP, iConfLevel); |
80 | 80 | end |
81 | 81 |
|
82 | 82 | % Confidence intervals computed using sample |
83 | 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]); |
| 84 | + pStruct.CI.S(iP,:,iConfLevel) = prctile(pStruct.S.(type)(iP,:,1),50 + 100*[-confLevels(iConfLevel)/2, confLevels(iConfLevel)/2]); |
85 | 85 | end |
86 | 86 | end |
87 | 87 | end |
88 | 88 |
|
89 | 89 | %% Output |
90 | 90 | switch options.mode |
91 | 91 | case 'visual' |
92 | | - plotConfidenceIntervals(pStruct, alpha, [], options); |
| 92 | + plotConfidenceIntervals(pStruct, [], options); |
93 | 93 | disp('-> Calculation of confidence intervals for parameters FINISHED.'); |
94 | 94 | case 'text' |
95 | 95 | disp('-> Calculation of confidence intervals for parameters FINISHED.'); |
|
99 | 99 | end |
100 | 100 |
|
101 | 101 |
|
102 | | -function pStruct = getCIfromOptimization(pStruct, alpha, type, iMAP, iP, iConfLevel, options) |
| 102 | +function pStruct = getCIfromOptimization(pStruct, confLevel, type, iMAP, iP, iConfLevel, options) |
103 | 103 |
|
104 | 104 | if strcmp(type, 'par') |
105 | 105 | % Inversion of Hessian |
|
116 | 116 |
|
117 | 117 | % Confidence intervals computed using local approximation and a |
118 | 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)); |
| 119 | + pStruct.CI.local_PL(iP,1,iConfLevel) = pStruct.MS.(type)(iP,iMAP) - sqrt(icdf('chi2',confLevel,1)*Sigma(iP,iP)); |
| 120 | + pStruct.CI.local_PL(iP,2,iConfLevel) = pStruct.MS.(type)(iP,iMAP) + sqrt(icdf('chi2',confLevel,1)*Sigma(iP,iP)); |
121 | 121 |
|
122 | 122 | % Confidence intervals computed using local approximation and the |
123 | 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))); |
| 124 | + pStruct.CI.local_B(iP,1,iConfLevel) = icdf('norm', (1-confLevel)/2,pStruct.MS.(type)(iP,iMAP),sqrt(Sigma(iP,iP))); |
| 125 | + pStruct.CI.local_B(iP,2,iConfLevel) = icdf('norm',1-(1-confLevel)/2,pStruct.MS.(type)(iP,iMAP),sqrt(Sigma(iP,iP))); |
126 | 126 |
|
127 | 127 | end |
128 | 128 |
|
129 | 129 |
|
130 | | -function pStruct = getCIfromProfiles(pStruct, alpha, type, iMAP, iP, iConfLevel) |
| 130 | +function pStruct = getCIfromProfiles(pStruct, confLevel, type, iMAP, iP, iConfLevel) |
131 | 131 |
|
132 | 132 | if ~isempty(pStruct.P(iP).(type)) |
133 | 133 | % left bound |
134 | 134 | if strcmp(type, 'par') |
135 | 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'); |
| 136 | + j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',confLevel,1)/2),1,'last'); |
137 | 137 | if ~isempty(j) |
138 | 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)); |
| 139 | + pStruct.P(iP).(type)(iP,ind([j,j+1])),exp(-icdf('chi2',confLevel,1)/2)); |
140 | 140 | else |
141 | 141 | pStruct.CI.PL(iP,1,iConfLevel) = -inf; |
142 | 142 | end |
143 | 143 | else |
144 | 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'); |
| 145 | + j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',confLevel,1)/2),1,'last'); |
146 | 146 | if ~isempty(j) |
147 | 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)); |
| 148 | + pStruct.P(iP).(type)(ind([j,j+1])),exp(-icdf('chi2',confLevel,1)/2)); |
149 | 149 | else |
150 | 150 | pStruct.CI.PL(iP,1,iConfLevel) = -inf; |
151 | 151 | end |
|
154 | 154 | % right bound |
155 | 155 | if strcmp(type, 'par') |
156 | 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'); |
| 157 | + j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',confLevel,1)/2),1,'first'); |
158 | 158 | if ~isempty(j) |
159 | 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)); |
| 160 | + pStruct.P(iP).(type)(iP,ind([j-1,j])),exp(-icdf('chi2',confLevel,1)/2)); |
161 | 161 | else |
162 | 162 | pStruct.CI.PL(iP,2,iConfLevel) = inf; |
163 | 163 | end |
164 | 164 | else |
165 | 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'); |
| 166 | + j = find(pStruct.P(iP).R(ind) <= exp(-icdf('chi2',confLevel,1)/2),1,'first'); |
167 | 167 | if ~isempty(j) |
168 | 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)); |
| 169 | + pStruct.P(iP).(type)(ind([j-1,j])),exp(-icdf('chi2',confLevel,1)/2)); |
170 | 170 | else |
171 | 171 | pStruct.CI.PL(iP,2,iConfLevel) = inf; |
172 | 172 | end |
|
0 commit comments