ACT-CV - Machine Vision for Cognitive Modeling
distrib.h
Go to the documentation of this file.
1 // -*- mode: c++; indent-tabs-mode: nil; c-basic-offset: 4; coding: iso-8859-1; -*-
2 
3 /*
4 ACT-CV - Machine Vision for Cognitive Modeling
5 Copyright (c) 2008 Marc Halbruegge
6 
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11 
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21 
29 #ifndef _DISTRIB_H
30 #define _DISTRIB_H
31 
32 #include <cmath>
33 #include <stdexcept>
34 
35 
36 #define ITERMAX 100
37 #define STARTFRACTION 20
38 
39 // 1/sqrt(2*PI)
40 #define SQRTI2PI 0.398942280401432677939946059934382
41 
42 // log(2*PI)/2
43 #define HL2PI 0.918938533204672741780329736405618
44 
45 #define MIN_WAHRSCHEINLICHKEIT 1e-10
46 
48 class Logit {
50  Logit(){}
51 public:
52 
53  static double expit(double l) {
54  return std::exp(l)/(1.0+std::exp(l));
55  }
56 
57  static double logit(double p) {
58  return std::log(p/(1.0-p));
59  }
60  static double risiko(double p) {
61  return p/(1.0-p);
62  }
63  static double inv_risiko(double r) {
64  return r/(1.0+r);
65  }
66 
67 };
68 
69 
71 class Chi {
73  Chi(){};
74 protected:
75  static double log_gamma (double x) {
76  double a,g;
77  int i,n;
78 
79  a=x;
80  n = (int) std::floor(6.0 - a + 0.0001);
81  if (n>0) a = a + n;
82 
83  g = 1.0/(a*a);
84  g = (1.0 - g*(1.0/30.0 - g*(1.0/105.0 - g*(1.0/140.0 - g/99.0))))/(12.0*a);
85  g = g + ( (a - 0.5)*log(a) - a + HL2PI );
86  for (i=0; i<n; i++) {
87  a = a - 1.0;
88  g = g - log(a);
89  }
90  return g;
91  }
92 
93  static double dgamma(double x1, double x2) {
94  // Schnelle und genaue Berechnung der
95  // Gamma-Verteilung mit zwei Parametern
96  double alpha,y, CF, t, gamma;
97  int k=STARTFRACTION;
98 
99  alpha=x1;
100  y=x2;
101 
102  if(y<=alpha) {
103  t = y*(1.0-k-alpha)/(alpha+2.0*k-1.0+k*y/(alpha+2*k));
104  for (k=STARTFRACTION-1; k>=2; k--)
105  t = y*(1.0-k-alpha)/(alpha+2.0*k-1.0+k*y/(alpha+2*k+t));
106  CF = 1.0 - y/(alpha+1.0+y/(alpha+2.0+t));
107  gamma = std::exp(alpha*log(y)-y-log_gamma(alpha+1.0)-log(CF));
108  } else {
109  t = (k-alpha)/(y+k);
110  for (k=STARTFRACTION-1; k>=2; k--)
111  t = (k-alpha)/(y+k/(1.0+t));
112  CF = 1.0 + (1-alpha)/(y+1.0/(1.0+t));
113  gamma = (alpha-1.0)*std::log(y)-y-log_gamma(alpha)-log(CF);
114  if(gamma>0) gamma=0;
115  if(gamma<-100) gamma=1.0;
116  else gamma = 1.0 - std::exp(gamma);
117  }
118  return gamma;
119  }
120 
121 public:
122 
123  static double ChiQVerteilung(double chi, double df) {
124  //Schwanz der Verteilung, d.h. von "rechts" gelesen
125  double p;
126  if (df<=0) throw std::range_error("df too low");
127  if (chi<0.00001) return 1.0-MIN_WAHRSCHEINLICHKEIT;
128  if (df>4 && chi>10*df) return MIN_WAHRSCHEINLICHKEIT;
129  p=dgamma(df/2.0, chi/2.0);
130  p=1-p;
131  return p > MIN_WAHRSCHEINLICHKEIT ? p : MIN_WAHRSCHEINLICHKEIT;
132  }
133 
134  static double ChiSqTail(double chi, double df) {
135  //Schwanz der Verteilung, d.h. von "rechts" gelesen
136  double p;
137  if (df<=0) throw std::range_error("df too low");
138  if (chi<0.00001) return 1.0-MIN_WAHRSCHEINLICHKEIT;
139  if (df>4 && chi>10*df) return MIN_WAHRSCHEINLICHKEIT;
140  p=dgamma(df/2.0, chi/2.0);
141  p=1-p;
142  return p > MIN_WAHRSCHEINLICHKEIT ? p : MIN_WAHRSCHEINLICHKEIT;
143  }
144 
145  static double ChiSqHead(double chi, double df) {
146  double p;
147  if (df<=0) throw std::range_error("df too low");
148  if (chi<0.00001) return MIN_WAHRSCHEINLICHKEIT;
149  if (df>4 && chi>10*df) return 1.0-MIN_WAHRSCHEINLICHKEIT;
150  p=dgamma(df/2.0, chi/2.0);
151  return p > MIN_WAHRSCHEINLICHKEIT ? p : MIN_WAHRSCHEINLICHKEIT;
152  }
153 
154 };
155 
156 
158 class Phi {
159 public:
160 
161  static double NormalVerteilung(double z) {
162  double t,t1,a,z1,w,p1,zz,p;
163  int i,weiter;
164 
165  a=z;
166  if (a<0) a=-a;
167 
168  if (a>6.8) {
169  zz = a*a;
170  w = 1.0-1.0/(zz+3.0-1.0/(0.22*zz+0.704));
171  p = 0.5-SQRTI2PI*exp(-zz/2.0)*w/a;
172  } else {
173  i = 2;
174  t = 1.0;
175  p = 1.0;
176  t1 = a*a/4.0;
177  z1 = t1;
178  p1 = 0.0;
179  weiter = 1;
180  while (weiter) {
181  t = z1*(t1-t)/i;
182  i++;
183  p += (t/i);
184  if (std::abs(p1-p)<MIN_WAHRSCHEINLICHKEIT) weiter=0;
185  if (weiter) {
186  t1 = z1*(t-t1)/i;
187  i++;
188  p1 = p;
189  }
190  }
191  p = p * SQRTI2PI * a * std::exp(-z1/2.0);
192  }
193  if (z < 0.0) return (0.5 + p);
194  else return (0.5 - p);
195  }
196 
197 
198  static double InvNormalVerteilung(double x) {
199  double p, unten, oben, mitte, pu, po, pm;
200  int iter;
201 
202  p=x;
203  if (p<0.005) p = 0.005;
204  if (p>0.995) p = 0.995;
205 
206  if (p<0.5) {
207  oben = 0.0;
208  unten = -2.0;
209  } else {
210  unten = 0.0;
211  oben =2.0;
212  }
213 
214  pu=1-NormalVerteilung(unten);
215  po=1-NormalVerteilung(oben);
216 
217  for (iter=0; iter<ITERMAX; iter++) {
218  mitte = (unten+oben)/2.0;
219  if (pu>p) {
220  oben= unten;
221  unten *= 2;
222  pu=1-NormalVerteilung(unten);
223  } else if (po<p) {
224  unten = oben;
225  oben *= 2;
226  po=1-NormalVerteilung(oben);
227  } else {
228  pm=1-NormalVerteilung(mitte);
229  if (pm>=p) {
230  oben=mitte;
231  po = pm;
232  } else {
233  unten=mitte;
234  pu = pm;
235  }
236  }
237 
238  }
239 
240  return mitte;
241 
242  }
243 
244 };
245 
247 class FisherZ {
248 public:
249  static double Z(double r) {
250  return .5 * std::log( (1.0+r)/(1.0-r) );
251  }
252  static double InvZ(double Z) {
253  return (std::exp(2*Z)-1) / (std::exp(2*Z)+1);
254  }
255 };
256 
257 //abgeschrieben bei G. Gediga
258 
259 // References:
260 //
261 // F-Verteilung:
262 // Emerson, P.L., FISHTAIL: Practical F integration on small
263 // computers in BASIC and C. Behavioral Research Methods, Instru-
264 // mentations, & Computers, 1988, 20, 65-69.
265 //
266 // Chi-Quadrat-Verteilung / Gamma-Verteilung:
267 // Posten, H.O., An effective algorithm for the noncentral chi-
268 // squared distribution function. American Statistician, 1989,
269 // 43, 261-263.
270 //
271 // Normalverteilung:
272 // Kerridge, D.F. & Cook, G.W., Yet another series for the normal
273 // integral. Biometrika, 1976, 63, 401-403.
274 //
275 // Nicht-zentrale Chi-Quadrat- und F-Verteilung
276 // Johnson,N.L. & Kotz,S. (1970) Distributions in statistics.
277 // Continous univariate distributions - 2. New York: Wiley,
278 // pp 139 (Formel 22), pp 195 (3-Momenten-Approx).
279 //
280 
281 #endif
282 

ACT-CV - Machine Vision for Cognitive Modeling
© 2015 Marc Halbruegge (actcvlibrary@googlemail.com)