Commit 25d43a4b authored by Andreas J. Guelzow 's avatar Andreas J. Guelzow

implement test_random_randnorm

2012-02-11  Andreas J. Guelzow <aguelzow@pyrshep.ca>

	* src/sstest.c (test_random_randnorm): new
	(test_random): enable test_random_randnorm
	(test_random_normality): new
	(test_random): use sample size of 20000
	* src/rangefunc.h (gnm_range_adtest): new
	* src/rangefunc.c (gnm_range_adtest): new, extracted from
	plugins/fn-stat/functions.c

2012-02-11  Andreas J. Guelzow <aguelzow@pyrshep.ca>

	* functions.c (gnumeric_adtest): use gnm_range_adtest
parent 608872a5
2012-02-11 Andreas J. Guelzow <aguelzow@pyrshep.ca>
* src/sstest.c (test_random_randnorm): new
(test_random): enable test_random_randnorm
(test_random_normality): new
(test_random): use sample size of 20000
* src/rangefunc.h (gnm_range_adtest): new
* src/rangefunc.c (gnm_range_adtest): new, extracted from
plugins/fn-stat/functions.c
2012-02-08 Jean Brefort <jean.brefort@normalesup.org>
* src/wbc-gtk-edit.c (wbcg_insert_object): don't destroy the object
......
2012-02-11 Andreas J. Guelzow <aguelzow@pyrshep.ca>
* functions.c (gnumeric_adtest): use gnm_range_adtest
2011-11-27 Morten Welinder <terra@gnome.org>
* Release 1.11.1
......
......@@ -4851,8 +4851,8 @@ gnumeric_adtest (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
gnm_float *xs;
int n;
GnmValue *result = NULL;
gnm_float mu = 0.;
gnm_float sigma = 1.;
gnm_float statistics = 0.;
gnm_float p = 0.;
xs = collect_floats_value (argv[0], ei->pos,
COLLECT_IGNORE_STRINGS |
......@@ -4868,43 +4868,16 @@ gnumeric_adtest (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
value_array_set (result, 0, 2,
value_new_int (n));
if ((n < 8) || gnm_range_average (xs, n, &mu)
|| gnm_range_stddev_est (xs, n, &sigma)) {
if ((n < 8) || gnm_range_adtest (xs, n, &p, &statistics)) {
value_array_set (result, 0, 0,
value_new_error_VALUE (ei->pos));
value_array_set (result, 0, 1,
value_new_error_VALUE (ei->pos));
value_new_error_VALUE (ei->pos));
} else {
int i;
gnm_float total = 0.;
gnm_float p;
gnm_float *ys;
ys = range_sort (xs, n);
for (i = 0; i < n; i++) {
gnm_float val = (pnorm (ys[i], mu, sigma, TRUE, TRUE) + pnorm (ys[n - i - 1], mu, sigma, FALSE, TRUE));
total += ((2*i+1)* val);
}
total = - n - total/n;
value_array_set (result, 0, 1,
value_new_float (total));
g_free (ys);
total *= (1 + 0.75 / n + 2.25 / (n * n));
if (total < 0.2)
p = 1. - gnm_exp (-13.436 + 101.14 * total - 223.73 * total * total);
else if (total < 0.34)
p = 1. - gnm_exp (-8.318 + 42.796 * total - 59.938 * total * total);
else if (total < 0.6)
p = gnm_exp (0.9177 - 4.279 * total - 1.38 * total * total);
else
p = gnm_exp (1.2937 - 5.709 * total + 0.0186 * total * total);
value_array_set (result, 0, 0,
value_new_float (p));
value_array_set (result, 0, 1,
value_new_float (statistics));
}
out:
......
/* vim: set sw=8: -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
/*
* rangefunc.c: Functions on ranges (data sets).
*
......@@ -14,9 +15,9 @@
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "tools/analysis-tools.h"
int
gnm_range_count (gnm_float const *xs, int n, gnm_float *res)
gnm_range_count (G_GNUC_UNUSED gnm_float const *xs, int n, gnm_float *res)
{
*res = n;
return 0;
......@@ -431,3 +432,48 @@ gnm_range_mode (gnm_float const *xs, int n, gnm_float *res)
*res = mode;
return 0;
}
int
gnm_range_adtest (gnm_float const *xs, int n, gnm_float *pvalue,
gnm_float *statistics)
{
gnm_float mu = 0.;
gnm_float sigma = 1.;
if ((n < 8) || gnm_range_average (xs, n, &mu)
|| gnm_range_stddev_est (xs, n, &sigma))
return 1;
else {
int i;
gnm_float total = 0.;
gnm_float p;
gnm_float *ys;
ys = range_sort (xs, n);
for (i = 0; i < n; i++) {
gnm_float val = (pnorm (ys[i], mu, sigma, TRUE, TRUE) +
pnorm (ys[n - i - 1],
mu, sigma, FALSE, TRUE));
total += ((2*i+1)* val);
}
total = - n - total/n;
g_free (ys);
total *= (1 + 0.75 / n + 2.25 / (n * n));
if (total < 0.2)
p = 1. - gnm_exp (-13.436 + 101.14 * total - 223.73 * total * total);
else if (total < 0.34)
p = 1. - gnm_exp (-8.318 + 42.796 * total - 59.938 * total * total);
else if (total < 0.6)
p = gnm_exp (0.9177 - 4.279 * total - 1.38 * total * total);
else
p = gnm_exp (1.2937 - 5.709 * total + 0.0186 * total * total);
if (statistics)
*statistics = total;
if (pvalue)
*pvalue = p;
return 0;
}
}
......@@ -59,6 +59,9 @@ int gnm_range_rsq_pop (gnm_float const *xs, const gnm_float *ys, int n, gnm_floa
int gnm_range_mode (gnm_float const *xs, int n, gnm_float *res);
int gnm_range_adtest (gnm_float const *xs, int n, gnm_float *p,
gnm_float *statistics);
G_END_DECLS
#endif /* _GNM_RANGEFUNC_H_ */
......@@ -56,7 +56,7 @@ mark_test_end (const char *name)
}
static void
cb_collect_names (const char *name, GnmNamedExpr *nexpr, GSList **names)
cb_collect_names (G_GNUC_UNUSED const char *name, GnmNamedExpr *nexpr, GSList **names)
{
*names = g_slist_prepend (*names, nexpr);
}
......@@ -471,6 +471,73 @@ test_random_1 (int N, const char *expr,
return res;
}
static gnm_float *
test_random_normality (int N, const char *expr,
gnm_float *mean, gnm_float *var,
gnm_float *adtest, gnm_float *cvmtest,
gnm_float *lkstest, gnm_float *sftest)
{
Workbook *wb = workbook_new ();
Sheet *sheet = workbook_sheet_add
(wb, -1, GNM_DEFAULT_COLS, GNM_DEFAULT_ROWS);
gnm_float *res = g_new (gnm_float, N);
int i;
char *s;
g_printerr ("Testing %s\n", expr);
for (i = 0; i < N; i++)
define_cell (sheet, 0, i, expr);
s = g_strdup_printf ("=average(a1:a%d)", N);
define_cell (sheet, 1, 0, s);
g_free (s);
s = g_strdup_printf ("=var(a1:a%d)", N);
define_cell (sheet, 1, 1, s);
g_free (s);
s = g_strdup_printf ("=adtest(a1:a%d)", N);
define_cell (sheet, 1, 2, s);
g_free (s);
s = g_strdup_printf ("=cvmtest(a1:a%d)", N);
define_cell (sheet, 1, 3, s);
g_free (s);
s = g_strdup_printf ("=lkstest(a1:a%d)", N);
define_cell (sheet, 1, 4, s);
g_free (s);
s = g_strdup_printf ("=sftest(a1:a%d)", N > 5000 ? 5000 : N);
define_cell (sheet, 1, 5, s);
g_free (s);
workbook_recalc (sheet->workbook);
for (i = 0; i < N; i++)
res[i] = value_get_as_float (sheet_cell_get (sheet, 0, i)->value);
*mean = value_get_as_float (sheet_cell_get (sheet, 1, 0)->value);
g_printerr ("Mean: %.10" GNM_FORMAT_g "\n", *mean);
*var = value_get_as_float (sheet_cell_get (sheet, 1, 1)->value);
g_printerr ("Var: %.10" GNM_FORMAT_g "\n", *var);
*adtest = value_get_as_float (sheet_cell_get (sheet, 1, 2)->value);
g_printerr ("ADTest: %.10" GNM_FORMAT_g "\n", *adtest);
*cvmtest = value_get_as_float (sheet_cell_get (sheet, 1, 3)->value);
g_printerr ("CVMTest: %.10" GNM_FORMAT_g "\n", *cvmtest);
*lkstest = value_get_as_float (sheet_cell_get (sheet, 1, 4)->value);
g_printerr ("LKSTest: %.10" GNM_FORMAT_g "\n", *lkstest);
*sftest = value_get_as_float (sheet_cell_get (sheet, 1, 5)->value);
g_printerr ("CVMTest: %.10" GNM_FORMAT_g "\n", *sftest);
g_object_unref (wb);
return res;
}
static void
test_random_rand (int N)
{
......@@ -569,6 +636,59 @@ test_random_randbernoulli (int N)
g_printerr ("\n");
}
static void
test_random_randnorm (int N)
{
gnm_float mean, var, adtest, cvmtest, lkstest, sftest;
gnm_float mean_target = 0, var_target = 1;
gnm_float *vals;
gboolean ok;
char *expr;
gnm_float T;
expr = g_strdup_printf ("=RANDNORM(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")",
mean_target, var_target);
vals = test_random_normality (N, expr, &mean, &var, &adtest, &cvmtest, &lkstest, &sftest);
g_free (expr);
g_free (vals);
ok = TRUE;
T = mean_target;
if (gnm_abs (mean - T) > 0.02) {
g_printerr ("Mean failure [%.10" GNM_FORMAT_g "]\n", T);
ok = FALSE;
}
T = var_target;
if (gnm_abs (var - T) > 0.02) {
g_printerr ("Var failure [%.10" GNM_FORMAT_g "]\n", T);
ok = FALSE;
}
if (adtest < 0.05) {
g_printerr ("Anderson Darling Test rejected [%.10" GNM_FORMAT_g "]\n", adtest);
ok = FALSE;
}
if (cvmtest < 0.05) {
g_printerr ("Cramér-von Mises Test rejected [%.10" GNM_FORMAT_g "]\n", cvmtest);
ok = FALSE;
}
if (lkstest < 0.05) {
g_printerr ("Lilliefors (Kolmogorov-Smirnov) Test rejected [%.10" GNM_FORMAT_g "]\n",
lkstest);
ok = FALSE;
}
if (sftest < 0.05) {
g_printerr ("Shapiro-Francia Test rejected [%.10" GNM_FORMAT_g "]\n", sftest);
ok = FALSE;
}
if (ok)
g_printerr ("OK\n");
g_printerr ("\n");
}
static void
test_random_randsnorm (int N)
{
......@@ -619,11 +739,12 @@ static void
test_random (void)
{
const char *test_name = "test_random";
const int N = 10000;
const int N = 20000;
mark_test_start (test_name);
test_random_rand (N);
test_random_randbernoulli (N);
test_random_randnorm (N);
test_random_randsnorm (N);
#if 0
test_random_randbeta (N);
......@@ -647,7 +768,6 @@ test_random (void)
test_random_randlogistic (N);
test_random_randlognorm (N);
test_random_randnegbinom (N);
test_random_randnorm (N);
test_random_randpareto (N);
test_random_randpoisson (N);
test_random_randrayleigh (N);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment