X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fts%2Finnovations.c;h=3ab2f3edce165ed0a709688e711867986683eee3;hb=5c3291dc396b795696e94f47780308fd7ace6fc4;hp=107d8ba47897e2eafde521b1ec02e1d7bfc4ad8b;hpb=755ecfd2e8d86bc134fe7202c46fee354ec166d0;p=pspp-builds.git diff --git a/src/math/ts/innovations.c b/src/math/ts/innovations.c index 107d8ba4..3ab2f3ed 100644 --- a/src/math/ts/innovations.c +++ b/src/math/ts/innovations.c @@ -1,22 +1,19 @@ -/* - src/math/ts/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. - */ +/* PSPP - a program for statistical analysis. + Copyright (C) 2006 Free Software Foundation, Inc. + + 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 3 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, see . */ + /* Find preliminary ARMA coefficients via the innovations algorithm. Also compute the sample mean and covariance matrix for each series. @@ -28,19 +25,24 @@ 0-387-97429-6. Sections 5.2, 8.3 and 8.4. */ +#include + +#include +#include + #include #include -#include -#include -#include #include +#include #include #include +#include "xalloc.h" + static void get_mean (const gsl_matrix *data, struct innovations_estimate **est) - + { size_t n; size_t i; @@ -57,7 +59,7 @@ get_mean (const gsl_matrix *data, for (n = 0; n < data->size2; n++) { tmp = gsl_matrix_get (data, i, n); - if (!gsl_isnan (tmp)) + if (!isnan (tmp)) { est[n]->n_obs += 1.0; d = (tmp - est[n]->mean) / est[n]->n_obs; @@ -66,7 +68,7 @@ get_mean (const gsl_matrix *data, } } } -static void +static void update_cov (struct innovations_estimate **est, gsl_vector_const_view x, gsl_vector_const_view y, size_t lag) { @@ -78,9 +80,9 @@ update_cov (struct innovations_estimate **est, gsl_vector_const_view x, { xj = gsl_vector_get (&x.vector, j); yj = gsl_vector_get (&y.vector, j); - if (!gsl_isnan (xj)) + if (!isnan (xj)) { - if (!gsl_isnan (yj)) + if (!isnan (yj)) { xj -= est[j]->mean; yj -= est[j]->mean; @@ -90,7 +92,7 @@ update_cov (struct innovations_estimate **est, gsl_vector_const_view x, } } static int -get_covariance (const gsl_matrix *data, +get_covariance (const gsl_matrix *data, struct innovations_estimate **est, size_t max_lag) { size_t lag; @@ -116,7 +118,7 @@ get_covariance (const gsl_matrix *data, { for (lag = 0; lag <= max_lag && lag < data->size1 - i; lag++) { - update_cov (est, gsl_matrix_const_row (data, i), + update_cov (est, gsl_matrix_const_row (data, i), gsl_matrix_const_row (data, i + lag), lag); } } @@ -162,7 +164,7 @@ innovations_update_scale (struct innovations_estimate *est, double *theta, for (j = 0; j < i; j++) { k = i - j - 1; - result -= theta[k] * theta[k] * est->scale[j]; + result -= pow2 (theta[k]) * est->scale[j]; } est->scale[i] = result; } @@ -195,12 +197,12 @@ innovations_update_coeff (double **theta, struct innovations_estimate *est, for (j = 1; j <= i; j++) { k = i - j; - theta[i][k] = (est->cov[k+1] - + theta[i][k] = (est->cov[k+1] - innovations_convolve (theta[i] + k + 1, theta[j - 1], est, j)) / est->scale[j]; } innovations_update_scale (est, theta[i], i + 1); - } + } } static void get_coef (const gsl_matrix *data, @@ -231,7 +233,7 @@ get_coef (const gsl_matrix *data, Let X[m], X[m-1],... denote the original series. Let X_hat[0] denote the best predicted value of X[0], X_hat[1] denote the projection of X[1] onto the subspace - spanned by {X[0] - X_hat[0]}. Let X_hat[m] denote the + spanned by {X[0] - X_hat[0]}. Let X_hat[m] denote the projection of X[m] onto the subspace spanned by {X[m-1] - X_hat[m-1], X[m-2] - X_hat[m-2],...,X[0] - X_hat[0]}. @@ -252,8 +254,8 @@ get_coef (const gsl_matrix *data, } static void -innovations_struct_init (struct innovations_estimate *est, - const struct design_matrix *dm, +innovations_struct_init (struct innovations_estimate *est, + const struct design_matrix *dm, size_t lag) { size_t j; @@ -286,7 +288,7 @@ innovations_struct_init (struct innovations_estimate *est, */ static void subtract_mean (gsl_matrix *m, struct innovations_estimate **est) -{ +{ size_t i; size_t j; double tmp; @@ -300,7 +302,7 @@ subtract_mean (gsl_matrix *m, struct innovations_estimate **est) } } } -struct innovations_estimate ** +struct innovations_estimate ** pspp_innovations (const struct design_matrix *dm, size_t lag) { struct innovations_estimate **est; @@ -318,11 +320,11 @@ pspp_innovations (const struct design_matrix *dm, size_t lag) subtract_mean (dm->m, est); get_covariance (dm->m, est, lag); get_coef (dm->m, est, lag); - + return est; } -static void +static void pspp_innovations_free_one (struct innovations_estimate *est) { size_t i;