1 /******************************************************************************\
2 * *
3 * SimpleCUDD library (www.cs.kuleuven.be/~theo/tools/simplecudd.html) *
4 * SimpleCUDD was developed at Katholieke Universiteit Leuven(www.kuleuven.be) *
5 * *
6 * Copyright Katholieke Universiteit Leuven 2008, 2009, 2010 *
7 * *
8 * Author: Bernd Gutmann *
9 * File: problogmath.c *
10 * $Date:: 2010-12-17 12:21:58 +0100 (Fri, 17 Dec 2010) $ *
11 * $Revision:: 5159 $ *
12 * *
13 ********************************************************************************
14 * *
15 * Artistic License 2.0 *
16 * *
17 * Copyright (c) 2000-2006, The Perl Foundation. *
18 * *
19 * Everyone is permitted to copy and distribute verbatim copies of this license *
20 * document, but changing it is not allowed. *
21 * *
22 * Preamble *
23 * *
24 * This license establishes the terms under which a given free software Package *
25 * may be copied, modified, distributed, and/or redistributed. The intent is *
26 * that the Copyright Holder maintains some artistic control over the *
27 * development of that Package while still keeping the Package available as *
28 * open source and free software. *
29 * *
30 * You are always permitted to make arrangements wholly outside of this license *
31 * directly with the Copyright Holder of a given Package. If the terms of this *
32 * license do not permit the full use that you propose to make of the Package, *
33 * you should contact the Copyright Holder and seek a different licensing *
34 * arrangement. *
35 * Definitions *
36 * *
37 * "Copyright Holder" means the individual(s) or organization(s) named in the *
38 * copyright notice for the entire Package. *
39 * *
40 * "Contributor" means any party that has contributed code or other material to *
41 * the Package, in accordance with the Copyright Holder's procedures. *
42 * *
43 * "You" and "your" means any person who would like to copy, distribute, or *
44 * modify the Package. *
45 * *
46 * "Package" means the collection of files distributed by the Copyright Holder, *
47 * and derivatives of that collection and/or of those files. A given Package *
48 * may consist of either the Standard Version, or a Modified Version. *
49 * *
50 * "Distribute" means providing a copy of the Package or making it accessible *
51 * to anyone else, or in the case of a company or organization, to others *
52 * outside of your company or organization. *
53 * *
54 * "Distributor Fee" means any fee that you charge for Distributing this *
55 * Package or providing support for this Package to another party. It does not *
56 * mean licensing fees. *
57 * *
58 * "Standard Version" refers to the Package if it has not been modified, or has *
59 * been modified only in ways explicitly requested by the Copyright Holder. *
60 * *
61 * "Modified Version" means the Package, if it has been changed, and such *
62 * changes were not explicitly requested by the Copyright Holder. *
63 * *
64 * "Original License" means this Artistic License as Distributed with the *
65 * Standard Version of the Package, in its current version or as it may be *
66 * modified by The Perl Foundation in the future. *
67 * *
68 * "Source" form means the source code, documentation source, and configuration *
69 * files for the Package. *
70 * *
71 * "Compiled" form means the compiled bytecode, object code, binary, or any *
72 * other form resulting from mechanical transformation or translation of the *
73 * Source form. *
74 * Permission for Use and Modification Without Distribution *
75 * *
76 * (1) You are permitted to use the Standard Version and create and use *
77 * Modified Versions for any purpose without restriction, provided that you do *
78 * not Distribute the Modified Version. *
79 * Permissions for Redistribution of the Standard Version *
80 * *
81 * (2) You may Distribute verbatim copies of the Source form of the Standard *
82 * Version of this Package in any medium without restriction, either gratis or *
83 * for a Distributor Fee, provided that you duplicate all of the original *
84 * copyright notices and associated disclaimers. At your discretion, such *
85 * verbatim copies may or may not include a Compiled form of the Package. *
86 * *
87 * (3) You may apply any bug fixes, portability changes, and other *
88 * modifications made available from the Copyright Holder. The resulting *
89 * Package will still be considered the Standard Version, and as such will be *
90 * subject to the Original License. *
91 * Distribution of Modified Versions of the Package as Source *
92 * *
93 * (4) You may Distribute your Modified Version as Source (either gratis or for *
94 * a Distributor Fee, and with or without a Compiled form of the Modified *
95 * Version) provided that you clearly document how it differs from the Standard *
96 * Version, including, but not limited to, documenting any non-standard *
97 * features, executables, or modules, and provided that you do at least ONE of *
98 * the following: *
99 * *
100 * (a) make the Modified Version available to the Copyright Holder of the *
101 * Standard Version, under the Original License, so that the Copyright Holder *
102 * may include your modifications in the Standard Version. *
103 * (b) ensure that installation of your Modified Version does not prevent the *
104 * user installing or running the Standard Version. In addition, the Modified *
105 * Version must bear a name that is different from the name of the Standard *
106 * Version. *
107 * (c) allow anyone who receives a copy of the Modified Version to make the *
108 * Source form of the Modified Version available to others under *
109 * (i) the Original License or *
110 * (ii) a license that permits the licensee to freely copy, modify and *
111 * redistribute the Modified Version using the same licensing terms that apply *
112 * to the copy that the licensee received, and requires that the Source form of *
113 * the Modified Version, and of any works derived from it, be made freely *
114 * available in that license fees are prohibited but Distributor Fees are *
115 * allowed. *
116 * Distribution of Compiled Forms of the Standard Version or Modified Versions *
117 * without the Source *
118 * *
119 * (5) You may Distribute Compiled forms of the Standard Version without the *
120 * Source, provided that you include complete instructions on how to get the *
121 * Source of the Standard Version. Such instructions must be valid at the time *
122 * of your distribution. If these instructions, at any time while you are *
123 * carrying out such distribution, become invalid, you must provide new *
124 * instructions on demand or cease further distribution. If you provide valid *
125 * instructions or cease distribution within thirty days after you become aware *
126 * that the instructions are invalid, then you do not forfeit any of your *
127 * rights under this license. *
128 * *
129 * (6) You may Distribute a Modified Version in Compiled form without the *
130 * Source, provided that you comply with Section 4 with respect to the Source *
131 * of the Modified Version. *
132 * Aggregating or Linking the Package *
133 * *
134 * (7) You may aggregate the Package (either the Standard Version or Modified *
135 * Version) with other packages and Distribute the resulting aggregation *
136 * provided that you do not charge a licensing fee for the Package. Distributor *
137 * Fees are permitted, and licensing fees for other components in the *
138 * aggregation are permitted. The terms of this license apply to the use and *
139 * Distribution of the Standard or Modified Versions as included in the *
140 * aggregation. *
141 * *
142 * (8) You are permitted to link Modified and Standard Versions with other *
143 * works, to embed the Package in a larger work of your own, or to build *
144 * stand-alone binary or bytecode versions of applications that include the *
145 * Package, and Distribute the result without restriction, provided the result *
146 * does not expose a direct interface to the Package. *
147 * Items That are Not Considered Part of a Modified Version *
148 * *
149 * (9) Works (including, but not limited to, modules and scripts) that merely *
150 * extend or make use of the Package, do not, by themselves, cause the Package *
151 * to be a Modified Version. In addition, such works are not considered parts *
152 * of the Package itself, and are not subject to the terms of this license. *
153 * General Provisions *
154 * *
155 * (10) Any use, modification, and distribution of the Standard or Modified *
156 * Versions is governed by this Artistic License. By using, modifying or *
157 * distributing the Package, you accept this license. Do not use, modify, or *
158 * distribute the Package, if you do not accept this license. *
159 * *
160 * (11) If your Modified Version has been derived from a Modified Version made *
161 * by someone other than you, you are nevertheless required to ensure that your *
162 * Modified Version complies with the requirements of this license. *
163 * *
164 * (12) This license does not grant you the right to use any trademark, service *
165 * mark, tradename, or logo of the Copyright Holder. *
166 * *
167 * (13) This license includes the non-exclusive, worldwide, free-of-charge *
168 * patent license to make, have made, use, offer to sell, sell, import and *
169 * otherwise transfer the Package with respect to any patent claims licensable *
170 * by the Copyright Holder that are necessarily infringed by the Package. If *
171 * you institute patent litigation (including a cross-claim or counterclaim) *
172 * against any party alleging that the Package constitutes direct or *
173 * contributory patent infringement, then this Artistic License to you shall *
174 * terminate on the date that such litigation is filed. *
175 * *
176 * (14) Disclaimer of Warranty: THE PACKAGE IS PROVIDED BY THE COPYRIGHT HOLDER *
177 * AND CONTRIBUTORS "AS IS' AND WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES. THE *
178 * IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR *
179 * NON-INFRINGEMENT ARE DISCLAIMED TO THE EXTENT PERMITTED BY YOUR LOCAL LAW. *
180 * UNLESS REQUIRED BY LAW, NO COPYRIGHT HOLDER OR CONTRIBUTOR WILL BE LIABLE *
181 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING IN *
182 * ANY WAY OUT OF THE USE OF THE PACKAGE, EVEN IF ADVISED OF THE POSSIBILITY OF *
183 * SUCH DAMAGE. *
184 * *
185 * The End *
186 * *
187 \******************************************************************************/
188
189 #include "problogmath.h"
190 #include "general.h"
191
sigmoid(double x,double slope)192 double sigmoid(double x, double slope) {
193 return 1.0 / (1.0 + exp(-x * slope));
194 }
195
196 // This function calculates the accumulated density of the normal distribution
197 // For details see G. Marsaglia, Evaluating the Normal Distribution, Journal of Statistical Software, 2004:11(4).
Phi(double x)198 double Phi(double x) {
199 double s=x;
200 double t=0.0;
201 double b=x;
202 double q=x*x;
203 double i=1;
204
205 // if the value is too small or too big, return
206 // 0/1 to avoid long computations
207 if (x < -10.0) {
208 return 0.0;
209 }
210
211 if (x > 10.0) {
212 return 1.0;
213 }
214
215 // t is the value from last iteration
216 // s is the value from the current iteration
217 // iterate until they are equal
218 while(fabs(s-t) >= DBL_MIN) {
219 t=s;
220 i+=2;
221 b*=q/i;
222 s+=b;
223 }
224
225 return 0.5+s*exp(-0.5*q-0.91893853320467274178);
226 }
227
228 // integrates the normal distribution over [low,high]
cumulative_normal(double low,double high,double mu,double sigma)229 double cumulative_normal(double low, double high, double mu, double sigma) {
230 return Phi((high-mu)/sigma) - Phi((low-mu)/sigma);
231 }
232
233 // integrates the normal distribution over [-oo,high]
cumulative_normal_upper(double high,double mu,double sigma)234 double cumulative_normal_upper(double high, double mu, double sigma) {
235 return Phi((high-mu)/sigma);
236 }
237
238
239 // evaluates the density of the normal distribution
normal(double x,double mu,double sigma)240 double normal(double x, double mu,double sigma) {
241 double inner=(x-mu)/sigma;
242 double denom=sigma*sqrt(2*3.14159265358979323846);
243 return exp(-inner*inner/2)/denom;
244 }
245
cumulative_normal_dmu(double low,double high,double mu,double sigma)246 double cumulative_normal_dmu(double low, double high,double mu,double sigma) {
247 return normal(low,mu,sigma) - normal(high,mu,sigma);
248 }
249
cumulative_normal_upper_dmu(double high,double mu,double sigma)250 double cumulative_normal_upper_dmu(double high,double mu,double sigma) {
251 return - normal(high,mu,sigma);
252 }
253
254
cumulative_normal_dsigma(double low,double high,double mu,double sigma)255 double cumulative_normal_dsigma(double low, double high,double mu,double sigma) {
256 return (((mu-high)*normal(high,mu,sigma) - (mu-low)*normal(low,mu,sigma))/sigma);
257 }
258
cumulative_normal_upper_dsigma(double high,double mu,double sigma)259 double cumulative_normal_upper_dsigma(double high,double mu,double sigma) {
260 return (mu-high)*normal(high,mu,sigma);
261 }
262
263
264 // this function parses two strings "$a;$b" and "???_???l$ch$d" where $a-$d are (real) numbers
265 // it is used to parse in the parameters of continues variables from the input file
parse_density_integral_string(char * input,char * variablename)266 density_integral parse_density_integral_string(char *input, char *variablename) {
267 density_integral result;
268 double sigma;
269 int i;
270 char garbage[64], s1[64],s2[64],s3[64],s4[64];
271
272 if(sscanf(input, "%64[^;];%64[^;]", s1,s2) != 2) {
273 fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
274 fprintf(stderr, "The string should contain 2 fields seperated by ; characters.\n");
275 exit(EXIT_FAILURE);
276 }
277
278 if (!getRealNumber(s1, &result.mu)) {
279 fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
280 fprintf(stderr, "%s is not a number\n",s1);
281 exit(EXIT_FAILURE);
282 }
283
284 if (!getRealNumber(s2, &sigma) || sigma<=0.0) {
285 fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
286 fprintf(stderr, "%s is not a number\n",s2);
287 exit(EXIT_FAILURE);
288 }
289 result.log_sigma=log(sigma);
290
291 /* if (result.sigma<=0) { */
292 /* fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string",input); */
293 /* fprintf(stderr, "The value for sigma has to be larger than 0.\n"); */
294
295 /* exit(EXIT_FAILURE); */
296 /* } */
297
298 if (sscanf(variablename,"%64[^lh]l%64[^lh]h%64[^lh]",garbage,s3,s4) != 3) {
299 fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",variablename);
300 fprintf(stderr, "The string should contain 2 fields seperated by ; characters.\n");
301 exit(EXIT_FAILURE);
302 }
303
304 // replace the d by . in s1 and s2
305 for(i=0; s3[i]!='\0' ; i++) {
306 if (s3[i]=='d') {
307 s3[i]='.';
308 }
309 if (s3[i]=='m') {
310 s3[i]='-';
311 }
312 }
313 for(i=0; s4[i]!='\0' ; i++) {
314 if (s4[i]=='d') {
315 s4[i]='.';
316 }
317 if (s4[i]=='m') {
318 s4[i]='-';
319 }
320 }
321
322 if (!getRealNumber(s3, &result.low)) {
323 fprintf(stderr, "Error at parsing the string %s in the function parse_density_integral_string\n",input);
324 fprintf(stderr, "%s is not a number\n",s1);
325 exit(EXIT_FAILURE);
326 }
327
328 if (!getRealNumber(s4, &result.high)) {
329 fprintf(stderr, "Error ar parsing the string %s in the function parse_density_integral_string\n",input);
330 fprintf(stderr, "%s is not a number\n",s1);
331 exit(EXIT_FAILURE);
332 }
333
334
335 if (result.low>result.high) {
336 fprintf(stderr, "Error ar parsing the string %s in the function parse_density_integral_string\n",input);
337 fprintf(stderr, "The value for low has to be larger than then value for high.\n");
338 fprintf(stderr, " was [%f, %f]\n",result.low, result.high);
339 fprintf(stderr, " input %s \n",input);
340 fprintf(stderr, " variablename %s \n",variablename);
341
342 exit(EXIT_FAILURE);
343 }
344
345
346 return result;
347 }
348