Line data Source code
1 : /*******************************************************************************
2 : * *
3 : * Cosmos:(C)oncept et (O)utils (S)tatistique pour les (Mo)deles *
4 : * (S)tochastiques *
5 : * *
6 : * Copyright (C) 2009-2012 LSV & LACL *
7 : * Authors: Paolo Ballarini Benoît Barbot & Hilal Djafri *
8 : * Website: http://www.lsv.ens-cachan.fr/Software/cosmos *
9 : * *
10 : * This program is free software; you can redistribute it and/or modify *
11 : * it under the terms of the GNU General Public License as published by *
12 : * the Free Software Foundation; either version 3 of the License, or *
13 : * (at your option) any later version. *
14 : * *
15 : * This program is distributed in the hope that it will be useful, *
16 : * but WITHOUT ANY WARRANTY; without even the implied warranty of *
17 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18 : * GNU General Public License for more details. *
19 : * *
20 : * You should have received a copy of the GNU General Public License along *
21 : * with this program; if not, write to the Free Software Foundation, Inc., *
22 : * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. *
23 : * file parameters.hpp *
24 : * Created by Benoit Barbot on 12/12/12. *
25 : *******************************************************************************
26 : */
27 :
28 : #include <iostream>
29 : #include <limits>
30 : #include <boost/math/distributions/normal.hpp>
31 : #include <boost/math/distributions/binomial.hpp>
32 :
33 : #include "HaslFormula.hpp"
34 : #include "parameters.hpp"
35 :
36 : /**
37 : * Trivial confidence interval containing all R.
38 : */
39 215 : ConfInt::ConfInt(){
40 215 : mean = 0;
41 215 : low = - std::numeric_limits<double>::infinity();
42 215 : up = std::numeric_limits<double>::infinity();
43 215 : min = - std::numeric_limits<double>::infinity();
44 215 : max = std::numeric_limits<double>::infinity();
45 215 : conf = 1;
46 215 : }
47 :
48 : /**
49 : * Symetric confidence interval
50 : * @param meanArg is the middle of the confidence interval
51 : * @param width is the width of the confidence interval
52 : * @param min the smallest value observed
53 : * @param max the maximal value observed
54 : */
55 1916 : ConfInt::ConfInt(double meanArg,double width,double minval,double maxval,double c){
56 1916 : mean = meanArg;
57 1916 : low = meanArg - width/2;
58 1916 : up = meanArg + width/2;
59 1916 : min= minval;
60 1916 : max= maxval;
61 1916 : conf=c;
62 1916 : }
63 :
64 : /**
65 : * Asymetric confidence interval.
66 : * @param meanArg is the mean of the confidence interval it is not necessary
67 : * the middle of the confidence interval if the distribution is not symetric
68 : * and is not approximated by the normal distribution wich is symetric.
69 : * @param lowArg the Lower bound of the confidence interval.
70 : * @param upArg the Upper bound of the confidence interval.
71 : * @param min the smallest value observed
72 : * @param max the maximal value observed
73 : */
74 87 : ConfInt::ConfInt(double meanArg,double lowArg,double upArg,double minval,double maxval,double c){
75 87 : mean = meanArg;
76 87 : low = lowArg;
77 87 : up = upArg;
78 87 : min= minval;
79 87 : max= maxval;
80 87 : conf =c;
81 87 : }
82 :
83 2235 : ConfInt::~ConfInt(){}
84 :
85 : /**
86 : * Return the width of a confidence interval
87 : */
88 4349 : double ConfInt::width(){
89 4349 : return (up-low);
90 : }
91 :
92 93 : ConfInt & ConfInt::operator+=(const ConfInt& rhs){
93 93 : low += rhs.low;
94 93 : mean += rhs.mean;
95 93 : up += rhs.up;
96 93 : min += rhs.min;
97 93 : max += rhs.max;
98 93 : conf = conf +rhs.conf -1;
99 93 : return *this;
100 : }
101 :
102 0 : ConfInt & ConfInt::operator-=(const ConfInt& rhs){
103 0 : low -= rhs.up;
104 0 : mean -= rhs.mean;
105 0 : up -= rhs.low;
106 0 : min -= rhs.max;
107 0 : max -= rhs.min;
108 0 : conf = conf +rhs.conf -1;
109 0 : return *this;
110 : }
111 :
112 173 : ConfInt & ConfInt::operator*=(const ConfInt& rhs){
113 173 : mean = mean*rhs.mean;
114 346 : low = fmin(fmin(low*rhs.low,up*rhs.up),
115 346 : fmin(low*rhs.up,up*rhs.low));
116 346 : up = fmax(fmax(low*rhs.low,up*rhs.up),
117 346 : fmax(low*rhs.up,up*rhs.low));
118 346 : min = fmin(fmin(min*rhs.min,max*rhs.max),
119 346 : fmin(min*rhs.max,max*rhs.min));
120 346 : max = fmax(fmax(min*rhs.min,max*rhs.max),
121 346 : fmax(min*rhs.max,max*rhs.min));
122 173 : conf = conf +rhs.conf -1;
123 173 : return *this;
124 : }
125 :
126 5 : ConfInt & ConfInt::operator/=(const ConfInt& rhs){
127 5 : mean /= rhs.mean;
128 :
129 5 : if(rhs.low * rhs.up<=0){
130 0 : low = - std::numeric_limits<double>::infinity();
131 0 : up = std::numeric_limits<double>::infinity();
132 : } else {
133 5 : if(rhs.low > 0){
134 5 : low /= rhs.up;
135 5 : up /= rhs.low;
136 : }else{
137 0 : low /= rhs.low;
138 0 : up /= rhs.up;
139 : }
140 : }
141 5 : if(rhs.min * rhs.max<=0){
142 0 : min = - std::numeric_limits<double>::infinity();
143 0 : max = std::numeric_limits<double>::infinity();
144 : } else {
145 5 : if(rhs.min > 0){
146 5 : min /= rhs.max;
147 5 : max /= rhs.min;
148 : }else{
149 0 : min /= rhs.min;
150 0 : max /= rhs.max;
151 : }
152 : }
153 5 : conf = conf +rhs.conf -1;
154 5 : return *this;
155 : }
156 :
157 : /**
158 : * Build an HASL formula with the PROB,EXISTS,ALLS operators.
159 : * Compute exact confidence interval.
160 : */
161 9 : HaslFormulasTop::HaslFormulasTop(const HaslType t){
162 9 : assert( t==PROBABILITY || t==EXISTS || t== NOTALLS);
163 9 : TypeOp = t;
164 9 : Level =0;
165 9 : Value =0;
166 9 : Value2=0;
167 9 : Algebraic = 0;
168 9 : left = NULL;
169 9 : right = NULL;
170 9 : }
171 :
172 : /**
173 : * Build an HASL formula with the PROB,EXISTS,ALLS operators.
174 : * Compute exact confidence interval.
175 : */
176 0 : HaslFormulasTop::HaslFormulasTop(const HaslType t,size_t al){
177 0 : assert( t==PROBCOND);
178 0 : TypeOp = t;
179 0 : Level =0;
180 0 : Value =0;
181 0 : Value2=0;
182 0 : Algebraic = al;
183 0 : left = NULL;
184 0 : right = NULL;
185 0 : }
186 :
187 : /**
188 : * Build a CONSTANT Hasl formula.
189 : * @param v the value of the constant
190 : */
191 3 : HaslFormulasTop::HaslFormulasTop(double v1,double v2,double l){
192 3 : TypeOp = CONSTANT;
193 3 : Level = l;
194 3 : Value = v1;
195 3 : Value2= v2;
196 3 : Algebraic = 0;
197 3 : left = NULL;
198 3 : right = NULL;
199 3 : }
200 :
201 : /**
202 : * Build an HASL formula computing the expected value of a subformula.
203 : * Use the approximation to the normal low to compute the confidence interval
204 : * @param al the index of the sub formula return by the simulator.
205 : * @param l the confidence level of the confidence interval.
206 : */
207 203 : HaslFormulasTop::HaslFormulasTop(size_t al){
208 203 : TypeOp = EXPECTANCY;
209 203 : Level = 0;
210 203 : Value = 0;
211 203 : Value2 = 0;
212 203 : Algebraic = al;
213 203 : left = NULL;
214 203 : right = NULL;
215 203 : }
216 :
217 : /**
218 : * Build an HASL formula testing if the probability of accepting a run is above some threshold.
219 : * In this setting the confidence level is the probability of type
220 : * I and type II errors.
221 : * @param p is the threshold.
222 : * @param delta is the half width of the indiference region.
223 : */
224 0 : HaslFormulasTop::HaslFormulasTop(double p, double delta){
225 0 : TypeOp = HYPOTHESIS ;
226 0 : Level = 0;
227 0 : Value = p -delta;
228 0 : Value2= p + delta;
229 0 : Algebraic = 0;
230 0 : left = NULL;
231 0 : right = NULL;
232 0 : }
233 :
234 :
235 : /*
236 : * Build an Hasl formula of a given type as a tree with two descendant.
237 : * @param t is the type of the node either HASL_PLUS, HASL_TIME or HASL_DIV.
238 : * @param l The left Hasl Formula.
239 : * @param r the right Hasl formula.
240 : */
241 3 : HaslFormulasTop::HaslFormulasTop(HaslType t, HaslFormulasTop* l,HaslFormulasTop* r){
242 3 : TypeOp = t;
243 3 : Algebraic = 0;
244 3 : Level = 0;
245 3 : Value = 0;
246 3 : Value2= 0;
247 3 : left = l;
248 3 : right = r;
249 3 : }
250 :
251 : /**
252 : * Destructor call destructor recusively on descendant if neaded.
253 : */
254 0 : HaslFormulasTop::~HaslFormulasTop(){
255 0 : if(!left) delete left;
256 0 : if(!right) delete right;
257 0 : }
258 :
259 : /**
260 : * This method return true if a haslformula is constant.
261 : * This is usfull to know if a subformula have a confidence level
262 : * different than 1.
263 : */
264 6 : bool HaslFormulasTop::isConstant(){
265 6 : switch (TypeOp) {
266 : case CONSTANT:
267 4 : return true;
268 :
269 : case HASL_PLUS:
270 : case HASL_TIME:
271 : case HASL_DIV:
272 1 : return (left->isConstant() && right->isConstant());
273 : default:
274 1 : return false;
275 : }
276 : }
277 :
278 : /**
279 : * This function set the confidence level of an HASL formula.
280 : * When the fomula contain
281 : * @param l The confidence level
282 : */
283 217 : void HaslFormulasTop::setLevel(double l){
284 217 : switch (TypeOp) {
285 : case HYPOTHESIS:
286 : case PROBABILITY:
287 : case PROBCOND:
288 7 : Level = l;
289 7 : break;
290 :
291 : case EXPECTANCY:
292 : case CDF_PART:
293 : case PDF_PART:
294 : case RE_Likelyhood:
295 201 : Level = l;
296 201 : Value = quantile(boost::math::normal() , 0.5 + l / 2.0);
297 201 : break;
298 :
299 : case RE_AVG://temporary
300 4 : Level = l;
301 : //Here we divide the error by two because half the error come from the bernoulli evaluation.
302 4 : Value = quantile(boost::math::normal() , 0.75 + l / 4.0);
303 4 : break;
304 :
305 : case RE_Continuous:
306 2 : Level = l;
307 2 : Value = quantile(boost::math::normal() , 0.75 + l / 4.0);
308 2 : break;
309 :
310 : case CONSTANT:
311 0 : break;
312 :
313 : case HASL_PLUS:
314 : case HASL_MINUS:
315 : case HASL_TIME:
316 : case HASL_DIV:
317 3 : Level=l;
318 3 : if (left->isConstant()){right->setLevel(l);}
319 1 : else if (right->isConstant()){left->setLevel(l);}
320 : else{
321 : //We need to increase the condidence level on each subformula
322 0 : left->setLevel((1+l)/2);
323 0 : right->setLevel((1+l)/2);
324 : }
325 3 : break;
326 : case NOTALLS:
327 : case EXISTS:
328 0 : Level =1;
329 0 : break;
330 :
331 : default:
332 0 : std::cerr << "Fail to parse Hasl Formula"<< std::endl;
333 0 : exit(EXIT_FAILURE);
334 : }
335 217 : }
336 :
337 :
338 0 : void HaslFormulasTop::computeChernoff(){
339 0 : double b = 0.0;
340 0 : for(let it : P.HaslFormulas) b = fmax(b,it->bound());
341 0 : if(b== INFINITY){
342 0 : std::cerr << "Cannot use Chernoff-Hoeffding bounds: no bounds on the computed value" << std::endl;
343 0 : exit(EXIT_FAILURE);
344 : }
345 0 : if(P.MaxRuns== (unsigned long)-1){
346 0 : P.MaxRuns = (int)(2.0*2.0*2.0*b*b/(P.Width*P.Width) * log(2/(1-P.Level)));
347 0 : }else if(P.Width == 0){
348 0 : P.Width = 2.0*b * sqrt( (2.0/P.MaxRuns) * log(2.0/(1.0-P.Level)));
349 0 : }else if(P.Level ==0){
350 0 : P.Level = (1.0 - (2.0* fmin(0.5 ,exp( (double)P.MaxRuns * P.Width*P.Width / (-2.0*2.0*2.0*b*b)))));
351 : }
352 0 : }
353 :
354 :
355 : /**
356 : * The function eval compute confidence interval for the HASL formula
357 : * from the simualation result return by the simulators.
358 : * It always return a confdience interval containing a confidence level
359 : * Depending of the context it should be interpreted differently.
360 : * In sequential setting, one should check that the intervals satisfy the
361 : * stopping criteria.
362 : * In static setting, one should check independantly the number of sample to
363 : * decide whether the simulation should be stopped.
364 : *
365 : * Some remarks about the estimation of the confidence interval adopted here
366 : * Let \f$ l \f$ =ConfLevel, the confidence level
367 : * \f$ l=1-alpha \f$
368 : * Let \f$ w \f$ = ConfWidth, the size of the confidence interval
369 : */
370 2009 : ConfInt HaslFormulasTop::eval(const BatchR &batch)const{
371 2009 : switch (TypeOp) {
372 : case PROBABILITY:
373 :
374 : //if (P.Width ==0)// In static setting use exact binomial computation
375 : {
376 : /*
377 : * Here we used the boost library for computing the binomial
378 : * confidence interval.
379 : * According to boost documentation the algorithme used is
380 : * Clopper-person:
381 : * Clopper, C. J. and Pearson, E. S. (1934). The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika 26 404-413.
382 : */
383 77 : double mean = (double)batch.Isucc / (double)batch.I;
384 77 : double l = boost::math::binomial_distribution<>::find_lower_bound_on_p(batch.I,batch.Isucc, (1-Level)/2);
385 77 : double u = boost::math::binomial_distribution<>::find_upper_bound_on_p(batch.I,batch.Isucc, (1-Level)/2);
386 :
387 77 : if(batch.Isucc ==0)return ConfInt(mean,l,u,0.0,0.0,Level);
388 75 : if(batch.Isucc == batch.I)return ConfInt(mean,l,u,1.0,1.0,Level);
389 :
390 75 : return ConfInt(mean,l,u,0.0,1.0,Level);
391 : } // In Sequential case fall back to show robin algorithm
392 :
393 : // NO BREAK on purpose here !
394 : case PROBCOND:
395 : {
396 : //Same as PROB but with a more specific constraint
397 0 : double mean = (double)batch.bernVar[Algebraic] / (double)batch.Isucc;
398 0 : double l = boost::math::binomial_distribution<>::find_lower_bound_on_p(batch.Isucc,batch.bernVar[Algebraic], (1-Level)/2);
399 0 : double u = boost::math::binomial_distribution<>::find_upper_bound_on_p(batch.Isucc,batch.bernVar[Algebraic], (1-Level)/2);
400 :
401 0 : if(batch.bernVar[Algebraic] ==0)return ConfInt(mean,l,u,0.0,0.0,Level);
402 0 : if(batch.bernVar[Algebraic] == batch.Isucc)return ConfInt(mean,l,u,1.0,1.0,Level);
403 :
404 0 : return ConfInt(mean,l,u,0.0,1.0,Level);
405 :
406 : }
407 :
408 : /*
409 : * Here Gaussian confidence interval are computed.
410 : * Let \f$ \mu \f$ the value to estimate, and \f$ x \f$ the
411 : * estimation of \f$ \mu \f$
412 : * then \f$ \mathbb{P}( \mu \in [x-\frac{w}{2} , x+\frac{w}{2}]) = 1-alpha \f$
413 : *
414 : * The gaussian confidence interval is given by :
415 : * \f$ [x-z(1-\frac{alpha}{2}) * \frac{StandardDeviation~ }
416 : * { \sqrt{NumberOfObservations}} ,
417 : * x+z(1-\frac{alpha}{2}) * \frac{StandardDeviation~}
418 : * { \sqrt{NumberOfObservations}}] \f$
419 : * with : \f$ z(1-\frac{alpha}{2})=z(1-\frac{1-l}{2}) = z(0.5+\frac{l}{2}) \f$
420 : *
421 : * \f$ StandartDeviation~ = \sqrt{ Variance +\frac{1}{n} } \f$
422 : * This correction come from the Chows and Robbin algorithm to ensure
423 : * The correctness of the stopping condition.
424 : *
425 : * We use the boost library: quantile(boost::math::normal() , level)
426 : * to compute quantile of the normal low. Its value have been stored
427 : * in Value by the function setLevel.
428 : */
429 : case EXPECTANCY:
430 : case CDF_PART:
431 : case PDF_PART:
432 : {
433 : //No accepting trajectory the condidence interval is R.
434 1884 : if(batch.Isucc==0)return ConfInt();
435 :
436 1883 : double mean = batch.Mean[Algebraic]/batch.Isucc;
437 1883 : double m2 = batch.M2[Algebraic]/batch.Isucc;
438 1883 : double variance = m2 - mean * mean;
439 :
440 1883 : if(P.Width!=0){ //P.width =0 indicate static setting
441 1883 : variance += 1.0/batch.Isucc;
442 : //Here the +1 come from the Chows and Robbin algorithm in Sequential setting
443 : }
444 :
445 1883 : double width = 2 * Value * sqrt(variance/(batch.Isucc -1));
446 :
447 1883 : return ConfInt(mean,width,batch.Min[Algebraic],batch.Max[Algebraic],Level);
448 : }
449 :
450 : case RE_Likelyhood:
451 : {
452 : //No accepting trajectory the condidence interval is R.
453 8 : if(batch.Isucc==0)return ConfInt();
454 8 : double mean = batch.Mean[Algebraic]/batch.Isucc;
455 8 : double m2 = batch.M2[Algebraic]/batch.Isucc;
456 :
457 : //If the variance is equal to 0 numerical instability can make it a negative number.
458 8 : double variance = fmax(0,m2 - mean * mean);
459 :
460 8 : double width = 2 * Value * sqrt(variance/batch.Isucc);
461 :
462 8 : return ConfInt(mean,width,batch.Min[Algebraic],batch.Max[Algebraic],Level);
463 : }
464 :
465 : case RE_AVG://temporary
466 : {
467 : //No accepting trajectory => the condidence interval is R.
468 8 : if(batch.Isucc==0)return ConfInt();
469 :
470 8 : double mean = (double)batch.Isucc / (double)batch.I;
471 8 : double l = boost::math::binomial_distribution<>::find_lower_bound_on_p(batch.I,batch.Isucc, (1-Level)/2);
472 8 : double u = boost::math::binomial_distribution<>::find_upper_bound_on_p(batch.I,batch.Isucc, (1-Level)/2);
473 8 : double mean2 = batch.Mean[Algebraic]/batch.Isucc;
474 8 : double m2 = batch.M2[Algebraic]/batch.Isucc;
475 8 : double variance = fmax(0,m2 - mean2 * mean2);
476 8 : double width = Value * sqrt(variance/batch.Isucc);
477 :
478 : return ConfInt(mean*mean2,
479 8 : (mean2 - width)*l,
480 16 : (mean2 + width)*u,0.0,1.0 * batch.Max[Algebraic],Level );
481 : }
482 :
483 : case RE_Continuous:
484 : {
485 : /*
486 : * Only used by the rare event engine in for countinuous bounded
487 : * Until property. This case apply the uniformization method.
488 : * The simulator returns the estimate of each \mu_u as well
489 : * as the Poisson probabilities computed using Fox-glynn algorithm.
490 : *
491 : */
492 : //batch.print();
493 2 : size_t N = batch.Mean.size()/2;
494 2 : std::cout << "tablelength = "<< batch.Mean.size() << std::endl;
495 2 : double error = 1- Level;
496 2 : double LevelN = 1- error/(N+2);
497 2 : double ValueN = quantile(boost::math::normal() , (1+LevelN)/2);
498 :
499 4 : ConfInt totalInt = ConfInt(0.0,0.0,0.0,0.0,Level);
500 4 : ConfInt likelyhood = ConfInt(0.0,0.0,0.0,0.0,Level);
501 4 : ConfInt reacheabilityprob = ConfInt(0.0,0.0,0.0,0.0,Level);
502 4 : ConfInt poisson = ConfInt(0.0,0.0,0.0,0.0,Level);
503 :
504 2 : if(P.verbose>1)std::cerr << "i,\tsucc,\tbatch,\tMean,\tM2,\tMin,\tMax,\tPoisson" << std::endl;
505 2 : if(P.verbose>1)std::cerr << "continuous step:\t" << P.continuousStep << std::endl;
506 86 : for(size_t i=0; i< N; i++){
507 84 : int i2 = (i>0 && P.continuousStep>1.0 )? i-1 : i;
508 :
509 : //evaluate the likelyhood:
510 84 : likelyhood.mean = batch.Mean[2*i+1]/batch.Mean[2*i];
511 84 : double var = (batch.M2[2*i+1]/batch.Mean[2*i]) - pow((batch.Mean[2*i+1]/batch.Mean[2*i]),2);
512 84 : double widthN = 0.0;
513 84 : if(var>0)widthN = 2* ValueN * sqrt(var/batch.Mean[2*i]);
514 84 : likelyhood.low = batch.Min[2*i2+1];
515 84 : likelyhood.up = batch.Max[2*i+1]; //likelyhood->mean + widthN/2.0;
516 84 : likelyhood.min = batch.Min[2*i2+1];
517 84 : likelyhood.max = batch.Max[2*i+1];
518 :
519 : //evaluate probability to reach final state:
520 84 : reacheabilityprob.mean = batch.Mean[2*i]/batch.M2[2*i];
521 84 : reacheabilityprob.low = boost::math::binomial_distribution<>::find_lower_bound_on_p(batch.M2[2*i2],batch.Mean[2*i2], (1-LevelN)/2);
522 84 : reacheabilityprob.up = boost::math::binomial_distribution<>::find_upper_bound_on_p(batch.M2[2*i],batch.Mean[2*i], (1-LevelN)/2);
523 84 : reacheabilityprob.min = 0.0;
524 84 : reacheabilityprob.max = 1.0;
525 : //evaluate poisson:
526 84 : poisson.mean = batch.Min[2*i];
527 84 : poisson.low = (1-Value2)*batch.Min[2*i];
528 84 : poisson.up = batch.Min[2*i];
529 84 : poisson.min=poisson.low;
530 84 : poisson.max=poisson.max;
531 :
532 84 : if(P.verbose>1){
533 84 : std::cout << "i: " << i<< "\tsucc: "<< batch.Mean[2*i]<< "\tLikelyhood: " << likelyhood.mean << "\t[" << likelyhood.low<<";"<<likelyhood.up << "]\twidth: "<< widthN <<"\tpoisson: " << batch.Min[2*i] << "\tconfint: ["<< poisson.low <<";"<< poisson.up << "]";
534 84 : std::cerr << i << ",\t" << batch.Mean[2*i] << ",\t" << batch.M2[2*i] << ",\t" << batch.Mean[2*i+1]/batch.Mean[2*i] << ",\t" << batch.M2[2*i+1]/batch.Mean[2*i] << ",\t" << batch.Min[2*i+1] << ",\t" << batch.Max[2*i+1] << ",\t" << batch.Mean[2*i] << std::endl;
535 : }
536 :
537 : //Multiply the three confidence interval:
538 84 : poisson *= likelyhood;
539 84 : poisson *= reacheabilityprob;
540 :
541 84 : if(P.verbose>1){
542 84 : std::cout << "\tresult: ["<< poisson.low <<";"<< poisson.up << "]" << std::endl;
543 : }
544 :
545 :
546 : //Add the confidence interval to the total one.
547 : //if(poisson->width() < 1.0)
548 84 : totalInt += poisson;
549 :
550 : }
551 :
552 : //Add the error made on the left of the truncation point.
553 : ConfInt leftTruncationError =
554 4 : ConfInt(0.0,0.0,batch.Mean[1]/batch.M2[0]* Value2/2,0,batch.Mean[1]/batch.M2[0]* Value2/2,Level);
555 :
556 : //Add the error made on the rigth trucation point TODO!
557 4 : ConfInt rightTruncationError = ConfInt(0.0,0.0,0.0,0.0,Level);
558 :
559 2 : totalInt += leftTruncationError;
560 2 : totalInt += rightTruncationError;
561 :
562 2 : return totalInt;
563 : }
564 :
565 : case HYPOTHESIS:
566 : {
567 : /*
568 : * Implementation of the SPRT.
569 : * Value < Value2
570 : * Value <= p <= Value2 => Indiferrence region
571 : * H0 = { Value2 < p <= 1 } => p0=Value2
572 : * H1 = { 0 <= p < Value } => p1=Value
573 : *
574 : * To comply with the architecture of cosmos requiring confidence
575 : * interval this function also return confidence interval.
576 : * Let l be the minimal confidence level for wich either H0 or H1 is
577 : * accepted
578 : * return [0,Value2] with confidence level l when H0 is accepted
579 : * return [Value, 1] with confidence level l when H1 is accepted
580 : *
581 : * The stopping criterion is reached when l is greater than the expected
582 : * level. This stopping criterion is exaclty the one of
583 : * A. Wald. Sequential tests of statistical hypotheses.
584 : * The Annals of Mathematical Statistics, 16(2):117–186, 06 1945.
585 : */
586 :
587 0 : const double alog = log((1-Value)/(1-Value2));
588 0 : const double blog = log(Value/Value2);
589 :
590 0 : double loglikelyhoodRatio = alog * (double)(batch.I - batch.Isucc) + blog * (double)batch.Isucc;
591 :
592 :
593 : //std::cerr << "[" << blog << " -- " << loglikelyhoodRatio << " -- " << alog << "]" << std::endl;
594 0 : double l = 1- 1/(1+ exp(fabs(loglikelyhoodRatio)));
595 :
596 0 : if(loglikelyhoodRatio <= 0){ //Accept H0
597 0 : return ConfInt((double)batch.Isucc/(double)batch.I, Value ,1,0.0,1.0,l);
598 : }else{
599 0 : return ConfInt((double)batch.Isucc/(double)batch.I, 0 ,Value2,0.0,1.0,l);
600 : }
601 : }
602 :
603 : case CONSTANT:
604 15 : return ConfInt(Value,Value2,Value-Value2/2.0,Value+Value2/2.0,1.0);
605 :
606 : case HASL_PLUS:
607 5 : return (left->eval(batch) += right->eval(batch));
608 :
609 : case HASL_MINUS:
610 0 : return (left->eval(batch) -= right->eval(batch));
611 :
612 : case HASL_TIME:
613 5 : return (left->eval(batch) *= right->eval(batch));
614 :
615 : case HASL_DIV:
616 5 : return (left->eval(batch) /= right->eval(batch));
617 :
618 : case EXISTS:
619 0 : if ( batch.Isucc > 0) return ConfInt(1.0,1.0,1.0,0.0,1.0,1.0);
620 0 : else return ConfInt(0.0,0.0,1.0,0.0,0.0,1.0);
621 : case NOTALLS:
622 0 : if ( batch.Isucc < batch.I) return ConfInt(1.0,1.0,1.0,0.0,1.0,1.0);
623 0 : else return ConfInt(0.0,0.0,1.0,0.0,0.0,1.0);
624 :
625 :
626 : default:
627 0 : std::cerr << "Fail to parse Hasl Formula"<< std::endl;
628 0 : exit(EXIT_FAILURE);
629 : }
630 : }
631 :
632 0 : double HaslFormulasTop::bound()const{
633 0 : switch (TypeOp) {
634 : case PROBABILITY:
635 : case PROBCOND:
636 0 : return 0.5;
637 :
638 : case EXPECTANCY:
639 : case CDF_PART:
640 : case PDF_PART:
641 : case RE_AVG:
642 : case RE_Continuous:
643 0 : return std::numeric_limits<double>::infinity();
644 :
645 : case CONSTANT:
646 0 : return fabs(Value);
647 :
648 : case HASL_PLUS:
649 : {
650 0 : double v1 = left->bound();
651 0 : double v2 = right->bound();
652 0 : return (fmax(v1,v2));
653 : }
654 :
655 : case HASL_TIME:
656 : {
657 0 : double v1 = left->bound();
658 0 : double v2 = right->bound();
659 0 : return (v1*v2);
660 : }
661 :
662 : case HASL_DIV:
663 : {
664 0 : return (std::numeric_limits<double>::infinity());
665 : }
666 :
667 : default:
668 0 : std::cerr << "Fail to parse Hasl Formula"<< std::endl;
669 0 : exit(EXIT_FAILURE);
670 : }
671 105 : }
|
