--- /dev/null
+/*******************************************************************************\r
+ * Copyright (c) 2013 Association for Decentralized Information Management in\r
+ * Industry THTH ry.\r
+ * All rights reserved. This program and the accompanying materials\r
+ * are made available under the terms of the Eclipse Public License v1.0\r
+ * which accompanies this distribution, and is available at\r
+ * http://www.eclipse.org/legal/epl-v10.html\r
+ *\r
+ * Contributors:\r
+ * Semantum Oy - initial API and implementation\r
+ * VTT Technical Research Centre of Finland\r
+ *******************************************************************************/\r
+package fi.semantum.sysdyn.solver;\r
+\r
+public class NormalDistribution {\r
+\r
+ public static double uniformVariableToNormal(double min, double max, double mean, double stdev, double prob) {\r
+ double minProb = cdf(min, mean, stdev, false);\r
+ double maxProb = cdf(max, mean, stdev, false);\r
+ return quantile(minProb + (maxProb - minProb) * prob, mean, stdev);\r
+ }\r
+\r
+ //The following is from BEAST-MCMC (https://code.google.com/p/beast-mcmc/)\r
+ \r
+ /*\r
+ * NormalDistribution.java\r
+ *\r
+ * Copyright (c) 2002-2011 Alexei Drummond, Andrew Rambaut and Marc Suchard\r
+ *\r
+ * This file is part of BEAST.\r
+ * See the NOTICE file distributed with this work for additional\r
+ * information regarding copyright ownership and licensing.\r
+ *\r
+ * BEAST is free software; you can redistribute it and/or modify\r
+ * it under the terms of the GNU Lesser General Public License as\r
+ * published by the Free Software Foundation; either version 2\r
+ * of the License, or (at your option) any later version.\r
+ *\r
+ * BEAST is distributed in the hope that it will be useful,\r
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\r
+ * GNU Lesser General Public License for more details.\r
+ *\r
+ * You should have received a copy of the GNU Lesser General Public\r
+ * License along with BEAST; if not, write to the\r
+ * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,\r
+ * Boston, MA 02110-1301 USA\r
+ */\r
+ \r
+ /**\r
+ * quantiles (=inverse cumulative density function)\r
+ *\r
+ * @param z argument\r
+ * @param m mean\r
+ * @param sd standard deviation\r
+ * @return icdf at z\r
+ */\r
+ public static double quantile(double z, double m, double sd) {\r
+ return m + Math.sqrt(2.0) * sd * inverseErf(2.0 * z - 1.0);\r
+ }\r
+ \r
+ /** A more accurate and faster implementation of the cdf (taken from function pnorm in the R statistical language)\r
+ * This implementation has discrepancies depending on the programming language and system architecture\r
+ * In Java, returned values become zero once z reaches -37.5193 exactly on the machine tested\r
+ * In the other implementation, the returned value 0 at about z = -8\r
+ * In C, this 0 value is reached approximately z = -37.51938\r
+ *\r
+ * Will later need to be optimised for BEAST\r
+ *\r
+ * @param x argument\r
+ * @param mu mean\r
+ * @param sigma standard deviation\r
+ * @param log_p is p logged\r
+ * @return cdf at x\r
+ */\r
+ public static double cdf(double x, double mu, double sigma, boolean log_p) {\r
+ boolean i_tail=false;\r
+ double p, cp = Double.NaN;\r
+\r
+ if(Double.isNaN(x) || Double.isNaN(mu) || Double.isNaN(sigma)) {\r
+ return Double.NaN;\r
+ }\r
+ if(Double.isInfinite(x) && mu == x) { /* x-mu is NaN */\r
+ return Double.NaN;\r
+ }\r
+ if (sigma <= 0) {\r
+ if(sigma < 0) {\r
+ return Double.NaN;\r
+ }\r
+ return (x < mu) ? 0.0 : 1.0;\r
+ }\r
+ p = (x - mu) / sigma;\r
+ if(Double.isInfinite(p)) {\r
+ return (x < mu) ? 0.0 : 1.0;\r
+ }\r
+ x = p;\r
+ if(Double.isNaN(x)) {\r
+ return Double.NaN;\r
+ }\r
+\r
+ double xden, xnum, temp, del, eps, xsq, y;\r
+ int i;\r
+ boolean lower, upper;\r
+ eps = DBL_EPSILON * 0.5;\r
+ lower = !i_tail;\r
+ upper = i_tail;\r
+\r
+ y = Math.abs(x);\r
+ if (y <= 0.67448975) { /* Normal.quantile(3/4, 1, 0) = 0.67448975 */\r
+ if (y > eps) {\r
+ xsq = x * x;\r
+ xnum = a[4] * xsq;\r
+ xden = xsq;\r
+ for (i = 0; i < 3; i++) {\r
+ xnum = (xnum + a[i]) * xsq;\r
+ xden = (xden + b[i]) * xsq;\r
+ }\r
+ }\r
+ else {\r
+ xnum = xden = 0.0;\r
+ }\r
+ temp = x * (xnum + a[3]) / (xden + b[3]);\r
+ if(lower) {\r
+ p = 0.5 + temp;\r
+ }\r
+ if(upper) {\r
+ cp = 0.5 - temp;\r
+ }\r
+ if(log_p) {\r
+ if(lower) {\r
+ p = Math.log(p);\r
+ }\r
+ if(upper) {\r
+ cp = Math.log(cp);\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ else if (y <= M_SQRT_32) {\r
+ /* Evaluate pnorm for 0.67448975 = Normal.quantile(3/4, 1, 0) < |x| <= sqrt(32) ~= 5.657 */\r
+\r
+ xnum = c[8] * y;\r
+ xden = y;\r
+ for (i = 0; i < 7; i++) {\r
+ xnum = (xnum + c[i]) * y;\r
+ xden = (xden + d[i]) * y;\r
+ }\r
+ temp = (xnum + c[7]) / (xden + d[7]);\r
+\r
+ //do_del(y);\r
+ //swap_tail;\r
+ //#define do_del(X) \\r
+ xsq = ((int) (y * CUTOFF)) * 1.0 / CUTOFF;\r
+ del = (y - xsq) * (y + xsq);\r
+ if(log_p) {\r
+ p = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);\r
+ if((lower && x > 0.0) || (upper && x <= 0.0)) {\r
+ cp = Math.log(1.0-Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp);\r
+ }\r
+ }\r
+ else {\r
+ p = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;\r
+ cp = 1.0 - p;\r
+ }\r
+ //#define swap_tail \\r
+ if (x > 0.0) {\r
+ temp = p;\r
+ if(lower) {\r
+ p = cp;\r
+ }\r
+ cp = temp;\r
+ }\r
+ }\r
+ /* else |x| > sqrt(32) = 5.657 :\r
+ * the next two case differentiations were really for lower=T, log=F\r
+ * Particularly *not* for log_p !\r
+ * Cody had (-37.5193 < x && x < 8.2924) ; R originally had y < 50\r
+ * Note that we do want symmetry(0), lower/upper -> hence use y\r
+ */\r
+ else if(log_p || (lower && -37.5193 < x && x < 8.2924)\r
+ || (upper && -8.2924 < x && x < 37.5193)) {\r
+\r
+ /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */\r
+ xsq = 1.0 / (x * x);\r
+ xnum = p_[5] * xsq;\r
+ xden = xsq;\r
+ for (i = 0; i < 4; i++) {\r
+ xnum = (xnum + p_[i]) * xsq;\r
+ xden = (xden + q[i]) * xsq;\r
+ }\r
+ temp = xsq * (xnum + p_[4]) / (xden + q[4]);\r
+ temp = (M_1_SQRT_2PI - temp) / y;\r
+\r
+ //do_del(x);\r
+ xsq = ((int) (x * CUTOFF)) * 1.0 / CUTOFF;\r
+ del = (x - xsq) * (x + xsq);\r
+ if(log_p) {\r
+ p = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);\r
+ if((lower && x > 0.0) || (upper && x <= 0.0)) {\r
+ cp = Math.log(1.0-Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp);\r
+ }\r
+ }\r
+ else {\r
+ p = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;\r
+ cp = 1.0 - p;\r
+ }\r
+ //swap_tail;\r
+ if (x > 0.0) {\r
+ temp = p;\r
+ if(lower) {\r
+ p = cp;\r
+ }\r
+ cp = temp;\r
+ }\r
+ }\r
+ else { /* no log_p , large x such that probs are 0 or 1 */\r
+ if(x > 0) {\r
+ p = 1.0;\r
+ cp = 0.0;\r
+ }\r
+ else {\r
+ p = 0.0;\r
+ cp = 1.0;\r
+ }\r
+ }\r
+ return p;\r
+\r
+ }\r
+\r
+ // Private\r
+\r
+ protected double m, sd;\r
+\r
+ private static final double[] a = {\r
+ 2.2352520354606839287,\r
+ 161.02823106855587881,\r
+ 1067.6894854603709582,\r
+ 18154.981253343561249,\r
+ 0.065682337918207449113\r
+ };\r
+ private static final double[] b = {\r
+ 47.20258190468824187,\r
+ 976.09855173777669322,\r
+ 10260.932208618978205,\r
+ 45507.789335026729956\r
+ };\r
+ private static final double[] c = {\r
+ 0.39894151208813466764,\r
+ 8.8831497943883759412,\r
+ 93.506656132177855979,\r
+ 597.27027639480026226,\r
+ 2494.5375852903726711,\r
+ 6848.1904505362823326,\r
+ 11602.651437647350124,\r
+ 9842.7148383839780218,\r
+ 1.0765576773720192317e-8\r
+ };\r
+ private static final double[] d = {\r
+ 22.266688044328115691,\r
+ 235.38790178262499861,\r
+ 1519.377599407554805,\r
+ 6485.558298266760755,\r
+ 18615.571640885098091,\r
+ 34900.952721145977266,\r
+ 38912.003286093271411,\r
+ 19685.429676859990727\r
+ };\r
+ private static final double[] p_ = {\r
+ 0.21589853405795699,\r
+ 0.1274011611602473639,\r
+ 0.022235277870649807,\r
+ 0.001421619193227893466,\r
+ 2.9112874951168792e-5,\r
+ 0.02307344176494017303\r
+ };\r
+ private static final double[] q = {\r
+ 1.28426009614491121,\r
+ 0.468238212480865118,\r
+ 0.0659881378689285515,\r
+ 0.00378239633202758244,\r
+ 7.29751555083966205e-5\r
+ };\r
+\r
+ private static final int CUTOFF = 16; /* Cutoff allowing exact "*" and "/" */\r
+\r
+ private static final double M_SQRT_32 = 5.656854249492380195206754896838; /* The square root of 32 */\r
+ private static final double M_1_SQRT_2PI = 0.398942280401432677939946059934;\r
+ private static final double DBL_EPSILON = 2.2204460492503131e-016;\r
+ \r
+ /*\r
+ * ErrorFunction.java\r
+ *\r
+ * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut\r
+ *\r
+ * This file is part of BEAST.\r
+ * See the NOTICE file distributed with this work for additional\r
+ * information regarding copyright ownership and licensing.\r
+ *\r
+ * BEAST is free software; you can redistribute it and/or modify\r
+ * it under the terms of the GNU Lesser General Public License as\r
+ * published by the Free Software Foundation; either version 2\r
+ * of the License, or (at your option) any later version.\r
+ *\r
+ * BEAST is distributed in the hope that it will be useful,\r
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\r
+ * GNU Lesser General Public License for more details.\r
+ *\r
+ * You should have received a copy of the GNU Lesser General Public\r
+ * License along with BEAST; if not, write to the\r
+ * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,\r
+ * Boston, MA 02110-1301 USA\r
+ */\r
+ \r
+ /**\r
+ * inverse error function\r
+ *\r
+ * @param z argument\r
+ *\r
+ * @return function value\r
+ */\r
+ public static double inverseErf(double z)\r
+ {\r
+ return pointNormal(0.5*z+0.5)/Math.sqrt(2.0);\r
+ }\r
+ \r
+ // Private\r
+ \r
+ // Returns z so that Prob{x<z}=prob where x ~ N(0,1) and (1e-12) < prob<1-(1e-12)\r
+ private static double pointNormal(double prob)\r
+ {\r
+ // Odeh RE & Evans JO (1974) The percentage points of the normal distribution.\r
+ // Applied Statistics 22: 96-97 (AS70)\r
+ \r
+ // Newer methods:\r
+ // Wichura MJ (1988) Algorithm AS 241: the percentage points of the\r
+ // normal distribution. 37: 477-484.\r
+ // Beasley JD & Springer SG (1977). Algorithm AS 111: the percentage \r
+ // points of the normal distribution. 26: 118-121.\r
+ \r
+ double a0 = -0.322232431088, a1 = -1, a2 = -0.342242088547, a3 = -0.0204231210245;\r
+ double a4 = -0.453642210148e-4, b0 = 0.0993484626060, b1 = 0.588581570495;\r
+ double b2 = 0.531103462366, b3 = 0.103537752850, b4 = 0.0038560700634;\r
+ double y, z = 0, p = prob, p1;\r
+ \r
+ p1 = (p < 0.5 ? p : 1-p);\r
+ if (p1 < 1e-20)\r
+ {\r
+ new IllegalArgumentException("Argument prob out of range");\r
+ }\r
+ \r
+ y = Math.sqrt(Math.log(1/(p1*p1))); \r
+ z = y + ((((y*a4+a3)*y+a2)*y+a1)*y+a0)/((((y*b4+b3)*y+b2)*y+b1)*y+b0);\r
+ return (p < 0.5 ? -z : z);\r
+ }\r
+\r
+}\r
Code Review / simantics / sysdyn.git / commitdiff