[gnumeric] fn-r: implement skew-t cdf and inverse.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] fn-r: implement skew-t cdf and inverse.
- Date: Thu, 4 Apr 2013 18:47:04 +0000 (UTC)
commit c7582a2644a5bb8b1ec6120eeb9726473d6729e3
Author: Morten Welinder <terra gnome org>
Date: Thu Apr 4 14:46:01 2013 -0400
fn-r: implement skew-t cdf and inverse.
NEWS | 2 +-
plugins/fn-r/ChangeLog | 4 ++
plugins/fn-r/extra.c | 77 +++++++++++++++++++++++++++++--
plugins/fn-r/functions.c | 110 ++++++++++++++++++++++++++++++++++++++++++-
plugins/fn-r/generate | 6 +--
plugins/fn-r/plugin.xml.in | 3 +
src/mathfunc.c | 2 +-
src/mathfunc.h | 1 +
tools/import-R | 2 +-
9 files changed, 193 insertions(+), 14 deletions(-)
---
diff --git a/NEWS b/NEWS
index 1ebceda..bcd7aa0 100644
--- a/NEWS
+++ b/NEWS
@@ -20,7 +20,7 @@ Morten:
* Fix focus problem in cell comment dialog. [#696826]
* Move sheet renaming into a dialog.
* Make it easier to see what sheet is selected. [#659317]
-
+ * Implement R.PST, R.QST, and R.QSNORM.
--------------------------------------------------------------------------
Gnumeric 1.12.1
diff --git a/plugins/fn-r/ChangeLog b/plugins/fn-r/ChangeLog
index 64651ef..f18192f 100644
--- a/plugins/fn-r/ChangeLog
+++ b/plugins/fn-r/ChangeLog
@@ -1,3 +1,7 @@
+2013-04-04 Morten Welinder <terra gnome org>
+
+ * extra.c (pst): Implement.
+
2013-03-09 Morten Welinder <terra gnome org>
* Release 1.12.1
diff --git a/plugins/fn-r/extra.c b/plugins/fn-r/extra.c
index 6366b79..9a50c3f 100644
--- a/plugins/fn-r/extra.c
+++ b/plugins/fn-r/extra.c
@@ -210,13 +210,22 @@ dst (gnm_float x, gnm_float n, gnm_float shape, gboolean give_log)
}
}
-
gnm_float
pst (gnm_float x, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean log_p)
{
+ gnm_float p;
+
+ if (n <= 0)
+ return gnm_nan;
+
if (shape == 0.)
return pt (x, n, lower_tail, log_p);
+ if (n > 100) {
+ /* Approximation */
+ return psnorm (x, shape, 0.0, 1.0, lower_tail, log_p);
+ }
+
/* Generic fallback. */
if (!lower_tail)
return log_p
@@ -225,10 +234,70 @@ pst (gnm_float x, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean lo
if (log_p)
gnm_log (pst (x, n, shape, TRUE, FALSE));
- /* FIXME: Numerical integration of dst from -inf to x. */
+ if (n != gnm_floor (n)) {
+ /* We would need numerical integration for this. */
+ return gnm_nan;
+ }
+
+ /*
+ * Use recurrence formula from "Recurrent relations for
+ * distributions of a skew-t and a linear combination of order
+ * statistics form a bivariate-t", Computational Statistics
+ * and Data Analysis volume 52, 2009 by Jamallizadeh,
+ * Khosravi, Balakrishnan.
+ *
+ * This brings us down to n==1 or n==2 for which explicit formulas
+ * are available.
+ */
+
+ p = 0;
+ while (n > 2) {
+ double a, lb, c, d, pv, v = n - 1;
+
+ d = v == 2
+ ? M_LN2gnum - gnm_log (M_PIgnum) + gnm_log (3) / 2
+ : (0.5 + M_LN2gnum / 2 - gnm_log (M_PIgnum) / 2 +
+ v / 2 * (gnm_log1p (-1 / (v - 1)) + gnm_log (v + 1)) -
+ 0.5 * (gnm_log (v - 2) + gnm_log (v + 1)) +
+ stirlerr (v / 2 - 1) -
+ stirlerr ((v - 1) / 2));
+
+ a = v + 1 + x * x;
+ lb = (d - gnm_log (a) * v / 2);
+ c = pt (gnm_sqrt (v) * shape * x / gnm_sqrt (a), v, TRUE, FALSE);
+ pv = x * gnm_exp (lb) * c;
+ p += pv;
+
+ n -= 2;
+ x *= gnm_sqrt ((v - 1) / (v + 1));
+ }
+
+ g_return_val_if_fail (n == 1 || n == 2, gnm_nan);
+ if (n == 1) {
+ gnm_float p1;
+
+ p1 = (gnm_atan (x) + gnm_acos (shape / gnm_sqrt ((1 + shape * shape) * (1 + x * x)))) /
M_PIgnum;
+ p += p1;
+ } else if (n == 2) {
+ gnm_float p2, f;
+
+ f = x / gnm_sqrt (2 + x * x);
+
+ p2 = (0.5 - gnm_atan (shape) / M_PIgnum) +
+ f * (0.5 + gnm_atan (shape * f) / M_PIgnum);
+
+ p += p2;
+ } else {
+ return gnm_nan;
+ }
+
+ /*
+ * Negatives can occur due to rounding errors and hopefully for no
+ * other reason.
+ */
+ p = CLAMP (p, 0.0, 1.0);
- /* Give up. */
- return gnm_nan;
+ return p;
}
diff --git a/plugins/fn-r/functions.c b/plugins/fn-r/functions.c
index 45ef4ce..066bd4b 100644
--- a/plugins/fn-r/functions.c
+++ b/plugins/fn-r/functions.c
@@ -1187,7 +1187,7 @@ static GnmFuncHelp const help_r_dsnorm[] = {
{ GNM_FUNC_HELP_ARG, F_("scale:the scale parameter of the distribution") },
{ GNM_FUNC_HELP_ARG, F_("give_log:if true, log of the result will be returned instead") },
{ GNM_FUNC_HELP_DESCRIPTION, F_("This function returns the probability density function of the
skew-normal distribution.") },
- { GNM_FUNC_HELP_SEEALSO, "R.PSNORM" },
+ { GNM_FUNC_HELP_SEEALSO, "R.PSNORM,R.QSNORM" },
{ GNM_FUNC_HELP_END }
};
@@ -1215,7 +1215,7 @@ static GnmFuncHelp const help_r_psnorm[] = {
{ 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
skew-normal distribution.") },
- { GNM_FUNC_HELP_SEEALSO, "R.DSNORM" },
+ { GNM_FUNC_HELP_SEEALSO, "R.DSNORM,R.QSNORM" },
{ GNM_FUNC_HELP_END }
};
@@ -1235,6 +1235,35 @@ gnumeric_r_psnorm (GnmFuncEvalInfo *ei, GnmValue const * const *args)
/* ------------------------------------------------------------------------- */
+static GnmFuncHelp const help_r_qsnorm[] = {
+ { GNM_FUNC_HELP_NAME, F_("R.QSNORM:probability quantile function of the skew-normal distribution") },
+ { GNM_FUNC_HELP_ARG, F_("p:probability") },
+ { GNM_FUNC_HELP_ARG, F_("shape:the shape parameter of the distribution") },
+ { GNM_FUNC_HELP_ARG, F_("location:the location parameter of the distribution") },
+ { GNM_FUNC_HELP_ARG, F_("scale:the scale parameter of the distribution") },
+ { 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 probability quantile function, i.e., the
inverse of the cumulative distribution function, of the skew-normal distribution.") },
+ { GNM_FUNC_HELP_SEEALSO, "R.DSNORM,R.PSNORM" },
+ { GNM_FUNC_HELP_END }
+};
+
+static GnmValue *
+gnumeric_r_qsnorm (GnmFuncEvalInfo *ei, GnmValue const * const *args)
+{
+ gnm_float p = value_get_as_float (args[0]);
+ gnm_float shape = value_get_as_float (args[1]);
+ gnm_float location = value_get_as_float (args[2]);
+ gnm_float scale = value_get_as_float (args[3]);
+ 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 (qsnorm (p, shape, location, scale, lower_tail, log_p));
+}
+
+/* ------------------------------------------------------------------------- */
+
+
static GnmFuncHelp const help_r_dst[] = {
{ GNM_FUNC_HELP_NAME, F_("R.DST:probability density function of the skew-t distribution") },
{ GNM_FUNC_HELP_ARG, F_("x:observation") },
@@ -1242,7 +1271,7 @@ static GnmFuncHelp const help_r_dst[] = {
{ GNM_FUNC_HELP_ARG, F_("shape:the shape parameter of the distribution") },
{ GNM_FUNC_HELP_ARG, F_("give_log:if true, log of the result will be returned instead") },
{ GNM_FUNC_HELP_DESCRIPTION, F_("This function returns the probability density function of the skew-t
distribution.") },
- { GNM_FUNC_HELP_SEEALSO, "" },
+ { GNM_FUNC_HELP_SEEALSO, "R.PST,R.QST" },
{ GNM_FUNC_HELP_END }
};
@@ -1259,6 +1288,60 @@ gnumeric_r_dst (GnmFuncEvalInfo *ei, GnmValue const * const *args)
/* ------------------------------------------------------------------------- */
+
+static GnmFuncHelp const help_r_pst[] = {
+ { GNM_FUNC_HELP_NAME, F_("R.PST:cumulative distribution function of the skew-t distribution") },
+ { GNM_FUNC_HELP_ARG, F_("x:observation") },
+ { GNM_FUNC_HELP_ARG, F_("n:the number of degrees of freedom of the distribution") },
+ { GNM_FUNC_HELP_ARG, F_("shape:the shape parameter of the distribution") },
+ { 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
skew-t distribution.") },
+ { GNM_FUNC_HELP_SEEALSO, "R.DST,R.QST" },
+ { GNM_FUNC_HELP_END }
+};
+
+static GnmValue *
+gnumeric_r_pst (GnmFuncEvalInfo *ei, GnmValue const * const *args)
+{
+ gnm_float x = value_get_as_float (args[0]);
+ gnm_float n = value_get_as_float (args[1]);
+ gnm_float shape = value_get_as_float (args[2]);
+ gboolean lower_tail = args[3] ? value_get_as_checked_bool (args[3]) : TRUE;
+ gboolean log_p = args[4] ? value_get_as_checked_bool (args[4]) : FALSE;
+
+ return value_new_float (pst (x, n, shape, lower_tail, log_p));
+}
+
+/* ------------------------------------------------------------------------- */
+
+
+static GnmFuncHelp const help_r_qst[] = {
+ { GNM_FUNC_HELP_NAME, F_("R.QST:probability quantile function of the skew-t distribution") },
+ { GNM_FUNC_HELP_ARG, F_("p:probability") },
+ { GNM_FUNC_HELP_ARG, F_("n:the number of degrees of freedom of the distribution") },
+ { GNM_FUNC_HELP_ARG, F_("shape:the shape parameter of the distribution") },
+ { 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 probability quantile function, i.e., the
inverse of the cumulative distribution function, of the skew-t distribution.") },
+ { GNM_FUNC_HELP_SEEALSO, "R.DST,R.PST" },
+ { GNM_FUNC_HELP_END }
+};
+
+static GnmValue *
+gnumeric_r_qst (GnmFuncEvalInfo *ei, GnmValue const * const *args)
+{
+ gnm_float p = value_get_as_float (args[0]);
+ gnm_float n = value_get_as_float (args[1]);
+ gnm_float shape = value_get_as_float (args[2]);
+ gboolean lower_tail = args[3] ? value_get_as_checked_bool (args[3]) : TRUE;
+ gboolean log_p = args[4] ? value_get_as_checked_bool (args[4]) : FALSE;
+
+ return value_new_float (qst (p, n, shape, lower_tail, log_p));
+}
+
+/* ------------------------------------------------------------------------- */
+
G_MODULE_EXPORT void
go_plugin_init (GOPlugin *plugin, GOCmdContext *cc)
{
@@ -1602,11 +1685,32 @@ GnmFuncDescriptor const stat_functions[] = {
GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
},
{
+ "r.qsnorm",
+ "ffff|bb",
+ help_r_qsnorm,
+ gnumeric_r_qsnorm, NULL, NULL, NULL,
+ GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
+ },
+ {
"r.dst",
"fff|b",
help_r_dst,
gnumeric_r_dst, NULL, NULL, NULL,
GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
},
+ {
+ "r.pst",
+ "fff|bb",
+ help_r_pst,
+ gnumeric_r_pst, NULL, NULL, NULL,
+ GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
+ },
+ {
+ "r.qst",
+ "fff|bb",
+ help_r_qst,
+ gnumeric_r_qst, NULL, NULL, NULL,
+ GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
+ },
{ NULL }
};
diff --git a/plugins/fn-r/generate b/plugins/fn-r/generate
index 53e18f7..c136a8d 100644
--- a/plugins/fn-r/generate
+++ b/plugins/fn-r/generate
@@ -121,15 +121,13 @@ my %argoverride = ();
({ 'psuc' => "the probability of success in each trial",
@common })];
- $funcs{'dsnorm'} = $funcs{'psnorm'} =
-#$funcs{'qsnorm'} =
+ $funcs{'dsnorm'} = $funcs{'psnorm'} = $funcs{'qsnorm'} =
[\&distribution,
'skew-normal',
({ 'location' => "the location parameter $of",
@common })];
- $funcs{'dst'} =
-#$funcs{'pst'} = $funcs{'qst'} =
+ $funcs{'dst'} = $funcs{'pst'} = $funcs{'qst'} =
[\&distribution,
'skew-t',
({ 'n' => "the number of degrees of freedom $of",
diff --git a/plugins/fn-r/plugin.xml.in b/plugins/fn-r/plugin.xml.in
index b5ee8c3..901f254 100644
--- a/plugins/fn-r/plugin.xml.in
+++ b/plugins/fn-r/plugin.xml.in
@@ -42,6 +42,7 @@
<function name="r.pnorm"/>
<function name="r.ppois"/>
<function name="r.psnorm"/>
+ <function name="r.pst"/>
<function name="r.pt"/>
<function name="r.pweibull"/>
<function name="r.qbeta"/>
@@ -57,6 +58,8 @@
<function name="r.qnbinom"/>
<function name="r.qnorm"/>
<function name="r.qpois"/>
+ <function name="r.qsnorm"/>
+ <function name="r.qst"/>
<function name="r.qt"/>
<function name="r.qweibull"/>
</functions>
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 18680dc..bd8e827 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -905,7 +905,7 @@ gnm_float ppois(gnm_float x, gnm_float lambda, gboolean lower_tail, gboolean log
* see also lgammacor() in ./lgammacor.c which computes almost the same!
*/
-static gnm_float stirlerr(gnm_float n)
+gnm_float stirlerr(gnm_float n)
{
#define S0 GNM_const(0.083333333333333333333) /* 1/12 */
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 73304f8..67acb33 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -34,6 +34,7 @@ gnm_float gnm_trunc (gnm_float x);
gnm_float logfbit (gnm_float x);
gnm_float logspace_add (gnm_float logx, gnm_float logy);
gnm_float logspace_sub (gnm_float logx, gnm_float logy);
+gnm_float stirlerr(gnm_float n);
gnm_float gnm_cot (gnm_float x);
gnm_float gnm_acot (gnm_float x);
diff --git a/tools/import-R b/tools/import-R
index 6b3f1a0..7b8b8ee 100644
--- a/tools/import-R
+++ b/tools/import-R
@@ -152,7 +152,7 @@ sub import_file {
s~([-+]?\d+\.\d*)(\s*/\s*)([-+e0-9.]+)~GNM_const\($1\)$2$3~;
# These are made static.
-
s/^gnm_float\s+(pbeta_both|stirlerr|bd0|dpois_raw|chebyshev_eval|lgammacor|lbeta|pbeta_raw|dbinom_raw)/static
gnm_float $1/;
+
s/^gnm_float\s+(pbeta_both|bd0|dpois_raw|chebyshev_eval|lgammacor|lbeta|pbeta_raw|dbinom_raw)/static
gnm_float $1/;
s/^int\s+(chebyshev_init)/static int $1/;
# Fix a bug...
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
