Файл betaDF.h /***********************************************************/ /* Beta distribution */ /***********************************************************/ #ifndef __BETA_H__ /* To prevent redefinition */ #define ENTRY extern #define LOCAL static class BetaDF { public: BetaDF(double u, double w); double value(double x); // Функция распределения Beta(x|a,b) double inv(double p); // Обратная функция: Beta(x|a,b)=p private: double a,b, logBeta; double fraction(double a, double b, double x); }; #define __BETA_H__ /* Prevents redefinition */ #endif /* Ends #ifndef __BETA_H__ */ Файл betaDF.сpp /***********************************************************/ /* Бета-распределение */ /***********************************************************/ #include #include #include "BetaDF.h" #include "logGamma.h" BetaDF::BetaDF(double u, double w): a(u), b(w), logBeta(logGamma(a) + logGamma(b) - logGamma(a + b)) { assert(a > 0 && b > 0); } double BetaDF::fraction(double a, double b, double x) // // См. Abramowitz & Stegun, // Handbook of Mathematical Functions, 1964 [26.5.8] // М.Абрамовиц, И.Стиган // Справочник по специальным функциям (М: Мир, 1979) // // Неполная бета-функция вычисляется с помощью разложения в цепную дробь // // i_beta(a,b,x) = x^{a}*(1-x)^{b}*fraction / a * beta(a,b), // где // // 1 d1 d2 d3 d4 // fraction = --- ---- ---- ---- ---- .... // 1+ 1+ 1+ 1+ 1+ // // Подходящие дроби: A(n) / B(n) // // где // A(n) = (s(n) * A(n-1) + r(n) * A(n-2)) * factor // B(n) = (s(n) * B(n-1) + r(n) * B(n-2)) * factor // и // A(-1) = 1, B(-1) = 0, A(0) = s(0), B(0) = 1. // // Здесь s(0) = 0 и при n >= 1 s(n) = 1, // а r(1) = 1 и при i >= 2 // // r(i) = m(b-m)x / (a+i-1)(a+i) когда i = 2m, // r(i) = -(a+m)(a+b+m)x / (a+i-1)(a+i) когда i = 2m+1, // // factor - шкалирующий множитель, позволяющий избежать переполнения. // // Итак, A(0) = 0 , B(0) = 1, // r(1) = -(a+b)*x / (a+1) // A(1) = A(0) + r(1)*A(-1) = r(1) = 1 // B(1) = B(0) + r(1)*B(-1) = 1 // { double old_bta = 0, factor = 1; double A0 = 0, A1 = 1, B0 = 1, B1 = 1; double bta = 1, am = a, ai = a; double iter = 0, r; do { // часть цикла, вычисляющая нечетные подходящие дроби // начинаем с i = 1, iter = 0 ai += 1; // = a+i r = -am * (am + b) * x / ((ai - 1) * ai); /* пересчет A и B в два шага */ A0 = (A1 + r * A0) * factor; /* i НЕчетно */ B0 = (B1 + r * B0) * factor; // часть цикла, вычисляющая нечетные подходящие дроби // начинаем с i = 2, iter = 1 am += 1; iter += 1; ai += 1; r = iter * (b - iter) * x * factor / ((ai - 1) * ai); A1 = A0 + r * A1; /* i четно, A0 и B0 уже шкалированы */ B1 = B0 + r * B1; old_bta = bta; factor = 1 / B1; bta = A1 * factor; } while (fabs(old_bta) != fabs(bta)); return bta * exp(a * log(x) + b * log(1 - x) - logBeta) / a; }/*incBeta_fraction*/ double BetaDF::value(double x) // // Вычисляет Beta(x|a,b): // вероятность того, что случайная величина, // подчиняющаяся бета-распределению с параметрами 'a' и 'b', // меньше или равна 'x'. // { if (x <= 0) return 0; /* НЕ ошибка! */ else if (x >= 1) return 1; /* НЕ ошибка! */ if (x < (a + 1) / (a + b + 2)) return fraction(a, b, x); else return 1 - fraction(b, a, 1 - x); }/*value*/ double BetaDF::inv(double p) // // Ищет такое значение 'x', для которого Beta(x|a,b) = p, // т.е. равна 'p' вероятность того, что случайная величина, // подчиняющаяся бета-распределению с параметрами 'a' и 'b', // меньше или равна 'x'. // { double fx, l = 0, r = 1, x = 0.5; assert(p >= 0 && p <= 1); if (p == 0 || p == 1) return p; do { fx = value(x); if (fx > p) r = x; else if (fx < p) l = x; else return x; x = (l + r)* 0.5; } while ((l!=x) && (r!=x)); return x; }/*inv*/ <<<  Назад  1  2   8  Комментарии к статье  8 8  Обсудить в чате