source: sasview/src/sas/models/c_extension/cephes/kolmogorov.c @ 8f46df7

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 8f46df7 was 79492222, checked in by krzywon, 10 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 4.9 KB
Line 
1
2/* Re Kolmogorov statistics, here is Birnbaum and Tingey's formula for the
3   distribution of D+, the maximum of all positive deviations between a
4   theoretical distribution function P(x) and an empirical one Sn(x)
5   from n samples.
6
7     +
8    D  =         sup     [P(x) - S (x)]
9     n     -inf < x < inf         n
10
11
12                  [n(1-e)]
13        +            -                    v-1              n-v
14    Pr{D   > e} =    >    C    e (e + v/n)    (1 - e - v/n)
15        n            -   n v
16                    v=0
17
18    [n(1-e)] is the largest integer not exceeding n(1-e).
19    nCv is the number of combinations of n things taken v at a time.  */
20
21
22#include "mconf.h"
23#ifdef ANSIPROT
24extern double pow ( double, double );
25extern double floor ( double );
26extern double lgam ( double );
27extern double exp ( double );
28extern double sqrt ( double );
29extern double log ( double );
30extern double fabs ( double );
31double smirnov ( int, double );
32double kolmogorov ( double );
33#else
34double pow (), floor (), lgam (), exp (), sqrt (), log (), fabs ();
35double smirnov (), kolmogorov ();
36#endif
37extern double MAXLOG;
38
39/* Exact Smirnov statistic, for one-sided test.  */
40double
41smirnov (n, e)
42     int n;
43     double e;
44{
45  int v, nn;
46  double evn, omevn, p, t, c, lgamnp1;
47
48  if (n <= 0 || e < 0.0 || e > 1.0)
49    return (-1.0);
50  nn = floor ((double) n * (1.0 - e));
51  p = 0.0;
52  if (n < 1013)
53    {
54      c = 1.0;
55      for (v = 0; v <= nn; v++)
56        {
57          evn = e + ((double) v) / n;
58          p += c * pow (evn, (double) (v - 1))
59            * pow (1.0 - evn, (double) (n - v));
60          /* Next combinatorial term; worst case error = 4e-15.  */
61          c *= ((double) (n - v)) / (v + 1);
62        }
63    }
64  else
65    {
66      lgamnp1 = lgam ((double) (n + 1));
67      for (v = 0; v <= nn; v++)
68        {
69          evn = e + ((double) v) / n;
70          omevn = 1.0 - evn;
71          if (fabs (omevn) > 0.0)
72            {
73              t = lgamnp1
74                - lgam ((double) (v + 1))
75                - lgam ((double) (n - v + 1))
76                + (v - 1) * log (evn)
77                + (n - v) * log (omevn);
78              if (t > -MAXLOG)
79                p += exp (t);
80            }
81        }
82    }
83  return (p * e);
84}
85
86
87/* Kolmogorov's limiting distribution of two-sided test, returns
88   probability that sqrt(n) * max deviation > y,
89   or that max deviation > y/sqrt(n).
90   The approximation is useful for the tail of the distribution
91   when n is large.  */
92double
93kolmogorov (y)
94     double y;
95{
96  double p, t, r, sign, x;
97
98  x = -2.0 * y * y;
99  sign = 1.0;
100  p = 0.0;
101  r = 1.0;
102  do
103    {
104      t = exp (x * r * r);
105      p += sign * t;
106      if (t == 0.0)
107        break;
108      r += 1.0;
109      sign = -sign;
110    }
111  while ((t / p) > 1.1e-16);
112  return (p + p);
113}
114
115/* Functional inverse of Smirnov distribution
116   finds e such that smirnov(n,e) = p.  */
117double
118smirnovi (n, p)
119     int n;
120     double p;
121{
122  double e, t, dpde;
123
124  if (p <= 0.0 || p > 1.0)
125    {
126      mtherr ("smirnovi", DOMAIN);
127      return 0.0;
128    }
129  /* Start with approximation p = exp(-2 n e^2).  */
130  e = sqrt (-log (p) / (2.0 * n));
131  do
132    {
133      /* Use approximate derivative in Newton iteration. */
134      t = -2.0 * n * e;
135      dpde = 2.0 * t * exp (t * e);
136      if (fabs (dpde) > 0.0)
137        t = (p - smirnov (n, e)) / dpde;
138      else
139        {
140          mtherr ("smirnovi", UNDERFLOW);
141          return 0.0;
142        }
143      e = e + t;
144      if (e >= 1.0 || e <= 0.0)
145        {
146          mtherr ("smirnovi", OVERFLOW);
147          return 0.0;
148        }
149    }
150  while (fabs (t / e) > 1e-10);
151  return (e);
152}
153
154
155/* Functional inverse of Kolmogorov statistic for two-sided test.
156   Finds y such that kolmogorov(y) = p.
157   If e = smirnovi (n,p), then kolmogi(2 * p) / sqrt(n) should
158   be close to e.  */
159double
160kolmogi (p)
161     double p;
162{
163  double y, t, dpdy;
164
165  if (p <= 0.0 || p > 1.0)
166    {
167      mtherr ("kolmogi", DOMAIN);
168      return 0.0;
169    }
170  /* Start with approximation p = 2 exp(-2 y^2).  */
171  y = sqrt (-0.5 * log (0.5 * p));
172  do
173    {
174      /* Use approximate derivative in Newton iteration. */
175      t = -2.0 * y;
176      dpdy = 4.0 * t * exp (t * y);
177      if (fabs (dpdy) > 0.0)
178        t = (p - kolmogorov (y)) / dpdy;
179      else
180        {
181          mtherr ("kolmogi", UNDERFLOW);
182          return 0.0;
183        }
184      y = y + t;
185    }
186  while (fabs (t / y) > 1e-10);
187  return (y);
188}
189
190
191#ifdef SALONE
192/* Type in a number.  */
193void
194getnum (s, px)
195     char *s;
196     double *px;
197{
198  char str[30];
199
200  printf (" %s (%.15e) ? ", s, *px);
201  gets (str);
202  if (str[0] == '\0' || str[0] == '\n')
203    return;
204  sscanf (str, "%lf", px);
205  printf ("%.15e\n", *px);
206}
207
208/* Type in values, get answers.  */
209void
210main ()
211{
212  int n;
213  double e, p, ps, pk, ek, y;
214
215  n = 5;
216  e = 0.0;
217  p = 0.1;
218loop:
219  ps = n;
220  getnum ("n", &ps);
221  n = ps;
222  if (n <= 0)
223    {
224      printf ("? Operator error.\n");
225      goto loop;
226    }
227  /*
228  getnum ("e", &e);
229  ps = smirnov (n, e);
230  y = sqrt ((double) n) * e;
231  printf ("y = %.4e\n", y);
232  pk = kolmogorov (y);
233  printf ("Smirnov = %.15e, Kolmogorov/2 = %.15e\n", ps, pk / 2.0);
234*/
235  getnum ("p", &p);
236  e = smirnovi (n, p);
237  printf ("Smirnov e = %.15e\n", e);
238  y = kolmogi (2.0 * p);
239  ek = y / sqrt ((double) n);
240  printf ("Kolmogorov e = %.15e\n", ek);
241  goto loop;
242}
243#endif
Note: See TracBrowser for help on using the repository browser.