Fix loop (by switching method).

2002-07-09  Morten Welinder  <terra@diku.dk>

	* src/mathfunc.c (random_poisson): Fix loop (by switching method).

ChangeLog

+2002-07-09  Morten Welinder  <terra@diku.dk>
+
+	* src/mathfunc.c (random_poisson): Fix loop (by switching method).
+
 2002-07-09  Jody Goldberg <jody@gnome.org>
 
 	* plugins/Makefile.am : excel is no longer conditional

NEWS

@@ -56,6 +56,7 @@ Morten:
 	* Update mathfunc.c to R 1.5.1.
 	* Improve precision of HYPGEOMDIST and NEGBINOMDIST.
 	* Fix PERMUT crash.
+	* Fix RANDPOISSON loop and performance.
 
 Jon Kre:
 	* Use GsfInput in the Python plugin loader

OChangeLog-2003-12-23

+2002-07-09  Morten Welinder  <terra@diku.dk>
+
+	* src/mathfunc.c (random_poisson): Fix loop (by switching method).
+
 2002-07-09  Jody Goldberg <jody@gnome.org>
 
 	* plugins/Makefile.am : excel is no longer conditional

OChangeLog-2005-11-14

+2002-07-09  Morten Welinder  <terra@diku.dk>
+
+	* src/mathfunc.c (random_poisson): Fix loop (by switching method).
+
 2002-07-09  Jody Goldberg <jody@gnome.org>
 
 	* plugins/Makefile.am : excel is no longer conditional

src/mathfunc.c

@@ -854,6 +854,106 @@ gnum_float ppois(gnum_float x, gnum_float lambda, gboolean lower_tail, gboolean
     return pgamma(lambda, x + 1, 1., !lower_tail, log_p);
 }
 
+/* ------------------------------------------------------------------------ */
+/* Imported src/nmath/qpois.c from R.  */
+/*
+ *  Mathlib : A C Library of Special Functions
+ *  Copyright (C) 1998 Ross Ihaka
+ *  Copyright (C) 2000 The R Development Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
+ *
+ *  DESCRIPTION
+ *
+ *	The quantile function of the Poisson distribution.
+ *
+ *  METHOD
+ *
+ *	Uses the Cornish-Fisher Expansion to include a skewness
+ *	correction to a normal approximation.  This gives an
+ *	initial value which never seems to be off by more than
+ *	1 or 2.	 A search is then conducted of values close to
+ *	this initial start point.
+ */
+
+
+gnum_float qpois(gnum_float p, gnum_float lambda, gboolean lower_tail, gboolean log_p)
+{
+    gnum_float mu, sigma, gamma, z, y;
+#ifdef IEEE_754
+    if (isnangnum(p) || isnangnum(lambda))
+	return p + lambda;
+#endif
+    if(!finitegnum(lambda))
+	ML_ERR_return_NAN;
+    R_Q_P01_check(p);
+    if(lambda < 0) ML_ERR_return_NAN;
+
+    if (p == R_DT_0) return 0;
+    if (p == R_DT_1) return ML_POSINF;
+
+    if(lambda == 0) return 0;
+
+    mu = lambda;
+    sigma = sqrtgnum(lambda);
+    gamma = sigma;
+
+    /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
+     * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
+    if(!lower_tail || log_p) {
+	p = R_DT_qIv(p); /* need check again (cancellation!): */
+	if (p == 0.) return 0;
+	if (p == 1.) return ML_POSINF;
+    }
+    /* temporary hack --- FIXME --- */
+    if (p + 1.01*GNUM_EPSILON >= 1.) return ML_POSINF;
+
+    /* y := approx.value (Cornish-Fisher expansion) :  */
+    z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
+    y = floorgnum(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);
+
+    z = ppois(y, lambda, /*lower_tail*/TRUE, /*log_p*/FALSE);
+
+    /* fuzz to ensure left continuity; 1 - 1e-7 may lose too much : */
+    p *= 1 - 64*GNUM_EPSILON;
+
+/*-- Fixme, here y can be way off --
+  should use interval search instead of primitive stepping down or up */
+
+#ifdef maybe_future
+    if((lower_tail && z >= p) || (!lower_tail && z <= p)) {
+#else
+    if(z >= p) {
+#endif
+			/* search to the left */
+	for(;;) {
+	    if(y == 0 ||
+	       (z = ppois(y - 1, lambda, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
+		return y;
+	    y = y - 1;
+	}
+    }
+    else {		/* search to the right */
+	for(;;) {
+	    y = y + 1;
+	    if((z = ppois(y, lambda, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p)
+		return y;
+	}
+    }
+}
+
 /* ------------------------------------------------------------------------ */
 /* Imported src/nmath/stirlerr.c from R.  */
 /*
@@ -4275,19 +4375,11 @@ random_normal (void)
 gnum_float
 random_poisson (gnum_float lambda)
 {
-        gnum_float x = expgnum (-lambda);
-	gnum_float r = random_01 ();
-	gnum_float t = x;
-	gnum_float i = 0;
-
-	/* FIXME: Looks dubious, performanc-wise.  */
-	while (r > t) {
-	      x *= lambda / (i + 1);
-	      i += 1;
-	      t += x;
-	}
-
-	return i;
+	/*
+	 * This may not be optimal code, but it sure is easy to
+	 * understand compared to R's code.
+	 */
+	return qpois (random_01 (), lambda, TRUE, FALSE);
 }
 
 /*

src/mathfunc.h

@@ -87,6 +87,7 @@ gnum_float pweibull (gnum_float x, gnum_float shape, gnum_float scale, gboolean
 /* The Poisson distribution.  */
 gnum_float dpois (gnum_float x, gnum_float lambda, gboolean give_log);
 gnum_float ppois (gnum_float x, gnum_float lambda, gboolean lower_tail, gboolean log_p);
+gnum_float qpois (gnum_float p, gnum_float lambda, gboolean lower_tail, gboolean log_p);
 
 /* The exponential distribution.  */
 gnum_float dexp (gnum_float x, gnum_float scale, gboolean give_log);

tools/ChangeLog

+2002-07-09  Morten Welinder  <terra@diku.dk>
+
+	* import-R (import_file): Import qpois.c also.
+
 2002-06-21  Morten Welinder  <terra@diku.dk>
 
 	* import-R: Import also dnbinom.c and dhyper.c

tools/import-R

@@ -11,6 +11,7 @@ my @files =
      "plnorm.c",
      "qlnorm.c",
      "ppois.c",
+     "qpois.c",
      "stirlerr.c",
      "bd0.c",
      "dpois.c",
@@ -121,6 +122,10 @@ sub import_file {
 	    s/int log_p/gboolean log_p/g;
 	    s/int lower_tail/gboolean lower_tail/g;
 
+	    if (/^.*char\s+\*vmax;\s*$/) {
+		$_ = "#ifndef MATHLIB_STANDALONE\n$_\n#endif";
+	    }
+
 	    $incomment = 1 if m|/\*$|;
 	}