From: Jason Stover Date: Wed, 7 Jun 2006 01:12:37 +0000 (+0000) Subject: Renamed time-series directory ts X-Git-Tag: v0.6.0~829 X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?p=pspp-builds.git;a=commitdiff_plain;h=0a490aaca0b5a6350d3428101b8c496f396ac26b Renamed time-series directory ts --- diff --git a/src/math/time-series/innovations.c b/src/math/time-series/innovations.c deleted file mode 100644 index b9a9dde1..00000000 --- a/src/math/time-series/innovations.c +++ /dev/null @@ -1,214 +0,0 @@ -/* - src/math/time-series/arma/innovations.c - - Copyright (C) 2006 Free Software Foundation, Inc. Written by Jason H. Stover. - - 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., 51 - Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA. - */ -/* - Find preliminary ARMA coefficients via the innovations algorithm. - Also compute the sample mean and covariance matrix for each series. - - Reference: - - P. J. Brockwell and R. A. Davis. Time Series: Theory and - Methods. Second edition. Springer. New York. 1991. ISBN - 0-387-97429-6. Sections 5.2, 8.3 and 8.4. - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#define _(msgid) gettext (msgid) - -static void -get_mean_variance (size_t n_vars, const struct casefile *cf, - struct innovations_estimate **est) - -{ - struct casereader *r; - struct ccase *c; - struct ccase *c2; - size_t n; - double *x; - double d; - double tmp; - double variance; - - for (n = 0; n < n_vars; n++) - { - est[n]->n_obs = 2.0; - est[n]->mean = 0.0; - est[n]->variance = 0.0; - } - for (r = casefile_get_reader (cf); casereader_read (r, &c); - case_destroy (&c)) - { - for (n = 0; n < n_vars; n++) - { - if (!mv_is_value_missing (&v->miss, val)) - { - tmp = case_data (&c, est[n]->variable->fv); - d = (tmp - est[n]->mean) / est[n]->n_obs; - est[n]->mean += d; - est[n]->variance += est[n]->n_obs * est[n]->n_obs * d * d; - est[n]->n_obs += 1.0; - } - } - } - for (n = 0; n < n_vars; n++) - { - /* Maximum likelihood estimate of the variance. */ - est[n]->variance /= est[n]->n_obs; - } -} - -/* - Read the first MAX_LAG cases. - */ -static bool -innovations_init_cases (struct casereader *r, struct ccase **c, size_t max_lag) -{ - bool value = true; - size_t lag = 0; - - while (value) - { - lag++; - value = casereader_read (r, c + lag); - } - return value; -} - -/* - Read one case and update C, which contains the last MAX_LAG cases. - */ -static bool -innovations_update_cases (struct casereader *r, struct ccase **c, size_t max_lag) -{ - size_t lag; - bool value = false; - - for (lag = 0; lag < max_lag - 1; lag++) - { - c[lag] = c[lag+1]; - } - value = casereader_read (r, c + lag); - return value; -} -static void -get_covariance (size_t n_vars, const struct casefile *cf, - struct innovations **est, size_t max_lag) -{ - struct casereader *r; - struct ccase **c; - struct ccase *cur_case; - size_t lag; - size_t n_vars; - bool read_case = false; - double d; - double tmp; - - c = xnmalloc (max_lag, sizeof (*c)); - - for (lag = 0; lag < max_lag; lag++) - { - c[lag] = xmalloc (sizeof *c[i]); - } - - r = casefile_get_reader (cf); - read_case = innovations_init_cases (r, c, max_lag); - - while (read_case) - { - for (n = 0; n < n_vars; n++) - { - cur_case = case_data (c[0], est[n]->variable->fv); - if (!mv_is_value_missing (&est[n]->variable->miss, cur_case)) - { - cur_case -= est[n]->mean; - for (lag = 1; lag <= max_lag; lag++) - { - tmp = case_data (c[lag], est[n]->variable->fv); - if (!mv_is_value_missing (&est[n]->variable->miss, tmp)) - { - d = (tmp - est[n]->mean); - *(est[n]->cov + lag) += d * cur_case; - } - } - } - } - read_case = innovations_update_cases (r, c, max_lag); - } - for (lag = 0; lag <= max_lag; lag++) - { - for (n = 0; n < n_vars; n++) - { - *(est[n]->cov + lag) /= (est[n]->n_obs - lag); - } - } - for (lag = 0; lag < max_lag; lag++) - { - free (c[lag]); - } - free (c); -} - -struct innovations_estimate ** pspp_innovations (const struct variable **vars, size_t *n_vars, - size_t lag, const struct casefile *cf) -{ - struct innovations_estimate **est; - struct casereader *r; - struct ccase *c; - size_t i; - - est = xnmalloc (*n_vars, sizeof *est); - for (i = 0; i < *n_vars; i++) - { - if (vars[i]->type == NUMERIC) - { - est[i] = xmalloc (sizeof **est); - est[i]->variable = vars[i]; - est[i]->mean = 0.0; - est[i]->variance = 0.0; - est[i]->cov = gsl_matrix_calloc (lag, lag); - est[i]->coeff = xnmalloc (lag, sizeof (*est[i]->coeff)); - for (j = 0; j < lag; j++) - { - est[i]->coeff + j = xmalloc (sizeof (*(est[i]->coeff + j))); - } - } - else - { - *n_vars--; - msg (MW, _("Cannot compute autocovariance for a non-numeric variable %s"), - var_to_string (vars[i])); - } - } - - /* - First data pass to get the mean and variance. - */ - get_mean_variance (*n_vars, cf, est); - get_covariance (*n_vars, cf, est, lag); -} diff --git a/src/math/ts/innovations.c b/src/math/ts/innovations.c new file mode 100644 index 00000000..dba44849 --- /dev/null +++ b/src/math/ts/innovations.c @@ -0,0 +1,215 @@ +/* + src/math/time-series/arma/innovations.c + + Copyright (C) 2006 Free Software Foundation, Inc. Written by Jason H. Stover. + + 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., 51 + Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA. + */ +/* + Find preliminary ARMA coefficients via the innovations algorithm. + Also compute the sample mean and covariance matrix for each series. + + Reference: + + P. J. Brockwell and R. A. Davis. Time Series: Theory and + Methods. Second edition. Springer. New York. 1991. ISBN + 0-387-97429-6. Sections 5.2, 8.3 and 8.4. + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#define _(msgid) gettext (msgid) + +static void +get_mean_variance (size_t n_vars, const struct casefile *cf, + struct innovations_estimate **est) + +{ + struct casereader *r; + struct ccase *c; + struct ccase *c2; + size_t n; + double *x; + double d; + double tmp; + double variance; + + for (n = 0; n < n_vars; n++) + { + est[n]->n_obs = 2.0; + est[n]->mean = 0.0; + est[n]->variance = 0.0; + } + for (r = casefile_get_reader (cf); casereader_read (r, &c); + case_destroy (&c)) + { + for (n = 0; n < n_vars; n++) + { + if (!mv_is_value_missing (&v->miss, val)) + { + tmp = case_data (&c, est[n]->variable->fv); + d = (tmp - est[n]->mean) / est[n]->n_obs; + est[n]->mean += d; + est[n]->variance += est[n]->n_obs * est[n]->n_obs * d * d; + est[n]->n_obs += 1.0; + } + } + } + for (n = 0; n < n_vars; n++) + { + /* Maximum likelihood estimate of the variance. */ + est[n]->variance /= est[n]->n_obs; + } +} + +/* + Read the first MAX_LAG cases. + */ +static bool +innovations_init_cases (struct casereader *r, struct ccase **c, size_t max_lag) +{ + bool value = true; + size_t lag = 0; + + while (value) + { + lag++; + value = casereader_read (r, c + lag); + } + return value; +} + +/* + Read one case and update C, which contains the last MAX_LAG cases. + */ +static bool +innovations_update_cases (struct casereader *r, struct ccase **c, size_t max_lag) +{ + size_t lag; + bool value = false; + + for (lag = 0; lag < max_lag - 1; lag++) + { + c[lag] = c[lag+1]; + } + value = casereader_read (r, c + lag); + return value; +} +static void +get_covariance (size_t n_vars, const struct casefile *cf, + struct innovations **est, size_t max_lag) +{ + struct casereader *r; + struct ccase **c; + struct ccase *cur_case; + size_t lag; + size_t n_vars; + bool read_case = false; + double d; + double tmp; + + c = xnmalloc (max_lag, sizeof (*c)); + + for (lag = 0; lag < max_lag; lag++) + { + c[lag] = xmalloc (sizeof *c[i]); + } + + r = casefile_get_reader (cf); + read_case = innovations_init_cases (r, c, max_lag); + + while (read_case) + { + for (n = 0; n < n_vars; n++) + { + cur_case = case_data (c[0], est[n]->variable->fv); + if (!mv_is_value_missing (&est[n]->variable->miss, cur_case)) + { + cur_case -= est[n]->mean; + for (lag = 1; lag <= max_lag; lag++) + { + tmp = case_data (c[lag], est[n]->variable->fv); + if (!mv_is_value_missing (&est[n]->variable->miss, tmp)) + { + d = (tmp - est[n]->mean); + *(est[n]->cov + lag) += d * cur_case; + } + } + } + } + read_case = innovations_update_cases (r, c, max_lag); + } + for (lag = 0; lag <= max_lag; lag++) + { + for (n = 0; n < n_vars; n++) + { + *(est[n]->cov + lag) /= (est[n]->n_obs - lag); + } + } + for (lag = 0; lag < max_lag; lag++) + { + free (c[lag]); + } + free (c); +} + +struct innovations_estimate ** pspp_innovations (const struct variable **vars, size_t *n_vars, + size_t lag, const struct casefile *cf) +{ + struct innovations_estimate **est; + struct casereader *r; + struct ccase *c; + size_t i; + + est = xnmalloc (*n_vars, sizeof *est); + for (i = 0; i < *n_vars; i++) + { + if (vars[i]->type == NUMERIC) + { + est[i] = xmalloc (sizeof **est); + est[i]->variable = vars[i]; + est[i]->mean = 0.0; + est[i]->variance = 0.0; + est[i]->cov = xnmalloc (lag, sizeof (est[i]->cov)); + est[i]->coeff = xnmalloc (lag, sizeof (*est[i]->coeff)); + for (j = 0; j < lag; j++) + { + est[i]->coeff + j = xmalloc (sizeof (*(est[i]->coeff + j))); + } + } + else + { + *n_vars--; + msg (MW, _("Cannot compute autocovariance for a non-numeric variable %s"), + var_to_string (vars[i])); + } + } + + /* + First data pass to get the mean and variance. + */ + get_mean_variance (*n_vars, cf, est); + get_covariance (*n_vars, cf, est, lag); +}