GNOME.org
 

[gnumeric] Add R.PSNORM and R.DSNORM.

commit a634a18f72d5f807e1121514a41306bf024733e9
Author: Andreas J. Guelzow <aguelzow pyrshep ca>
Date:   Mon Oct 12 17:51:27 2009 -0600

    Add R.PSNORM and R.DSNORM.
    
    2009-10-12  Andreas J. Guelzow <aguelzow pyrshep ca>
    
    	* src/mathfunc.c (random_skew_normal): simplify
    	(random_skew_tdist): simplify
    
    2009-10-12  Andreas J. Guelzow <aguelzow pyrshep ca>
    
    	* generate: handle *snorm and *st
    	* extra.h (dsnorm): new
    	(psnorm): new
    	(qsnorm): new
    	(dst): new
    	(pst): new
    	(qst): new
    	* extra.c (dsnorm): new
    	(gnm_owent): new
    	(psnorm): new
    	(qsnorm): new stub
    	(dst): new stub
    	(pst): new stub
    	(qst): new stub
    	(plugin.xml.in): updated due to changes in generate
    	(functions.c): ditto

 ChangeLog                  |    7 ++-
 NEWS                       |    1 +
 plugins/fn-r/ChangeLog     |   19 +++++
 plugins/fn-r/extra.c       |  177 ++++++++++++++++++++++++++++++++++++++++++++
 plugins/fn-r/extra.h       |   11 +++
 plugins/fn-r/functions.c   |   70 +++++++++++++++++
 plugins/fn-r/generate      |   14 ++++
 plugins/fn-r/plugin.xml.in |    2 +
 src/mathfunc.c             |    7 +-
 9 files changed, 302 insertions(+), 6 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index fbba355..7771fc9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,10 +1,15 @@
+2009-10-11  Andreas J. Guelzow <aguelzow pyrshep ca>
+
+	* src/mathfunc.c (random_skew_normal): simplify
+	(random_skew_tdist): simplify
+
 2009-10-11  Morten Welinder  <terra gnome org>
 
 	* src/sheet.c (gnm_sheet_resize_main): Reduce ->cols.max_used and
 	->rows.max_used as appropriate.
 	(gnm_sheet_resize): Check for merges.  Add new perr argument.  All
 	callers changed.
-
+	
 2009-10-11  Andreas J. Guelzow <aguelzow pyrshep ca>
 
 	* src/mathfunc.c (random_skew_normal): new
diff --git a/NEWS b/NEWS
index aa2b500..57f6ad9 100644
--- a/NEWS
+++ b/NEWS
@@ -2,6 +2,7 @@ Gnumeric 1.9.15
 
 Andreas:
 	* Add RANDSNORM and RANDSTDIST. [#144717]
+	* Add R.PSNORM and R.DSNORM.
 
 Jody:
 	* First steps towards a turnkey win32 build.
diff --git a/plugins/fn-r/ChangeLog b/plugins/fn-r/ChangeLog
index 5c2533e..9f75215 100644
--- a/plugins/fn-r/ChangeLog
+++ b/plugins/fn-r/ChangeLog
@@ -1,3 +1,22 @@
+2009-10-12  Andreas J. Guelzow <aguelzow pyrshep ca>
+
+	* generate: handle *snorm and *st
+	* extra.h (dsnorm): new
+	(psnorm): new	
+	(qsnorm): new
+	(dst): new
+	(pst): new
+	(qst): new
+	* extra.c (dsnorm): new
+	(gnm_owent): new
+	(psnorm): new	
+	(qsnorm): new stub
+	(dst): new stub
+	(pst): new stub
+	(qst): new stub
+	(plugin.xml.in): updated due to changes in generate
+	(functions.c): ditto
+
 2009-10-11  Morten Welinder <terra gnome org>
 
 	* Release 1.9.14
diff --git a/plugins/fn-r/extra.c b/plugins/fn-r/extra.c
index 955af00..0d42bf6 100644
--- a/plugins/fn-r/extra.c
+++ b/plugins/fn-r/extra.c
@@ -38,4 +38,181 @@ qcauchy (gnm_float p, gnm_float location, gnm_float scale,
 	return location + scale / gnm_tan(M_PIgnum * p);
 }
 
+
 /* ------------------------------------------------------------------------- */
+
+/* This implementation of Owen's T function is based on code licensed under GPL v.2: */
+
+/*  GNU General Public License Agreement */
+/*  Copyright (C) 2004-2007 CodeCogs, Zyba Ltd, Broadwood, Holford, TA5 1DU, England. */
+
+/*  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 CodeCogs.  */
+/*  You must retain a copy of this licence in all copies.  */
+
+/*  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. */
+/*  --------------------------------------------------------------------------------- */
+/* ! Evaluates the Owen&#039;s T function. */
+
+#define LIM1 1E-35
+#define LIM2 15.0
+#define LIM3 15.0
+#define LIM4 1E-5
+ 
+#define TWOPI_INVERSE 1/(2*M_PIgnum)
+
+static gnm_float
+gnm_owent (gnm_float h, gnm_float a)
+{
+	gnm_float weight[10] = { 0.0666713443086881375935688098933,
+				 0.149451349150580593145776339658,
+				 0.219086362515982043995534934228,
+				 0.269266719309996355091226921569,
+				 0.295524224714752870173892994651,
+				 0.295524224714752870173892994651,
+				 0.269266719309996355091226921569,
+				 0.219086362515982043995534934228,
+				 0.149451349150580593145776339658,
+				 0.0666713443086881375935688098933 };
+	gnm_float xtab[10] = {0.026093471482828279922035987916,
+			      0.134936633311015489267903311577,
+			      0.320590431700975593765672634885,
+			      0.566604605870752809200734056834,
+			      0.85112566101836878911517399887,
+			      1.148874338981631210884826001130,
+			      1.433395394129247190799265943166,
+			      1.679409568299024406234327365115,
+			      1.865063366688984510732096688423,
+			      1.973906528517171720077964012084};
+	gnm_float hs, h2, as, rt;
+	int i;
+
+	if (fabs(h) < LIM1) return atan(a) * TWOPI_INVERSE;
+	if (fabs(h) > LIM2 || fabs(a) < LIM1) return 0.0;
+	
+	hs = -0.5 * h * h;
+	h2 = a;
+	as = a * a;
+	
+	if (log(1.0 + as) - hs * as >= LIM3)
+	{
+		gnm_float h1 = 0.5 * a;
+		as *= 0.25;
+		for (;;)
+		{
+			gnm_float rt = as + 1.0;
+			h2 = h1 + (hs * as + LIM3 - log(rt)) 
+				/ (2.0 * h1 * (1.0 / rt - hs));
+			as = h2 * h2;
+			if (fabs(h2 - h1) < LIM4) break;
+			h1 = h2;
+		}
+	}
+		
+	rt = 0.0;
+	for (i = 0; i < 10; i++)
+	{
+		gnm_float x = 0.5 * h2 * xtab[i], tmp = 1.0 + x * x;
+		rt += weight[i] * gnm_exp (hs * tmp) / tmp;
+	}
+	return 0.5 * rt * h2 * TWOPI_INVERSE;
+}
+
+
+#undef LIM1
+#undef LIM2
+#undef LIM3
+#undef LIM4
+#undef TWOPI_INVERSE
+
+/* ------------------------------------------------------------------------- */
+
+/* The skew-normal distribution.  */
+
+gnm_float 
+dsnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean give_log)
+{
+	if (shape == 0.) 
+		return dnorm (x, location, scale, give_log);
+	else if (give_log)
+		return gnm_log (2.) + dnorm (x, location, scale, TRUE) + pnorm (shape * x, shape * location, scale, TRUE, TRUE);
+	else
+		return 2 * dnorm (x, location, scale, FALSE) * pnorm (shape * x, location/shape, scale, TRUE, FALSE);
+}
+
+gnm_float 
+psnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p)
+{
+	gnm_float result;
+
+	if (shape == 0.) 
+		return pnorm (x, location, scale, lower_tail, log_p);
+	
+	result = pnorm (x, location, scale, TRUE, FALSE) - 2 * gnm_owent ((x - location)/scale, shape);
+
+	if (!lower_tail)
+		result = 1. - result;
+
+	if (log_p)
+		return gnm_log (result);
+	else
+		return result;
+}
+
+
+gnm_float 
+qsnorm (gnm_float p, gnm_float shape, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p)
+{
+	if (shape == 0.) 
+		return qnorm (p, location, scale, lower_tail, log_p);
+	else if (log_p)
+		return 0.;
+	else
+		return 0;
+}
+
+/* ------------------------------------------------------------------------- */
+
+/* The skew-t distribution.  */
+
+gnm_float 
+dst (gnm_float x, gnm_float n, gnm_float shape, gboolean give_log)
+{
+	if (shape == 0.) 
+		return dt (x, n, give_log);
+	else if (give_log)
+		return 0.;
+	else
+		return 0.;
+}
+
+
+gnm_float 
+pst (gnm_float x, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean log_p)
+{
+	if (shape == 0.) 
+		return pt (x, n, lower_tail, log_p);
+	else if (log_p)
+		return 0.;
+	else
+		return 0.;
+}
+
+
+gnm_float 
+qst (gnm_float p, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean log_p)
+{
+	if (shape == 0.) 
+		return qt (p, n, lower_tail, log_p);
+	else if (log_p)
+		return 0.;
+	else
+		return 0.;
+}
+
+
+
+/* ------------------------------------------------------------------------- */
+
diff --git a/plugins/fn-r/extra.h b/plugins/fn-r/extra.h
index 9d84a71..4a5ded3 100644
--- a/plugins/fn-r/extra.h
+++ b/plugins/fn-r/extra.h
@@ -6,4 +6,15 @@
 
 gnm_float qcauchy (gnm_float p, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p);
 
+/* The skew-normal distribution.  */
+gnm_float dsnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean give_log);
+gnm_float psnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p);
+gnm_float qsnorm (gnm_float p, gnm_float shape, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p);
+
+/* The skew-t distribution.  */
+gnm_float 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 qst (gnm_float p, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean log_p);
+
+
 #endif
diff --git a/plugins/fn-r/functions.c b/plugins/fn-r/functions.c
index 90e5e00..5bf1630 100644
--- a/plugins/fn-r/functions.c
+++ b/plugins/fn-r/functions.c
@@ -1178,6 +1178,62 @@ gnumeric_r_qcauchy (GnmFuncEvalInfo *ei, GnmValue const * const *args)
 
 /* ------------------------------------------------------------------------- */
 
+
+static GnmFuncHelp const help_r_dsnorm[] = {
+	{ GNM_FUNC_HELP_NAME, F_("R.DSNORM:probability density function of the skew-normal distribution") },
+	{ GNM_FUNC_HELP_ARG, F_("x:observation") },
+	{ 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_("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_END }
+};
+
+static GnmValue *
+gnumeric_r_dsnorm (GnmFuncEvalInfo *ei, GnmValue const * const *args)
+{
+	gnm_float x = 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 give_log = args[4] ? value_get_as_checked_bool (args[4]) : FALSE;
+
+	return value_new_float (dsnorm (x, shape, location, scale, give_log));
+}
+
+/* ------------------------------------------------------------------------- */
+
+
+static GnmFuncHelp const help_r_psnorm[] = {
+	{ GNM_FUNC_HELP_NAME, F_("R.PSNORM:cumulative distribution function of the skew-normal distribution") },
+	{ GNM_FUNC_HELP_ARG, F_("x:observation") },
+	{ 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 cumulative distribution function of the skew-normal distribution.") },
+	{ GNM_FUNC_HELP_SEEALSO, "R.DSNORM" },
+	{ GNM_FUNC_HELP_END }
+};
+
+static GnmValue *
+gnumeric_r_psnorm (GnmFuncEvalInfo *ei, GnmValue const * const *args)
+{
+	gnm_float x = 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 (psnorm (x, shape, location, scale, lower_tail, log_p));
+}
+
+/* ------------------------------------------------------------------------- */
+
 G_MODULE_EXPORT void
 go_plugin_init (GOPlugin *plugin, GOCmdContext *cc)
 {
@@ -1506,5 +1562,19 @@ GnmFuncDescriptor const stat_functions[] = {
 		gnumeric_r_qcauchy, NULL, NULL, NULL, NULL,
 		GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
 	},
+	{
+		"r.dsnorm",
+		"ffff|b",
+		help_r_dsnorm,
+		gnumeric_r_dsnorm, NULL, NULL, NULL, NULL,
+		GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
+	},
+	{
+		"r.psnorm",
+		"ffff|bb",
+		help_r_psnorm,
+		gnumeric_r_psnorm, NULL, 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 86bb584..a8cc3bf 100644
--- a/plugins/fn-r/generate
+++ b/plugins/fn-r/generate
@@ -138,8 +138,22 @@ my %argoverride = ();
 	    @common })];
     $argoverride{"dgeom:p"} = $argoverride{"pgeom:p"} = $argoverride{"qgeom:prob"} =
 	"psuc";
+
+    $funcs{'dsnorm'} =  $funcs{'psnorm'} =  
+#$funcs{'qsnorm'} =
+	[\&distribution,
+	 'skew-normal',
+	 ({ 'location' => "the location parameter $of",
+	    @common })];
+
+#    $funcs{'dst'} = $funcs{'pst'} = $funcs{'qst'} =
+#	[\&distribution,
+#	 'skew-t',
+#	 ({ 'n' => "the number of degrees of freedom $of",
+#	    @common })];
 }
 
+
 my %odf_note =
     ('qchisq' => 'A two argument invocation R.QCHISQ(@{p},@{df}) is exported to OpenFormula as CHISQINV(@{p},@{df}).',
      'pchisq' => 'A two argument invocation R.PCHISQ(@{x},@{df}) is exported to OpenFormula as CHISQDIST(@{x},@{df}).',
diff --git a/plugins/fn-r/plugin.xml.in b/plugins/fn-r/plugin.xml.in
index 8345670..79662be 100644
--- a/plugins/fn-r/plugin.xml.in
+++ b/plugins/fn-r/plugin.xml.in
@@ -24,6 +24,7 @@
 				<function name="r.dnbinom"/>
 				<function name="r.dnorm"/>
 				<function name="r.dpois"/>
+				<function name="r.dsnorm"/>
 				<function name="r.dt"/>
 				<function name="r.dweibull"/>
 				<function name="r.pbeta"/>
@@ -39,6 +40,7 @@
 				<function name="r.pnbinom"/>
 				<function name="r.pnorm"/>
 				<function name="r.ppois"/>
+				<function name="r.psnorm"/>
 				<function name="r.pt"/>
 				<function name="r.pweibull"/>
 				<function name="r.qbeta"/>
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 42bf902..6881fa4 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -7629,13 +7629,10 @@ gnm_float
 random_skew_normal (gnm_float a)
 {
 	gnm_float result;
-	gnm_float asq = a * a;
-	gnm_float delta = gnm_sqrt (asq / (1 + asq));
+	gnm_float delta = a / gnm_sqrt(1 + a * a);
 	gnm_float u = random_normal ();
 	gnm_float v = random_normal ();
 	
-	if (a < 0.)
-		delta *= -1.;
 	result = delta * u + gnm_sqrt (1-delta*delta) * v;
 	
 	return ((u < 0.) ? -result : result);
@@ -7658,7 +7655,7 @@ random_skew_tdist (gnm_float nu, gnm_float a)
 	gnm_float chi = random_chisq (nu);
 	gnm_float z = random_skew_normal (a);;
 
-	return (z / gnm_sqrt(chi/nu));
+	return (z * gnm_sqrt(nu/chi));
 }
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]