2 * mathfunc.c: Mathematical functions.
5 * Ross Ihaka. (See note 1.)
6 * The R Development Core Team. (See note 1.)
7 * Morten Welinder <terra@gnome.org>
8 * Miguel de Icaza (miguel@gnu.org)
9 * Jukka-Pekka Iivonen (iivonen@iki.fi)
10 * Ian Smith (iandjmsmith@aol.com). (See note 2.)
14 * NOTE 1: most of this file comes from the "R" package, notably version 2
15 * or newer (we re-sync from time to time).
16 * "R" is distributed under GPL licence, see file COPYING.
17 * The relevant parts are copyright (C) 1998 Ross Ihaka and
18 * 2000-2004 The R Development Core Team.
24 * NOTE 2: the pbeta (and support) code comes from Ian Smith. (Translated
25 * into C, adapted to Gnumeric naming convensions, and R's API conventions
26 * by Morten Welinder. Blame me for problems.)
28 * Copyright © Ian Smith 2002-2003
30 * Thanks to Jerry W. Lewis for help with testing of and improvements to the code.
36 #include <gnumeric-config.h>
41 #include <glib/gi18n-lib.h>
48 #include <goffice/goffice.h>
49 #include <glib/gstdio.h>
53 /* R code wants this, so provide it. */
58 #define M_SQRT_32 GNM_const(5.656854249492380195206754896838) /* sqrt(32) */
59 #define M_1_SQRT_2PI GNM_const(0.398942280401432677939946059934) /* 1/sqrt(2pi) */
60 #define M_2PIgnum (2 * M_PIgnum)
62 #define ML_ERROR(cause) do { } while(0)
63 #define MATHLIB_WARNING g_warning
64 #define REprintf g_warning
66 static inline gnm_float
fmin2 (gnm_float x
, gnm_float y
) { return MIN (x
, y
); }
67 static inline gnm_float
fmax2 (gnm_float x
, gnm_float y
) { return MAX (x
, y
); }
69 #define MATHLIB_STANDALONE
70 #define ML_ERR_return_NAN { return gnm_nan; }
71 static void pnorm_both (gnm_float x
, gnm_float
*cum
, gnm_float
*ccum
, int i_tail
, gboolean log_p
);
73 /* MW ---------------------------------------------------------------------- */
78 /* Nothing, for the time being. */
82 * A table of logs for scaling purposes. Each value has four parts with
83 * 23 bits in each. That means each part can be multiplied by a double
84 * with at most 30 bits set and not have any rounding error. Note, that
85 * the first entry is log(2).
87 * Entry i is associated with the value r = 0.5 + i / 256.0. The
88 * argument to log is p/q where q=1024 and p=floor(q / r + 0.5).
89 * Thus r*p/q is close to 1.
91 static const float bd0_scale
[128 + 1][4] = {
92 { +0x1.62e430p
-1, -0x1.05c610p
-29, -0x1.950d88p
-54, +0x1.d9cc02p
-79 }, /* 128: log(2048/1024.) */
93 { +0x1.5ee02cp
-1, -0x1.6dbe98p
-25, -0x1.51e540p
-50, +0x1.2bfa48p
-74 }, /* 129: log(2032/1024.) */
94 { +0x1.5ad404p
-1, +0x1.86b3e4p
-26, +0x1.9f6534p
-50, +0x1.54be04p
-74 }, /* 130: log(2016/1024.) */
95 { +0x1.570124p
-1, -0x1.9ed750p
-25, -0x1.f37dd0p
-51, +0x1.10b770p
-77 }, /* 131: log(2001/1024.) */
96 { +0x1.5326e4p
-1, -0x1.9b9874p
-25, -0x1.378194p
-49, +0x1.56feb2p
-74 }, /* 132: log(1986/1024.) */
97 { +0x1.4f4528p
-1, +0x1.aca70cp
-28, +0x1.103e74p
-53, +0x1.9c410ap
-81 }, /* 133: log(1971/1024.) */
98 { +0x1.4b5bd8p
-1, -0x1.6a91d8p
-25, -0x1.8e43d0p
-50, -0x1.afba9ep
-77 }, /* 134: log(1956/1024.) */
99 { +0x1.47ae54p
-1, -0x1.abb51cp
-25, +0x1.19b798p
-51, +0x1.45e09cp
-76 }, /* 135: log(1942/1024.) */
100 { +0x1.43fa00p
-1, -0x1.d06318p
-25, -0x1.8858d8p
-49, -0x1.1927c4p
-75 }, /* 136: log(1928/1024.) */
101 { +0x1.3ffa40p
-1, +0x1.1a427cp
-25, +0x1.151640p
-53, -0x1.4f5606p
-77 }, /* 137: log(1913/1024.) */
102 { +0x1.3c7c80p
-1, -0x1.19bf48p
-34, +0x1.05fc94p
-58, -0x1.c096fcp
-82 }, /* 138: log(1900/1024.) */
103 { +0x1.38b320p
-1, +0x1.6b5778p
-25, +0x1.be38d0p
-50, -0x1.075e96p
-74 }, /* 139: log(1886/1024.) */
104 { +0x1.34e288p
-1, +0x1.d9ce1cp
-25, +0x1.316eb8p
-49, +0x1.2d885cp
-73 }, /* 140: log(1872/1024.) */
105 { +0x1.315124p
-1, +0x1.c2fc60p
-29, -0x1.4396fcp
-53, +0x1.acf376p
-78 }, /* 141: log(1859/1024.) */
106 { +0x1.2db954p
-1, +0x1.720de4p
-25, -0x1.d39b04p
-49, -0x1.f11176p
-76 }, /* 142: log(1846/1024.) */
107 { +0x1.2a1b08p
-1, -0x1.562494p
-25, +0x1.a7863cp
-49, +0x1.85dd64p
-73 }, /* 143: log(1833/1024.) */
108 { +0x1.267620p
-1, +0x1.3430e0p
-29, -0x1.96a958p
-56, +0x1.f8e636p
-82 }, /* 144: log(1820/1024.) */
109 { +0x1.23130cp
-1, +0x1.7bebf4p
-25, +0x1.416f1cp
-52, -0x1.78dd36p
-77 }, /* 145: log(1808/1024.) */
110 { +0x1.1faa34p
-1, +0x1.70e128p
-26, +0x1.81817cp
-50, -0x1.c2179cp
-76 }, /* 146: log(1796/1024.) */
111 { +0x1.1bf204p
-1, +0x1.3a9620p
-28, +0x1.2f94c0p
-52, +0x1.9096c0p
-76 }, /* 147: log(1783/1024.) */
112 { +0x1.187ce4p
-1, -0x1.077870p
-27, +0x1.655a80p
-51, +0x1.eaafd6p
-78 }, /* 148: log(1771/1024.) */
113 { +0x1.1501c0p
-1, -0x1.406cacp
-25, -0x1.e72290p
-49, +0x1.5dd800p
-73 }, /* 149: log(1759/1024.) */
114 { +0x1.11cb80p
-1, +0x1.787cd0p
-25, -0x1.efdc78p
-51, -0x1.5380cep
-77 }, /* 150: log(1748/1024.) */
115 { +0x1.0e4498p
-1, +0x1.747324p
-27, -0x1.024548p
-51, +0x1.77a5a6p
-75 }, /* 151: log(1736/1024.) */
116 { +0x1.0b036cp
-1, +0x1.690c74p
-25, +0x1.5d0cc4p
-50, -0x1.c0e23cp
-76 }, /* 152: log(1725/1024.) */
117 { +0x1.077070p
-1, -0x1.a769bcp
-27, +0x1.452234p
-52, +0x1.6ba668p
-76 }, /* 153: log(1713/1024.) */
118 { +0x1.04240cp
-1, -0x1.a686acp
-27, -0x1.ef46b0p
-52, -0x1.5ce10cp
-76 }, /* 154: log(1702/1024.) */
119 { +0x1.00d22cp
-1, +0x1.fc0e10p
-25, +0x1.6ee034p
-50, -0x1.19a2ccp
-74 }, /* 155: log(1691/1024.) */
120 { +0x1.faf588p
-2, +0x1.ef1e64p
-27, -0x1.26504cp
-54, -0x1.b15792p
-82 }, /* 156: log(1680/1024.) */
121 { +0x1.f4d87cp
-2, +0x1.d7b980p
-26, -0x1.a114d8p
-50, +0x1.9758c6p
-75 }, /* 157: log(1670/1024.) */
122 { +0x1.ee1414p
-2, +0x1.2ec060p
-26, +0x1.dc00fcp
-52, +0x1.f8833cp
-76 }, /* 158: log(1659/1024.) */
123 { +0x1.e7e32cp
-2, -0x1.ac796cp
-27, -0x1.a68818p
-54, +0x1.235d02p
-78 }, /* 159: log(1649/1024.) */
124 { +0x1.e108a0p
-2, -0x1.768ba4p
-28, -0x1.f050a8p
-52, +0x1.00d632p
-82 }, /* 160: log(1638/1024.) */
125 { +0x1.dac354p
-2, -0x1.d3a6acp
-30, +0x1.18734cp
-57, -0x1.f97902p
-83 }, /* 161: log(1628/1024.) */
126 { +0x1.d47424p
-2, +0x1.7dbbacp
-31, -0x1.d5ada4p
-56, +0x1.56fcaap
-81 }, /* 162: log(1618/1024.) */
127 { +0x1.ce1af0p
-2, +0x1.70be7cp
-27, +0x1.6f6fa4p
-51, +0x1.7955a2p
-75 }, /* 163: log(1608/1024.) */
128 { +0x1.c7b798p
-2, +0x1.ec36ecp
-26, -0x1.07e294p
-50, -0x1.ca183cp
-75 }, /* 164: log(1598/1024.) */
129 { +0x1.c1ef04p
-2, +0x1.c1dfd4p
-26, +0x1.888eecp
-50, -0x1.fd6b86p
-75 }, /* 165: log(1589/1024.) */
130 { +0x1.bb7810p
-2, +0x1.478bfcp
-26, +0x1.245b8cp
-50, +0x1.ea9d52p
-74 }, /* 166: log(1579/1024.) */
131 { +0x1.b59da0p
-2, -0x1.882b08p
-27, +0x1.31573cp
-53, -0x1.8c249ap
-77 }, /* 167: log(1570/1024.) */
132 { +0x1.af1294p
-2, -0x1.b710f4p
-27, +0x1.622670p
-51, +0x1.128578p
-76 }, /* 168: log(1560/1024.) */
133 { +0x1.a925d4p
-2, -0x1.0ae750p
-27, +0x1.574ed4p
-51, +0x1.084996p
-75 }, /* 169: log(1551/1024.) */
134 { +0x1.a33040p
-2, +0x1.027d30p
-29, +0x1.b9a550p
-53, -0x1.b2e38ap
-78 }, /* 170: log(1542/1024.) */
135 { +0x1.9d31c0p
-2, -0x1.5ec12cp
-26, -0x1.5245e0p
-52, +0x1.2522d0p
-79 }, /* 171: log(1533/1024.) */
136 { +0x1.972a34p
-2, +0x1.135158p
-30, +0x1.a5c09cp
-56, +0x1.24b70ep
-80 }, /* 172: log(1524/1024.) */
137 { +0x1.911984p
-2, +0x1.0995d4p
-26, +0x1.3bfb5cp
-50, +0x1.2c9dd6p
-75 }, /* 173: log(1515/1024.) */
138 { +0x1.8bad98p
-2, -0x1.1d6144p
-29, +0x1.5b9208p
-53, +0x1.1ec158p
-77 }, /* 174: log(1507/1024.) */
139 { +0x1.858b58p
-2, -0x1.1b4678p
-27, +0x1.56cab4p
-53, -0x1.2fdc0cp
-78 }, /* 175: log(1498/1024.) */
140 { +0x1.7f5fa0p
-2, +0x1.3aaf48p
-27, +0x1.461964p
-51, +0x1.4ae476p
-75 }, /* 176: log(1489/1024.) */
141 { +0x1.79db68p
-2, -0x1.7e5054p
-26, +0x1.673750p
-51, -0x1.a11f7ap
-76 }, /* 177: log(1481/1024.) */
142 { +0x1.744f88p
-2, -0x1.cc0e18p
-26, -0x1.1e9d18p
-50, -0x1.6c06bcp
-78 }, /* 178: log(1473/1024.) */
143 { +0x1.6e08ecp
-2, -0x1.5d45e0p
-26, -0x1.c73ec8p
-50, +0x1.318d72p
-74 }, /* 179: log(1464/1024.) */
144 { +0x1.686c80p
-2, +0x1.e9b14cp
-26, -0x1.13bbd4p
-50, -0x1.efeb1cp
-78 }, /* 180: log(1456/1024.) */
145 { +0x1.62c830p
-2, -0x1.a8c70cp
-27, -0x1.5a1214p
-51, -0x1.bab3fcp
-79 }, /* 181: log(1448/1024.) */
146 { +0x1.5d1bdcp
-2, -0x1.4fec6cp
-31, +0x1.423638p
-56, +0x1.ee3feep
-83 }, /* 182: log(1440/1024.) */
147 { +0x1.576770p
-2, +0x1.7455a8p
-26, -0x1.3ab654p
-50, -0x1.26be4cp
-75 }, /* 183: log(1432/1024.) */
148 { +0x1.5262e0p
-2, -0x1.146778p
-26, -0x1.b9f708p
-52, -0x1.294018p
-77 }, /* 184: log(1425/1024.) */
149 { +0x1.4c9f08p
-2, +0x1.e152c4p
-26, -0x1.dde710p
-53, +0x1.fd2208p
-77 }, /* 185: log(1417/1024.) */
150 { +0x1.46d2d8p
-2, +0x1.c28058p
-26, -0x1.936284p
-50, +0x1.9fdd68p
-74 }, /* 186: log(1409/1024.) */
151 { +0x1.41b940p
-2, +0x1.cce0c0p
-26, -0x1.1a4050p
-50, +0x1.bc0376p
-76 }, /* 187: log(1402/1024.) */
152 { +0x1.3bdd24p
-2, +0x1.d6296cp
-27, +0x1.425b48p
-51, -0x1.cddb2cp
-77 }, /* 188: log(1394/1024.) */
153 { +0x1.36b578p
-2, -0x1.287ddcp
-27, -0x1.2d0f4cp
-51, +0x1.38447ep
-75 }, /* 189: log(1387/1024.) */
154 { +0x1.31871cp
-2, +0x1.2a8830p
-27, +0x1.3eae54p
-52, -0x1.898136p
-77 }, /* 190: log(1380/1024.) */
155 { +0x1.2b9304p
-2, -0x1.51d8b8p
-28, +0x1.27694cp
-52, -0x1.fd852ap
-76 }, /* 191: log(1372/1024.) */
156 { +0x1.265620p
-2, -0x1.d98f3cp
-27, +0x1.a44338p
-51, -0x1.56e85ep
-78 }, /* 192: log(1365/1024.) */
157 { +0x1.211254p
-2, +0x1.986160p
-26, +0x1.73c5d0p
-51, +0x1.4a861ep
-75 }, /* 193: log(1358/1024.) */
158 { +0x1.1bc794p
-2, +0x1.fa3918p
-27, +0x1.879c5cp
-51, +0x1.16107cp
-78 }, /* 194: log(1351/1024.) */
159 { +0x1.1675ccp
-2, -0x1.4545a0p
-26, +0x1.c07398p
-51, +0x1.f55c42p
-76 }, /* 195: log(1344/1024.) */
160 { +0x1.111ce4p
-2, +0x1.f72670p
-37, -0x1.b84b5cp
-61, +0x1.a4a4dcp
-85 }, /* 196: log(1337/1024.) */
161 { +0x1.0c81d4p
-2, +0x1.0c150cp
-27, +0x1.218600p
-51, -0x1.d17312p
-76 }, /* 197: log(1331/1024.) */
162 { +0x1.071b84p
-2, +0x1.fcd590p
-26, +0x1.a3a2e0p
-51, +0x1.fe5ef8p
-76 }, /* 198: log(1324/1024.) */
163 { +0x1.01ade4p
-2, -0x1.bb1844p
-28, +0x1.db3cccp
-52, +0x1.1f56fcp
-77 }, /* 199: log(1317/1024.) */
164 { +0x1.fa01c4p
-3, -0x1.12a0d0p
-29, -0x1.f71fb0p
-54, +0x1.e287a4p
-78 }, /* 200: log(1311/1024.) */
165 { +0x1.ef0adcp
-3, +0x1.7b8b28p
-28, -0x1.35bce4p
-52, -0x1.abc8f8p
-79 }, /* 201: log(1304/1024.) */
166 { +0x1.e598ecp
-3, +0x1.5a87e4p
-27, -0x1.134bd0p
-51, +0x1.c2cebep
-76 }, /* 202: log(1298/1024.) */
167 { +0x1.da85d8p
-3, -0x1.df31b0p
-27, +0x1.94c16cp
-57, +0x1.8fd7eap
-82 }, /* 203: log(1291/1024.) */
168 { +0x1.d0fb80p
-3, -0x1.bb5434p
-28, -0x1.ea5640p
-52, -0x1.8ceca4p
-77 }, /* 204: log(1285/1024.) */
169 { +0x1.c765b8p
-3, +0x1.e4d68cp
-27, +0x1.5b59b4p
-51, +0x1.76f6c4p
-76 }, /* 205: log(1279/1024.) */
170 { +0x1.bdc46cp
-3, -0x1.1cbb50p
-27, +0x1.2da010p
-51, +0x1.eb282cp
-75 }, /* 206: log(1273/1024.) */
171 { +0x1.b27980p
-3, -0x1.1b9ce0p
-27, +0x1.7756f8p
-52, +0x1.2ff572p
-76 }, /* 207: log(1266/1024.) */
172 { +0x1.a8bed0p
-3, -0x1.bbe874p
-30, +0x1.85cf20p
-56, +0x1.b9cf18p
-80 }, /* 208: log(1260/1024.) */
173 { +0x1.9ef83cp
-3, +0x1.2769a4p
-27, -0x1.85bda0p
-52, +0x1.8c8018p
-79 }, /* 209: log(1254/1024.) */
174 { +0x1.9525a8p
-3, +0x1.cf456cp
-27, -0x1.7137d8p
-52, -0x1.f158e8p
-76 }, /* 210: log(1248/1024.) */
175 { +0x1.8b46f8p
-3, +0x1.11b12cp
-30, +0x1.9f2104p
-54, -0x1.22836ep
-78 }, /* 211: log(1242/1024.) */
176 { +0x1.83040cp
-3, +0x1.2379e4p
-28, +0x1.b71c70p
-52, -0x1.990cdep
-76 }, /* 212: log(1237/1024.) */
177 { +0x1.790ed4p
-3, +0x1.dc4c68p
-28, -0x1.910ac8p
-52, +0x1.dd1bd6p
-76 }, /* 213: log(1231/1024.) */
178 { +0x1.6f0d28p
-3, +0x1.5cad68p
-28, +0x1.737c94p
-52, -0x1.9184bap
-77 }, /* 214: log(1225/1024.) */
179 { +0x1.64fee8p
-3, +0x1.04bf88p
-28, +0x1.6fca28p
-52, +0x1.8884a8p
-76 }, /* 215: log(1219/1024.) */
180 { +0x1.5c9400p
-3, +0x1.d65cb0p
-29, -0x1.b2919cp
-53, +0x1.b99bcep
-77 }, /* 216: log(1214/1024.) */
181 { +0x1.526e60p
-3, -0x1.c5e4bcp
-27, -0x1.0ba380p
-52, +0x1.d6e3ccp
-79 }, /* 217: log(1208/1024.) */
182 { +0x1.483bccp
-3, +0x1.9cdc7cp
-28, -0x1.5ad8dcp
-54, -0x1.392d3cp
-83 }, /* 218: log(1202/1024.) */
183 { +0x1.3fb25cp
-3, -0x1.a6ad74p
-27, +0x1.5be6b4p
-52, -0x1.4e0114p
-77 }, /* 219: log(1197/1024.) */
184 { +0x1.371fc4p
-3, -0x1.fe1708p
-27, -0x1.78864cp
-52, -0x1.27543ap
-76 }, /* 220: log(1192/1024.) */
185 { +0x1.2cca10p
-3, -0x1.4141b4p
-28, -0x1.ef191cp
-52, +0x1.00ee08p
-76 }, /* 221: log(1186/1024.) */
186 { +0x1.242310p
-3, +0x1.3ba510p
-27, -0x1.d003c8p
-51, +0x1.162640p
-76 }, /* 222: log(1181/1024.) */
187 { +0x1.1b72acp
-3, +0x1.52f67cp
-27, -0x1.fd6fa0p
-51, +0x1.1a3966p
-77 }, /* 223: log(1176/1024.) */
188 { +0x1.10f8e4p
-3, +0x1.129cd8p
-30, +0x1.31ef30p
-55, +0x1.a73e38p
-79 }, /* 224: log(1170/1024.) */
189 { +0x1.08338cp
-3, -0x1.005d7cp
-27, -0x1.661a9cp
-51, +0x1.1f138ap
-79 }, /* 225: log(1165/1024.) */
190 { +0x1.fec914p
-4, -0x1.c482a8p
-29, -0x1.55746cp
-54, +0x1.99f932p
-80 }, /* 226: log(1160/1024.) */
191 { +0x1.ed1794p
-4, +0x1.d06f00p
-29, +0x1.75e45cp
-53, -0x1.d0483ep
-78 }, /* 227: log(1155/1024.) */
192 { +0x1.db5270p
-4, +0x1.87d928p
-32, -0x1.0f52a4p
-57, +0x1.81f4a6p
-84 }, /* 228: log(1150/1024.) */
193 { +0x1.c97978p
-4, +0x1.af1d24p
-29, -0x1.0977d0p
-60, -0x1.8839d0p
-84 }, /* 229: log(1145/1024.) */
194 { +0x1.b78c84p
-4, -0x1.44f124p
-28, -0x1.ef7bc4p
-52, +0x1.9e0650p
-78 }, /* 230: log(1140/1024.) */
195 { +0x1.a58b60p
-4, +0x1.856464p
-29, +0x1.c651d0p
-55, +0x1.b06b0cp
-79 }, /* 231: log(1135/1024.) */
196 { +0x1.9375e4p
-4, +0x1.5595ecp
-28, +0x1.dc3738p
-52, +0x1.86c89ap
-81 }, /* 232: log(1130/1024.) */
197 { +0x1.814be4p
-4, -0x1.c073fcp
-28, -0x1.371f88p
-53, -0x1.5f4080p
-77 }, /* 233: log(1125/1024.) */
198 { +0x1.6f0d28p
-4, +0x1.5cad68p
-29, +0x1.737c94p
-53, -0x1.9184bap
-78 }, /* 234: log(1120/1024.) */
199 { +0x1.60658cp
-4, -0x1.6c8af4p
-28, +0x1.d8ef74p
-55, +0x1.c4f792p
-80 }, /* 235: log(1116/1024.) */
200 { +0x1.4e0110p
-4, +0x1.146b5cp
-29, +0x1.73f7ccp
-54, -0x1.d28db8p
-79 }, /* 236: log(1111/1024.) */
201 { +0x1.3b8758p
-4, +0x1.8b1b70p
-28, -0x1.20aca4p
-52, -0x1.651894p
-76 }, /* 237: log(1106/1024.) */
202 { +0x1.28f834p
-4, +0x1.43b6a4p
-30, -0x1.452af8p
-55, +0x1.976892p
-80 }, /* 238: log(1101/1024.) */
203 { +0x1.1a0fbcp
-4, -0x1.e4075cp
-28, +0x1.1fe618p
-52, +0x1.9d6dc2p
-77 }, /* 239: log(1097/1024.) */
204 { +0x1.075984p
-4, -0x1.4ce370p
-29, -0x1.d9fc98p
-53, +0x1.4ccf12p
-77 }, /* 240: log(1092/1024.) */
205 { +0x1.f0a30cp
-5, +0x1.162a68p
-37, -0x1.e83368p
-61, -0x1.d222a6p
-86 }, /* 241: log(1088/1024.) */
206 { +0x1.cae730p
-5, -0x1.1a8f7cp
-31, -0x1.5f9014p
-55, +0x1.2720c0p
-79 }, /* 242: log(1083/1024.) */
207 { +0x1.ac9724p
-5, -0x1.e8ee08p
-29, +0x1.a7de04p
-54, -0x1.9bba74p
-78 }, /* 243: log(1079/1024.) */
208 { +0x1.868a84p
-5, -0x1.ef8128p
-30, +0x1.dc5eccp
-54, -0x1.58d250p
-79 }, /* 244: log(1074/1024.) */
209 { +0x1.67f950p
-5, -0x1.ed684cp
-30, -0x1.f060c0p
-55, -0x1.b1294cp
-80 }, /* 245: log(1070/1024.) */
210 { +0x1.494accp
-5, +0x1.a6c890p
-32, -0x1.c3ad48p
-56, -0x1.6dc66cp
-84 }, /* 246: log(1066/1024.) */
211 { +0x1.22c71cp
-5, -0x1.8abe2cp
-32, -0x1.7e7078p
-56, -0x1.ddc3dcp
-86 }, /* 247: log(1061/1024.) */
212 { +0x1.03d5d8p
-5, +0x1.79cfbcp
-31, -0x1.da7c4cp
-58, +0x1.4e7582p
-83 }, /* 248: log(1057/1024.) */
213 { +0x1.c98d18p
-6, +0x1.a01904p
-31, -0x1.854164p
-55, +0x1.883c36p
-79 }, /* 249: log(1053/1024.) */
214 { +0x1.8b31fcp
-6, -0x1.356500p
-30, +0x1.c3ab48p
-55, +0x1.b69bdap
-80 }, /* 250: log(1049/1024.) */
215 { +0x1.3cea44p
-6, +0x1.a352bcp
-33, -0x1.8865acp
-57, -0x1.48159cp
-81 }, /* 251: log(1044/1024.) */
216 { +0x1.fc0a8cp
-7, -0x1.e07f84p
-32, +0x1.e7cf6cp
-58, +0x1.3a69c0p
-82 }, /* 252: log(1040/1024.) */
217 { +0x1.7dc474p
-7, +0x1.f810a8p
-31, -0x1.245b5cp
-56, -0x1.a1f4f8p
-80 }, /* 253: log(1036/1024.) */
218 { +0x1.fe02a8p
-8, -0x1.4ef988p
-32, +0x1.1f86ecp
-57, +0x1.20723cp
-81 }, /* 254: log(1032/1024.) */
219 { +0x1.ff00acp
-9, -0x1.d4ef44p
-33, +0x1.2821acp
-63, +0x1.5a6d32p
-87 }, /* 255: log(1028/1024.) */
223 #define PAIR_ADD(d, H, L) do { \
224 gnm_float d_ = (d); \
225 gnm_float dh_ = gnm_floor (d_ * 65536 + 0.5) / 65536; \
226 gnm_float dl_ = d_ - dh_; \
227 if (0) g_printerr ("Adding %.50g (%a)\n", d_, d_); \
232 #define ADD1(d) PAIR_ADD(d,*yh,*yl)
235 * Compute x * gnm_log (x / M) + (M - x)
236 * aka -x * log1pmx ((M - x) / x)
238 * Deliver the result back in two parts, *yh and *yl.
242 ebd0(gnm_float x
, gnm_float M
, gnm_float
*yh
, gnm_float
*yl
)
244 gnm_float r
, f
, fg
, M1
;
248 const double S
= 1u << Sb
;
249 const int N
= G_N_ELEMENTS(bd0_scale
) - 1;
250 const double e4
= GNM_EPSILON
* GNM_EPSILON
* GNM_EPSILON
* GNM_EPSILON
;
252 if (gnm_isnan (x
) || gnm_isnan (M
)) {
260 if (x
< M
* e4
) { ADD1(M
); return; }
261 if (M
== 0) { *yh
= gnm_pinf
; return; }
264 ADD1(x
* (gnm_log(x
) - 1));
265 ADD1(-x
* gnm_log(M
));
270 g_printerr ("x=%.20g M=%.20g\n", x
, M
);
273 r
= gnm_frexp (M
/ x
, &e
);
274 i
= gnm_floor ((r
- 0.5) * (2 * N
) + 0.5);
275 g_assert (i
>= 0 && i
<= N
);
276 f
= gnm_floor (S
/ (0.5 + i
/ (2.0 * N
)) + 0.5);
277 fg
= gnm_ldexp (f
, -(e
+ Sb
));
279 /* We now have (M * fg / x) close to 1. */
282 * We need to compute this:
283 * (x/M)^x * exp(M-x) =
284 * (M/x)^-x * exp(M-x) =
285 * (M*fg/x)^-x * (fg)^x * exp(M-x) =
286 * (M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg)
289 * log((x/M)^x * exp(M-x)) =
290 * log((M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg)) =
291 * log((M*fg/x)^-x * exp(M*fg-x)) + x*log(fg) + (M-M*fg) =
292 * -x*log1pmx((M*fg-x)/x) + x*log(fg) + M - M*fg =
294 * Note, that fg has at most 10 bits. If M and x are suitably
295 * "nice" -- such as being integers or half-integers -- then
296 * we can compute M*fg as well as x * bd0_scale[.][.] without
300 for (j
= G_N_ELEMENTS(bd0_scale
[i
]) - 1; j
>= 0; j
--) {
301 ADD1(x
* bd0_scale
[i
][j
]); /* Handles x*log(fg*2^e) */
302 ADD1(-x
* e
* bd0_scale
[0][j
]); /* Handles x*log(1/2^e) */
306 M1
= gnm_floor (M
+ 0.5);
310 ADD1(-x
* log1pmx ((M
* fg
- x
) / x
));
313 g_printerr ("res=%.20g + %.20g\n", *yh
, *yl
);
319 /* Legacy function. */
321 bd0(gnm_float x
, gnm_float M
)
324 ebd0 (x
, M
, &yh
, &yl
);
328 /* ------------------------------------------------------------------------- */
329 /* --- BEGIN MAGIC R SOURCE MARKER --- */
331 /* The following source code was imported from the R project. */
332 /* It was automatically transformed by tools/import-R. */
334 /* Imported src/nmath/dpq.h from R. */
335 /* Utilities for `dpq' handling (density/probability/quantile) */
337 /* give_log in "d"; log_p in "p" & "q" : */
338 #define give_log log_p
341 #define R_D__0 (log_p ? gnm_ninf : 0.) /* 0 */
342 #define R_D__1 (log_p ? 0. : 1.) /* 1 */
343 #define R_DT_0 (lower_tail ? R_D__0 : R_D__1) /* 0 */
344 #define R_DT_1 (lower_tail ? R_D__1 : R_D__0) /* 1 */
346 #define R_D_Lval(p) (lower_tail ? (p) : (1 - (p))) /* p */
347 #define R_D_Cval(p) (lower_tail ? (1 - (p)) : (p)) /* 1 - p */
349 #define R_D_val(x) (log_p ? gnm_log(x) : (x)) /* x in pF(x,..) */
350 #define R_D_qIv(p) (log_p ? gnm_exp(p) : (p)) /* p in qF(p,..) */
351 #define R_D_exp(x) (log_p ? (x) : gnm_exp(x)) /* gnm_exp(x) */
352 #define R_D_log(p) (log_p ? (p) : gnm_log(p)) /* gnm_log(p) */
353 #define R_D_Clog(p) (log_p ? gnm_log1p(-(p)) : (1 - (p)))/* [log](1-p) */
355 /* gnm_log(1-gnm_exp(x)): R_D_LExp(x) == (gnm_log1p(- R_D_qIv(x))) but even more stable:*/
356 #define R_D_LExp(x) (log_p ? swap_log_tail(x) : gnm_log1p(-x))
359 * #define R_DT_val(x) R_D_val(R_D_Lval(x))
360 * #define R_DT_Cval(x) R_D_val(R_D_Cval(x)) */
361 #define R_DT_val(x) (lower_tail ? R_D_val(x) : R_D_Clog(x))
362 #define R_DT_Cval(x) (lower_tail ? R_D_Clog(x) : R_D_val(x))
364 /*#define R_DT_qIv(p) R_D_Lval(R_D_qIv(p)) * p in qF ! */
365 #define R_DT_qIv(p) (log_p ? (lower_tail ? gnm_exp(p) : - gnm_expm1(p)) \
368 /*#define R_DT_CIv(p) R_D_Cval(R_D_qIv(p)) * 1 - p in qF */
369 #define R_DT_CIv(p) (log_p ? (lower_tail ? -gnm_expm1(p) : gnm_exp(p)) \
372 #define R_DT_exp(x) R_D_exp(R_D_Lval(x)) /* gnm_exp(x) */
373 #define R_DT_Cexp(x) R_D_exp(R_D_Cval(x)) /* gnm_exp(1 - x) */
375 #define R_DT_log(p) (lower_tail? R_D_log(p) : R_D_LExp(p))/* gnm_log(p) in qF */
376 #define R_DT_Clog(p) (lower_tail? R_D_LExp(p): R_D_log(p))/* gnm_log1p (-p) in qF*/
377 #define R_DT_Log(p) (lower_tail? (p) : swap_log_tail(p))
378 /* == R_DT_log when we already "know" log_p == TRUE :*/
381 #define R_Q_P01_check(p) \
382 if ((log_p && p > 0) || \
383 (!log_p && (p < 0 || p > 1)) ) \
387 /* additions for density functions (C.Loader) */
388 #define R_D_fexp(f,x) (give_log ? -0.5*gnm_log(f)+(x) : gnm_exp(x)/gnm_sqrt(f))
389 #define R_D_forceint(x) gnm_floor((x) + 0.5)
390 #define R_D_nonint(x) (gnm_abs((x) - gnm_floor((x)+0.25)) > 1e-7)
391 /* [neg]ative or [non int]eger : */
392 #define R_D_negInonint(x) (x < 0. || R_D_nonint(x))
394 #define R_D_nonint_check(x) \
395 if(R_D_nonint(x)) { \
396 MATHLIB_WARNING("non-integer x = %" GNM_FORMAT_f "", x); \
400 /* ------------------------------------------------------------------------ */
401 /* Imported src/nmath/ftrunc.c from R. */
403 * Mathlib : A C Library of Special Functions
404 * Copyright (C) 1998 Ross Ihaka
406 * This program is free software; you can redistribute it and/or modify
407 * it under the terms of the GNU General Public License as published by
408 * the Free Software Foundation; either version 2 of the License, or
409 * (at your option) any later version.
411 * This program is distributed in the hope that it will be useful,
412 * but WITHOUT ANY WARRANTY; without even the implied warranty of
413 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
414 * GNU General Public License for more details.
416 * You should have received a copy of the GNU General Public License
417 * along with this program; if not, write to the Free Software
418 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
424 * double ftrunc(double x);
428 * Truncation toward zero.
432 gnm_float
gnm_trunc(gnm_float x
)
434 if(x
>= 0) return gnm_floor(x
);
435 else return gnm_ceil(x
);
438 /* ------------------------------------------------------------------------ */
439 /* Imported src/nmath/pnorm.c from R. */
441 * Mathlib : A C Library of Special Functions
442 * Copyright (C) 1998 Ross Ihaka
443 * Copyright (C) 2000-2002 The R Development Core Team
444 * Copyright (C) 2003 The R Foundation
446 * This program is free software; you can redistribute it and/or modify
447 * it under the terms of the GNU General Public License as published by
448 * the Free Software Foundation; either version 2 of the License, or
449 * (at your option) any later version.
451 * This program is distributed in the hope that it will be useful,
452 * but WITHOUT ANY WARRANTY; without even the implied warranty of
453 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
454 * GNU General Public License for more details.
456 * You should have received a copy of the GNU General Public License
457 * along with this program; if not, write to the Free Software
458 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
465 * double pnorm5(double x, double mu, double sigma, int lower_tail,int log_p);
466 * {pnorm (..) is synonymous and preferred inside R}
468 * void pnorm_both(double x, double *cum, double *ccum,
469 * int i_tail, int log_p);
473 * The main computation evaluates near-minimax approximations derived
474 * from those in "Rational Chebyshev approximations for the error
475 * function" by W. J. Cody, Math. Comp., 1969, 631-637. This
476 * transportable program uses rational functions that theoretically
477 * approximate the normal distribution function to at least 18
478 * significant decimal digits. The accuracy achieved depends on the
479 * arithmetic system, the compiler, the intrinsic functions, and
480 * proper selection of the machine-dependent constants.
484 * Cody, W. D. (1993).
485 * ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
486 * Special Function Routines and Test Drivers".
487 * ACM Transactions on Mathematical Software. 19, 22-32.
491 * The "_both" , lower, upper, and log_p variants were added by
492 * Martin Maechler, Jan.2000;
493 * as well as log1p() and similar improvements later on.
495 * James M. Rath contributed bug report PR#699 and patches correcting SIXTEN
496 * and if() clauses {with a bug: "|| instead of &&" -> PR #2883) more in line
497 * with the original Cody code.
500 gnm_float
pnorm(gnm_float x
, gnm_float mu
, gnm_float sigma
, gboolean lower_tail
, gboolean log_p
)
502 gnm_float p
, cp
= 0.;
504 /* Note: The structure of these checks has been carefully thought through.
505 * For example, if x == mu and sigma == 0, we get the correct answer 1.
508 if(gnm_isnan(x
) || gnm_isnan(mu
) || gnm_isnan(sigma
))
509 return x
+ mu
+ sigma
;
511 if(!gnm_finite(x
) && mu
== x
) return gnm_nan
;/* x-mu is NaN */
513 if(sigma
< 0) ML_ERR_return_NAN
;
515 return (x
< mu
) ? R_DT_0
: R_DT_1
;
517 p
= (x
- mu
) / sigma
;
519 return (x
< mu
) ? R_DT_0
: R_DT_1
;
522 pnorm_both(x
, &p
, &cp
, (lower_tail
? 0 : 1), log_p
);
524 return(lower_tail
? p
: cp
);
527 #define SIXTEN 16 /* Cutoff allowing exact "*" and "/" */
529 void pnorm_both(gnm_float x
, gnm_float
*cum
, gnm_float
*ccum
, int i_tail
, gboolean log_p
)
531 /* i_tail in {0,1,2} means: "lower", "upper", or "both" :
532 if(lower) return *cum := P[X <= x]
533 if(upper) return *ccum := P[X > x] = 1 - P[X <= x]
535 static const gnm_float a
[5] = {
536 GNM_const(2.2352520354606839287),
537 GNM_const(161.02823106855587881),
538 GNM_const(1067.6894854603709582),
539 GNM_const(18154.981253343561249),
540 GNM_const(0.065682337918207449113)
542 static const gnm_float b
[4] = {
543 GNM_const(47.20258190468824187),
544 GNM_const(976.09855173777669322),
545 GNM_const(10260.932208618978205),
546 GNM_const(45507.789335026729956)
548 static const gnm_float c
[9] = {
549 GNM_const(0.39894151208813466764),
550 GNM_const(8.8831497943883759412),
551 GNM_const(93.506656132177855979),
552 GNM_const(597.27027639480026226),
553 GNM_const(2494.5375852903726711),
554 GNM_const(6848.1904505362823326),
555 GNM_const(11602.651437647350124),
556 GNM_const(9842.7148383839780218),
557 GNM_const(1.0765576773720192317e-8)
559 static const gnm_float d
[8] = {
560 GNM_const(22.266688044328115691),
561 GNM_const(235.38790178262499861),
562 GNM_const(1519.377599407554805),
563 GNM_const(6485.558298266760755),
564 GNM_const(18615.571640885098091),
565 GNM_const(34900.952721145977266),
566 GNM_const(38912.003286093271411),
567 GNM_const(19685.429676859990727)
569 static const gnm_float p
[6] = {
570 GNM_const(0.21589853405795699),
571 GNM_const(0.1274011611602473639),
572 GNM_const(0.022235277870649807),
573 GNM_const(0.001421619193227893466),
574 GNM_const(2.9112874951168792e-5),
575 GNM_const(0.02307344176494017303)
577 static const gnm_float q
[5] = {
578 GNM_const(1.28426009614491121),
579 GNM_const(0.468238212480865118),
580 GNM_const(0.0659881378689285515),
581 GNM_const(0.00378239633202758244),
582 GNM_const(7.29751555083966205e-5)
585 gnm_float xden
, xnum
, temp
, del
, eps
, xsq
, y
;
587 gnm_float min
= GNM_MIN
;
592 if(gnm_isnan(x
)) { *cum
= *ccum
= x
; return; }
595 /* Consider changing these : */
596 eps
= GNM_EPSILON
* 0.5;
598 /* i_tail in {0,1,2} =^= {lower, upper, both} */
603 if (y
<= GNM_const(0.67448975)) { /* qnorm(3/4) = .6744.... -- earlier had GNM_const(0.66291) */
608 for (i
= 0; i
< 3; ++i
) {
609 xnum
= (xnum
+ a
[i
]) * xsq
;
610 xden
= (xden
+ b
[i
]) * xsq
;
612 } else xnum
= xden
= 0.0;
614 temp
= x
* (xnum
+ a
[3]) / (xden
+ b
[3]);
615 if(lower
) *cum
= 0.5 + temp
;
616 if(upper
) *ccum
= 0.5 - temp
;
618 if(lower
) *cum
= gnm_log(*cum
);
619 if(upper
) *ccum
= gnm_log(*ccum
);
622 else if (y
<= M_SQRT_32
) {
624 /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= gnm_sqrt(32) ~= 5.657 */
628 for (i
= 0; i
< 7; ++i
) {
629 xnum
= (xnum
+ c
[i
]) * y
;
630 xden
= (xden
+ d
[i
]) * y
;
632 temp
= (xnum
+ c
[7]) / (xden
+ d
[7]);
635 xsq = gnm_trunc(X * SIXTEN) / SIXTEN; \
636 del = (X - xsq) * (X + xsq); \
638 *cum = (-xsq * xsq * 0.5) + (-del * 0.5) + gnm_log(temp); \
639 if((lower && x > 0.) || (upper && x <= 0.)) \
640 *ccum = gnm_log1p(-gnm_exp(-xsq * xsq * 0.5) * \
641 gnm_exp(-del * 0.5) * temp); \
644 *cum = gnm_exp(-xsq * xsq * 0.5) * gnm_exp(-del * 0.5) * temp; \
645 *ccum = 1.0 - *cum; \
649 if (x > 0.) {/* swap ccum <--> cum */ \
650 temp = *cum; if(lower) *cum = *ccum; *ccum = temp; \
657 /* else |x| > gnm_sqrt(32) = 5.657 :
658 * the next two case differentiations were really for lower=T, log=F
659 * Particularly *not* for log_p !
661 * Cody had (-37.5193 < x && x < 8.2924) ; R originally had y < 50
663 * Note that we do want symmetry(0), lower/upper -> hence use y
666 /* ^^^^^ MM FIXME: can speedup for log_p and much larger |x| !
667 * Then, make use of Abramowitz & Stegun, 26.2.13, something like
671 if(xsq * GNM_EPSILON < 1.)
672 del = (1. - (1. - 5./(xsq+6.)) / (xsq+4.)) / (xsq+2.);
675 *cum = -.5*xsq - M_LN_SQRT_2PI - gnm_log(x) + gnm_log1p(-del);
676 *ccum = gnm_log1p(-gnm_exp(*cum)); /.* ~ gnm_log(1) = 0 *./
681 || (lower
&& -37.5193 < x
&& x
< 8.2924)
682 || (upper
&& -8.2924 < x
&& x
< 37.5193)
685 /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
689 for (i
= 0; i
< 4; ++i
) {
690 xnum
= (xnum
+ p
[i
]) * xsq
;
691 xden
= (xden
+ q
[i
]) * xsq
;
693 temp
= xsq
* (xnum
+ p
[4]) / (xden
+ q
[4]);
694 temp
= (M_1_SQRT_2PI
- temp
) / y
;
699 else { /* no log_p , large x such that probs are 0 or 1 */
700 if(x
> 0) { *cum
= 1.; *ccum
= 0.; }
701 else { *cum
= 0.; *ccum
= 1.; }
706 /* do not return "denormalized" -- we do in R */
708 if(*cum
> -min
) *cum
= -0.;
709 if(*ccum
> -min
)*ccum
= -0.;
712 if(*cum
< min
) *cum
= 0.;
713 if(*ccum
< min
) *ccum
= 0.;
718 /* Cleaning up done by tools/import-R: */
723 /* ------------------------------------------------------------------------ */
724 /* Imported src/nmath/qnorm.c from R. */
726 * Mathlib : A C Library of Special Functions
727 * Copyright (C) 1998 Ross Ihaka
728 * Copyright (C) 2000 The R Development Core Team
729 * based on AS 111 (C) 1977 Royal Statistical Society
730 * and on AS 241 (C) 1988 Royal Statistical Society
732 * This program is free software; you can redistribute it and/or modify
733 * it under the terms of the GNU General Public License as published by
734 * the Free Software Foundation; either version 2 of the License, or
735 * (at your option) any later version.
737 * This program is distributed in the hope that it will be useful,
738 * but WITHOUT ANY WARRANTY; without even the implied warranty of
739 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
740 * GNU General Public License for more details.
742 * You should have received a copy of the GNU General Public License
743 * along with this program; if not, write to the Free Software
744 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
749 * double qnorm5(double p, double mu, double sigma,
750 * int lower_tail, int log_p)
751 * {qnorm (..) is synonymous and preferred inside R}
755 * Compute the quantile function for the normal distribution.
757 * For small to moderate probabilities, algorithm referenced
758 * below is used to obtain an initial approximation which is
759 * polished with a final Newton step.
761 * For very large arguments, an algorithm of Wichura is used.
765 * Beasley, J. D. and S. G. Springer (1977).
766 * Algorithm AS 111: The percentage points of the normal distribution,
767 * Applied Statistics, 26, 118-121.
769 * Wichura, M.J. (1988).
770 * Algorithm AS 241: The Percentage Points of the Normal Distribution.
771 * Applied Statistics, 37, 477-484.
775 gnm_float
qnorm(gnm_float p
, gnm_float mu
, gnm_float sigma
, gboolean lower_tail
, gboolean log_p
)
777 gnm_float p_
, q
, r
, val
;
780 if (gnm_isnan(p
) || gnm_isnan(mu
) || gnm_isnan(sigma
))
781 return p
+ mu
+ sigma
;
783 if (p
== R_DT_0
) return gnm_ninf
;
784 if (p
== R_DT_1
) return gnm_pinf
;
787 if(sigma
< 0) ML_ERR_return_NAN
;
788 if(sigma
== 0) return mu
;
790 p_
= R_DT_qIv(p
);/* real lower_tail prob. p */
794 REprintf("qnorm(p=%10.7" GNM_FORMAT_g
", m=%" GNM_FORMAT_g
", s=%" GNM_FORMAT_g
", l.t.= %d, log= %d): q = %" GNM_FORMAT_g
"\n",
795 p
,mu
,sigma
, lower_tail
, log_p
, q
);
799 /*-- use AS 241 --- */
800 /* gnm_float ppnd16_(gnm_float *p, long *ifault)*/
801 /* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
803 Produces the normal deviate Z corresponding to a given lower
804 tail area of P; Z is accurate to about 1 part in 10**16.
806 (original fortran code used PARAMETER(..) for the coefficients
807 and provided hash codes for checking them...)
809 if (gnm_abs(q
) <= .425) {/* 0.075 <= p <= 0.925 */
810 r
= GNM_const(.180625) - q
* q
;
812 q
* (((((((r
* GNM_const(2509.0809287301226727) +
813 GNM_const(33430.575583588128105)) * r
+ GNM_const(67265.770927008700853)) * r
+
814 GNM_const(45921.953931549871457)) * r
+ GNM_const(13731.693765509461125)) * r
+
815 GNM_const(1971.5909503065514427)) * r
+ GNM_const(133.14166789178437745)) * r
+
816 GNM_const(3.387132872796366608))
817 / (((((((r
* GNM_const(5226.495278852854561) +
818 GNM_const(28729.085735721942674)) * r
+ GNM_const(39307.89580009271061)) * r
+
819 GNM_const(21213.794301586595867)) * r
+ GNM_const(5394.1960214247511077)) * r
+
820 GNM_const(687.1870074920579083)) * r
+ GNM_const(42.313330701600911252)) * r
+ 1.);
822 else { /* closer than 0.075 from {0,1} boundary */
824 /* r = min(p, 1-p) < 0.075 */
826 r
= R_DT_CIv(p
);/* 1-p */
828 r
= p_
;/* = R_DT_Iv(p) ^= p */
830 r
= gnm_sqrt(- ((log_p
&&
831 ((lower_tail
&& q
<= 0) || (!lower_tail
&& q
> 0))) ?
832 p
: /* else */ gnm_log(r
)));
833 /* r = gnm_sqrt(-gnm_log(r)) <==> min(p, 1-p) = gnm_exp( - r^2 ) */
835 REprintf("\t close to 0 or 1: r = %7" GNM_FORMAT_g
"\n", r
);
838 if (r
<= 5.) { /* <==> min(p,1-p) >= gnm_exp(-25) ~= 1.3888e-11 */
840 val
= (((((((r
* GNM_const(7.7454501427834140764e-4) +
841 GNM_const(.0227238449892691845833)) * r
+ GNM_const(.24178072517745061177)) *
842 r
+ GNM_const(1.27045825245236838258)) * r
+
843 GNM_const(3.64784832476320460504)) * r
+ GNM_const(5.7694972214606914055)) *
844 r
+ GNM_const(4.6303378461565452959)) * r
+
845 GNM_const(1.42343711074968357734))
847 GNM_const(1.05075007164441684324e-9) + GNM_const(5.475938084995344946e-4)) *
848 r
+ GNM_const(.0151986665636164571966)) * r
+
849 GNM_const(.14810397642748007459)) * r
+ GNM_const(.68976733498510000455)) *
850 r
+ GNM_const(1.6763848301838038494)) * r
+
851 GNM_const(2.05319162663775882187)) * r
+ 1.);
853 else { /* very close to 0 or 1 */
855 val
= (((((((r
* GNM_const(2.01033439929228813265e-7) +
856 GNM_const(2.71155556874348757815e-5)) * r
+
857 GNM_const(.0012426609473880784386)) * r
+ GNM_const(.026532189526576123093)) *
858 r
+ GNM_const(.29656057182850489123)) * r
+
859 GNM_const(1.7848265399172913358)) * r
+ GNM_const(5.4637849111641143699)) *
860 r
+ GNM_const(6.6579046435011037772))
862 GNM_const(2.04426310338993978564e-15) + GNM_const(1.4215117583164458887e-7))*
863 r
+ GNM_const(1.8463183175100546818e-5)) * r
+
864 GNM_const(7.868691311456132591e-4)) * r
+ GNM_const(.0148753612908506148525))
865 * r
+ GNM_const(.13692988092273580531)) * r
+
866 GNM_const(.59983220655588793769)) * r
+ 1.);
871 /* return (q >= 0.)? r : -r ;*/
873 return mu
+ sigma
* val
;
877 /* ------------------------------------------------------------------------ */
878 /* Imported src/nmath/ppois.c from R. */
880 * Mathlib : A C Library of Special Functions
881 * Copyright (C) 1998 Ross Ihaka
882 * Copyright (C) 2000 The R Development Core Team
884 * This program is free software; you can redistribute it and/or modify
885 * it under the terms of the GNU General Public License as published by
886 * the Free Software Foundation; either version 2 of the License, or
887 * (at your option) any later version.
889 * This program is distributed in the hope that it will be useful,
890 * but WITHOUT ANY WARRANTY; without even the implied warranty of
891 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
892 * GNU General Public License for more details.
894 * You should have received a copy of the GNU General Public License
895 * along with this program; if not, write to the Free Software
896 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
901 * The distribution function of the Poisson distribution.
905 gnm_float
ppois(gnm_float x
, gnm_float lambda
, gboolean lower_tail
, gboolean log_p
)
908 if (gnm_isnan(x
) || gnm_isnan(lambda
))
911 if(lambda
< 0.) ML_ERR_return_NAN
;
913 x
= gnm_fake_floor(x
);
914 if (x
< 0) return R_DT_0
;
915 if (lambda
== 0.) return R_DT_1
;
916 if (!gnm_finite(x
)) return R_DT_1
;
918 return pgamma(lambda
, x
+ 1, 1., !lower_tail
, log_p
);
921 /* ------------------------------------------------------------------------ */
922 /* Imported src/nmath/dpois.c from R. */
925 * Catherine Loader, catherine@research.bell-labs.com.
929 * Copyright (C) 2000, The R Core Development Team
931 * This program is free software; you can redistribute it and/or modify
932 * it under the terms of the GNU General Public License as published by
933 * the Free Software Foundation; either version 2 of the License, or
934 * (at your option) any later version.
936 * This program is distributed in the hope that it will be useful,
937 * but WITHOUT ANY WARRANTY; without even the implied warranty of
938 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
939 * GNU General Public License for more details.
941 * You should have received a copy of the GNU General Public License
942 * along with this program; if not, write to the Free Software
943 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
949 * dpois() checks argument validity and calls dpois_raw().
951 * dpois_raw() computes the Poisson probability lb^x exp(-lb) / x!.
952 * This does not check that x is an integer, since dgamma() may
953 * call this with a fractional x argument. Any necessary argument
954 * checks should be done in the calling function.
959 static gnm_float
dpois_raw(gnm_float x
, gnm_float lambda
, gboolean give_log
)
961 gnm_float yh
, yl
, f2
;
963 /* x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
966 if (lambda
== 0) return( (x
== 0) ? R_D__1
: R_D__0
);
967 if (!gnm_finite(lambda
)) return R_D__0
;
968 if (x
< 0) return( R_D__0
);
969 if (x
<= lambda
* GNM_MIN
) return(R_D_exp(-lambda
) );
970 if (lambda
< x
* GNM_MIN
) return(R_D_exp(-lambda
+ x
*gnm_log(lambda
) -lgamma1p (x
)));
972 ebd0 (x
, lambda
, &yh
, &yl
);
973 PAIR_ADD (stirlerr(x
), yh
, yl
);
977 ? -yl
- yh
- 0.5 * gnm_log (f2
)
978 : gnm_exp (-yl
) * gnm_exp (-yh
) / gnm_sqrt (f2
);
981 gnm_float
dpois(gnm_float x
, gnm_float lambda
, gboolean give_log
)
984 if(gnm_isnan(x
) || gnm_isnan(lambda
))
988 if (lambda
< 0) ML_ERR_return_NAN
;
990 if (x
< 0 || !gnm_finite(x
))
995 return( dpois_raw(x
,lambda
,give_log
) );
998 /* ------------------------------------------------------------------------ */
999 /* Imported src/nmath/dgamma.c from R. */
1002 * Catherine Loader, catherine@research.bell-labs.com.
1006 * Copyright (C) 2000 The R Core Development Team
1007 * Copyright (C) 2004 The R Foundation
1009 * This program is free software; you can redistribute it and/or modify
1010 * it under the terms of the GNU General Public License as published by
1011 * the Free Software Foundation; either version 2 of the License, or
1012 * (at your option) any later version.
1014 * This program is distributed in the hope that it will be useful,
1015 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1016 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1017 * GNU General Public License for more details.
1019 * You should have received a copy of the GNU General Public License
1020 * along with this program; if not, write to the Free Software
1021 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1027 * Computes the density of the gamma distribution,
1029 * 1/s (x/s)^{a-1} exp(-x/s)
1030 * p(x;a,s) = -----------------------
1033 * where `s' is the scale (= 1/lambda in other parametrizations)
1034 * and `a' is the shape parameter ( = alpha in other contexts).
1036 * The old (R 1.1.1) version of the code is available via `#define D_non_pois'
1040 gnm_float
dgamma(gnm_float x
, gnm_float shape
, gnm_float scale
, gboolean give_log
)
1044 if (gnm_isnan(x
) || gnm_isnan(shape
) || gnm_isnan(scale
))
1045 return x
+ shape
+ scale
;
1047 if (shape
<= 0 || scale
<= 0) ML_ERR_return_NAN
;
1051 if (shape
< 1) return gnm_pinf
;
1052 if (shape
> 1) return R_D__0
;
1054 return give_log
? -gnm_log(scale
) : 1 / scale
;
1058 pr
= dpois_raw(shape
, x
/scale
, give_log
);
1059 return give_log
? pr
+ gnm_log(shape
/x
) : pr
*shape
/x
;
1061 /* else shape >= 1 */
1062 pr
= dpois_raw(shape
-1, x
/scale
, give_log
);
1063 return give_log
? pr
- gnm_log(scale
) : pr
/scale
;
1066 /* ------------------------------------------------------------------------ */
1067 /* Imported src/nmath/pgamma.c from R. */
1069 * Mathlib : A C Library of Special Functions
1070 * Copyright (C) 2005 Morten Welinder <terra@gnome.org>
1071 * Copyright (C) 2005 The R Foundation
1073 * This program is free software; you can redistribute it and/or modify
1074 * it under the terms of the GNU General Public License as published by
1075 * the Free Software Foundation; either version 2 of the License, or
1076 * (at your option) any later version.
1078 * This program is distributed in the hope that it will be useful,
1079 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1080 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1081 * GNU General Public License for more details.
1083 * A copy of the GNU General Public License is available via WWW at
1084 * http://www.gnu.org/copyleft/gpl.html. You can also obtain it by
1085 * writing to the Free Software Foundation, Inc., 51 Franklin St, Fifth
1086 * Floor, Boston, MA 02110-1301 USA.
1090 * #include <Rmath.h>
1092 * double pgamma (double x, double alph, double scale,
1093 * int lower_tail, int log_p)
1095 * double log1pmx (double x)
1096 * double lgamma1p (double a)
1098 * double logspace_add (double logx, double logy)
1099 * double logspace_sub (double logx, double logy)
1104 * This function computes the distribution function for the
1105 * gamma distribution with shape parameter alph and scale parameter
1106 * scale. This is also known as the incomplete gamma function.
1107 * See Abramowitz and Stegun (6.5.1) for example.
1111 * Complete redesign by Morten Welinder, originally for Gnumeric.
1112 * Improvements (e.g. "while NEEDED_SCALE") by Martin Maechler
1113 * The old version can be activated by compiling with -DR_USE_OLD_PGAMMA
1119 /*----------- DEBUGGING -------------
1120 * make CFLAGS='-DDEBUG_p -g -I/usr/local/include -I../include'
1124 /* Scalefactor:= (2^32)^8 = 2^256 = GNM_const(1.157921e+77) */
1125 #define SQR(x) ((x)*(x))
1126 static const gnm_float scalefactor
= SQR(SQR(SQR(4294967296.0)));
1129 /* If |x| > |k| * M_cutoff, then log[ gnm_exp(-x) * k^x ] =~= -x */
1130 static const gnm_float M_cutoff
= M_LN2gnum
* GNM_MAX_EXP
/ GNM_EPSILON
;/*=GNM_const(3.196577e18)*/
1132 /* Continued fraction for calculation of
1133 * 1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ...
1135 * auxilary in log1pmx() and lgamma1p()
1138 gnm_logcf (gnm_float x
, gnm_float i
, gnm_float d
)
1140 gnm_float c1
= 2 * d
;
1141 gnm_float c2
= i
+ d
;
1142 gnm_float c4
= c2
+ d
;
1144 gnm_float b1
= i
* (c2
- i
* x
);
1145 gnm_float b2
= d
* d
* x
;
1146 gnm_float a2
= c4
* c2
- b2
;
1147 const gnm_float cfVSmall
= 1.0e-14;/* ~= relative tolerance */
1154 b2
= c4
* b1
- i
* b2
;
1156 while (gnm_abs (a2
* b1
- a1
* b2
) > gnm_abs (cfVSmall
* b1
* b2
)) {
1157 gnm_float c3
= c2
*c2
*x
;
1160 a1
= c4
* a2
- c3
* a1
;
1161 b1
= c4
* b2
- c3
* b1
;
1166 a2
= c4
* a1
- c3
* a2
;
1167 b2
= c4
* b1
- c3
* b2
;
1169 if (gnm_abs (b2
) > scalefactor
) {
1174 } else if (gnm_abs (b2
) < 1 / scalefactor
) {
1185 /* Accurate calculation of gnm_log1p (x)-x, particularly for small x. */
1186 gnm_float
log1pmx (gnm_float x
)
1188 static const gnm_float minLog1Value
= GNM_const(-0.79149064);
1189 static const gnm_float two
= 2;
1191 if (x
> 1 || x
< minLog1Value
)
1192 return gnm_log1p(x
) - x
;
1193 else { /* expand in [x/(2+x)]^2 */
1194 gnm_float term
= x
/ (2 + x
);
1195 gnm_float y
= term
* term
;
1196 if (gnm_abs(x
) < 1e-2)
1197 return term
* ((((two
/ 9 * y
+ two
/ 7) * y
+ two
/ 5) * y
+
1200 return term
* (2 * y
* gnm_logcf (y
, 3, 2) - x
);
1206 * Compute the log of a sum from logs of terms, i.e.,
1208 * log (exp (logx) + exp (logy))
1210 * without causing overflows and without throwing away large handfuls
1213 gnm_float
logspace_add (gnm_float logx
, gnm_float logy
)
1215 return fmax2 (logx
, logy
) + gnm_log1p (gnm_exp (-gnm_abs (logx
- logy
)));
1220 * Compute the log of a difference from logs of terms, i.e.,
1222 * log (exp (logx) - exp (logy))
1224 * without causing overflows and without throwing away large handfuls
1227 gnm_float
logspace_sub (gnm_float logx
, gnm_float logy
)
1229 return logx
+ gnm_log1p (-gnm_exp (logy
- logx
));
1233 #ifndef R_USE_OLD_PGAMMA
1235 /* dpois_wrap (x_P_1, lambda, g_log) ==
1236 * dpois (x_P_1 - 1, lambda, g_log)
1239 dpois_wrap (gnm_float x_plus_1
, gnm_float lambda
, gboolean give_log
)
1242 REprintf (" dpois_wrap(x+1=%.14" GNM_FORMAT_g
", lambda=%.14" GNM_FORMAT_g
", log=%d)\n",
1243 x_plus_1
, lambda
, give_log
);
1245 if (!gnm_finite(lambda
))
1248 return dpois_raw (x_plus_1
- 1, lambda
, give_log
);
1249 if (lambda
> gnm_abs(x_plus_1
- 1) * M_cutoff
)
1250 return R_D_exp(-lambda
- gnm_lgamma(x_plus_1
));
1252 gnm_float d
= dpois_raw (x_plus_1
, lambda
, give_log
);
1254 REprintf (" -> d=dpois_raw(..)=%.14" GNM_FORMAT_g
"\n", d
);
1257 ? d
+ gnm_log (x_plus_1
/ lambda
)
1258 : d
* (x_plus_1
/ lambda
);
1263 * Abramowitz and Stegun 6.5.29 [right]
1266 pgamma_smallx (gnm_float x
, gnm_float alph
, gboolean lower_tail
, gboolean log_p
)
1268 gnm_float sum
= 0, c
= alph
, n
= 0, term
;
1271 REprintf (" pg_smallx(x=%.12" GNM_FORMAT_g
", alph=%.12" GNM_FORMAT_g
"): ", x
, alph
);
1275 * Relative to 6.5.29 all terms have been multiplied by alph
1276 * and the first, thus being 1, is omitted.
1282 term
= c
/ (alph
+ n
);
1284 } while (gnm_abs (term
) > GNM_EPSILON
* gnm_abs (sum
));
1287 REprintf (" conv.sum=%" GNM_FORMAT_g
";", sum
);
1290 gnm_float f1
= log_p
? gnm_log1p (sum
) : 1 + sum
;
1293 f2
= dpois_raw (alph
, x
, log_p
);
1294 f2
= log_p
? f2
+ x
: f2
* gnm_exp (x
);
1296 f2
= alph
* gnm_log (x
) - lgamma1p (alph
);
1298 f2
= gnm_pow (x
, alph
) / gnm_exp (lgamma1p (alph
));
1300 REprintf (" (f1,f2)= (%" GNM_FORMAT_g
",%" GNM_FORMAT_g
")\n", f1
,f2
);
1302 return log_p
? f1
+ f2
: f1
* f2
;
1304 gnm_float lf2
= alph
* gnm_log (x
) - lgamma1p (alph
);
1306 REprintf (" 1:%.14" GNM_FORMAT_g
" 2:%.14" GNM_FORMAT_g
"\n", alph
* gnm_log (x
), lgamma1p (alph
));
1307 REprintf (" sum=%.14" GNM_FORMAT_g
" gnm_log1p (sum)=%.14" GNM_FORMAT_g
" lf2=%.14" GNM_FORMAT_g
"\n",
1308 sum
, gnm_log1p (sum
), lf2
);
1311 return swap_log_tail (gnm_log1p (sum
) + lf2
);
1313 gnm_float f1m1
= sum
;
1314 gnm_float f2m1
= gnm_expm1 (lf2
);
1315 return -(f1m1
+ f2m1
+ f1m1
* f2m1
);
1318 } /* pgamma_smallx() */
1321 pd_upper_series (gnm_float x
, gnm_float y
, gboolean log_p
)
1323 gnm_float term
= x
/ y
;
1324 gnm_float sum
= term
;
1330 } while (term
> sum
* GNM_EPSILON
);
1332 /* sum = \sum_{n=1}^ oo x^n / (y*(y+1)*...*(y+n-1))
1333 * = \sum_{n=0}^ oo x^(n+1) / (y*(y+1)*...*(y+n))
1334 * = x/y * (1 + \sum_{n=1}^oo x^n / ((y+1)*...*(y+n)))
1335 * ~ x/y + o(x/y) {which happens when alph -> Inf}
1337 return log_p
? gnm_log (sum
) : sum
;
1340 /* Continued fraction for calculation of
1342 * = (i / d) + o(i/d)
1345 pd_lower_cf (gnm_float i
, gnm_float d
)
1347 gnm_float f
= -42, of
, rf
= i
/ d
;
1349 gnm_float c1
= 0, c2
, c3
, c4
;
1350 gnm_float a1
= 0, b1
= 1;
1351 gnm_float a2
= i
, b2
= d
;
1353 #define NEEDED_SCALE \
1354 (b2 > scalefactor) { \
1355 a1 /= scalefactor; \
1356 b1 /= scalefactor; \
1357 a2 /= scalefactor; \
1358 b2 /= scalefactor; \
1361 #define max_it 200000
1364 REprintf("pd_lower_cf(i=%.14" GNM_FORMAT_g
", d=%.14" GNM_FORMAT_g
")\n", i
, d
);
1370 return 0;/* when d >>>> i originally */
1375 while (c1
< max_it
) {
1380 a1
= c4
* a2
+ c3
* a1
;
1381 b1
= c4
* b2
+ c3
* b1
;
1387 a2
= c4
* a1
+ c3
* a2
;
1388 b2
= c4
* b1
+ c3
* b2
;
1395 /* convergence check: relative; absolute for small f : */
1396 if (gnm_abs (f
- of
) <= GNM_EPSILON
* fmax2(rf
, gnm_abs(f
)))
1401 REprintf(" ** NON-convergence in pgamma()'s pd_lower_cf() f= %" GNM_FORMAT_g
".\n", f
);
1402 return f
;/* should not happen ... */
1403 } /* pd_lower_cf() */
1408 pd_lower_series (gnm_float lambda
, gnm_float y
)
1410 gnm_float term
= 1, sum
= 0;
1413 REprintf("pd_lower_series(lam=%.14" GNM_FORMAT_g
", y=%.14" GNM_FORMAT_g
") ...", lambda
, y
);
1415 while (y
>= 1 && term
> sum
* GNM_EPSILON
) {
1420 /* sum = \sum_{n=0}^ oo y*(y-1)*...*(y - n) / lambda^(n+1)
1421 * = y/lambda * (1 + \sum_{n=1}^Inf (y-1)*...*(y-n) / lambda^n
1422 * ~ y/lambda + o(y/lambda)
1425 REprintf(" done: term=%" GNM_FORMAT_g
", sum=%" GNM_FORMAT_g
", y= %" GNM_FORMAT_g
"\n", term
, sum
, y
);
1428 if (y
!= gnm_floor (y
)) {
1430 * The series does not converge as the terms start getting
1431 * bigger (besides flipping sign) for y < -lambda.
1435 REprintf(" y not int: add another term ");
1437 f
= pd_lower_cf (y
, lambda
+ 1 - y
);
1439 REprintf(" (= %.14" GNM_FORMAT_g
") * term = %.14" GNM_FORMAT_g
" to sum %" GNM_FORMAT_g
"\n", f
, term
* f
, sum
);
1445 } /* pd_lower_series() */
1448 * Asymptotic expansion to calculate the probability that poisson variate
1452 ppois_asymp (gnm_float x
, gnm_float lambda
, gboolean lower_tail
, gboolean log_p
)
1454 static const gnm_float coef15
= 1 / GNM_const(12.);
1455 static const gnm_float coef25
= 1 / GNM_const(288.);
1456 static const gnm_float coef35
= -139 / GNM_const(51840.);
1457 static const gnm_float coef45
= -571 / GNM_const(2488320.);
1458 static const gnm_float coef55
= 163879 / GNM_const(209018880.);
1459 static const gnm_float coef65
= 5246819 / GNM_const(75246796800.);
1460 static const gnm_float coef75
= -534703531 / GNM_const(902961561600.);
1461 static const gnm_float coef1
= 2 / GNM_const(3.);
1462 static const gnm_float coef2
= -4 / GNM_const(135.);
1463 static const gnm_float coef3
= 8 / GNM_const(2835.);
1464 static const gnm_float coef4
= 16 / GNM_const(8505.);
1465 static const gnm_float coef5
= -8992 / GNM_const(12629925.);
1466 static const gnm_float coef6
= -334144 / GNM_const(492567075.);
1467 static const gnm_float coef7
= 698752 / GNM_const(1477701225.);
1468 static const gnm_float two
= 2;
1470 gnm_float dfm
, pt_
,s2pt
,res1
,res2
,elfb
,term
;
1471 gnm_float ig2
,ig3
,ig4
,ig5
,ig6
,ig7
,ig25
,ig35
,ig45
,ig55
,ig65
,ig75
;
1472 gnm_float f
, np
, nd
;
1475 pt_
= -x
* log1pmx (dfm
/ x
);
1476 s2pt
= gnm_sqrt (2 * pt_
);
1477 if (dfm
< 0) s2pt
= -s2pt
;
1480 term
= pt_
* pt_
* 0.5;
1491 term
= pt_
* (two
/ 3);
1493 term
*= pt_
* (two
/ 5);
1495 term
*= pt_
* (two
/ 7);
1497 term
*= pt_
* (two
/ 9);
1499 term
*= pt_
* (two
/ 11);
1501 term
*= pt_
* (two
/ 13);
1504 elfb
= ((((((coef75
/x
+ coef65
)/x
+ coef55
)/x
+ coef45
)/x
+ coef35
)/x
+
1505 coef25
)/x
+ coef15
) + x
;
1506 res1
= ((((((ig7
*coef7
/x
+ ig6
*coef6
)/x
+ ig5
*coef5
)/x
+ ig4
*coef4
)/x
+
1507 ig3
*coef3
)/x
+ ig2
*coef2
)/x
+ coef1
)*gnm_sqrt(x
);
1508 res2
= ((((((ig75
*coef75
/x
+ ig65
*coef65
)/x
+ ig55
*coef55
)/x
+ ig45
*coef45
)/
1509 x
+ ig35
*coef35
)/x
+ ig25
*coef25
)/x
+ coef15
)*s2pt
;
1511 if (!lower_tail
) elfb
= -elfb
;
1512 f
= (res1
+ res2
) / elfb
;
1514 np
= pnorm (s2pt
, 0.0, 1.0, !lower_tail
, log_p
);
1515 nd
= dnorm (s2pt
, 0.0, 1.0, log_p
);
1518 REprintf ("pp*_asymp(): f=%.14" GNM_FORMAT_g
" np=%.14" GNM_FORMAT_g
" nd=%.14" GNM_FORMAT_g
" f*nd=%.14" GNM_FORMAT_g
"\n",
1524 ? logspace_add (np
, gnm_log (gnm_abs (f
)) + nd
)
1525 : logspace_sub (np
, gnm_log (gnm_abs (f
)) + nd
);
1528 } /* ppois_asymp() */
1532 pgamma_raw (gnm_float x
, gnm_float alph
, gboolean lower_tail
, gboolean log_p
)
1537 REprintf("pgamma_raw(x=%.14" GNM_FORMAT_g
", alph=%.14" GNM_FORMAT_g
", low=%d, log=%d)\n",
1538 x
, alph
, lower_tail
, log_p
);
1541 res
= pgamma_smallx (x
, alph
, lower_tail
, log_p
);
1542 } else if (x
<= alph
- 1 && x
< 0.8 * (alph
+ 50)) {/* incl. large alph */
1543 gnm_float sum
= pd_upper_series (x
, alph
, log_p
);/* = x/alph + o(x/alph) */
1544 gnm_float d
= dpois_wrap (alph
, x
, log_p
);
1546 REprintf(" alph `large': sum=pd_upper*()= %.12" GNM_FORMAT_g
", d=dpois_w(*)= %.12" GNM_FORMAT_g
" ",
1551 ? swap_log_tail (d
+ sum
)
1554 res
= log_p
? sum
+ d
: sum
* d
;
1555 } else if (alph
- 1 < x
&& alph
< 0.8 * (x
+ 50)) {/* incl. large x */
1557 gnm_float d
= dpois_wrap (alph
, x
, log_p
);
1559 REprintf(" x `large': d=dpois_w(*)= %.14" GNM_FORMAT_g
" ", d
);
1563 if (x
* GNM_EPSILON
> 1 - alph
)
1566 gnm_float f
= pd_lower_cf (alph
, x
- (alph
- 1)) * x
/ alph
;
1567 /* = [alph/(x - alph+1) + o(alph/(x-alph+1))] * x/alph = 1 + o(1) */
1568 sum
= log_p
? gnm_log (f
) : f
;
1571 sum
= pd_lower_series (x
, alph
- 1);/* = (alph-1)/x + o((alph-1)/x) */
1572 sum
= log_p
? gnm_log1p (sum
) : 1 + sum
;
1575 REprintf(", sum= %.14" GNM_FORMAT_g
"\n", sum
);
1578 res
= log_p
? sum
+ d
: sum
* d
;
1581 ? swap_log_tail (d
+ sum
)
1585 REprintf(" using ppois_asymp()\n");
1587 res
= ppois_asymp (alph
- 1, x
, !lower_tail
, log_p
);
1591 * We lose a fair amount of accuracy to underflow in the cases
1592 * where the final result is very close to DBL_MIN. In those
1593 * cases, simply redo via log space.
1595 if (!log_p
&& res
< GNM_MIN
/ GNM_EPSILON
) {
1596 /* with(.Machine, gnm_float.xmin / gnm_float.eps) #|-> GNM_const(1.002084e-292) */
1598 REprintf(" very small res=%.14" GNM_FORMAT_g
"; -> recompute via log\n", res
);
1600 return gnm_exp (pgamma_raw (x
, alph
, lower_tail
, 1));
1606 gnm_float
pgamma(gnm_float x
, gnm_float alph
, gnm_float scale
, gboolean lower_tail
, gboolean log_p
)
1609 if (gnm_isnan(x
) || gnm_isnan(alph
) || gnm_isnan(scale
))
1610 return x
+ alph
+ scale
;
1612 if(alph
<= 0. || scale
<= 0.)
1616 if (gnm_isnan(x
)) /* eg. original x = scale = +Inf */
1619 if (x
<= 0.) /* also for scale=Inf and finite x */
1622 return pgamma_raw (x
, alph
, lower_tail
, log_p
);
1624 /* From: terra@gnome.org (Morten Welinder)
1625 * To: R-bugs@biostat.ku.dk
1626 * Cc: maechler@stat.math.ethz.ch
1627 * Subject: Re: [Rd] pgamma discontinuity (PR#7307)
1628 * Date: Tue, 11 Jan 2005 13:57:26 -0500 (EST)
1630 * this version of pgamma appears to be quite good and certainly a vast
1631 * improvement over current R code. (I last looked at 2.0.1) Apart from
1632 * type naming, this is what I have been using for Gnumeric 1.4.1.
1634 * This could be included into R as-is, but you might want to benefit from
1635 * making gnm_logcf, log1pmx, lgamma1p, and possibly logspace_add/logspace_sub
1636 * available to other parts of R.
1638 * MM: I've not (yet?) taken gnm_logcf(), but the other four
1643 /* R_USE_OLD_PGAMMA */
1645 * Copyright (C) 1998 Ross Ihaka
1646 * Copyright (C) 1999-2000 The R Development Core Team
1647 * Copyright (C) 2003-2004 The R Foundation
1648 * based on AS 239 (C) 1988 Royal Statistical Society
1654 * This function is an adaptation of Algorithm 239 from the
1655 * Applied Statistics Series. The algorithm is faster than
1656 * those by W. Fullerton in the FNLIB library and also the
1657 * TOMS 542 alorithm of W. Gautschi. It provides comparable
1658 * accuracy to those algorithms and is considerably simpler.
1662 * Algorithm AS 239, Incomplete Gamma Function
1663 * Applied Statistics 37, 1988.
1666 gnm_float
pgamma(gnm_float x
, gnm_float alph
, gnm_float scale
, gboolean lower_tail
, gboolean log_p
)
1672 /* normal approx. for alph > alphlimit */
1673 alphlimit
= 1e5
;/* was 1000. till R.1.8.x */
1675 gnm_float pn1
, pn2
, pn3
, pn4
, pn5
, pn6
, arg
, a
, b
, c
, an
, osum
, sum
;
1679 /* check that we have valid values for x and alph */
1682 if (gnm_isnan(x
) || gnm_isnan(alph
) || gnm_isnan(scale
))
1683 return x
+ alph
+ scale
;
1686 REprintf("pgamma(x=%4" GNM_FORMAT_g
", alph=%4" GNM_FORMAT_g
", scale=%4" GNM_FORMAT_g
"): ",x
,alph
,scale
);
1688 if(alph
<= 0. || scale
<= 0.)
1693 REprintf("-> x=%4" GNM_FORMAT_g
"; ",x
);
1696 if (gnm_isnan(x
)) /* eg. original x = scale = Inf */
1703 pn1 = gnm_sqrt(alph) * 3. * (gnm_pow(x/alph, GNM_const(1.)/3.) + 1. / (9. * alph) - 1.); \
1704 return pnorm(pn1, 0., 1., lower_tail, log_p);
1706 if (alph
> alphlimit
) { /* use a normal approximation */
1710 if (x
> xbig
* alph
) {
1711 if (x
> GNM_MAX
* alph
)
1712 /* if x is extremely large __compared to alph__ then return 1 */
1714 else { /* this only "helps" when log_p = TRUE */
1719 if (x
<= 1. || x
< alph
) {
1721 pearson
= 1;/* use pearson's series expansion. */
1723 arg
= alph
* gnm_log(x
) - x
- gnm_lgamma(alph
+ 1.);
1725 REprintf("Pearson arg=%" GNM_FORMAT_g
" ", arg
);
1734 } while (c
> GNM_EPSILON
* sum
);
1736 else { /* x >= max( 1, alph) */
1738 pearson
= 0;/* use a continued fraction expansion */
1740 arg
= alph
* gnm_log(x
) - x
- gnm_lgamma(alph
);
1742 REprintf("Cont.Fract. arg=%" GNM_FORMAT_g
" ", arg
);
1751 for (n
= 1; ; n
++) {
1752 a
+= 1.;/* = n+1 -alph */
1753 b
+= 2.;/* = 2(n+1)-alph+x */
1755 pn5
= b
* pn3
- an
* pn1
;
1756 pn6
= b
* pn4
- an
* pn2
;
1757 if (gnm_abs(pn6
) > 0.) {
1760 if (gnm_abs(osum
- sum
) <= GNM_EPSILON
* fmin2(1., sum
))
1767 if (gnm_abs(pn5
) >= xlarge
) {
1768 /* re-scale the terms in continued fraction if they are large */
1780 arg
+= gnm_log(sum
);
1782 lower_tail
= (lower_tail
== pearson
);
1784 if (log_p
&& lower_tail
)
1787 /* sum = gnm_exp(arg); and return if(lower_tail) sum else 1-sum : */
1788 return (lower_tail
) ? gnm_exp(arg
) : (log_p
? swap_log_tail(arg
) : -gnm_expm1(arg
));
1792 /* R_USE_OLD_PGAMMA */
1793 /* Cleaning up done by tools/import-R: */
1797 /* ------------------------------------------------------------------------ */
1798 /* Imported src/nmath/dt.c from R. */
1801 * Catherine Loader, catherine@research.bell-labs.com.
1805 * Copyright (C) 2000, The R Core Development Team
1807 * This program is free software; you can redistribute it and/or modify
1808 * it under the terms of the GNU General Public License as published by
1809 * the Free Software Foundation; either version 2 of the License, or
1810 * (at your option) any later version.
1812 * This program is distributed in the hope that it will be useful,
1813 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1814 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1815 * GNU General Public License for more details.
1817 * You should have received a copy of the GNU General Public License
1818 * along with this program; if not, write to the Free Software
1819 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1825 * The t density is evaluated as
1826 * sqrt(n/2) / ((n+1)/2) * Gamma((n+3)/2) / Gamma((n+2)/2).
1827 * * (1+x^2/n)^(-n/2)
1828 * / sqrt( 2 pi (1+x^2/n) )
1830 * This form leads to a stable computation for all
1831 * values of n, including n -> 0 and n -> infinity.
1835 gnm_float
dt(gnm_float x
, gnm_float n
, gboolean give_log
)
1839 if (gnm_isnan(x
) || gnm_isnan(n
))
1843 if (n
<= 0) ML_ERR_return_NAN
;
1847 return dnorm(x
, 0., 1., give_log
);
1849 t
= -bd0(n
/2.,(n
+1)/2.) + stirlerr((n
+1)/2.) - stirlerr(n
/2.);
1851 u
= gnm_log1p (x
*x
/n
) * n
/2;
1853 u
= -bd0(n
/2.,(n
+x
*x
)/2.) + x
*x
/2.;
1855 return R_D_fexp(M_2PIgnum
*(1+x
*x
/n
), t
-u
);
1858 /* ------------------------------------------------------------------------ */
1859 /* Imported src/nmath/pt.c from R. */
1861 * R : A Computer Language for Statistical Data Analysis
1862 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
1863 * Copyright (C) 2000 The R Development Core Team
1865 * This program is free software; you can redistribute it and/or modify
1866 * it under the terms of the GNU General Public License as published by
1867 * the Free Software Foundation; either version 2 of the License, or
1868 * (at your option) any later version.
1870 * This program is distributed in the hope that it will be useful,
1871 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1872 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1873 * GNU General Public License for more details.
1875 * You should have received a copy of the GNU General Public License
1876 * along with this program; if not, write to the Free Software
1877 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1882 gnm_float
pt(gnm_float x
, gnm_float n
, gboolean lower_tail
, gboolean log_p
)
1884 /* return P[ T <= x ] where
1885 * T ~ t_{n} (t distrib. with n degrees of freedom).
1887 * --> ./pnt.c for NON-central
1891 if (gnm_isnan(x
) || gnm_isnan(n
))
1894 if (n
<= 0.0) ML_ERR_return_NAN
;
1897 return (x
< 0) ? R_DT_0
: R_DT_1
;
1899 return pnorm(x
, 0.0, 1.0, lower_tail
, log_p
);
1900 if (0 && n
> 4e5
) { /*-- Fixme(?): test should depend on `n' AND `x' ! */
1901 /* Approx. from Abramowitz & Stegun 26.7.8 (p.949) */
1903 return pnorm(x
*(1. - val
)/gnm_sqrt(1. + x
*x
*2.*val
), 0.0, 1.0,
1908 ? pbeta (x
* x
/ (n
+ x
* x
), 0.5, n
/ 2, /*lower_tail*/0, log_p
)
1909 : pbeta (n
/ (n
+ x
* x
), n
/ 2.0, 0.5, /*lower_tail*/1, log_p
);
1911 /* Use "1 - v" if lower_tail and x > 0 (but not both):*/
1913 lower_tail
= !lower_tail
;
1916 if(lower_tail
) return gnm_log1p(-0.5*gnm_exp(val
));
1917 else return val
- M_LN2gnum
; /* = gnm_log(.5* pbeta(....)) */
1921 return R_D_Cval(val
);
1925 /* ------------------------------------------------------------------------ */
1926 /* Imported src/nmath/qt.c from R. */
1928 * Mathlib : A C Library of Special Functions
1929 * Copyright (C) 1998 Ross Ihaka
1930 * Copyright (C) 2000-2002 The R Development Core Team
1931 * Copyright (C) 2003 The R Foundation
1933 * This program is free software; you can redistribute it and/or modify
1934 * it under the terms of the GNU General Public License as published by
1935 * the Free Software Foundation; either version 2 of the License, or
1936 * (at your option) any later version.
1938 * This program is distributed in the hope that it will be useful,
1939 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1940 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1941 * GNU General Public License for more details.
1943 * You should have received a copy of the GNU General Public License
1944 * along with this program; if not, write to the Free Software
1945 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1950 * The "Student" t distribution quantile function.
1954 * This is a C translation of the Fortran routine given in:
1955 * Hill, G.W (1970) "Algorithm 396: Student's t-quantiles"
1956 * CACM 13(10), 619-620.
1959 * - lower_tail, log_p
1960 * - using expm1() : takes care of Lozy (1979) "Remark on Algo.", TOMS
1961 * - Apply 2-term Taylor expansion as in
1962 * Hill, G.W (1981) "Remark on Algo.396", ACM TOMS 7, 250-1
1963 * - Improve the formula decision for 1 < df < 2
1967 gnm_float
qt(gnm_float p
, gnm_float ndf
, gboolean lower_tail
, gboolean log_p
)
1969 const gnm_float eps
= 1.e
-12;
1971 gnm_float a
, b
, c
, d
, p_
, P
, q
, x
, y
;
1975 if (gnm_isnan(p
) || gnm_isnan(ndf
))
1978 if (p
== R_DT_0
) return gnm_ninf
;
1979 if (p
== R_DT_1
) return gnm_pinf
;
1982 if (ndf
< 1) /* FIXME: not yet treated here */
1985 /* FIXME: This test should depend on ndf AND p !!
1986 * ----- and in fact should be replaced by
1987 * something like Abramowitz & Stegun 26.7.5 (p.949)
1989 if (ndf
> 1e20
) return qnorm(p
, 0., 1., lower_tail
, log_p
);
1991 p_
= R_D_qIv(p
); /* note: gnm_exp(p) may underflow to 0; fix later */
1993 if((lower_tail
&& p_
> 0.5) || (!lower_tail
&& p_
< 0.5)) {
1994 neg
= FALSE
; P
= 2 * R_D_Cval(p_
);
1996 neg
= TRUE
; P
= 2 * R_D_Lval(p_
);
1997 } /* 0 <= P <= 1 ; P = 2*min(p_, 1 - p_) in all cases */
1999 if (gnm_abs(ndf
- 2) < eps
) { /* df ~= 2 */
2001 q
= gnm_sqrt(2 / (P
* (2 - P
)) - 2);
2002 else { /* P = 0, but maybe = gnm_exp(p) ! */
2003 if(log_p
) q
= M_SQRT2gnum
* gnm_exp(- .5 * R_D_Lval(p
));
2007 else if (ndf
< 1 + eps
) { /* df ~= 1 (df < 1 excluded above): Cauchy */
2009 q
= gnm_cotpi(P
/ 2);
2011 else { /* P = 0, but maybe p_ = gnm_exp(p) ! */
2012 if(log_p
) q
= M_1_PI
* gnm_exp(-R_D_Lval(p
));/* cot(e) ~ 1/e */
2016 else { /*-- usual case; including, e.g., df = 1.1 */
2017 a
= 1 / (ndf
- 0.5);
2019 c
= ((20700 * a
/ b
- 98) * a
- 16) * a
+ 96.36;
2020 d
= ((94.5 / (b
+ c
) - 3) / b
+ 1) * gnm_sqrt(a
* M_PI_2gnum
) * ndf
;
2022 y
= gnm_pow(d
* P
, 2 / ndf
);
2023 else /* P = 0 && log_p; P = 2*gnm_exp(p) */
2024 y
= gnm_exp(2 / ndf
* (gnm_log(d
) + M_LN2gnum
+ R_D_Lval(p
)));
2026 if ((ndf
< 2.1 && P
> 0.5) || y
> 0.05 + a
) { /* P > P0(df) */
2027 /* Asymptotic inverse expansion about normal */
2029 x
= qnorm(0.5 * P
, 0., 1., /*lower_tail*/TRUE
, /*log_p*/FALSE
);
2030 else /* P = 0 && log_p; P = 2*gnm_exp(p') */
2031 x
= qnorm( p
, 0., 1., lower_tail
, /*log_p*/TRUE
);
2035 c
+= 0.3 * (ndf
- 4.5) * (x
+ 0.6);
2036 c
= (((0.05 * d
* x
- 5) * x
- 7) * x
- 2) * x
+ b
+ c
;
2037 y
= (((((0.4 * y
+ 6.3) * y
+ 36) * y
+ 94.5) / c
2038 - y
- 3) / b
+ 1) * x
;
2039 y
= gnm_expm1(a
* y
* y
);
2041 y
= ((1 / (((ndf
+ 6) / (ndf
* y
) - 0.089 * d
- 0.822)
2042 * (ndf
+ 2) * 3) + 0.5 / (ndf
+ 4))
2043 * y
- 1) * (ndf
+ 1) / (ndf
+ 2) + 1 / y
;
2045 q
= gnm_sqrt(ndf
* y
);
2047 /* Now apply 2-term Taylor expansion improvement (1-term = Newton):
2048 * as by Hill (1981) [ref.above] */
2050 /* FIXME: This is can be far from optimal when log_p = TRUE !
2051 * and probably also improvable when lower_tail = FALSE */
2052 x
= (pt(q
, ndf
, /*lower_tail = */FALSE
, /*log_p = */FALSE
) - P
/2) /
2053 dt(q
, ndf
, /* give_log = */FALSE
);
2054 /* Newton (=Taylor 1 term):
2056 * Taylor 2-term : */
2057 q
+= x
* (1. + x
* q
* (ndf
+ 1) / (2 * (q
* q
+ ndf
)));
2063 /* ------------------------------------------------------------------------ */
2064 /* Imported src/nmath/pchisq.c from R. */
2066 * Mathlib : A C Library of Special Functions
2067 * Copyright (C) 1998 Ross Ihaka
2068 * Copyright (C) 2000 The R Development Core Team
2070 * This program is free software; you can redistribute it and/or modify
2071 * it under the terms of the GNU General Public License as published by
2072 * the Free Software Foundation; either version 2 of the License, or
2073 * (at your option) any later version.
2075 * This program is distributed in the hope that it will be useful, but
2076 * WITHOUT ANY WARRANTY; without even the implied warranty of
2077 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2078 * General Public License for more details.
2080 * You should have received a copy of the GNU General Public License
2081 * along with this program; if not, write to the Free Software
2082 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2087 * The distribution function of the chi-squared distribution.
2091 gnm_float
pchisq(gnm_float x
, gnm_float df
, gboolean lower_tail
, gboolean log_p
)
2093 return pgamma(x
, df
/2., 2., lower_tail
, log_p
);
2096 /* ------------------------------------------------------------------------ */
2097 /* Imported src/nmath/qchisq.c from R. */
2099 * Mathlib : A C Library of Special Functions
2100 * Copyright (C) 1998 Ross Ihaka
2101 * Copyright (C) 2000 The R Development Core Team
2103 * This program is free software; you can redistribute it and/or modify
2104 * it under the terms of the GNU General Public License as published by
2105 * the Free Software Foundation; either version 2 of the License, or
2106 * (at your option) any later version.
2108 * This program is distributed in the hope that it will be useful,
2109 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2110 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2111 * GNU General Public License for more details.
2113 * You should have received a copy of the GNU General Public License
2114 * along with this program; if not, write to the Free Software
2115 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2120 * The quantile function of the chi-squared distribution.
2124 gnm_float
qchisq(gnm_float p
, gnm_float df
, gboolean lower_tail
, gboolean log_p
)
2126 return qgamma(p
, 0.5 * df
, 2.0, lower_tail
, log_p
);
2129 /* ------------------------------------------------------------------------ */
2130 /* Imported src/nmath/dweibull.c from R. */
2132 * Mathlib : A C Library of Special Functions
2133 * Copyright (C) 1998 Ross Ihaka
2134 * Copyright (C) 2000 The R Development Core Team
2136 * This program is free software; you can redistribute it and/or modify
2137 * it under the terms of the GNU General Public License as published by
2138 * the Free Software Foundation; either version 2 of the License, or
2139 * (at your option) any later version.
2141 * This program is distributed in the hope that it will be useful,
2142 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2143 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2144 * GNU General Public License for more details.
2146 * You should have received a copy of the GNU General Public License
2147 * along with this program; if not, write to the Free Software
2148 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2153 * The density function of the Weibull distribution.
2157 gnm_float
dweibull(gnm_float x
, gnm_float shape
, gnm_float scale
, gboolean give_log
)
2159 gnm_float tmp1
, tmp2
;
2161 if (gnm_isnan(x
) || gnm_isnan(shape
) || gnm_isnan(scale
))
2162 return x
+ shape
+ scale
;
2164 if (shape
<= 0 || scale
<= 0) ML_ERR_return_NAN
;
2166 if (x
< 0) return R_D__0
;
2167 if (!gnm_finite(x
)) return R_D__0
;
2168 tmp1
= gnm_pow(x
/ scale
, shape
- 1);
2169 tmp2
= tmp1
* (x
/ scale
);
2171 -tmp2
+ gnm_log(shape
* tmp1
/ scale
) :
2172 shape
* tmp1
* gnm_exp(-tmp2
) / scale
;
2175 /* ------------------------------------------------------------------------ */
2176 /* Imported src/nmath/pweibull.c from R. */
2178 * Mathlib : A C Library of Special Functions
2179 * Copyright (C) 1998 Ross Ihaka
2180 * Copyright (C) 2000-2002 The R Development Core Team
2182 * This program is free software; you can redistribute it and/or modify
2183 * it under the terms of the GNU General Public License as published by
2184 * the Free Software Foundation; either version 2 of the License, or
2185 * (at your option) any later version.
2187 * This program is distributed in the hope that it will be useful,
2188 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2189 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2190 * GNU General Public License for more details.
2192 * You should have received a copy of the GNU General Public License
2193 * along with this program; if not, write to the Free Software
2194 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2199 * The distribution function of the Weibull distribution.
2203 gnm_float
pweibull(gnm_float x
, gnm_float shape
, gnm_float scale
, gboolean lower_tail
, gboolean log_p
)
2206 if (gnm_isnan(x
) || gnm_isnan(shape
) || gnm_isnan(scale
))
2207 return x
+ shape
+ scale
;
2209 if(shape
<= 0 || scale
<= 0) ML_ERR_return_NAN
;
2213 x
= -gnm_pow(x
/ scale
, shape
);
2216 /* gnm_log(1 - gnm_exp(x)) for x < 0 : */
2217 ? (x
> -M_LN2gnum
? gnm_log(-gnm_expm1(x
)) : gnm_log1p(-gnm_exp(x
)))
2219 /* else: !lower_tail */
2223 /* ------------------------------------------------------------------------ */
2224 /* Imported src/nmath/pbinom.c from R. */
2226 * Mathlib : A C Library of Special Functions
2227 * Copyright (C) 1998 Ross Ihaka
2228 * Copyright (C) 2000, 2002 The R Development Core Team
2229 * Copyright (C) 2004 The R Foundation
2231 * This program is free software; you can redistribute it and/or modify
2232 * it under the terms of the GNU General Public License as published by
2233 * the Free Software Foundation; either version 2 of the License, or
2234 * (at your option) any later version.
2236 * This program is distributed in the hope that it will be useful,
2237 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2238 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2239 * GNU General Public License for more details.
2241 * You should have received a copy of the GNU General Public License
2242 * along with this program; if not, write to the Free Software
2243 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2248 * The distribution function of the binomial distribution.
2251 gnm_float
pbinom(gnm_float x
, gnm_float n
, gnm_float p
, gboolean lower_tail
, gboolean log_p
)
2254 if (gnm_isnan(x
) || gnm_isnan(n
) || gnm_isnan(p
))
2256 if (!gnm_finite(n
) || !gnm_finite(p
)) ML_ERR_return_NAN
;
2259 if(R_D_nonint(n
)) ML_ERR_return_NAN
;
2260 n
= R_D_forceint(n
);
2261 if(n
<= 0 || p
< 0 || p
> 1) ML_ERR_return_NAN
;
2263 x
= gnm_fake_floor(x
);
2264 if (x
< 0.0) return R_DT_0
;
2265 if (n
<= x
) return R_DT_1
;
2266 return pbeta(p
, x
+ 1, n
- x
, !lower_tail
, log_p
);
2269 /* ------------------------------------------------------------------------ */
2270 /* Imported src/nmath/dbinom.c from R. */
2273 * Catherine Loader, catherine@research.bell-labs.com.
2277 * Copyright (C) 2000, The R Core Development Team
2279 * This program is free software; you can redistribute it and/or modify
2280 * it under the terms of the GNU General Public License as published by
2281 * the Free Software Foundation; either version 2 of the License, or
2282 * (at your option) any later version.
2284 * This program is distributed in the hope that it will be useful,
2285 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2286 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2287 * GNU General Public License for more details.
2289 * You should have received a copy of the GNU General Public License
2290 * along with this program; if not, write to the Free Software
2291 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2297 * To compute the binomial probability, call dbinom(x,n,p).
2298 * This checks for argument validity, and calls dbinom_raw().
2300 * dbinom_raw() does the actual computation; note this is called by
2301 * other functions in addition to dbinom()).
2302 * (1) dbinom_raw() has both p and q arguments, when one may be represented
2303 * more accurately than the other (in particular, in df()).
2304 * (2) dbinom_raw() does NOT check that inputs x and n are integers. This
2305 * should be done in the calling function, where necessary.
2306 * (3) Also does not check for 0 <= p <= 1 and 0 <= q <= 1 or NaN's.
2307 * Do this in the calling function.
2311 static gnm_float
dbinom_raw(gnm_float x
, gnm_float n
, gnm_float p
, gnm_float q
, gboolean give_log
)
2313 gnm_float f2
, xh
, xl
, yh
, yl
;
2316 * We (ought to) have p+q = 1 with the assumption that the smaller is
2317 * more accurate that 1-other, i.e., that the bigger may already have
2318 * suffered some cancellation.
2321 if (p
== 0) return((x
== 0) ? R_D__1
: R_D__0
);
2322 if (q
== 0) return((x
== n
) ? R_D__1
: R_D__0
);
2325 /* Switch p and q to reduce code duplication. */
2333 /* Probability is p^n. */
2335 return give_log
? n
* gnm_log (p
) : gnm_pow (p
, n
);
2337 return give_log
? n
* gnm_log1p (-q
) : pow1p (-q
, n
);
2339 if (x
< 0 || x
> n
) return( R_D__0
);
2341 ebd0 (x
, n
*p
, &xh
, &xl
);
2342 PAIR_ADD(stirlerr(x
), xh
, xl
);
2344 ebd0 (n
-x
, n
*q
, &yh
, &yl
);
2345 PAIR_ADD(stirlerr(n
-x
), yh
, yl
);
2347 PAIR_ADD(yl
, xh
, xl
);
2348 PAIR_ADD(yh
, xh
, xl
);
2349 PAIR_ADD(-stirlerr(n
), xh
, xl
);
2351 f2
= (M_2PIgnum
*x
*(n
-x
))/n
;
2354 ? -xl
- xh
- 0.5 * gnm_log (f2
)
2355 : gnm_exp (-xl
) * gnm_exp (-xh
) / gnm_sqrt (f2
);
2358 gnm_float
dbinom(gnm_float x
, gnm_float n
, gnm_float p
, gboolean give_log
)
2361 /* NaNs propagated correctly */
2362 if (gnm_isnan(x
) || gnm_isnan(n
) || gnm_isnan(p
)) return x
+ n
+ p
;
2365 if (p
< 0 || p
> 1 || R_D_negInonint(n
))
2367 R_D_nonint_check(x
);
2369 n
= R_D_forceint(n
);
2370 x
= R_D_forceint(x
);
2372 return dbinom_raw(x
,n
,p
,1-p
,give_log
);
2375 /* ------------------------------------------------------------------------ */
2376 /* Imported src/nmath/qbinom.c from R. */
2378 * Mathlib : A C Library of Special Functions
2379 * Copyright (C) 1998 Ross Ihaka
2380 * Copyright (C) 2000, 2002 The R Development Core Team
2381 * Copyright (C) 2003--2004 The R Foundation
2383 * This program is free software; you can redistribute it and/or modify
2384 * it under the terms of the GNU General Public License as published by
2385 * the Free Software Foundation; either version 2 of the License, or
2386 * (at your option) any later version.
2388 * This program is distributed in the hope that it will be useful,
2389 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2390 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2391 * GNU General Public License for more details.
2393 * You should have received a copy of the GNU General Public License
2394 * along with this program; if not, write to the Free Software
2395 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2400 * The quantile function of the binomial distribution.
2404 * Uses the Cornish-Fisher Expansion to include a skewness
2405 * correction to a normal approximation. This gives an
2406 * initial value which never seems to be off by more than
2407 * 1 or 2. A search is then conducted of values close to
2408 * this initial start point.
2412 gnm_float
qbinom(gnm_float p
, gnm_float n
, gnm_float pr
, gboolean lower_tail
, gboolean log_p
)
2414 gnm_float q
, mu
, sigma
, gamma
, z
, y
;
2417 if (gnm_isnan(p
) || gnm_isnan(n
) || gnm_isnan(pr
))
2420 if(!gnm_finite(p
) || !gnm_finite(n
) || !gnm_finite(pr
))
2424 if(n
!= gnm_floor(n
+ 0.5)) ML_ERR_return_NAN
;
2425 if (pr
< 0 || pr
> 1 || n
< 0)
2428 if (pr
== 0. || n
== 0) return 0.;
2429 if (p
== R_DT_0
) return 0.;
2430 if (p
== R_DT_1
) return n
;
2433 if(q
== 0.) return n
; /* covers the full range of the distribution */
2435 sigma
= gnm_sqrt(n
* pr
* q
);
2436 gamma
= (q
- pr
) / sigma
;
2439 REprintf("qbinom(p=%7" GNM_FORMAT_g
", n=%" GNM_FORMAT_g
", pr=%7" GNM_FORMAT_g
", l.t.=%d, log=%d): sigm=%" GNM_FORMAT_g
", gam=%" GNM_FORMAT_g
"\n",
2440 p
,n
,pr
, lower_tail
, log_p
, sigma
, gamma
);
2442 /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
2443 * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
2444 if(!lower_tail
|| log_p
) {
2445 p
= R_DT_qIv(p
); /* need check again (cancellation!): */
2446 if (p
== 0.) return 0.;
2447 if (p
== 1.) return n
;
2449 /* temporary hack --- FIXME --- */
2450 if (p
+ 1.01*GNM_EPSILON
>= 1.) return n
;
2452 /* y := approx.value (Cornish-Fisher expansion) : */
2453 z
= qnorm(p
, 0., 1., /*lower_tail*/TRUE
, /*log_p*/FALSE
);
2454 y
= gnm_floor(mu
+ sigma
* (z
+ gamma
* (z
*z
- 1) / 6) + 0.5);
2455 if(y
> n
) /* way off */ y
= n
;
2458 REprintf(" new (p,1-p)=(%7" GNM_FORMAT_g
",%7" GNM_FORMAT_g
"), z=qnorm(..)=%7" GNM_FORMAT_g
", y=%5" GNM_FORMAT_g
"\n", p
, 1-p
, z
, y
);
2460 z
= pbinom(y
, n
, pr
, /*lower_tail*/TRUE
, /*log_p*/FALSE
);
2462 /* fuzz to ensure left continuity: */
2463 p
*= 1 - 64*GNM_EPSILON
;
2465 /*-- Fixme, here y can be way off --
2466 should use interval search instead of primitive stepping down or up */
2469 if((lower_tail
&& z
>= p
) || (!lower_tail
&& z
<= p
)) {
2473 /* search to the left */
2475 REprintf("\tnew z=%7" GNM_FORMAT_g
" >= p = %7" GNM_FORMAT_g
" --> search to left (y--) ..\n", z
,p
);
2479 (z
= pbinom(y
- 1, n
, pr
, /*l._t.*/TRUE
, /*log_p*/FALSE
)) < p
)
2484 else { /* search to the right */
2486 REprintf("\tnew z=%7" GNM_FORMAT_g
" < p = %7" GNM_FORMAT_g
" --> search to right (y++) ..\n", z
,p
);
2491 (z
= pbinom(y
, n
, pr
, /*l._t.*/TRUE
, /*log_p*/FALSE
)) >= p
)
2501 /* ------------------------------------------------------------------------ */
2502 /* Imported src/nmath/dnbinom.c from R. */
2505 * Catherine Loader, catherine@research.bell-labs.com.
2506 * October 23, 2000 and Feb, 2001.
2509 * Copyright (C) 2000--2001, The R Core Development Team
2511 * This program is free software; you can redistribute it and/or modify
2512 * it under the terms of the GNU General Public License as published by
2513 * the Free Software Foundation; either version 2 of the License, or
2514 * (at your option) any later version.
2516 * This program is distributed in the hope that it will be useful,
2517 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2518 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2519 * GNU General Public License for more details.
2521 * You should have received a copy of the GNU General Public License
2522 * along with this program; if not, write to the Free Software
2523 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2529 * Computes the negative binomial distribution. For integer n,
2530 * this is probability of x failures before the nth success in a
2531 * sequence of Bernoulli trials. We do not enforce integer n, since
2532 * the distribution is well defined for non-integers,
2533 * and this can be useful for e.g. overdispersed discrete survival times.
2537 gnm_float
dnbinom(gnm_float x
, gnm_float n
, gnm_float p
, gboolean give_log
)
2542 if (gnm_isnan(x
) || gnm_isnan(n
) || gnm_isnan(p
))
2546 if (p
< 0 || p
> 1 || n
<= 0) ML_ERR_return_NAN
;
2547 R_D_nonint_check(x
);
2548 if (x
< 0 || !gnm_finite(x
)) return R_D__0
;
2549 x
= R_D_forceint(x
);
2551 prob
= dbinom_raw(n
, x
+n
, p
, 1-p
, give_log
);
2552 p
= ((gnm_float
)n
)/(n
+x
);
2553 return((give_log
) ? gnm_log(p
) + prob
: p
* prob
);
2556 /* ------------------------------------------------------------------------ */
2557 /* Imported src/nmath/pnbinom.c from R. */
2559 * Mathlib : A C Library of Special Functions
2560 * Copyright (C) 1998 Ross Ihaka
2561 * Copyright (C) 2000 The R Development Core Team
2563 * This program is free software; you can redistribute it and/or modify
2564 * it under the terms of the GNU General Public License as published by
2565 * the Free Software Foundation; either version 2 of the License, or
2566 * (at your option) any later version.
2568 * This program is distributed in the hope that it will be useful,
2569 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2570 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2571 * GNU General Public License for more details.
2573 * You should have received a copy of the GNU General Public License
2574 * along with this program; if not, write to the Free Software
2575 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2580 * The distribution function of the negative binomial distribution.
2584 * x = the number of failures before the n-th success
2588 gnm_float
pnbinom(gnm_float x
, gnm_float n
, gnm_float p
, gboolean lower_tail
, gboolean log_p
)
2591 if (gnm_isnan(x
) || gnm_isnan(n
) || gnm_isnan(p
))
2593 if(!gnm_finite(n
) || !gnm_finite(p
)) ML_ERR_return_NAN
;
2595 if (n
< 0 || p
<= 0 || p
> 1) ML_ERR_return_NAN
;
2597 x
= gnm_fake_floor(x
);
2598 if (x
< 0) return R_DT_0
;
2599 if (!gnm_finite(x
)) return R_DT_1
;
2600 return pbeta(p
, n
, x
+ 1, lower_tail
, log_p
);
2603 /* ------------------------------------------------------------------------ */
2604 /* Imported src/nmath/qnbinom.c from R. */
2606 * Mathlib : A C Library of Special Functions
2607 * Copyright (C) 1998 Ross Ihaka
2608 * Copyright (C) 2000 The R Development Core Team
2610 * This program is free software; you can redistribute it and/or modify
2611 * it under the terms of the GNU General Public License as published by
2612 * the Free Software Foundation; either version 2 of the License, or
2613 * (at your option) any later version.
2615 * This program is distributed in the hope that it will be useful,
2616 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2617 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2618 * GNU General Public License for more details.
2620 * You should have received a copy of the GNU General Public License
2621 * along with this program; if not, write to the Free Software
2622 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2627 * #include <Rmath.h>
2628 * double qnbinom(double p, double n, double pr, int lower_tail, int log_p)
2632 * The quantile function of the negative binomial distribution.
2636 * x = the number of failures before the n-th success
2640 * Uses the Cornish-Fisher Expansion to include a skewness
2641 * correction to a normal approximation. This gives an
2642 * initial value which never seems to be off by more than
2643 * 1 or 2. A search is then conducted of values close to
2644 * this initial start point.
2648 gnm_float
qnbinom(gnm_float p
, gnm_float n
, gnm_float pr
, gboolean lower_tail
, gboolean log_p
)
2650 gnm_float P
, Q
, mu
, sigma
, gamma
, z
, y
;
2653 if (gnm_isnan(p
) || gnm_isnan(n
) || gnm_isnan(pr
))
2657 if (pr
<= 0 || pr
>= 1 || n
<= 0) ML_ERR_return_NAN
;
2659 if (p
== R_DT_0
) return 0;
2660 if (p
== R_DT_1
) return gnm_pinf
;
2664 sigma
= gnm_sqrt(n
* P
* Q
);
2665 gamma
= (Q
+ P
)/sigma
;
2667 /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
2668 * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
2669 if(!lower_tail
|| log_p
) {
2670 p
= R_DT_qIv(p
); /* need check again (cancellation!): */
2671 if (p
== R_DT_0
) return 0;
2672 if (p
== R_DT_1
) return gnm_pinf
;
2674 /* temporary hack --- FIXME --- */
2675 if (p
+ 1.01*GNM_EPSILON
>= 1.) return gnm_pinf
;
2677 /* y := approx.value (Cornish-Fisher expansion) : */
2678 z
= qnorm(p
, 0., 1., /*lower_tail*/TRUE
, /*log_p*/FALSE
);
2679 y
= gnm_floor(mu
+ sigma
* (z
+ gamma
* (z
*z
- 1) / 6) + 0.5);
2681 z
= pnbinom(y
, n
, pr
, /*lower_tail*/TRUE
, /*log_p*/FALSE
);
2683 /* fuzz to ensure left continuity: */
2684 p
*= 1 - 64*GNM_EPSILON
;
2686 /*-- Fixme, here y can be way off --
2687 should use interval search instead of primitive stepping down or up */
2690 if((lower_tail
&& z
>= p
) || (!lower_tail
&& z
<= p
)) {
2694 /* search to the left */
2697 (z
= pnbinom(y
- 1, n
, pr
, /*l._t.*/TRUE
, /*log_p*/FALSE
)) < p
)
2702 else { /* search to the right */
2706 if((z
= pnbinom(y
, n
, pr
, /*l._t.*/TRUE
, /*log_p*/FALSE
)) >= p
)
2715 /* ------------------------------------------------------------------------ */
2716 /* Imported src/nmath/dbeta.c from R. */
2719 * Catherine Loader, catherine@research.bell-labs.com.
2723 * Copyright (C) 2000, The R Core Development Team
2725 * This program is free software; you can redistribute it and/or modify
2726 * it under the terms of the GNU General Public License as published by
2727 * the Free Software Foundation; either version 2 of the License, or
2728 * (at your option) any later version.
2730 * This program is distributed in the hope that it will be useful,
2731 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2732 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2733 * GNU General Public License for more details.
2735 * You should have received a copy of the GNU General Public License
2736 * along with this program; if not, write to the Free Software
2737 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2744 * p(x;a,b) = ------------ x (1-x)
2747 * = (a+b-1) dbinom(a-1; a+b-2,x)
2749 * We must modify this when a<1 or b<1, to avoid passing negative
2750 * arguments to dbinom_raw. Note that the modifications require
2751 * division by x and/or 1-x, so cannot be used except where necessary.
2755 gnm_float
dbeta(gnm_float x
, gnm_float a
, gnm_float b
, gboolean give_log
)
2758 volatile gnm_float am1
, bm1
; /* prevent roundoff trouble on some
2762 /* NaNs propagated correctly */
2763 if (gnm_isnan(x
) || gnm_isnan(a
) || gnm_isnan(b
)) return x
+ a
+ b
;
2766 if (a
<= 0 || b
<= 0) ML_ERR_return_NAN
;
2767 if (x
< 0 || x
> 1) return(R_D__0
);
2769 if(a
> 1) return(R_D__0
);
2770 if(a
< 1) return(gnm_pinf
);
2771 /* a == 1 : */ return(R_D_val(b
));
2774 if(b
> 1) return(R_D__0
);
2775 if(b
< 1) return(gnm_pinf
);
2776 /* b == 1 : */ return(R_D_val(a
));
2779 if (b
< 1) { /* a,b < 1 */
2780 f
= a
*b
/((a
+b
)*x
*(1-x
));
2781 p
= dbinom_raw(a
,a
+b
, x
,1-x
, give_log
);
2783 else { /* a < 1 <= b */
2786 p
= dbinom_raw(a
,a
+bm1
, x
,1-x
, give_log
);
2790 if (b
< 1) { /* a >= 1 > b */
2793 p
= dbinom_raw(am1
,am1
+b
, x
,1-x
, give_log
);
2795 else { /* a,b >= 1 */
2799 p
= dbinom_raw(am1
,am1
+bm1
, x
,1-x
, give_log
);
2802 return( (give_log
) ? p
+ gnm_log(f
) : p
*f
);
2805 /* ------------------------------------------------------------------------ */
2806 /* Imported src/nmath/dhyper.c from R. */
2809 * Catherine Loader, catherine@research.bell-labs.com.
2813 * Copyright (C) 2000, 2001 The R Core Development Team
2815 * This program is free software; you can redistribute it and/or modify
2816 * it under the terms of the GNU General Public License as published by
2817 * the Free Software Foundation; either version 2 of the License, or
2818 * (at your option) any later version.
2820 * This program is distributed in the hope that it will be useful,
2821 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2822 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2823 * GNU General Public License for more details.
2825 * You should have received a copy of the GNU General Public License
2826 * along with this program; if not, write to the Free Software
2827 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2833 * Given a sequence of r successes and b failures, we sample n (\le b+r)
2834 * items without replacement. The hypergeometric probability is the
2835 * probability of x successes:
2837 * choose(r, x) * choose(b, n-x)
2838 * p(x; r,b,n) = ----------------------------- =
2841 * dbinom(x,r,p) * dbinom(n-x,b,p)
2842 * = --------------------------------
2845 * for any p. For numerical stability, we take p=n/(r+b); with this choice,
2846 * the denominator is not exponentially small.
2850 gnm_float
dhyper(gnm_float x
, gnm_float r
, gnm_float b
, gnm_float n
, gboolean give_log
)
2852 gnm_float p
, q
, p1
, p2
, p3
;
2855 if (gnm_isnan(x
) || gnm_isnan(r
) || gnm_isnan(b
) || gnm_isnan(n
))
2856 return x
+ r
+ b
+ n
;
2859 if (R_D_negInonint(r
) || R_D_negInonint(b
) || R_D_negInonint(n
) || n
> r
+b
)
2861 if (R_D_negInonint(x
))
2864 x
= R_D_forceint(x
);
2865 r
= R_D_forceint(r
);
2866 b
= R_D_forceint(b
);
2867 n
= R_D_forceint(n
);
2869 if (n
< x
|| r
< x
|| n
- x
> b
) return(R_D__0
);
2870 if (n
== 0) return((x
== 0) ? R_D__1
: R_D__0
);
2873 * Round to float to make both p and q numbers with few (<26) bits set.
2874 * Unless p<2^-27 that also ensures that p+q=1 without rounding errors.
2876 p
= (float)(n
/ (gnm_float
)(r
+b
));
2879 p1
= dbinom_raw(x
, r
, p
,q
,give_log
);
2880 p2
= dbinom_raw(n
-x
,b
, p
,q
,give_log
);
2881 p3
= dbinom_raw(n
,r
+b
, p
,q
,give_log
);
2883 return( (give_log
) ? p1
+ p2
- p3
: p1
*p2
/p3
);
2887 /* ------------------------------------------------------------------------ */
2888 /* Imported src/nmath/phyper.c from R. */
2890 * Mathlib : A C Library of Special Functions
2891 * Copyright (C) 1998 Ross Ihaka
2892 * Copyright (C) 1999-2000 The R Development Core Team
2893 * Copyright (C) 2004 Morten Welinder
2894 * Copyright (C) 2004 The R Foundation
2896 * This program is free software; you can redistribute it and/or modify
2897 * it under the terms of the GNU General Public License as published by
2898 * the Free Software Foundation; either version 2 of the License, or
2899 * (at your option) any later version.
2901 * This program is distributed in the hope that it will be useful,
2902 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2903 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2904 * GNU General Public License for more details.
2906 * You should have received a copy of the GNU General Public License
2907 * along with this program; if not, write to the Free Software
2908 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2913 * The distribution function of the hypergeometric distribution.
2915 * Current implementation based on posting
2916 * From: Morten Welinder <terra@gnome.org>
2917 * Cc: R-bugs@biostat.ku.dk
2918 * Subject: [Rd] phyper accuracy and efficiency (PR#6772)
2919 * Date: Thu, 15 Apr 2004 18:06:37 +0200 (CEST)
2922 The current version has very serious cancellation issues. For example,
2923 if you ask for a small right-tail you are likely to get total cancellation.
2924 For example, phyper(59, 150, 150, 60, FALSE, FALSE) gives 6.372680161e-14.
2925 The right answer is dhyper(0, 150, 150, 60, FALSE) which is 5.111204798e-22.
2927 phyper is also really slow for large arguments.
2929 Therefore, I suggest using the code below. This is a sniplet from Gnumeric ...
2930 The code isn't perfect. In fact, if x*(NR+NB) is close to n*NR,
2931 then this code can take a while. Not longer than the old code, though.
2933 -- Thanks to Ian Smith for ideas.
2937 static gnm_float
pdhyper (gnm_float x
, gnm_float NR
, gnm_float NB
, gnm_float n
, gboolean log_p
)
2942 * phyper (x, NR, NB, n, TRUE, FALSE)
2943 * [log] ----------------------------------
2944 * dhyper (x, NR, NB, n, FALSE)
2946 * without actually calling phyper. This assumes that
2948 * x * (NR + NB) <= n * NR
2954 while (x
> 0 && term
> GNM_EPSILON
* sum
) {
2955 term
*= x
* (NB
- n
+ x
) / (n
+ 1 - x
) / (NR
+ 1 - x
);
2960 return log_p
? gnm_log1p(sum
) : 1 + sum
;
2964 /* FIXME: The old phyper() code was basically used in ./qhyper.c as well
2965 * ----- We need to sync this again!
2967 gnm_float
phyper (gnm_float x
, gnm_float NR
, gnm_float NB
, gnm_float n
,
2968 gboolean lower_tail
, gboolean log_p
)
2970 /* Sample of n balls from NR red and NB black ones; x are red */
2975 if(gnm_isnan(x
) || gnm_isnan(NR
) || gnm_isnan(NB
) || gnm_isnan(n
))
2976 return x
+ NR
+ NB
+ n
;
2979 x
= gnm_fake_floor(x
);
2980 NR
= R_D_forceint(NR
);
2981 NB
= R_D_forceint(NB
);
2982 n
= R_D_forceint(n
);
2984 if (NR
< 0 || NB
< 0 || !gnm_finite(NR
+ NB
) || n
< 0 || n
> NR
+ NB
)
2987 if (x
* (NR
+ NB
) > n
* NR
) {
2989 gnm_float oldNB
= NB
;
2993 lower_tail
= !lower_tail
;
2998 /* Warning: the following line is not in R: */
3002 d
= dhyper (x
, NR
, NB
, n
, log_p
);
3003 pd
= pdhyper(x
, NR
, NB
, n
, log_p
);
3005 return log_p
? R_DT_Log(d
+ pd
) : R_D_Lval(d
* pd
);
3008 /* ------------------------------------------------------------------------ */
3009 /* Imported src/nmath/dexp.c from R. */
3011 * Mathlib : A C Library of Special Functions
3012 * Copyright (C) 1998 Ross Ihaka
3013 * Copyright (C) 2000 The R Development Core Team
3015 * This program is free software; you can redistribute it and/or modify
3016 * it under the terms of the GNU General Public License as published by
3017 * the Free Software Foundation; either version 2 of the License, or
3018 * (at your option) any later version.
3020 * This program is distributed in the hope that it will be useful,
3021 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3023 * GNU General Public License for more details.
3025 * You should have received a copy of the GNU General Public License
3026 * along with this program; if not, write to the Free Software
3027 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3032 * The density of the exponential distribution.
3036 gnm_float
dexp(gnm_float x
, gnm_float scale
, gboolean give_log
)
3039 /* NaNs propagated correctly */
3040 if (gnm_isnan(x
) || gnm_isnan(scale
)) return x
+ scale
;
3042 if (scale
<= 0.0) ML_ERR_return_NAN
;
3047 (-x
/ scale
) - gnm_log(scale
) :
3048 gnm_exp(-x
/ scale
) / scale
);
3051 /* ------------------------------------------------------------------------ */
3052 /* Imported src/nmath/pexp.c from R. */
3054 * Mathlib : A C Library of Special Functions
3055 * Copyright (C) 1998 Ross Ihaka
3056 * Copyright (C) 2000-2002 The R Development Core Team
3058 * This program is free software; you can redistribute it and/or modify
3059 * it under the terms of the GNU General Public License as published by
3060 * the Free Software Foundation; either version 2 of the License, or
3061 * (at your option) any later version.
3063 * This program is distributed in the hope that it will be useful,
3064 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3065 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3066 * GNU General Public License for more details.
3068 * You should have received a copy of the GNU General Public License
3069 * along with this program; if not, write to the Free Software
3070 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3075 * The distribution function of the exponential distribution.
3078 gnm_float
pexp(gnm_float x
, gnm_float scale
, gboolean lower_tail
, gboolean log_p
)
3081 if (gnm_isnan(x
) || gnm_isnan(scale
))
3083 if (scale
< 0) ML_ERR_return_NAN
;
3085 if (scale
<= 0) ML_ERR_return_NAN
;
3090 /* same as weibull( shape = 1): */
3094 /* gnm_log(1 - gnm_exp(x)) for x < 0 : */
3095 ? (x
> -M_LN2gnum
? gnm_log(-gnm_expm1(x
)) : gnm_log1p(-gnm_exp(x
)))
3097 /* else: !lower_tail */
3101 /* ------------------------------------------------------------------------ */
3102 /* Imported src/nmath/dgeom.c from R. */
3105 * Catherine Loader, catherine@research.bell-labs.com.
3109 * Copyright (C) 2000, 2001 The R Core Development Team
3111 * This program is free software; you can redistribute it and/or modify
3112 * it under the terms of the GNU General Public License as published by
3113 * the Free Software Foundation; either version 2 of the License, or
3114 * (at your option) any later version.
3116 * This program is distributed in the hope that it will be useful,
3117 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3118 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3119 * GNU General Public License for more details.
3121 * You should have received a copy of the GNU General Public License
3122 * along with this program; if not, write to the Free Software
3123 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3129 * Computes the geometric probabilities, Pr(X=x) = p(1-p)^x.
3133 gnm_float
dgeom(gnm_float x
, gnm_float p
, gboolean give_log
)
3138 if (gnm_isnan(x
) || gnm_isnan(p
)) return x
+ p
;
3141 if (p
< 0 || p
> 1) ML_ERR_return_NAN
;
3143 R_D_nonint_check(x
);
3144 if (x
< 0 || !gnm_finite(x
) || p
== 0) return R_D__0
;
3145 x
= R_D_forceint(x
);
3147 /* prob = (1-p)^x, stable for small p */
3148 prob
= dbinom_raw(0.,x
, p
,1-p
, give_log
);
3150 return((give_log
) ? gnm_log(p
) + prob
: p
*prob
);
3153 /* ------------------------------------------------------------------------ */
3154 /* Imported src/nmath/pgeom.c from R. */
3156 * Mathlib : A C Library of Special Functions
3157 * Copyright (C) 1998 Ross Ihaka
3158 * Copyright (C) 2000, 2001 The R Development Core Team
3159 * Copyright (C) 2004 The R Foundation
3161 * This program is free software; you can redistribute it and/or modify
3162 * it under the terms of the GNU General Public License as published by
3163 * the Free Software Foundation; either version 2 of the License, or
3164 * (at your option) any later version.
3166 * This program is distributed in the hope that it will be useful,
3167 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3168 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3169 * GNU General Public License for more details.
3171 * You should have received a copy of the GNU General Public License
3172 * along with this program; if not, write to the Free Software
3173 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3178 * The distribution function of the geometric distribution.
3182 gnm_float
pgeom(gnm_float x
, gnm_float p
, gboolean lower_tail
, gboolean log_p
)
3185 if (gnm_isnan(x
) || gnm_isnan(p
))
3188 x
= gnm_fake_floor(x
);
3189 if(p
< 0 || p
> 1) ML_ERR_return_NAN
;
3191 if (x
< 0. || p
== 0.) return R_DT_0
;
3192 if (!gnm_finite(x
)) return R_DT_1
;
3194 if(p
== 1.) { /* we cannot assume IEEE */
3195 x
= lower_tail
? 1: 0;
3196 return log_p
? gnm_log(x
) : x
;
3198 x
= gnm_log1p(-p
) * (x
+ 1);
3200 return R_DT_Clog(x
);
3202 return lower_tail
? -gnm_expm1(x
) : gnm_exp(x
);
3205 /* ------------------------------------------------------------------------ */
3206 /* Imported src/nmath/dcauchy.c from R. */
3208 * Mathlib : A C Library of Special Functions
3209 * Copyright (C) 1998 Ross Ihaka
3210 * Copyright (C) 2000 The R Development Core Team
3212 * This program is free software; you can redistribute it and/or modify
3213 * it under the terms of the GNU General Public License as published by
3214 * the Free Software Foundation; either version 2 of the License, or
3215 * (at your option) any later version.
3217 * This program is distributed in the hope that it will be useful,
3218 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3219 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3220 * GNU General Public License for more details.
3222 * You should have received a copy of the GNU General Public License
3223 * along with this program; if not, write to the Free Software
3224 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
3229 * The density of the Cauchy distribution.
3233 gnm_float
dcauchy(gnm_float x
, gnm_float location
, gnm_float scale
, gboolean give_log
)
3237 /* NaNs propagated correctly */
3238 if (gnm_isnan(x
) || gnm_isnan(location
) || gnm_isnan(scale
))
3239 return x
+ location
+ scale
;
3241 if (scale
<= 0) ML_ERR_return_NAN
;
3243 y
= (x
- location
) / scale
;
3245 - gnm_log(M_PIgnum
* scale
* (1. + y
* y
)) :
3246 1. / (M_PIgnum
* scale
* (1. + y
* y
));
3249 /* ------------------------------------------------------------------------ */
3250 /* Imported src/nmath/pcauchy.c from R. */
3252 * Mathlib : A C Library of Special Functions
3253 * Copyright (C) 1998 Ross Ihaka
3254 * Copyright (C) 2000 The R Development Core Team
3255 * Copyright (C) 2004 The R Foundation
3257 * This program is free software; you can redistribute it and/or modify
3258 * it under the terms of the GNU General Public License as published by
3259 * the Free Software Foundation; either version 2 of the License, or
3260 * (at your option) any later version.
3262 * This program is distributed in the hope that it will be useful,
3263 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3264 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3265 * GNU General Public License for more details.
3267 * You should have received a copy of the GNU General Public License
3268 * along with this program; if not, write to the Free Software
3269 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3274 * The distribution function of the Cauchy distribution.
3278 gnm_float
pcauchy(gnm_float x
, gnm_float location
, gnm_float scale
,
3279 gboolean lower_tail
, gboolean log_p
)
3282 if (gnm_isnan(x
) || gnm_isnan(location
) || gnm_isnan(scale
))
3283 return x
+ location
+ scale
;
3285 if (scale
<= 0) ML_ERR_return_NAN
;
3287 x
= (x
- location
) / scale
;
3288 if (gnm_isnan(x
)) ML_ERR_return_NAN
;
3290 if(!gnm_finite(x
)) {
3291 if(x
< 0) return R_DT_0
;
3299 return R_D_Clog(gnm_atan2pi (1, x
));
3301 return R_D_val(gnm_atan2pi (1, -x
));
3304 /* ------------------------------------------------------------------------ */
3305 /* Imported src/nmath/df.c from R. */
3308 * Catherine Loader, catherine@research.bell-labs.com.
3312 * Copyright (C) 2000, The R Core Development Team
3314 * This program is free software; you can redistribute it and/or modify
3315 * it under the terms of the GNU General Public License as published by
3316 * the Free Software Foundation; either version 2 of the License, or
3317 * (at your option) any later version.
3319 * This program is distributed in the hope that it will be useful,
3320 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3321 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3322 * GNU General Public License for more details.
3324 * You should have received a copy of the GNU General Public License
3325 * along with this program; if not, write to the Free Software
3326 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3332 * The density function of the F distribution.
3333 * To evaluate it, write it as a Binomial probability with p = x*m/(n+x*m).
3334 * For m >= 2, we use the simplest conversion.
3335 * For m < 2, (m-2)/2 < 0 so the conversion will not work, and we must use
3336 * a second conversion.
3337 * Note the division by p; this seems unavoidable
3338 * for m < 2, since the F density has a singularity as x (or p) -> 0.
3342 gnm_float
df(gnm_float x
, gnm_float m
, gnm_float n
, gboolean give_log
)
3344 gnm_float p
, q
, f
, dens
;
3347 if (gnm_isnan(x
) || gnm_isnan(m
) || gnm_isnan(n
))
3350 if (m
<= 0 || n
<= 0) ML_ERR_return_NAN
;
3351 if (x
<= 0.) return(R_D__0
);
3359 dens
= dbinom_raw((m
-2)/2, (m
+n
-2)/2, p
, q
, give_log
);
3362 f
= m
*m
*q
/ (2*p
*(m
+n
));
3363 dens
= dbinom_raw(m
/2, (m
+n
)/2, p
, q
, give_log
);
3365 return(give_log
? gnm_log(f
)+dens
: f
*dens
);
3368 /* ------------------------------------------------------------------------ */
3369 /* Imported src/nmath/dchisq.c from R. */
3371 * Mathlib : A C Library of Special Functions
3372 * Copyright (C) 1998 Ross Ihaka
3373 * Copyright (C) 2000 The R Development Core Team
3375 * This program is free software; you can redistribute it and/or modify
3376 * it under the terms of the GNU General Public License as published by
3377 * the Free Software Foundation; either version 2 of the License, or
3378 * (at your option) any later version.
3380 * This program is distributed in the hope that it will be useful, but
3381 * WITHOUT ANY WARRANTY; without even the implied warranty of
3382 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3383 * General Public License for more details.
3385 * You should have received a copy of the GNU General Public License
3386 * along with this program; if not, write to the Free Software
3387 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3392 * The density of the chi-squared distribution.
3396 gnm_float
dchisq(gnm_float x
, gnm_float df
, gboolean give_log
)
3398 return dgamma(x
, df
/ 2., 2., give_log
);
3401 /* ------------------------------------------------------------------------ */
3402 /* Imported src/nmath/qweibull.c from R. */
3404 * Mathlib : A C Library of Special Functions
3405 * Copyright (C) 1998 Ross Ihaka
3406 * Copyright (C) 2000 The R Development Core Team
3408 * This program is free software; you can redistribute it and/or modify
3409 * it under the terms of the GNU General Public License as published by
3410 * the Free Software Foundation; either version 2 of the License, or
3411 * (at your option) any later version.
3413 * This program is distributed in the hope that it will be useful,
3414 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3415 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3416 * GNU General Public License for more details.
3418 * You should have received a copy of the GNU General Public License
3419 * along with this program; if not, write to the Free Software
3420 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3425 * The quantile function of the Weibull distribution.
3429 gnm_float
qweibull(gnm_float p
, gnm_float shape
, gnm_float scale
, gboolean lower_tail
, gboolean log_p
)
3432 if (gnm_isnan(p
) || gnm_isnan(shape
) || gnm_isnan(scale
))
3433 return p
+ shape
+ scale
;
3436 if (shape
<= 0 || scale
<= 0) ML_ERR_return_NAN
;
3438 if (p
== R_D__0
) return 0;
3439 if (p
== R_D__1
) return gnm_pinf
;
3440 return scale
* gnm_pow(- R_DT_Clog(p
), 1./shape
) ;
3443 /* ------------------------------------------------------------------------ */
3444 /* Imported src/nmath/qexp.c from R. */
3446 * Mathlib : A C Library of Special Functions
3447 * Copyright (C) 1998 Ross Ihaka
3448 * Copyright (C) 2000 The R Development Core Team
3450 * This program is free software; you can redistribute it and/or modify
3451 * it under the terms of the GNU General Public License as published by
3452 * the Free Software Foundation; either version 2 of the License, or
3453 * (at your option) any later version.
3455 * This program is distributed in the hope that it will be useful,
3456 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3457 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3458 * GNU General Public License for more details.
3460 * You should have received a copy of the GNU General Public License
3461 * along with this program; if not, write to the Free Software
3462 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3467 * The quantile function of the exponential distribution.
3472 gnm_float
qexp(gnm_float p
, gnm_float scale
, gboolean lower_tail
, gboolean log_p
)
3475 if (gnm_isnan(p
) || gnm_isnan(scale
))
3479 if (scale
< 0) ML_ERR_return_NAN
;
3484 return - scale
* R_DT_Clog(p
);
3487 /* ------------------------------------------------------------------------ */
3488 /* Imported src/nmath/qgeom.c from R. */
3490 * Mathlib : A C Library of Special Functions
3491 * Copyright (C) 1998 Ross Ihaka
3492 * Copyright (C) 2000 The R Development Core Team
3493 * Copyright (C) 2004 The R Foundation
3495 * This program is free software; you can redistribute it and/or modify
3496 * it under the terms of the GNU General Public License as published by
3497 * the Free Software Foundation; either version 2 of the License, or
3498 * (at your option) any later version.
3500 * This program is distributed in the hope that it will be useful,
3501 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3502 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3503 * GNU General Public License for more details.
3505 * You should have received a copy of the GNU General Public License
3506 * along with this program; if not, write to the Free Software
3507 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3512 * The quantile function of the geometric distribution.
3516 gnm_float
qgeom(gnm_float p
, gnm_float prob
, gboolean lower_tail
, gboolean log_p
)
3521 if (prob
<= 0 || prob
> 1) ML_ERR_return_NAN
;
3524 if (gnm_isnan(p
) || gnm_isnan(prob
))
3526 if (p
== R_DT_1
) return gnm_pinf
;
3528 if (p
== R_DT_1
) ML_ERR_return_NAN
;
3530 if (p
== R_DT_0
) return 0;
3532 /* add a fuzz to ensure left continuity */
3533 q
= gnm_ceil(R_DT_Clog(p
) / gnm_log1p(- prob
) - 1 - 1e-12);
3534 return MAX (q
, 0.0);
3537 /* ------------------------------------------------------------------------ */
3539 * Based on code imported from R by hand. Heavily modified to enhance
3540 * accuracy. See bug 700132.
3543 * Mathlib : A C Library of Special Functions
3544 * Copyright (C) 1998 Ross Ihaka
3545 * Copyright (C) 2000--2007 The R Core Team
3547 * This program is free software; you can redistribute it and/or modify
3548 * it under the terms of the GNU General Public License as published by
3549 * the Free Software Foundation; either version 2 of the License, or
3550 * (at your option) any later version.
3552 * This program is distributed in the hope that it will be useful,
3553 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3554 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3555 * GNU General Public License for more details.
3557 * You should have received a copy of the GNU General Public License
3558 * along with this program; if not, a copy is available at
3559 * http://www.r-project.org/Licenses/
3563 * #include <Rmath.h>
3564 * double ptukey(q, rr, cc, df, lower_tail, log_p);
3568 * Computes the probability that the maximum of rr studentized
3569 * ranges, each based on cc means and with df degrees of freedom
3570 * for the standard error, is less than q.
3572 * The algorithm is based on that of the reference.
3576 * Copenhaver, Margaret Diponzio & Holland, Burt S.
3577 * Multiple comparisons of simple effects in
3578 * the two-way analysis of variance with fixed effects.
3579 * Journal of Statistical Computation and Simulation,
3580 * Vol.30, pp.1-15, 1988.
3583 static gnm_float
ptukey_wprob(gnm_float w
, gnm_float rr
, gnm_float cc
)
3587 This function calculates probability integral of Hartley's
3591 rr = no. of rows or groups
3592 cc = no. of columns or treatments
3593 pr_w = returned probability integral
3595 program will not terminate if ir is raised.
3597 bb = upper limit of legendre integration
3598 nleg = order of legendre quadrature
3599 wlar = value of range above which wincr1 intervals are used to
3600 calculate second part of integral,
3601 else wincr2 intervals are used.
3602 or modifying a calculation.
3604 M_1_SQRT_2PI = 1 / sqrt(2 * pi); from abramowitz & stegun, p. 3.
3605 M_SQRT2gnum = sqrt(2)
3606 xleg = legendre 12-point nodes
3607 aleg = legendre 12-point coefficients
3610 static const gboolean debug
= FALSE
;
3611 static const gnm_float xleg
[] = {
3612 GNM_const(0.981560634246719250690549090149),
3613 GNM_const(0.904117256370474856678465866119),
3614 GNM_const(0.769902674194304687036893833213),
3615 GNM_const(0.587317954286617447296702418941),
3616 GNM_const(0.367831498998180193752691536644),
3617 GNM_const(0.125233408511468915472441369464)
3619 static const gnm_float aleg
[G_N_ELEMENTS (xleg
)] = {
3620 GNM_const(0.047175336386511827194615961485),
3621 GNM_const(0.106939325995318430960254718194),
3622 GNM_const(0.160078328543346226334652529543),
3623 GNM_const(0.203167426723065921749064455810),
3624 GNM_const(0.233492536538354808760849898925),
3625 GNM_const(0.249147045813402785000562436043)
3627 const int nleg
= G_N_ELEMENTS (xleg
) * 2;
3628 gnm_float pr_w
, binc
, qsqz
, blb
;
3633 /* find (2F(w/2) - 1) ^ cc */
3634 /* (first term in integral of hartley's form). */
3637 * Use different formulas for different size of qsqz to avoid
3638 * cancellation for pnorm or squeezing erf's result against 1.
3641 ? pow1p (-2.0 * pnorm (qsqz
, 0, 1, FALSE
, FALSE
), cc
)
3642 : gnm_pow (gnm_erf (qsqz
/ M_SQRT2gnum
), cc
);
3646 /* find the integral of second term of hartley's form */
3647 /* Limits of integration are (w/2, Inf). Large cc means that */
3648 /* that we need smaller step size. The formula for binc is */
3651 /* blb and bub are lower and upper limits of integration. */
3654 binc
= 3.0 / gnm_log1p (cc
);
3656 /* integrate over each interval */
3658 for (i
= 1; TRUE
; i
++) {
3659 gnm_float C
= blb
+ binc
* 0.5;
3660 gnm_float elsum
= 0.0;
3662 /* legendre quadrature with order = nleg */
3664 for (jj
= 0; jj
< nleg
; jj
++) {
3665 gnm_float xx
, aa
, v
, rinsum
;
3667 if (nleg
/ 2 <= jj
) {
3668 xx
= xleg
[nleg
- 1 - jj
];
3669 aa
= aleg
[nleg
- 1 - jj
];
3674 v
= C
+ binc
* 0.5 * xx
;
3676 rinsum
= pnorm2 (v
- w
, v
);
3677 elsum
+= gnm_pow (rinsum
, cc
- 1) * aa
* expmx2h(v
);
3679 elsum
*= (binc
* cc
* M_1_SQRT_2PI
);
3683 /* Nothing more will contribute. */
3688 if (elsum
<= pr_w
* (GNM_EPSILON
/ 2)) {
3689 /* This interval contributed nothing. We're done. */
3691 g_printerr ("End at i=%d for w=%g cc=%g dP=%g P=%g\n",
3692 i
, w
, cc
, elsum
, pr_w
);
3700 if (0) g_printerr ("w=%g pr_w=%.20g\n", w
, pr_w
);
3702 return gnm_pow (pr_w
, rr
);
3706 ptukey_otsum (gnm_float u0
, gnm_float u1
, gnm_float f2
, gnm_float f2lf
,
3707 gnm_float q
, gnm_float rr
, gnm_float cc
)
3709 gboolean debug
= FALSE
;
3710 static const gnm_float xlegq
[] = {
3711 GNM_const(0.989400934991649932596154173450),
3712 GNM_const(0.944575023073232576077988415535),
3713 GNM_const(0.865631202387831743880467897712),
3714 GNM_const(0.755404408355003033895101194847),
3715 GNM_const(0.617876244402643748446671764049),
3716 GNM_const(0.458016777657227386342419442984),
3717 GNM_const(0.281603550779258913230460501460),
3718 GNM_const(0.950125098376374401853193354250e-1)
3720 static const gnm_float alegq
[G_N_ELEMENTS (xlegq
)] = {
3721 GNM_const(0.271524594117540948517805724560e-1),
3722 GNM_const(0.622535239386478928628438369944e-1),
3723 GNM_const(0.951585116824927848099251076022e-1),
3724 GNM_const(0.124628971255533872052476282192),
3725 GNM_const(0.149595988816576732081501730547),
3726 GNM_const(0.169156519395002538189312079030),
3727 GNM_const(0.182603415044923588866763667969),
3728 GNM_const(0.189450610455068496285396723208)
3730 const int nlegq
= G_N_ELEMENTS (xlegq
) * 2;
3732 gnm_float C
= 0.5 * (u0
+ u1
);
3733 gnm_float L
= u1
- u0
;
3734 gnm_float otsum
= 0.0;
3736 for (jj
= 0; jj
< nlegq
; jj
++) {
3737 gnm_float xx
, aa
, u
, t1
, wprb
;
3739 if (nlegq
/ 2 <= jj
) {
3740 xx
= xlegq
[nlegq
- 1 - jj
];
3741 aa
= alegq
[nlegq
- 1 - jj
];
3747 u
= xx
* (0.5 * L
) + C
;
3749 t1
= f2lf
+ (f2
- 1) * gnm_log(u
) - u
* f2
;
3751 wprb
= ptukey_wprob(q
* gnm_sqrt(u
), rr
, cc
);
3752 otsum
+= aa
* (wprb
* (0.5 * L
) * gnm_exp(t1
));
3756 g_printerr ("Integral over [%g;%g] was %g\n",
3765 R_ptukey(gnm_float q
, gnm_float rr
, gnm_float cc
, gnm_float df
,
3766 gboolean lower_tail
, gboolean log_p
)
3768 /* function ptukey() [was qprob() ]:
3770 q = value of studentized range
3771 rr = no. of rows or groups
3772 cc = no. of columns or treatments
3773 df = degrees of freedom of error term
3774 ir[0] = error flag = 1 if wprob probability > 1
3775 ir[1] = error flag = 1 if qprob probability > 1
3777 qprob = returned probability integral over [0, q]
3779 The program will not terminate if ir[0] or ir[1] are raised.
3781 All references in wprob to Abramowitz and Stegun
3782 are from the following reference:
3784 Abramowitz, Milton and Stegun, Irene A.
3785 Handbook of Mathematical Functions.
3786 New York: Dover publications, Inc. (1970).
3788 All constants taken from this text are
3789 given to 25 significant digits.
3791 nlegq = order of legendre quadrature
3792 eps = max. allowable value of integral
3793 eps1 & eps2 = values below which there is
3794 no contribution to integral.
3796 d.f. <= dhaf: integral is divided into ulen1 length intervals. else
3797 d.f. <= dquar: integral is divided into ulen2 length intervals. else
3798 d.f. <= deigh: integral is divided into ulen3 length intervals. else
3799 d.f. <= dlarg: integral is divided into ulen4 length intervals.
3801 d.f. > dlarg: the range is used to calculate integral.
3803 xlegq = legendre 16-point nodes
3804 alegq = legendre 16-point coefficients
3806 The coefficients and nodes for the legendre quadrature used in
3807 qprob and wprob were calculated using the algorithms found in:
3809 Stroud, A. H. and Secrest, D.
3810 Gaussian Quadrature Formulas.
3812 New Jersey: Prentice-Hall, Inc, 1966.
3814 All values matched the tables (provided in same reference)
3815 to 30 significant digits.
3817 F(x) = .5 + erf(x / sqrt(2)) / 2 for x > 0
3819 F(x) = erfc( -x / sqrt(2)) / 2 for x < 0
3821 where F(x) is standard normal c. d. f.
3823 if degrees of freedom large, approximate integral
3824 with range distribution.
3827 static const gnm_float dhaf
= 100.0;
3828 static const gnm_float dquar
= 800.0;
3829 static const gnm_float deigh
= 5000.0;
3830 static const gnm_float dlarg
= 25000.0;
3831 static const gnm_float ulen1
= 1.0;
3832 static const gnm_float ulen2
= 0.5;
3833 static const gnm_float ulen3
= 0.25;
3834 static const gnm_float ulen4
= 0.125;
3835 gnm_float ans
, f2
, f2lf
, ulen
, u0
, u1
;
3837 gboolean fail
= FALSE
;
3838 gboolean debug
= TRUE
;
3841 if (gnm_isnan(q
) || gnm_isnan(rr
) || gnm_isnan(cc
) || gnm_isnan(df
))
3848 /* FIXME: Special case for cc==2&&rr=1: we have explicit formula */
3851 /* df must be > 1 */
3852 /* there must be at least two values */
3854 if (df
< 2 || rr
< 1 || cc
< 2) ML_ERR_return_NAN
;
3860 return R_DT_val(ptukey_wprob(q
, rr
, cc
));
3862 /* calculate leading constant */
3865 /* gnm_lgamma(u) = log(gamma(u)) */
3866 f2lf
= f2
* gnm_log(f2
) - gnm_lgamma(f2
);
3868 /* integral is divided into unit, half-unit, quarter-unit, or */
3869 /* eighth-unit length intervals depending on the value of the */
3870 /* degrees of freedom. */
3872 if (df
<= dhaf
) ulen
= ulen1
;
3873 else if (df
<= dquar
) ulen
= ulen2
;
3874 else if (df
<= deigh
) ulen
= ulen3
;
3877 /* integrate over each subinterval */
3881 * Integration for [0;ulen/2].
3883 * We use a sequence of intervals each covering an ever increasing
3884 * fraction of what is left. Breaking the interval takes care of
3885 * some very un-polynomial behaviour of the integrand:
3886 * (1) Infinite derivative at 0 for low df.
3887 * (2) Infinite roots (due to underflow) in the left part of an
3888 * interval with meaningful contributions from its right part
3891 for (i
= 1; TRUE
; i
++) {
3895 otsum
= ptukey_otsum (u0
, u1
, f2
, f2lf
, q
, rr
, cc
);
3898 if (otsum
<= ans
* (GNM_EPSILON
/ 2)) {
3899 /* This interval contributed nothing. We're done. */
3905 g_printerr ("PTUKEY FAIL LEFT: %d q=%g cc=%g df=%g otsum=%g ans=%g\n",
3906 i
, q
, cc
, df
, otsum
, ans
);
3915 * Integration for [ulen/2;Infinity].
3917 * We use a sequence of intervals starting with length ulen, but
3918 * when we think we're in the tail we ramp up the length.
3921 for (i
= 1; TRUE
; i
++) {
3925 otsum
= ptukey_otsum (u0
, u1
, f2
, f2lf
, q
, rr
, cc
);
3928 if (otsum
< ans
* GNM_EPSILON
) {
3930 * This interval contributed nothing. This can either be
3931 * because the integrand is so flat that we haven't seen
3932 * anyting yet or because we're done. Or both.
3934 * The maximum of the integrand falls not far from
3936 * but let's go a little further.
3938 if (ans
> 0 || u0
> 2) {
3945 g_printerr ("PTUKEY FAIL RIGHT: %i %g %g\n",
3952 * When we see a contribution that added less than 0.1%
3953 * to the result, start ramping up the interval size.
3954 * Note that this will not trigger unless ans>0.
3956 if (otsum
< ans
/ 1000)
3963 ML_ERROR(ME_PRECISION
);
3964 ans
= MIN (ans
, 1.0);
3965 return R_DT_val(ans
);
3968 /* ------------------------------------------------------------------------ */
3969 /* --- END MAGIC R SOURCE MARKER --- */
3971 /* Silly order-of-arguments wrapper. */
3973 ptukey(gnm_float x
, gnm_float nmeans
, gnm_float df
, gnm_float nranges
, gboolean lower_tail
, gboolean log_p
)
3975 return R_ptukey (x
, nranges
, nmeans
, df
, lower_tail
, log_p
);
3979 ptukey1 (gnm_float x
, const gnm_float shape
[],
3980 gboolean lower_tail
, gboolean log_p
)
3982 return ptukey (x
, shape
[0], shape
[1], shape
[2], lower_tail
, log_p
);
3986 qtukey (gnm_float p
, gnm_float nmeans
, gnm_float df
, gnm_float nranges
, gboolean lower_tail
, gboolean log_p
)
3988 gnm_float x0
, shape
[3];
3990 if (!log_p
&& p
> 0.9) {
3991 /* We're far into the tail. Flip. */
3993 lower_tail
= !lower_tail
;
3996 /* This is accurate for nmeans==2 and nranges==1. */
3997 x0
= M_SQRT2gnum
* qt ((1 + p
) / 2, df
, lower_tail
, log_p
);
4003 return pfuncinverter (p
, shape
, lower_tail
, log_p
, 0, gnm_pinf
, x0
, ptukey1
, NULL
);
4007 logspace_signed_add (gnm_float logx
, gnm_float logabsy
, gboolean ypos
)
4010 ? logspace_add (logx
, logabsy
)
4011 : logspace_sub (logx
, logabsy
);
4014 /* Accurate gnm_log (1 - gnm_exp (p)) for p <= 0. */
4016 swap_log_tail (gnm_float lp
)
4018 if (lp
> -1 / gnm_log (2))
4019 return gnm_log (-gnm_expm1 (lp
)); /* Good formula for lp near zero. */
4021 return gnm_log1p (-gnm_exp (lp
)); /* Good formula for small lp. */
4025 /* --- BEGIN IANDJMSMITH SOURCE MARKER --- */
4027 /* Calculation of logfbit(x)-logfbit(1+x). y2 must be < 1. */
4029 logfbitdif (gnm_float x
)
4031 gnm_float y
= 1 / (2 * x
+ 3);
4032 gnm_float y2
= y
* y
;
4033 return y2
* gnm_logcf (y2
, 3, 2);
4037 * lfbc{1-7} from this Mathematica program:
4039 * Table[Numerator[BernoulliB[2n]/(2n(2n - 1))], {n, 1, 22}]
4040 * Table[Denominator[BernoulliB[2n]/(2n(2n - 1))], {n, 1, 22}]
4042 static const gnm_float lfbc1
= GNM_const (1.0) / 12;
4043 static const gnm_float lfbc2
= GNM_const (1.0) / 30;
4044 static const gnm_float lfbc3
= GNM_const (1.0) / 105;
4045 static const gnm_float lfbc4
= GNM_const (1.0) / 140;
4046 static const gnm_float lfbc5
= GNM_const (1.0) / 99;
4047 static const gnm_float lfbc6
= GNM_const (691.0) / 30030;
4048 static const gnm_float lfbc7
= GNM_const (1.0) / 13;
4049 /* lfbc{8,9} to make logfbit(6) and logfbit(7) exact. */
4050 static const gnm_float lfbc8
= GNM_const (3.5068606896459316479e-01);
4051 static const gnm_float lfbc9
= GNM_const (1.6769998201671114808);
4053 /* This is also stirlerr(x+1). */
4055 logfbit (gnm_float x
)
4058 * Error part of Stirling's formula where
4059 * log(x!) = log(sqrt(twopi))+(x+0.5)*log(x+1)-(x+1)+logfbit(x).
4061 * Are we ever concerned about the relative error involved in this
4062 * function? I don't think so.
4064 if (x
>= 1e10
) return 1 / (12 * (x
+ 1));
4066 gnm_float x1
= x
+ 1;
4067 gnm_float x2
= 1 / (x1
* x1
);
4077 return lfbc1
* (1 - x3
) / x1
;
4079 else if (x
== 5) return GNM_const (0.13876128823070747998745727023762908562e-1);
4080 else if (x
== 4) return GNM_const (0.16644691189821192163194865373593391145e-1);
4081 else if (x
== 3) return GNM_const (0.20790672103765093111522771767848656333e-1);
4082 else if (x
== 2) return GNM_const (0.27677925684998339148789292746244666596e-1);
4083 else if (x
== 1) return GNM_const (0.41340695955409294093822081407117508025e-1);
4084 else if (x
== 0) return GNM_const (0.81061466795327258219670263594382360138e-1);
4089 x2
+= logfbitdif (x1
);
4092 return x2
+ logfbit (x1
);
4094 else return gnm_pinf
;
4097 /* Calculation of logfbit1(x)-logfbit1(1+x). */
4099 logfbit1dif (gnm_float x
)
4101 return (logfbitdif (x
) - 1 / (4 * (x
+ 1) * (x
+ 2))) / (x
+ 1.5);
4104 /* Derivative logfbit. */
4106 logfbit1 (gnm_float x
)
4108 if (x
>= 1e10
) return -lfbc1
* gnm_pow (x
+ 1, -2);
4110 gnm_float x1
= x
+ 1;
4111 gnm_float x2
= 1 / (x1
* x1
);
4113 x2
* (3 * lfbc2
- x2
*
4121 return -lfbc1
* (1 - x3
) * x2
;
4127 x2
+= logfbit1dif (x1
);
4130 return x2
+ logfbit1 (x1
);
4132 else return gnm_ninf
;
4136 /* Calculation of logfbit3(x)-logfbit3(1+x). */
4138 logfbit3dif (gnm_float x
)
4140 return -(2 * x
+ 3) * gnm_pow ((x
+ 1) * (x
+ 2), -3);
4144 /* Third derivative logfbit. */
4146 logfbit3 (gnm_float x
)
4148 if (x
>= 1e10
) return -0.5 * gnm_pow (x
+ 1, -4);
4150 gnm_float x1
= x
+ 1;
4151 gnm_float x2
= 1 / (x1
* x1
);
4153 x2
* (60 * lfbc2
- x2
*
4157 (1716 * lfbc6
- x2
*
4158 (2730 * lfbc7
- x2
*
4159 (4080 * lfbc8
- x2
*
4160 5814 * lfbc9
)))))));
4161 return -lfbc1
* (6 - x3
) * x2
* x2
;
4167 x2
+= logfbit3dif (x1
);
4170 return x2
+ logfbit3 (x1
);
4172 else return gnm_ninf
;
4175 /* Calculation of logfbit5(x)-logfbit5(1+x). */
4177 logfbit5dif (gnm_float x
)
4179 return -6 * (2 * x
+ 3) * ((5 * x
+ 15) * x
+ 12) *
4180 gnm_pow ((x
+ 1) * (x
+ 2), -5);
4183 /* Fifth derivative logfbit. */
4185 logfbit5 (gnm_float x
)
4187 if (x
>= 1e10
) return -10 * gnm_pow (x
+ 1, -6);
4189 gnm_float x1
= x
+ 1;
4190 gnm_float x2
= 1 / (x1
* x1
);
4192 x2
* (2520 * lfbc2
- x2
*
4193 (15120 * lfbc3
- x2
*
4194 (55440 * lfbc4
- x2
*
4195 (154440 * lfbc5
- x2
*
4196 (360360 * lfbc6
- x2
*
4197 (742560 * lfbc7
- x2
*
4198 (1395360 * lfbc8
- x2
*
4199 2441880 * lfbc9
)))))));
4200 return -lfbc1
* (120 - x3
) * x2
* x2
* x2
;
4201 } else if (x
> -1) {
4205 x2
+= logfbit5dif (x1
);
4208 return x2
+ logfbit5 (x1
);
4210 else return gnm_ninf
;
4213 /* Calculation of logfbit7(x)-logfbit7(1+x). */
4215 logfbit7dif (gnm_float x
)
4219 ((((14 * x
+ 84) * x
+ 196) * x
+ 210) * x
+ 87) *
4220 gnm_pow ((x
+ 1) * (x
+ 2), -7);
4223 /* Seventh derivative logfbit. */
4225 logfbit7 (gnm_float x
)
4227 if (x
>= 1e10
) return -420 * gnm_pow (x
+ 1, -8);
4229 gnm_float x1
= x
+ 1;
4230 gnm_float x2
= 1 / (x1
* x1
);
4232 x2
* (181440 * lfbc2
- x2
*
4233 (1663200 * lfbc3
- x2
*
4234 (8648640 * lfbc4
- x2
*
4235 (32432400 * lfbc5
- x2
*
4236 (98017920 * lfbc6
- x2
*
4237 (253955520 * lfbc7
- x2
*
4238 (586051200 * lfbc8
- x2
*
4239 1235591280 * lfbc9
)))))));
4240 return -lfbc1
* (5040 - x3
) * x2
* x2
* x2
* x2
;
4241 } else if (x
> -1) {
4245 x2
+= logfbit7dif (x1
);
4248 return x2
+ logfbit7 (x1
);
4250 else return gnm_ninf
;
4255 lfbaccdif (gnm_float a
, gnm_float b
)
4257 if (a
> 0.03 * (a
+ b
))
4258 return logfbit (a
+ b
) - logfbit (b
);
4260 gnm_float a2
= a
* a
;
4261 gnm_float ab2
= a
/ 2 + b
;
4262 return a
* (logfbit1 (ab2
) + a2
/ 24 *
4263 (logfbit3 (ab2
) + a2
/ 80 *
4264 (logfbit5 (ab2
) + a2
/ 168 *
4270 compbfunc (gnm_float x
, gnm_float a
, gnm_float b
)
4272 const gnm_float sumAcc
= 5E-16;
4275 gnm_float sum
= term
/ (a
+ 1);
4276 while (gnm_abs (term
) > gnm_abs (sum
* sumAcc
)) {
4277 term
*= (d
- b
) * x
/ d
;
4278 sum
+= term
/ (a
+ d
);
4281 return a
* (b
- 1) * sum
;
4285 pbeta_smalla (gnm_float x
, gnm_float a
, gnm_float b
, gboolean lower_tail
, gboolean log_p
)
4290 assert (a
>= 0 && b
>= 0);
4292 assert (b
< 1 || (1 + b
) * x
<= 1);
4300 lower_tail
= !lower_tail
;
4303 r
= (a
+ b
+ 0.5) * log1pmx (a
/ (1 + b
)) +
4304 a
* (a
- 0.5) / (1 + b
) +
4306 r
+= a
* gnm_log ((1 + b
) * x
) - lgamma1p (a
);
4309 return r
+ gnm_log1p (-compbfunc (x
, a
, b
)) + gnm_log (b
/ (a
+ b
));
4311 return gnm_exp (r
) * (1 - compbfunc (x
, a
, b
)) * (b
/ (a
+ b
));
4313 /* x=0.500000001 a=0.5 b=0.000001 ends up here [swapped]
4314 * with r=-7.94418987455065e-08 and cbf=-3.16694087508444e-07.
4316 * x=0.0000001 a=0.999999 b=0.02 end up here with
4317 * r=-16.098276918385 and cbf=-4.89999787339858e-08.
4320 return swap_log_tail (r
+ gnm_log1p (-compbfunc (x
, a
, b
)) + gnm_log (b
/ (a
+ b
)));
4323 r
+= compbfunc (x
, a
, b
) * (1 - r
);
4324 r
+= (a
/ (a
+ b
)) * (1 - r
);
4330 /* Cumulative Students t-distribution, with odd parameterisation for
4331 * use with binApprox.
4334 * logqk2 is LN(q)*k/2
4335 * approxtdistDens returns with approx density function, for use in
4339 tdistexp (gnm_float p
, gnm_float q
, gnm_float logqk2
, gnm_float k
,
4340 gboolean log_p
, gnm_float
*approxtdistDens
)
4342 const gnm_float sumAcc
= 5E-16;
4343 const gnm_float cfVSmall
= 1.0e-14;
4344 const gnm_float lstpi
= gnm_log (2 * M_PIgnum
) / 2;
4346 if (gnm_floor (k
/ 2) * 2 == k
) {
4347 gnm_float ldens
= logqk2
+ logfbit (k
- 1) - 2 * logfbit (k
* 0.5 - 1) - lstpi
;
4348 *approxtdistDens
= R_D_exp (ldens
);
4350 gnm_float ldens
= logqk2
+ k
* log1pmx (1 / k
) + 2 * logfbit ((k
- 1) * 0.5) - logfbit (k
- 1) - lstpi
;
4351 *approxtdistDens
= R_D_exp (ldens
);
4354 if (k
* p
< 4 * q
) {
4356 gnm_float aki
= k
+ 1;
4358 gnm_float term
= aki
* p
/ ai
;
4360 while (term
> sumAcc
* sum
) {
4364 term
*= aki
* p
/ ai
;
4369 ? logspace_sub (-M_LN2gnum
, *approxtdistDens
+ gnm_log1p (sum
) + gnm_log (k
* p
) / 2)
4370 : 0.5 - *approxtdistDens
* (sum
+ 1) * gnm_sqrt (k
* p
);
4372 gnm_float q1
= 2 * (1 + q
);
4373 gnm_float q8
= 8 * q
;
4376 gnm_float c1
= -6 * q
;
4378 gnm_float b2
= (k
- 1) * p
+ 3;
4379 gnm_float cadd
= -14 * q
;
4380 gnm_float c2
= b2
+ q1
;
4382 while (gnm_abs (a2
* b1
- a1
* b2
) > gnm_abs (cfVSmall
* b1
* b2
)) {
4383 a1
= c2
* a2
+ c1
* a1
;
4384 b1
= c2
* b2
+ c1
* b1
;
4388 a2
= c2
* a1
+ c1
* a2
;
4389 b2
= c2
* b1
+ c1
* b2
;
4394 if (gnm_abs (b2
) > scalefactor
) {
4395 a1
*= 1 / scalefactor
;
4396 b1
*= 1 / scalefactor
;
4397 a2
*= 1 / scalefactor
;
4398 b2
*= 1 / scalefactor
;
4399 } else if (gnm_abs (b2
) < 1 / scalefactor
) {
4408 ? *approxtdistDens
+ gnm_log1p (-q
* a2
/ b2
) - gnm_log (k
* p
) / 2
4409 : *approxtdistDens
* (1 - q
* a2
/ b2
) / gnm_sqrt (k
* p
);
4414 /* Asymptotic expansion to calculate the probability that binomial variate
4416 * (diffFromMean = (a+b)*p-a).
4419 binApprox (gnm_float a
, gnm_float b
, gnm_float diffFromMean
,
4420 gboolean lower_tail
, gboolean log_p
)
4422 gnm_float pq1
, res
, t
;
4423 gnm_float ib05
, ib15
, ib25
, ib35
, ib3
;
4424 gnm_float elfb
, coef15
, coef25
, coef35
;
4425 gnm_float approxtdistDens
;
4427 gnm_float n
= a
+ b
;
4428 gnm_float n1
= n
+ 1;
4429 gnm_float lvv
= b
- n
* diffFromMean
;
4430 gnm_float lval
= (a
* log1pmx (lvv
/ (a
* n1
)) +
4431 b
* log1pmx (-lvv
/ (b
* n1
))) / n
;
4432 gnm_float tp
= -gnm_expm1 (lval
);
4433 gnm_float mfac
= n1
* tp
;
4434 gnm_float ib2
= 1 + mfac
;
4435 gnm_float t1
= (n
+ 2) * tp
;
4439 ib3
= ib2
+ mfac
*t1
;
4440 ib05
= tdistexp (tp
, 1 - tp
, n1
* lval
, 2 * n1
, log_p
, &approxtdistDens
);
4442 ib15
= gnm_sqrt (mfac
);
4443 mfac
= t1
* (GNM_const (2.0) / 3);
4445 ib35
= ib25
+ mfac
* (n
+ 3) * tp
* (GNM_const (2.0) / 5);
4447 pq1
= (n
* n
) / (a
* b
);
4449 res
= (ib2
* (1 + 2 * pq1
) / 135 - 2 * ib3
* ((2 * pq1
- 43) * pq1
- 22) / (2835 * (n
+ 3))) / (n
+ 2);
4450 res
= (GNM_const (1.0) / 3 - res
) * 2 * gnm_sqrt (pq1
/ n1
) * (a
- b
) / n
;
4453 lower_tail
= !lower_tail
;
4456 n1
= (n
+ 1.5) * (n
+ 2.5);
4457 coef15
= (-17 + 2 * pq1
) / (24 * (n
+ 1.5));
4458 coef25
= (-503 + 4 * pq1
* (19 + pq1
)) / (1152 * n1
);
4459 coef35
= (-315733 + pq1
* (53310 + pq1
* (8196 - 1112 * pq1
))) /
4460 (414720 * n1
* (n
+ 3.5));
4461 elfb
= (coef35
+ coef25
) + coef15
;
4462 res
+= ib15
* ((coef35
* ib35
+ coef25
* ib25
) + coef15
);
4465 ? logspace_signed_add (ib05
,
4466 gnm_log (gnm_abs (res
)) + approxtdistDens
- gnm_log1p (elfb
),
4468 : ib05
+ res
* approxtdistDens
/ (1 + elfb
);
4472 : log_p
? swap_log_tail (t
) : (1 - t
);
4475 /* Probability that binomial variate with sample size i+j and
4476 * event prob p (=1-q) has value i (diffFromMean = (i+j)*p-i)
4479 binomialTerm (gnm_float i
, gnm_float j
, gnm_float p
, gnm_float q
,
4480 gnm_float diffFromMean
, gboolean log_p
)
4482 const gnm_float minLog1Value
= -0.79149064;
4484 gnm_float c4
,c5
,c6
,ps
,logbinomialTerm
,dfm
;
4487 if (i
== 0 && j
<= 0)
4490 if (i
<= -1 || j
< 0)
4503 dfm
= -diffFromMean
;
4506 c5
= (dfm
- (1 - ps
)) / (c2
+ 1);
4507 c6
= -(dfm
+ ps
) / (c3
+ 1);
4509 if (c5
< minLog1Value
) {
4511 logbinomialTerm
= c3
* gnm_log1p (-ps
);
4512 return log_p
? logbinomialTerm
: gnm_exp (logbinomialTerm
);
4513 } else if (ps
== 0 && c2
> 0) {
4516 t
= gnm_log ((ps
* c1
) / (c2
+ 1)) - c5
;
4522 c4
= logfbit (i
+ j
) - logfbit (i
) - logfbit (j
);
4523 logbinomialTerm
= c4
+ c2
* t
- c5
+ (c3
* log1pmx (c6
) - c6
);
4526 ? logbinomialTerm
+ 0.5 * gnm_log (c1
/ ((c2
+ 1) * (c3
+ 1) * 2 * M_PIgnum
))
4527 : gnm_exp (logbinomialTerm
) * gnm_sqrt (c1
/ ((c2
+ 1) * (c3
+ 1) * 2 * M_PIgnum
));
4532 * Probability that binomial variate with sample size ii+jj
4533 * and event prob pp (=1-qq) has value <=i.
4534 * (diffFromMean = (ii+jj)*pp-ii).
4537 binomialcf (gnm_float ii
, gnm_float jj
, gnm_float pp
, gnm_float qq
,
4538 gnm_float diffFromMean
, gboolean lower_tail
, gboolean log_p
)
4540 const gnm_float sumAlways
= 0;
4541 const gnm_float sumFactor
= 6;
4542 const gnm_float cfSmall
= 1.0e-15;
4544 gnm_float prob
,p
,q
,a1
,a2
,b1
,b2
,c1
,c2
,c3
,c4
,n1
,q1
,dfm
;
4545 gnm_float i
,j
,ni
,nj
,numb
,ip1
;
4549 if (ii
> -1 && (jj
<= 0 || pp
== 0)) {
4551 } else if (ii
> -1 && ii
< 0) {
4555 f
= ii
/ ((ii
+ jj
) * pp
);
4556 prob
= binomialTerm (ii
, jj
, pp
, qq
, (ii
+ jj
) * pp
- ii
, log_p
);
4557 prob
= log_p
? prob
+ gnm_log (f
) : prob
* f
;
4559 diffFromMean
= (ii
+ jj
) * pp
- ii
;
4561 prob
= binomialTerm (ii
, jj
, pp
, qq
, diffFromMean
, log_p
);
4567 swapped
= (n1
* qq
>= jj
+ 1);
4569 swapped
= (n1
* pp
<= ii
+ 2);
4571 if (prob
== R_D__0
) {
4572 if (swapped
== !lower_tail
)
4584 dfm
= 1 - diffFromMean
;
4593 if (i
> sumAlways
) {
4594 numb
= gnm_floor (sumFactor
* gnm_sqrt (p
+ 0.5) * gnm_exp (gnm_log (n1
* p
* q
) / 3));
4595 numb
= gnm_floor (numb
- dfm
);
4596 if (numb
> i
) numb
= gnm_floor (i
);
4598 numb
= gnm_floor (i
);
4599 if (numb
< 0) numb
= 0;
4604 a2
= (i
- numb
) * q
;
4605 b2
= dfm
+ numb
+ 1;
4611 while (gnm_abs (a2
* b1
- a1
* b2
) > gnm_abs (cfSmall
* b1
* b2
)) {
4616 a1
= c4
* a2
+ c3
* a1
;
4617 b1
= c4
* b2
+ c3
* b1
;
4622 a2
= c4
* a1
+ c3
* a2
;
4623 b2
= c4
* b1
+ c3
* b2
;
4625 if (gnm_abs (b2
) > scalefactor
) {
4626 a1
*= 1 / scalefactor
;
4627 b1
*= 1 / scalefactor
;
4628 a2
*= 1 / scalefactor
;
4629 b2
*= 1 / scalefactor
;
4630 } else if (gnm_abs (b2
) < 1 / scalefactor
) {
4640 ni
= (i
- numb
+ 1) * q
;
4641 nj
= (j
+ numb
) * p
;
4643 a1
= (1 + a1
) * (ni
/ nj
);
4649 prob
= log_p
? prob
+ gnm_log1p (a1
) : prob
* (1 + a1
);
4653 prob
+= gnm_log (ip1
* q
/ nj
);
4655 prob
*= ip1
* q
/ nj
;
4658 if (swapped
== !lower_tail
)
4661 return log_p
? swap_log_tail (prob
) : 1 - prob
;
4665 binomial (gnm_float ii
, gnm_float jj
, gnm_float pp
, gnm_float qq
,
4666 gnm_float diffFromMean
, gboolean lower_tail
, gboolean log_p
)
4668 gnm_float mij
= fmin2 (ii
, jj
);
4670 if (mij
> 500 && gnm_abs (diffFromMean
) < 0.01 * mij
)
4671 return binApprox (jj
- 1, ii
, diffFromMean
, lower_tail
, log_p
);
4673 return binomialcf (ii
, jj
, pp
, qq
, diffFromMean
, lower_tail
, log_p
);
4677 pbeta (gnm_float x
, gnm_float a
, gnm_float b
, gboolean lower_tail
, gboolean log_p
)
4681 if (gnm_isnan (x
) || gnm_isnan (a
) || gnm_isnan (b
))
4684 if (x
<= 0) return R_DT_0
;
4685 if (x
>= 1) return R_DT_1
;
4687 if (a
< 1 && (b
< 1 || (1 + b
) * x
<= 1))
4688 return pbeta_smalla (x
, a
, b
, lower_tail
, log_p
);
4690 if (b
< 1 && (1 + a
) * (1 - x
) <= 1)
4691 return pbeta_smalla (1 - x
, b
, a
, !lower_tail
, log_p
);
4694 return binomial (-a
, b
, x
, 1 - x
, 0, !lower_tail
, log_p
);
4697 return binomial (-b
, a
, 1 - x
, x
, 0, lower_tail
, log_p
);
4700 return binomial (am1
, b
, x
, 1 - x
, (am1
+ b
) * x
- am1
,
4701 !lower_tail
, log_p
);
4704 /* --- END IANDJMSMITH SOURCE MARKER --- */
4705 /* ------------------------------------------------------------------------ */
4708 pf (gnm_float x
, gnm_float n1
, gnm_float n2
, gboolean lower_tail
, gboolean log_p
)
4711 if (gnm_isnan (x
) || gnm_isnan (n1
) || gnm_isnan (n2
))
4714 if (n1
<= 0 || n2
<= 0) ML_ERR_return_NAN
;
4719 /* Avoid squeezing pbeta's first parameter against 1. */
4721 return pbeta (n2
/ (n2
+ n1
* x
), n2
/ 2, n1
/ 2,
4722 !lower_tail
, log_p
);
4724 return pbeta (n1
* x
/ (n2
+ n1
* x
), n1
/ 2, n2
/ 2,
4728 /* ------------------------------------------------------------------------ */
4731 ppois1 (gnm_float x
, const gnm_float
*plambda
,
4732 gboolean lower_tail
, gboolean log_p
)
4734 return ppois (x
, *plambda
, lower_tail
, log_p
);
4738 qpois (gnm_float p
, gnm_float lambda
, gboolean lower_tail
, gboolean log_p
)
4740 gnm_float mu
, sigma
, gamma
, y
, z
;
4746 sigma
= gnm_sqrt (lambda
);
4749 /* Cornish-Fisher expansion: */
4750 z
= qnorm (p
, 0., 1., lower_tail
, log_p
);
4751 y
= mu
+ sigma
* (z
+ gamma
* (z
* z
- 1) / 6);
4753 return discpfuncinverter (p
, &lambda
, lower_tail
, log_p
,
4758 /* ------------------------------------------------------------------------ */
4761 dgamma1 (gnm_float x
, const gnm_float
*palpha
, gboolean log_p
)
4763 return dgamma (x
, *palpha
, 1, log_p
);
4767 pgamma1 (gnm_float x
, const gnm_float
*palpha
,
4768 gboolean lower_tail
, gboolean log_p
)
4770 return pgamma (x
, *palpha
, 1, lower_tail
, log_p
);
4775 qgamma (gnm_float p
, gnm_float alpha
, gnm_float scale
,
4776 gboolean lower_tail
, gboolean log_p
)
4778 gnm_float res1
, x0
, v
;
4781 if (gnm_isnan(p
) || gnm_isnan(alpha
) || gnm_isnan(scale
))
4782 return p
+ alpha
+ scale
;
4785 if (alpha
<= 0) ML_ERR_return_NAN
;
4787 if (!log_p
&& p
> 0.9) {
4788 /* We're far into the tail. Flip. */
4790 lower_tail
= !lower_tail
;
4793 /* Make an initial guess, x0, assuming scale==1. */
4795 if (v
< -1.24 * R_DT_log (p
))
4796 x0
= gnm_pow (R_DT_qIv (p
) * alpha
* gnm_exp (gnm_lgamma (alpha
) + alpha
* M_LN2gnum
),
4799 gnm_float x1
= qnorm (p
, 0, 1, lower_tail
, log_p
);
4800 gnm_float p1
= 0.222222 / v
;
4801 x0
= v
* gnm_pow (x1
* gnm_sqrt (p1
) + 1 - p1
, 3) / 2;
4804 res1
= pfuncinverter (p
, &alpha
, lower_tail
, log_p
, 0, gnm_pinf
, x0
,
4807 return res1
* scale
;
4810 /* ------------------------------------------------------------------------ */
4813 dbeta1 (gnm_float x
, const gnm_float shape
[], gboolean log_p
)
4815 return dbeta (x
, shape
[0], shape
[1], log_p
);
4819 pbeta1 (gnm_float x
, const gnm_float shape
[],
4820 gboolean lower_tail
, gboolean log_p
)
4822 return pbeta (x
, shape
[0], shape
[1], lower_tail
, log_p
);
4826 abramowitz_stegun_26_5_22 (gnm_float p
, gnm_float a
, gnm_float b
,
4827 gboolean lower_tail
, gboolean log_p
)
4829 gnm_float yp
= qnorm (p
, 0, 1, !lower_tail
, log_p
);
4830 gnm_float ta
= 1 / (2 * a
- 1);
4831 gnm_float tb
= 1 / (2 * b
- 1);
4832 gnm_float h
= 2 / (ta
+ tb
);
4833 gnm_float l
= (yp
* yp
- 3) / 6;
4834 gnm_float w
= yp
* gnm_sqrt (h
+ l
) / h
-
4835 (tb
- ta
) * (l
+ (5 - 4 / h
) / 6);
4836 return a
/ (a
+ b
* gnm_exp (2 * w
));
4841 qbeta (gnm_float p
, gnm_float pin
, gnm_float qin
, gboolean lower_tail
, gboolean log_p
)
4843 gnm_float x0
, q
, shape
[2];
4844 gnm_float S
= pin
+ qin
;
4847 if (gnm_isnan (S
) || gnm_isnan (p
))
4852 if (pin
< 0. || qin
< 0.) ML_ERR_return_NAN
;
4854 if (!log_p
&& p
> 0.9) {
4855 /* We're far into the tail. Flip. */
4857 lower_tail
= !lower_tail
;
4860 if (pin
>= 1 && qin
>= 1)
4861 x0
= abramowitz_stegun_26_5_22 (p
, pin
, qin
, lower_tail
, log_p
);
4864 * The density function is U-shaped.
4866 gnm_float phalf
= pbeta (0.5, pin
, qin
, lower_tail
, log_p
);
4867 gnm_float lb
= gnm_lbeta (pin
, qin
);
4869 if (!lower_tail
== (p
> phalf
)) {
4871 * The following approximation follows from simply ignoring
4872 * the (1-t)^(qin-1) factor from the density. That works
4873 * fine as long as we stay far away from the right tail.
4875 x0
= gnm_exp ((gnm_log (pin
) + R_DT_log (p
) + lb
) / pin
);
4878 x0
= -gnm_expm1 ((gnm_log (qin
) + R_DT_Clog (p
) + lb
) / qin
);
4885 q
= pfuncinverter (p
, shape
, lower_tail
, log_p
, 0, 1, x0
,
4891 qf (gnm_float p
, gnm_float n1
, gnm_float n2
, gboolean lower_tail
, gboolean log_p
)
4895 if (gnm_isnan(p
) || gnm_isnan(n1
) || gnm_isnan(n2
))
4898 if (n1
<= 0. || n2
<= 0.) ML_ERR_return_NAN
;
4904 q
= qbeta (p
, n2
/ 2, n1
/ 2, !lower_tail
, log_p
);
4908 qc
= qbeta (p
, n1
/ 2, n2
/ 2, lower_tail
, log_p
);
4910 return (qc
/ q
) * (n2
/ n1
);
4915 pbinom2 (gnm_float x0
, gnm_float x1
, gnm_float n
, gnm_float p
)
4919 if (x0
> n
|| x1
< 0 || x0
> x1
)
4923 return dbinom (x0
, n
, p
, FALSE
);
4926 return pbinom (x1
, n
, p
, TRUE
, FALSE
);
4929 return pbinom (x0
- 1, n
, p
, FALSE
, FALSE
);
4931 Pl
= pbinom (x0
- 1, n
, p
, TRUE
, FALSE
);
4933 gnm_float PlC
= pbinom (x0
- 1, n
, p
, FALSE
, FALSE
);
4934 gnm_float PrC
= pbinom (x1
, n
, p
, FALSE
, FALSE
);
4937 gnm_float Pr
= pbinom (x1
, n
, p
, TRUE
, FALSE
);
4942 /* ------------------------------------------------------------------------- */
4947 * Returns: The result of exp(-0.5*@x*@x) with higher accuracy than the
4951 expmx2h (gnm_float x
)
4955 if (x
< 5 || gnm_isnan (x
))
4956 return gnm_exp (-0.5 * x
* x
);
4957 else if (x
>= GNM_MAX_EXP
* M_LN2gnum
+ 10)
4961 * Split x into two parts, x=x1+x2, such that |x2|<=2^-16.
4962 * Assuming that we are using IEEE doubles, that means that
4963 * x1*x1 is error free for x<1024 (above which we will underflow
4964 * anyway). If we are not using IEEE doubles then this is
4965 * still an improvement over the naive formula.
4967 gnm_float x1
= gnm_floor (x
* 65536 + 0.5) / 65536;
4968 gnm_float x2
= x
- x1
;
4969 return (gnm_exp (-0.5 * x1
* x1
) *
4970 gnm_exp ((-0.5 * x2
- x1
) * x2
));
4974 /* ------------------------------------------------------------------------- */
4981 * Returns: The arithmetic-geometric mean of @a and @b.
4984 agm (gnm_float a
, gnm_float b
)
4986 gnm_float ab
= a
* b
;
4987 gnm_float scale
= 1;
4990 if (a
< 0 || b
< 0 || gnm_isnan (ab
))
4993 if (a
== gnm_pinf
|| b
== gnm_pinf
)
4995 if (a
== 0 || b
== 0)
4998 if (ab
== 0 || ab
== gnm_pinf
) {
4999 // Underflow or overflow
5001 (void)gnm_frexp (a
, &ea
);
5002 (void)gnm_frexp (b
, &eb
);
5003 scale
= gnm_ldexp (1, -(ea
+ eb
) / 2);
5008 for (i
= 1; i
< 20; i
++) {
5009 gnm_float am
= (a
+ b
) / 2;
5010 gnm_float gm
= gnm_sqrt (a
* b
);
5013 if (gnm_abs (a
- b
) < a
* GNM_EPSILON
)
5017 g_warning ("AGM failed to converge.");
5028 * Returns: The result of (1+@x)^@y with less rounding error than the
5032 pow1p (gnm_float x
, gnm_float y
)
5035 * We defer to the naive algorithm in two cases: (1) when 1+x
5036 * is exact (let us hope the compiler does not mess this up),
5037 * and (2) when |x|>1/2 and we have no better algorithm.
5040 if ((x
+ 1) - x
== 1 || gnm_abs (x
) > 0.5 ||
5041 gnm_isnan (x
) || gnm_isnan (y
))
5042 return gnm_pow (1 + x
, y
);
5044 return 1 / pow1p (x
, -y
);
5046 gnm_float x1
= gnm_floor (x
* 65536 + 0.5) / 65536;
5047 gnm_float x2
= x
- x1
;
5049 ebd0 (y
, y
*(x
+1), &h
, &l
);
5050 PAIR_ADD (-y
*x1
, h
, l
);
5051 PAIR_ADD (-y
*x2
, h
, l
);
5054 g_printerr ("pow1p (%.20g, %.20g)\n", x
, y
);
5057 return gnm_exp (-l
) * gnm_exp (-h
);
5066 * Returns: The result of (1+@x)^@y-1 with less rounding error than the
5070 pow1pm1 (gnm_float x
, gnm_float y
)
5073 return gnm_pow (1 + x
, y
) - 1;
5075 return gnm_expm1 (y
* gnm_log1p (x
));
5080 ---------------------------------------------------------------------
5082 ---------------------------------------------------------------------
5086 * gnm_matrix_new: (skip)
5088 /* Note the order: y then x. */
5090 gnm_matrix_new (int rows
, int cols
)
5092 GnmMatrix
*m
= g_new (GnmMatrix
, 1);
5097 m
->data
= g_new (gnm_float
*, rows
);
5098 for (r
= 0; r
< rows
; r
++)
5099 m
->data
[r
] = g_new (gnm_float
, cols
);
5105 gnm_matrix_free (GnmMatrix
*m
)
5109 for (r
= 0; r
< m
->rows
; r
++)
5110 g_free (m
->data
[r
]);
5116 gnm_matrix_is_empty (GnmMatrix
const *m
)
5118 return m
== NULL
|| m
->rows
<= 0 || m
->cols
<= 0;
5122 * gnm_matrix_from_value: (skip)
5125 gnm_matrix_from_value (GnmValue
const *v
, GnmValue
**perr
, GnmEvalPos
const *ep
)
5129 GnmMatrix
*m
= NULL
;
5132 cols
= value_area_get_width (v
, ep
);
5133 rows
= value_area_get_height (v
, ep
);
5134 m
= gnm_matrix_new (rows
, cols
);
5135 for (r
= 0; r
< rows
; r
++) {
5136 for (c
= 0; c
< cols
; c
++) {
5137 GnmValue
const *v1
= value_area_fetch_x_y (v
, c
, r
, ep
);
5138 if (VALUE_IS_ERROR (v1
)) {
5139 *perr
= value_dup (v1
);
5140 gnm_matrix_free (m
);
5144 m
->data
[r
][c
] = value_get_as_float (v1
);
5151 gnm_matrix_to_value (GnmMatrix
const *m
)
5153 GnmValue
*res
= value_new_array_non_init (m
->cols
, m
->rows
);
5156 for (c
= 0; c
< m
->cols
; c
++) {
5157 res
->v_array
.vals
[c
] = g_new (GnmValue
*, m
->rows
);
5158 for (r
= 0; r
< m
->rows
; r
++)
5159 res
->v_array
.vals
[c
][r
] = value_new_float (m
->data
[r
][c
]);
5166 gnm_matrix_multiply (GnmMatrix
*C
, const GnmMatrix
*A
, const GnmMatrix
*B
)
5169 GnmAccumulator
*acc
;
5172 g_return_if_fail (C
!= NULL
);
5173 g_return_if_fail (A
!= NULL
);
5174 g_return_if_fail (B
!= NULL
);
5175 g_return_if_fail (C
->rows
== A
->rows
);
5176 g_return_if_fail (C
->cols
== B
->cols
);
5177 g_return_if_fail (A
->cols
== B
->rows
);
5179 state
= gnm_accumulator_start ();
5180 acc
= gnm_accumulator_new ();
5182 for (r
= 0; r
< C
->rows
; r
++) {
5183 for (c
= 0; c
< C
->cols
; c
++) {
5184 go_accumulator_clear (acc
);
5185 for (i
= 0; i
< A
->cols
; ++i
) {
5190 gnm_accumulator_add_quad (acc
, &p
);
5192 C
->data
[r
][c
] = gnm_accumulator_value (acc
);
5196 gnm_accumulator_free (acc
);
5197 gnm_accumulator_end (state
);
5200 /***************************************************************************/
5203 gnm_matrix_eigen_max_index (gnm_float
*row
, guint row_n
, guint size
)
5205 guint i
, res
= row_n
+ 1;
5211 max
= gnm_abs (row
[res
]);
5213 for (i
= res
+ 1; i
< size
; i
++)
5214 if (gnm_abs (row
[i
]) > max
) {
5216 max
= gnm_abs (row
[i
]);
5222 gnm_matrix_eigen_rotate (gnm_float
**matrix
, guint k
, guint l
, guint i
, guint j
, gnm_float c
, gnm_float s
)
5224 gnm_float x
= c
* matrix
[k
][l
] - s
* matrix
[i
][j
];
5225 gnm_float y
= s
* matrix
[k
][l
] + c
* matrix
[i
][j
];
5232 gnm_matrix_eigen_update (guint k
, gnm_float t
, gnm_float
*eigenvalues
, gboolean
*changed
, guint
*state
)
5234 gnm_float y
= eigenvalues
[k
];
5236 eigenvalues
[k
] += t
;
5237 unchanged
= (y
== eigenvalues
[k
]);
5238 if (changed
[k
] && unchanged
) {
5241 } else if (!changed
[k
] && !unchanged
) {
5248 * Calculates the eigenvalues and eigenvectors of a real symmetric matrix.
5250 * This is the Jacobi iterative process in which we use a sequence of
5251 * Jacobi rotations (two-sided Givens rotations) in order to reduce the
5252 * magnitude of off-diagonal elements while preserving eigenvalues.
5255 gnm_matrix_eigen (GnmMatrix
const *m
, GnmMatrix
*EIG
, gnm_float
*eigenvalues
)
5257 guint i
, state
, usize
, *ind
;
5261 gnm_float
**eigenvectors
;
5263 g_return_val_if_fail (m
!= NULL
, FALSE
);
5264 g_return_val_if_fail (m
->rows
== m
->cols
, FALSE
);
5265 g_return_val_if_fail (EIG
!= NULL
, FALSE
);
5266 g_return_val_if_fail (EIG
->rows
== EIG
->cols
, FALSE
);
5267 g_return_val_if_fail (EIG
->rows
== m
->rows
, FALSE
);
5270 eigenvectors
= EIG
->data
;
5275 ind
= g_new (guint
, usize
);
5276 changed
= g_new (gboolean
, usize
);
5278 for (i
= 0; i
< usize
; i
++) {
5280 for (j
= 0; j
< usize
; j
++)
5281 eigenvectors
[j
][i
] = 0.;
5282 eigenvectors
[i
][i
] = 1.;
5283 eigenvalues
[i
] = matrix
[i
][i
];
5284 ind
[i
] = gnm_matrix_eigen_max_index (matrix
[i
], i
, usize
);
5288 while (usize
> 1 && state
!= 0) {
5290 gnm_float c
, s
, y
, pivot
, t
;
5293 if (counter
> 400000) {
5296 g_print ("gnm_matrix_eigen exceeded iterations\n");
5299 for (k
= 1; k
< usize
- 1; k
++)
5300 if (gnm_abs (matrix
[k
][ind
[k
]]) > gnm_abs (matrix
[m
][ind
[m
]]))
5303 pivot
= matrix
[m
][l
];
5304 /* pivot is (m,l) */
5306 /* All remaining off-diagonal elements are zero. We're done. */
5310 y
= (eigenvalues
[l
] - eigenvalues
[m
]) / 2;
5311 t
= gnm_abs (y
) + gnm_hypot (pivot
, y
);
5312 s
= gnm_hypot (pivot
, t
);
5315 t
= pivot
* pivot
/ t
;
5321 gnm_matrix_eigen_update (m
, -t
, eigenvalues
, changed
, &state
);
5322 gnm_matrix_eigen_update (l
, t
, eigenvalues
, changed
, &state
);
5323 for (i
= 0; i
< m
; i
++)
5324 gnm_matrix_eigen_rotate (matrix
, i
, m
, i
, l
, c
, s
);
5325 for (i
= m
+ 1; i
< l
; i
++)
5326 gnm_matrix_eigen_rotate (matrix
, m
, i
, i
, l
, c
, s
);
5327 for (i
= l
+ 1; i
< usize
; i
++)
5328 gnm_matrix_eigen_rotate (matrix
, m
, i
, l
, i
, c
, s
);
5329 for (i
= 0; i
< usize
; i
++) {
5330 gnm_float x
= c
* eigenvectors
[i
][m
] - s
* eigenvectors
[i
][l
];
5331 gnm_float y
= s
* eigenvectors
[i
][m
] + c
* eigenvectors
[i
][l
];
5333 eigenvectors
[i
][m
] = x
;
5334 eigenvectors
[i
][l
] = y
;
5336 ind
[m
] = gnm_matrix_eigen_max_index (matrix
[m
], m
, usize
);
5337 ind
[l
] = gnm_matrix_eigen_max_index (matrix
[l
], l
, usize
);
5346 /* ------------------------------------------------------------------------- */
5349 swap_row_and_col (GnmMatrix
*M
, int a
, int b
)
5356 M
->data
[a
] = M
->data
[b
];
5360 for (i
= 0; i
< M
->rows
; i
++) {
5361 gnm_float d
= M
->data
[i
][a
];
5362 M
->data
[i
][a
] = M
->data
[i
][b
];
5369 gnm_matrix_modified_cholesky (GnmMatrix
const *A
,
5376 gnm_float nu
, xsi
, gam
, bsqr
, delta
;
5380 g_return_val_if_fail (A
->rows
== A
->cols
, FALSE
);
5381 g_return_val_if_fail (A
->rows
== L
->rows
, FALSE
);
5382 g_return_val_if_fail (A
->cols
== L
->cols
, FALSE
);
5384 // Copy A into L; Use G and C as aliases for L.
5385 for (i
= 0; i
< n
; i
++)
5386 for (j
= 0; j
< n
; j
++)
5387 L
->data
[i
][j
] = A
->data
[i
][j
];
5391 // Init permutation as identity.
5392 for (i
= 0; i
< n
; i
++)
5395 nu
= n
== 1 ? 1.0 : gnm_sqrt (n
* n
- 1);
5397 for (i
= 0; i
< n
; i
++) {
5398 gnm_float aii
= gnm_abs (G
->data
[i
][i
]);
5399 gam
= MAX (gam
, aii
);
5400 for (j
= i
+ 1; j
< n
; j
++) {
5401 gnm_float aij
= gnm_abs (G
->data
[i
][j
]);
5402 xsi
= MAX (xsi
, aij
);
5405 bsqr
= MAX (MAX (gam
, xsi
/ nu
), GNM_EPSILON
);
5406 delta
= MAX (gam
+ xsi
, 1.0) * GNM_EPSILON
;
5408 for (j
= 0; j
< n
; j
++) {
5410 gnm_float theta_j
= 0, dj
;
5413 for (i
= j
+ 1; i
< n
; i
++) {
5414 if (gnm_abs (C
->data
[i
][i
]) > gnm_abs (C
->data
[q
][q
]))
5422 swap_row_and_col (L
, j
, q
);
5423 a
= P
[j
]; P
[j
] = P
[q
]; P
[q
] = a
;
5424 b
= D
[j
]; D
[j
] = D
[q
]; D
[q
] = b
;
5425 if (E
) { b
= E
[j
]; E
[j
] = E
[q
]; E
[q
] = b
; }
5428 for (s
= 0; s
< j
; s
++)
5429 L
->data
[j
][s
] = C
->data
[j
][s
] / D
[s
];
5431 for (i
= j
+ 1; i
< n
; i
++) {
5433 gnm_float d
= G
->data
[i
][j
];
5435 for (s
= 0; s
< j
; s
++)
5436 d
-= L
->data
[j
][s
] * C
->data
[i
][s
];
5439 theta_j
= MAX (theta_j
, gnm_abs (d
));
5442 dj
= MAX (theta_j
* theta_j
/ bsqr
, delta
);
5443 dj
= MAX (dj
, gnm_abs (C
->data
[j
][j
]));
5445 if (E
) E
[j
] = dj
- C
->data
[j
][j
];
5447 for (i
= j
+ 1; i
< n
; i
++) {
5448 gnm_float cij
= C
->data
[i
][j
];
5449 C
->data
[i
][i
] -= cij
* cij
/ D
[j
];
5453 for (i
= 0; i
< n
; i
++) {
5454 for (j
= i
+ 1; j
< n
; j
++)
5463 gnm_linear_solve_posdef (GnmMatrix
const *A
, const gnm_float
*b
,
5470 GORegressionResult res
;
5473 g_return_val_if_fail (A
!= NULL
, GO_REG_invalid_dimensions
);
5474 g_return_val_if_fail (A
->rows
== A
->cols
, GO_REG_invalid_dimensions
);
5475 g_return_val_if_fail (b
!= NULL
, GO_REG_invalid_dimensions
);
5476 g_return_val_if_fail (x
!= NULL
, GO_REG_invalid_dimensions
);
5479 L
= gnm_matrix_new (n
, n
);
5480 D
= g_new (gnm_float
, n
);
5481 E
= g_new (gnm_float
, n
);
5484 ok
= gnm_matrix_modified_cholesky (A
, L
, D
, E
, P
);
5486 res
= GO_REG_invalid_data
;
5490 if (gnm_debug_flag ("posdef")) {
5491 for (i
= 0; i
< n
; i
++)
5492 g_printerr ("Posdef E[i] = %g\n", E
[P
[i
]]);
5495 // The only information from the above we use is E and P.
5496 // However, we reuse the memory for L for A+E
5497 for (i
= 0; i
< n
; i
++) {
5498 for (j
= 0; j
< n
; j
++)
5499 L
->data
[i
][j
] = A
->data
[i
][j
];
5500 L
->data
[i
][i
] += E
[P
[i
]];
5503 res
= gnm_linear_solve (L
, b
, x
);
5509 gnm_matrix_free (L
);
5514 /* ------------------------------------------------------------------------- */
5517 gnm_linear_solve (GnmMatrix
const *A
, const gnm_float
*b
,
5520 g_return_val_if_fail (A
!= NULL
, GO_REG_invalid_dimensions
);
5521 g_return_val_if_fail (A
->rows
== A
->cols
, GO_REG_invalid_dimensions
);
5522 g_return_val_if_fail (b
!= NULL
, GO_REG_invalid_dimensions
);
5523 g_return_val_if_fail (x
!= NULL
, GO_REG_invalid_dimensions
);
5526 #ifdef GNM_WITH_LONG_DOUBLE
5531 (A
->data
, b
, A
->rows
, x
);
5535 gnm_linear_solve_multiple (GnmMatrix
const *A
, GnmMatrix
*B
)
5537 g_return_val_if_fail (A
!= NULL
, GO_REG_invalid_dimensions
);
5538 g_return_val_if_fail (B
!= NULL
, GO_REG_invalid_dimensions
);
5539 g_return_val_if_fail (A
->rows
== A
->cols
, GO_REG_invalid_dimensions
);
5540 g_return_val_if_fail (A
->rows
== B
->rows
, GO_REG_invalid_dimensions
);
5543 #ifdef GNM_WITH_LONG_DOUBLE
5544 go_linear_solve_multiplel
5546 go_linear_solve_multiple
5548 (A
->data
, B
->data
, A
->rows
, B
->cols
);
5551 /* ------------------------------------------------------------------------- */
5553 #ifdef GNM_SUPPLIES_ERFL
5555 erfl (long double x
)
5557 if (fabsl (x
) < 0.125) {
5558 /* For small x the pnorm formula loses precision. */
5559 long double sum
= 0;
5560 long double term
= x
* 2 / sqrtl (M_PIgnum
);
5562 long double x2
= x
* x
;
5564 for (n
= 0; fabsl (term
) >= fabsl (sum
) * LDBL_EPSILON
; n
++) {
5565 sum
+= term
/ (2 * n
+ 1);
5566 term
*= -x2
/ (n
+ 1);
5571 return pnorm (x
* M_SQRT2gnum
, 0, 1, TRUE
, FALSE
) * 2 - 1;
5575 /* ------------------------------------------------------------------------- */
5577 #ifdef GNM_SUPPLIES_ERFCL
5579 erfcl (long double x
)
5581 return 2 * pnorm (x
* M_SQRT2gnum
, 0, 1, FALSE
, FALSE
);
5585 /* ------------------------------------------------------------------------- */
5588 gnm_owent_T1 (gnm_float h
, gnm_float a
, int order
)
5590 const gnm_float hs
= -0.5 * (h
* h
);
5591 const gnm_float dhs
= gnm_exp (hs
);
5592 const gnm_float as
= a
* a
;
5593 gnm_float aj
= a
/ (M_PIgnum
* 2);
5594 gnm_float dj
= gnm_expm1 (hs
);
5595 gnm_float gj
= hs
* dhs
;
5596 gnm_float res
= gnm_atan (a
) / (M_PIgnum
* 2);
5599 for (j
= 1; j
<= order
; j
++) {
5600 res
+= dj
* aj
/ (j
+ j
- 1);
5611 gnm_owent_T2 (gnm_float h
, gnm_float a
, int order
)
5613 const gnm_float ah
= a
* h
;
5614 const gnm_float as
= -a
* a
;
5615 const gnm_float y
= 1 / (h
* h
);
5617 gnm_float vi
= a
* dnorm (ah
, 0, 1, FALSE
);
5618 gnm_float z
= gnm_erf (ah
/ M_SQRT2gnum
) / (2 * h
);
5621 for (i
= 1; i
<= 2 * order
+ 1; i
+= 2) {
5623 z
= y
* (vi
- i
* z
);
5626 return val
* dnorm (h
, 0, 1, FALSE
);
5630 gnm_owent_T3 (gnm_float h
, gnm_float a
, int order
)
5632 static const gnm_float c2
[] = {
5633 GNM_const(+0.99999999999999987510),
5634 GNM_const(-0.99999999999988796462),
5635 GNM_const(+0.99999999998290743652),
5636 GNM_const(-0.99999999896282500134),
5637 GNM_const(+0.99999996660459362918),
5638 GNM_const(-0.99999933986272476760),
5639 GNM_const(+0.99999125611136965852),
5640 GNM_const(-0.99991777624463387686),
5641 GNM_const(+0.99942835555870132569),
5642 GNM_const(-0.99697311720723000295),
5643 GNM_const(+0.98751448037275303682),
5644 GNM_const(-0.95915857980572882813),
5645 GNM_const(+0.89246305511006708555),
5646 GNM_const(-0.76893425990463999675),
5647 GNM_const(+0.58893528468484693250),
5648 GNM_const(-0.38380345160440256652),
5649 GNM_const(+0.20317601701045299653),
5650 GNM_const(-0.82813631607004984866E-01),
5651 GNM_const(+0.24167984735759576523E-01),
5652 GNM_const(-0.44676566663971825242E-02),
5653 GNM_const(+0.39141169402373836468E-03)
5656 const gnm_float ah
= a
* h
;
5657 const gnm_float as
= a
* a
;
5658 const gnm_float y
= 1 / (h
* h
);
5659 gnm_float vi
= a
* dnorm (ah
, 0, 1, FALSE
);
5660 gnm_float zi
= gnm_erf (ah
/ M_SQRT2gnum
) / (2 * h
);
5664 g_return_val_if_fail (order
< (int)G_N_ELEMENTS(c2
), gnm_nan
);
5666 for (i
= 0; i
<= order
; i
++) {
5668 zi
= y
* ((i
+ i
+ 1) * zi
- vi
);
5671 return val
* dnorm (h
, 0, 1, FALSE
);
5675 gnm_owent_T4 (gnm_float h
, gnm_float a
, int order
)
5677 const gnm_float hs
= h
* h
;
5678 const gnm_float as
= -a
* a
;
5679 gnm_float ai
= a
* gnm_exp (-0.5 * hs
* (1 - as
)) / (2 * M_PIgnum
);
5684 for (i
= 1; i
<= 2 * order
+ 1; i
+= 2) {
5686 yi
= (1 - hs
* yi
) / (i
+ 2);
5693 gnm_owent_T5 (gnm_float h
, gnm_float a
, int order
)
5695 static const gnm_float pts
[] = {
5696 GNM_const(0.35082039676451715489E-02),
5697 GNM_const(0.31279042338030753740E-01),
5698 GNM_const(0.85266826283219451090E-01),
5699 GNM_const(0.16245071730812277011),
5700 GNM_const(0.25851196049125434828),
5701 GNM_const(0.36807553840697533536),
5702 GNM_const(0.48501092905604697475),
5703 GNM_const(0.60277514152618576821),
5704 GNM_const(0.71477884217753226516),
5705 GNM_const(0.81475510988760098605),
5706 GNM_const(0.89711029755948965867),
5707 GNM_const(0.95723808085944261843),
5708 GNM_const(0.99178832974629703586)
5710 static const gnm_float wts
[] = {
5711 0.18831438115323502887E-01,
5712 0.18567086243977649478E-01,
5713 0.18042093461223385584E-01,
5714 0.17263829606398753364E-01,
5715 0.16243219975989856730E-01,
5716 0.14994592034116704829E-01,
5717 0.13535474469662088392E-01,
5718 0.11886351605820165233E-01,
5719 0.10070377242777431897E-01,
5720 0.81130545742299586629E-02,
5721 0.60419009528470238773E-02,
5722 0.38862217010742057883E-02,
5723 0.16793031084546090448E-02
5725 const gnm_float as
= a
* a
;
5726 const gnm_float hs
= -0.5 * h
* h
;
5730 g_return_val_if_fail (order
<= (int)G_N_ELEMENTS(pts
), gnm_nan
);
5731 g_return_val_if_fail (order
<= (int)G_N_ELEMENTS(wts
), gnm_nan
);
5733 for (i
= 0; i
< order
; i
++) {
5734 gnm_float r
= 1 + as
* pts
[i
];
5735 val
+= wts
[i
] * gnm_exp (hs
* r
) / r
;
5742 gnm_owent_T6 (gnm_float h
, gnm_float a
, int order
)
5744 const gnm_float normh
= pnorm (h
, 0, 1, FALSE
, FALSE
);
5745 const gnm_float normhC
= 1 - normh
;
5746 const gnm_float y
= 1 - a
;
5747 const gnm_float r
= gnm_atan2 (y
, 1 + a
);
5748 gnm_float val
= 0.5 * normh
* normhC
;
5751 val
-= r
* gnm_exp (-0.5 * y
* h
* h
/ r
) / (2 * M_PIgnum
);
5757 gnm_owent_helper (gnm_float h
, gnm_float a
)
5759 static const double hrange
[] = {
5760 0.02, 0.06, 0.09, 0.125, 0.26, 0.4, 0.6,
5761 1.6, 1.7, 2.33, 2.4, 3.36, 3.4, 4.8
5763 static const double arange
[] = {
5764 0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999
5766 static const guint8 method
[] = {
5767 1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9,
5768 1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9,
5769 2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10,
5770 2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10,
5771 2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11,
5772 2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12,
5773 2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12,
5774 2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12
5778 g_return_val_if_fail (h
>= 0, gnm_nan
);
5779 g_return_val_if_fail (a
>= 0 && a
<= 1, gnm_nan
);
5781 for (ai
= 0; ai
< (int)G_N_ELEMENTS(arange
); ai
++)
5782 if (a
<= arange
[ai
])
5784 for (hi
= 0; hi
< (int)G_N_ELEMENTS(hrange
); hi
++)
5785 if (h
<= hrange
[hi
])
5788 switch (method
[ai
* (1 + G_N_ELEMENTS(hrange
)) + hi
]) {
5789 case 1: return gnm_owent_T1 (h
, a
, 2);
5790 case 2: return gnm_owent_T1 (h
, a
, 3);
5791 case 3: return gnm_owent_T1 (h
, a
, 4);
5792 case 4: return gnm_owent_T1 (h
, a
, 5);
5793 case 5: return gnm_owent_T1 (h
, a
, 7);
5794 case 6: return gnm_owent_T1 (h
, a
, 10);
5795 case 7: return gnm_owent_T1 (h
, a
, 12);
5796 case 8: return gnm_owent_T1 (h
, a
, 18);
5797 case 9: return gnm_owent_T2 (h
, a
, 10);
5798 case 10: return gnm_owent_T2 (h
, a
, 20);
5799 case 11: return gnm_owent_T2 (h
, a
, 30);
5800 case 12: return gnm_owent_T3 (h
, a
, 20);
5801 case 13: return gnm_owent_T4 (h
, a
, 4);
5802 case 14: return gnm_owent_T4 (h
, a
, 7);
5803 case 15: return gnm_owent_T4 (h
, a
, 8);
5804 case 16: return gnm_owent_T4 (h
, a
, 20);
5805 case 17: return gnm_owent_T5 (h
, a
, 13);
5806 case 18: return gnm_owent_T6 (h
, a
, 0);
5808 g_assert_not_reached ();
5813 * See "Fast and Accurate Calculation of Owen’s T-Function" by
5814 * Mike Patefield and David Tandy. Journal of Statistical Software,
5815 * Volume 5, Issue 5, July 2000.
5817 * Original code licensed under GPLv2+.
5820 gnm_owent (gnm_float h
, gnm_float a
)
5825 /* Even in the "h" argument. */
5828 /* Odd in the "a" argument. */
5835 res
= gnm_atan (a
) / (2 * M_PIgnum
);
5837 res
= 0.5 * pnorm (h
, 0, 1, TRUE
, FALSE
) *
5838 pnorm (h
, 0, 1, FALSE
, FALSE
);
5840 res
= gnm_owent_helper (h
, a
);
5842 gnm_float ah
= h
* a
;
5846 * T(h,a) = .5N(h) + .5N(ha) - N(h)N(ha) - T(ha,1/a)
5848 * with care to avoid cancellation.
5851 gnm_float nh
= 0.5 * gnm_erf (h
/ M_SQRT2gnum
);
5852 gnm_float nah
= 0.5 * gnm_erf (ah
/ M_SQRT2gnum
);
5853 res
= 0.25 - nh
* nah
-
5854 gnm_owent_helper (ah
, 1 / a
);
5856 gnm_float nh
= pnorm (h
, 0, 1, FALSE
, FALSE
);
5857 gnm_float nah
= pnorm (ah
, 0, 1, FALSE
, FALSE
);
5858 res
= 0.5 * (nh
+ nah
) - nh
* nah
-
5859 gnm_owent_helper (ah
, 1 / a
);
5863 /* Odd in the "a" argument. */
5870 /* ------------------------------------------------------------------------- */