[gnumeric] R.PTUKEY: New function.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] R.PTUKEY: New function.
- Date: Sun, 19 May 2013 00:43:01 +0000 (UTC)
commit e97df747d7a1e5eb80417814ddb234ed00b84048
Author: Morten Welinder <terra gnome org>
Date: Sat May 18 20:37:12 2013 -0400
R.PTUKEY: New function.
This is imported from R with a few improvements.
ChangeLog | 5 +
NEWS | 1 +
plugins/fn-r/functions.c | 35 ++++
plugins/fn-r/generate | 16 ++-
plugins/fn-r/plugin.xml.in | 1 +
samples/ptukey.gnumeric | Bin 0 -> 232075 bytes
src/mathfunc.c | 472 ++++++++++++++++++++++++++++++++++++++++++++
src/mathfunc.h | 5 +
8 files changed, 533 insertions(+), 2 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index e1b2a3d..0a3a237 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-05-18 Morten Welinder <terra gnome org>
+
+ * src/mathfunc.c (pnorm2): New function.
+ (R_ptukey): New function imported from R with local improvements.
+
2013-05-15 Andreas J. Guelzow <aguelzow pyrshep ca>
* src/print.c (gnm_create_widget_cb): by default observe
diff --git a/NEWS b/NEWS
index ded9624..837692d 100644
--- a/NEWS
+++ b/NEWS
@@ -18,6 +18,7 @@ Morten:
* Improve LINEST. [#551457]
* Improve statfuns tests. [#700294]
* Improve mathfuns tests. [#700295]
+ * Add new R.PTUKEY function. [#700132]
--------------------------------------------------------------------------
Gnumeric 1.12.2
diff --git a/plugins/fn-r/functions.c b/plugins/fn-r/functions.c
index 066bd4b..b25b8c7 100644
--- a/plugins/fn-r/functions.c
+++ b/plugins/fn-r/functions.c
@@ -1152,6 +1152,34 @@ gnumeric_r_pcauchy (GnmFuncEvalInfo *ei, GnmValue const * const *args)
/* ------------------------------------------------------------------------- */
+static GnmFuncHelp const help_r_ptukey[] = {
+ { GNM_FUNC_HELP_NAME, F_("R.PTUKEY:cumulative distribution function of the Studentized range
distribution") },
+ { GNM_FUNC_HELP_ARG, F_("x:observation") },
+ { GNM_FUNC_HELP_ARG, F_("nmeans:the number of means") },
+ { GNM_FUNC_HELP_ARG, F_("df:the number of degrees of freedom of the distribution") },
+ { GNM_FUNC_HELP_ARG, F_("nranges:the number of ranges; default is 1") },
+ { GNM_FUNC_HELP_ARG, F_("lower_tail:if true (the default), the lower tail of the distribution is
considered") },
+ { GNM_FUNC_HELP_ARG, F_("log_p:if true, log of the probability is used") },
+ { GNM_FUNC_HELP_DESCRIPTION, F_("This function returns the cumulative distribution function of the
Studentized range distribution.") },
+ { GNM_FUNC_HELP_END }
+};
+
+static GnmValue *
+gnumeric_r_ptukey (GnmFuncEvalInfo *ei, GnmValue const * const *args)
+{
+ gnm_float x = value_get_as_float (args[0]);
+ gnm_float nmeans = value_get_as_float (args[1]);
+ gnm_float df = value_get_as_float (args[2]);
+ gnm_float nranges = args[3] ? value_get_as_float (args[3]) : 1;
+ gboolean lower_tail = args[4] ? value_get_as_checked_bool (args[4]) : TRUE;
+ gboolean log_p = args[5] ? value_get_as_checked_bool (args[5]) : FALSE;
+
+ return value_new_float (ptukey (x, nmeans, df, nranges, lower_tail, log_p));
+}
+
+/* ------------------------------------------------------------------------- */
+
+
static GnmFuncHelp const help_r_qcauchy[] = {
{ GNM_FUNC_HELP_NAME, F_("R.QCAUCHY:probability quantile function of the Cauchy distribution") },
{ GNM_FUNC_HELP_ARG, F_("p:probability") },
@@ -1664,6 +1692,13 @@ GnmFuncDescriptor const stat_functions[] = {
GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
},
{
+ "r.ptukey",
+ "fff|fbb",
+ help_r_ptukey,
+ gnumeric_r_ptukey, NULL, NULL, NULL,
+ GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
+ },
+ {
"r.qcauchy",
"fff|bb",
help_r_qcauchy,
diff --git a/plugins/fn-r/generate b/plugins/fn-r/generate
index c136a8d..e01a1c8 100644
--- a/plugins/fn-r/generate
+++ b/plugins/fn-r/generate
@@ -8,6 +8,7 @@ sub distribution;
my %funcs = ();
my %argoverride = ();
+my %defaults;
{
my $of = "of the distribution";
@@ -132,6 +133,15 @@ my %argoverride = ();
'skew-t',
({ 'n' => "the number of degrees of freedom $of",
@common })];
+
+ $funcs{'ptukey'} =
+ [\&distribution,
+ 'Studentized range',
+ ({ 'nranges' => "the number of ranges; default is 1",
+ 'nmeans' => "the number of means",
+ 'df' => "the number of degrees of freedom $of",
+ @common })];
+ $defaults{'ptukey:nranges'} = 1;
}
@@ -281,7 +291,9 @@ sub distribution {
$seealso .= ",R.P$F1" if $func !~ /^p/ && exists $funcs{"p$f1"};
$seealso .= ",R.Q$F1" if $func !~ /^q/ && exists $funcs{"q$f1"};
$seealso =~ s/^,\s*//;
- &emit ("\t{ GNM_FUNC_HELP_SEEALSO, \"$seealso\" },\n");
+ if ($seealso) {
+ &emit ("\t{ GNM_FUNC_HELP_SEEALSO, \"$seealso\" },\n");
+ }
&emit ("\t{ GNM_FUNC_HELP_END }\n" .
"};\n\n");
@@ -297,7 +309,7 @@ sub distribution {
my $n = 0;
foreach (@args) {
my ($type,$name) = @{ $_ };
- my $def = undef;
+ my $def = $defaults{"$func:$name"};
$def = 'TRUE' if $name eq 'lower_tail';
$def = 'FALSE' if $name eq 'give_log' || $name eq 'log_p';
diff --git a/plugins/fn-r/plugin.xml.in b/plugins/fn-r/plugin.xml.in
index 901f254..f3eb79f 100644
--- a/plugins/fn-r/plugin.xml.in
+++ b/plugins/fn-r/plugin.xml.in
@@ -44,6 +44,7 @@
<function name="r.psnorm"/>
<function name="r.pst"/>
<function name="r.pt"/>
+ <function name="r.ptukey"/>
<function name="r.pweibull"/>
<function name="r.qbeta"/>
<function name="r.qbinom"/>
diff --git a/samples/ptukey.gnumeric b/samples/ptukey.gnumeric
new file mode 100644
index 0000000..69ba060
Binary files /dev/null and b/samples/ptukey.gnumeric differ
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 308fabb..7035487 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -5109,8 +5109,433 @@ gnm_float qgeom(gnm_float p, gnm_float prob, gboolean lower_tail, gboolean log_p
}
/* ------------------------------------------------------------------------ */
+/* Imported src/nmath/ptukey.c from R -- by hand for now. */
+/*
+ * Mathlib : A C Library of Special Functions
+ * Copyright (C) 1998 Ross Ihaka
+ * Copyright (C) 2000--2007 The R 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, a copy is available at
+ * http://www.r-project.org/Licenses/
+ *
+ * SYNOPSIS
+ *
+ * #include <Rmath.h>
+ * double ptukey(q, rr, cc, df, lower_tail, log_p);
+ *
+ * DESCRIPTION
+ *
+ * Computes the probability that the maximum of rr studentized
+ * ranges, each based on cc means and with df degrees of freedom
+ * for the standard error, is less than q.
+ *
+ * The algorithm is based on that of the reference.
+ *
+ * REFERENCE
+ *
+ * Copenhaver, Margaret Diponzio & Holland, Burt S.
+ * Multiple comparisons of simple effects in
+ * the two-way analysis of variance with fixed effects.
+ * Journal of Statistical Computation and Simulation,
+ * Vol.30, pp.1-15, 1988.
+ */
+
+static gnm_float ptukey_wprob(gnm_float w, gnm_float rr, gnm_float cc)
+{
+/* wprob() :
+
+ This function calculates probability integral of Hartley's
+ form of the range.
+
+ w = value of range
+ rr = no. of rows or groups
+ cc = no. of columns or treatments
+ ir = error flag = 1 if pr_w probability > 1
+ pr_w = returned probability integral from (0, w)
+
+ program will not terminate if ir is raised.
+
+ bb = upper limit of legendre integration
+ iMax = maximum acceptable value of integral
+ nleg = order of legendre quadrature
+ ihalf = int ((nleg + 1) / 2)
+ wlar = value of range above which wincr1 intervals are used to
+ calculate second part of integral,
+ else wincr2 intervals are used.
+ C1, C2, C3 = values which are used as cutoffs for terminating
+ or modifying a calculation.
+
+ M_1_SQRT_2PI = 1 / sqrt(2 * pi); from abramowitz & stegun, p. 3.
+ M_SQRT2 = sqrt(2)
+ xleg = legendre 12-point nodes
+ aleg = legendre 12-point coefficients
+ */
+#define nleg 12
+#define ihalf 6
+
+ /* looks like this is suboptimal for gnm_float precision.
+ (see how C1-C3 are used) <MM>
+ */
+ /* const gnm_float iMax = 1.; not used if = 1*/
+ const static gnm_float C1 = -30.;
+ const static gnm_float C2 = -50.;
+ const static gnm_float C3 = 60.;
+ const static gnm_float bb = 8.;
+ const static gnm_float wlar = 3.;
+ const static gnm_float wincr1 = 2.;
+ const static gnm_float wincr2 = 3.;
+ const static gnm_float xleg[ihalf] = {
+ GNM_const(0.981560634246719250690549090149),
+ GNM_const(0.904117256370474856678465866119),
+ GNM_const(0.769902674194304687036893833213),
+ GNM_const(0.587317954286617447296702418941),
+ GNM_const(0.367831498998180193752691536644),
+ GNM_const(0.125233408511468915472441369464)
+ };
+ const static gnm_float aleg[ihalf] = {
+ GNM_const(0.047175336386511827194615961485),
+ GNM_const(0.106939325995318430960254718194),
+ GNM_const(0.160078328543346226334652529543),
+ GNM_const(0.203167426723065921749064455810),
+ GNM_const(0.233492536538354808760849898925),
+ GNM_const(0.249147045813402785000562436043)
+ };
+ gnm_float a, ac, pr_w, b, binc, c, cc1,
+ qexpo, qsqz, rinsum, wi, wincr, xx;
+ /*LDOUBLE*/ gnm_float blb, bub, einsum, elsum;
+ int j, jj;
+
+
+ qsqz = w * 0.5;
+
+ /* if w >= 16 then the integral lower bound (occurs for c=20) */
+ /* is 0.99999999999995 so return a value of 1. */
+
+ if (qsqz >= bb)
+ return 1.0;
+
+ /* find (f(w/2) - 1) ^ cc */
+ /* (first term in integral of hartley's form). */
+
+ pr_w = gnm_erf (qsqz / M_SQRT2gnum);
+
+ /* if pr_w ^ cc < 2e-22 then set pr_w = 0 */
+ if (pr_w >= gnm_exp(C2 / cc))
+ pr_w = gnm_pow(pr_w, cc);
+ else
+ pr_w = 0.0;
+
+ /* if w is large then the second component of the */
+ /* integral is small, so fewer intervals are needed. */
+
+ if (w > wlar)
+ wincr = wincr1;
+ else
+ wincr = wincr2;
+
+ /* find the integral of second term of hartley's form */
+ /* for the integral of the range for equal-length */
+ /* intervals using legendre quadrature. limits of */
+ /* integration are from (w/2, 8). two or three */
+ /* equal-length intervals are used. */
+
+ /* blb and bub are lower and upper limits of integration. */
+
+ blb = qsqz;
+ binc = (bb - qsqz) / wincr;
+ bub = blb + binc;
+ einsum = 0.0;
+
+ /* integrate over each interval */
+
+ cc1 = cc - 1.0;
+ for (wi = 1; wi <= wincr; wi++) {
+ elsum = 0.0;
+ a = (gnm_float)(0.5 * (bub + blb));
+
+ /* legendre quadrature with order = nleg */
+
+ b = (gnm_float)(0.5 * (bub - blb));
+
+ for (jj = 1; jj <= nleg; jj++) {
+ if (ihalf < jj) {
+ j = (nleg - jj) + 1;
+ xx = xleg[j-1];
+ } else {
+ j = jj;
+ xx = -xleg[j-1];
+ }
+ c = b * xx;
+ ac = a + c;
+
+ /* if exp(-qexpo/2) < 9e-14, */
+ /* then doesn't contribute to integral */
+
+ qexpo = ac * ac;
+ if (qexpo > C3)
+ break;
+
+ /* if rinsum ^ (cc-1) < 9e-14, */
+ /* then doesn't contribute to integral */
+
+ // rinsum = (pplus * 0.5) - (pminus * 0.5);
+ rinsum = pnorm2 (ac - w, ac, 0.0, 1.0, FALSE);
+ if (rinsum >= exp(C1 / cc1)) {
+ rinsum = (aleg[j-1] * gnm_exp(-(0.5 * qexpo))) * gnm_pow(rinsum, cc1);
+ elsum += rinsum;
+ }
+ }
+ elsum *= (((2.0 * b) * cc) * M_1_SQRT_2PI);
+ einsum += elsum;
+ blb = bub;
+ bub += binc;
+ }
+
+ /* if pr_w ^ rr < 9e-14, then return 0 */
+ pr_w += (gnm_float) einsum;
+ if (pr_w <= gnm_exp(C1 / rr))
+ return 0.;
+
+ pr_w = gnm_pow(pr_w, rr);
+ if (pr_w >= 1.)/* 1 was iMax was eps */
+ return 1.;
+ return pr_w;
+} /* wprob() */
+
+
+static gnm_float
+R_ptukey(gnm_float q, gnm_float rr, gnm_float cc, gnm_float df,
+ gboolean lower_tail, gboolean log_p)
+{
+/* function ptukey() [was qprob() ]:
+
+ q = value of studentized range
+ rr = no. of rows or groups
+ cc = no. of columns or treatments
+ df = degrees of freedom of error term
+ ir[0] = error flag = 1 if wprob probability > 1
+ ir[1] = error flag = 1 if qprob probability > 1
+
+ qprob = returned probability integral over [0, q]
+
+ The program will not terminate if ir[0] or ir[1] are raised.
+
+ All references in wprob to Abramowitz and Stegun
+ are from the following reference:
+
+ Abramowitz, Milton and Stegun, Irene A.
+ Handbook of Mathematical Functions.
+ New York: Dover publications, Inc. (1970).
+
+ All constants taken from this text are
+ given to 25 significant digits.
+
+ nlegq = order of legendre quadrature
+ ihalfq = int ((nlegq + 1) / 2)
+ eps = max. allowable value of integral
+ eps1 & eps2 = values below which there is
+ no contribution to integral.
+
+ d.f. <= dhaf: integral is divided into ulen1 length intervals. else
+ d.f. <= dquar: integral is divided into ulen2 length intervals. else
+ d.f. <= deigh: integral is divided into ulen3 length intervals. else
+ d.f. <= dlarg: integral is divided into ulen4 length intervals.
+
+ d.f. > dlarg: the range is used to calculate integral.
+
+ M_LN2 = log(2)
+
+ xlegq = legendre 16-point nodes
+ alegq = legendre 16-point coefficients
+
+ The coefficients and nodes for the legendre quadrature used in
+ qprob and wprob were calculated using the algorithms found in:
+
+ Stroud, A. H. and Secrest, D.
+ Gaussian Quadrature Formulas.
+ Englewood Cliffs,
+ New Jersey: Prentice-Hall, Inc, 1966.
+
+ All values matched the tables (provided in same reference)
+ to 30 significant digits.
+
+ f(x) = .5 + erf(x / sqrt(2)) / 2 for x > 0
+
+ f(x) = erfc( -x / sqrt(2)) / 2 for x < 0
+
+ where f(x) is standard normal c. d. f.
+
+ if degrees of freedom large, approximate integral
+ with range distribution.
+ */
+#define nlegq 16
+#define ihalfq 8
+
+/* const gnm_float eps = 1.0; not used if = 1 */
+ const static gnm_float eps1 = -30.0;
+ const static gnm_float eps2 = GNM_const(1.0e-14);
+ const static gnm_float dhaf = 100.0;
+ const static gnm_float dquar = 800.0;
+ const static gnm_float deigh = 5000.0;
+ const static gnm_float dlarg = 25000.0;
+ const static gnm_float ulen1 = 1.0;
+ const static gnm_float ulen2 = 0.5;
+ const static gnm_float ulen3 = 0.25;
+ const static gnm_float ulen4 = 0.125;
+ const static gnm_float xlegq[ihalfq] = {
+ GNM_const(0.989400934991649932596154173450),
+ GNM_const(0.944575023073232576077988415535),
+ GNM_const(0.865631202387831743880467897712),
+ GNM_const(0.755404408355003033895101194847),
+ GNM_const(0.617876244402643748446671764049),
+ GNM_const(0.458016777657227386342419442984),
+ GNM_const(0.281603550779258913230460501460),
+ GNM_const(0.950125098376374401853193354250e-1)
+ };
+ const static gnm_float alegq[ihalfq] = {
+ GNM_const(0.271524594117540948517805724560e-1),
+ GNM_const(0.622535239386478928628438369944e-1),
+ GNM_const(0.951585116824927848099251076022e-1),
+ GNM_const(0.124628971255533872052476282192),
+ GNM_const(0.149595988816576732081501730547),
+ GNM_const(0.169156519395002538189312079030),
+ GNM_const(0.182603415044923588866763667969),
+ GNM_const(0.189450610455068496285396723208)
+ };
+ gnm_float ans, f2, f21, f2lf, ff4, otsum, qsqz, rotsum, t1, twa1, ulen, wprb;
+ int i, j, jj;
+
+#ifdef IEEE_754
+ if (gnm_isnan(q) || gnm_isnan(rr) || gnm_isnan(cc) || gnm_isnan(df))
+ ML_ERR_return_NAN;
+#endif
+
+ if (q <= 0)
+ return R_DT_0;
+
+ /* FIXME: Special case for cc==2: we have explicit formula */
+
+
+ /* df must be > 1 */
+ /* there must be at least two values */
+
+ if (df < 2 || rr < 1 || cc < 2) ML_ERR_return_NAN;
+
+ if(!gnm_finite(q))
+ return R_DT_1;
+
+ if (df > dlarg)
+ return R_DT_val(ptukey_wprob(q, rr, cc));
+
+ /* calculate leading constant */
+
+ f2 = df * 0.5;
+ /* gnm_lgamma(u) = log(gamma(u)) */
+ f2lf = ((f2 * gnm_log(df)) - (df * M_LN2gnum)) - gnm_lgamma(f2);
+ f21 = f2 - 1.0;
+
+ /* integral is divided into unit, half-unit, quarter-unit, or */
+ /* eighth-unit length intervals depending on the value of the */
+ /* degrees of freedom. */
+
+ ff4 = df * 0.25;
+ if (df <= dhaf) ulen = ulen1;
+ else if (df <= dquar) ulen = ulen2;
+ else if (df <= deigh) ulen = ulen3;
+ else ulen = ulen4;
+
+ f2lf += gnm_log(ulen);
+
+ /* integrate over each subinterval */
+
+ ans = 0.0;
+
+ for (i = 1; i <= 50; i++) {
+ otsum = 0.0;
+
+ /* legendre quadrature with order = nlegq */
+ /* nodes (stored in xlegq) are symmetric around zero. */
+
+ twa1 = (2 * i - 1) * ulen;
+
+ for (jj = 1; jj <= nlegq; jj++) {
+ if (ihalfq < jj) {
+ j = jj - ihalfq - 1;
+ t1 = (f2lf + (f21 * gnm_log(twa1 + (xlegq[j] * ulen))))
+ - (((xlegq[j] * ulen) + twa1) * ff4);
+ } else {
+ j = jj - 1;
+ t1 = (f2lf + (f21 * gnm_log(twa1 - (xlegq[j] * ulen))))
+ + (((xlegq[j] * ulen) - twa1) * ff4);
+
+ }
+
+ /* if exp(t1) < 9e-14, then doesn't contribute to integral */
+ if (t1 >= eps1) {
+ if (ihalfq < jj) {
+ qsqz = q * gnm_sqrt(((xlegq[j] * ulen) + twa1) * 0.5);
+ } else {
+ qsqz = q * gnm_sqrt(((-(xlegq[j] * ulen)) + twa1) * 0.5);
+ }
+
+ /* call wprob to find integral of range portion */
+
+ wprb = ptukey_wprob(qsqz, rr, cc);
+ rotsum = (wprb * alegq[j]) * gnm_exp(t1);
+ otsum += rotsum;
+ }
+ /* end legendre integral for interval i */
+ /* L200: */
+ }
+
+ /* if integral for interval i < 1e-14, then stop.
+ * However, in order to avoid small area under left tail,
+ * at least 1 / ulen intervals are calculated.
+ */
+ if (i * ulen >= 1.0 && otsum <= eps2)
+ break;
+
+ /* end of interval i */
+ /* L330: */
+
+ ans += otsum;
+ }
+
+ if(otsum > eps2) { /* not converged */
+ ML_ERROR(ME_PRECISION);
+ }
+ if (ans > 1.)
+ ans = 1.;
+ return R_DT_val(ans);
+}
+
+/* Cleaning up done by tools/import-R: */
+#undef nlegq
+#undef ihalfq
+#undef nleg
+#undef ihalf
+
+/* ------------------------------------------------------------------------ */
/* --- END MAGIC R SOURCE MARKER --- */
+/* Silly order-of-arguments wrapper. */
+gnm_float
+ptukey(gnm_float x, gnm_float nmeans, gnm_float df, gnm_float nranges, gboolean lower_tail, gboolean log_p)
+{
+ return R_ptukey (x, nranges, nmeans, df, lower_tail, log_p);
+}
+
static gnm_float
logspace_signed_add (gnm_float logx, gnm_float logabsy, gboolean ypos)
{
@@ -5129,6 +5554,53 @@ swap_log_tail (gnm_float lp)
return gnm_log1p (-gnm_exp (lp)); /* Good formula for small lp. */
}
+
+gnm_float
+pnorm2 (gnm_float x1, gnm_float x2, gnm_float mu, gnm_float sigma, gboolean log_p)
+{
+ gnm_float raw;
+ const gnm_float RAWOK = 0.1;
+ const gnm_float LOG_RAWOK = GNM_const(-2.3025850929940455);
+
+ if (gnm_isnan(x1) || gnm_isnan(x2) ||
+ gnm_isnan(mu) || gnm_isnan(sigma) ||
+ sigma <= 0)
+ return gnm_nan;
+
+ if (x1 > x2)
+ return log_p ? gnm_nan : 0 - pnorm2 (x2, x1, mu, sigma, log_p);
+
+ if (x1 == x2)
+ return R_D__0;
+
+ x1 = (x1 - mu) / sigma;
+ x2 = (x2 - mu) / sigma;
+
+ /* At this point x1<x2 and we can assume N(0,1). */
+
+ if (x1 < 0) {
+ gnm_float p1 = pnorm (x1, 0.0, 1.0, TRUE, log_p);
+ gnm_float p2 = pnorm (x2, 0.0, 1.0, TRUE, log_p);
+ raw = log_p ? logspace_sub (p2, p1) : p2 - p1;
+ } else {
+ gnm_float p1C = pnorm (x1, 0.0, 1.0, FALSE, log_p);
+ gnm_float p2C = pnorm (x2, 0.0, 1.0, FALSE, log_p);
+ raw = log_p ? logspace_sub (p1C, p2C) : p1C - p2C;
+ }
+
+ if (raw > (log_p ? LOG_RAWOK : RAWOK))
+ return raw;
+
+ if (0)
+ g_printerr ("x1=%.10g x2=%.10g p1=%.10g p2=%.10g\n",
+ x1, x2,
+ pnorm (x1, 0.0, 1.0, TRUE, FALSE),
+ pnorm (x2, 0.0, 1.0, TRUE, FALSE));
+ return raw;
+}
+
+
+
/* --- BEGIN IANDJMSMITH SOURCE MARKER --- */
/* Calculation of logfbit(x)-logfbit(1+x). y2 must be < 1. */
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 038850e..ccd7d72 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -56,6 +56,7 @@ gnm_float bessel_k (gnm_float x, gnm_float alpha, gnm_float expo);
/* The normal distribution. */
gnm_float dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log);
gnm_float pnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p);
+gnm_float pnorm2 (gnm_float x1, gnm_float x2, gnm_float mu, gnm_float sigma, gboolean log_p);
gnm_float qnorm (gnm_float p, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p);
/* The log-normal distribution. */
@@ -132,6 +133,10 @@ gnm_float pcauchy (gnm_float x, gnm_float location, gnm_float scale, gboolean lo
gnm_float random_exppow_pdf (gnm_float x, gnm_float a, gnm_float b);
gnm_float random_laplace_pdf (gnm_float x, gnm_float a);
+/* Studentized range distribution */
+/* Note: argument order differs from R. */
+gnm_float ptukey(gnm_float x, gnm_float nmeans, gnm_float df, gnm_float nranges, gboolean lower_tail,
gboolean log_p);
+
/* ------------------------------------------------------------------------- */
typedef gnm_float (*GnmPFunc) (gnm_float x, const gnm_float shape[],
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
