X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fts%2Finnovations.c;h=00d366cef7bbd9fe617209c949e9a655f1ff2407;hb=39929aa077830c708adcc5dfd224fd973428d0bc;hp=3921ea1c7e4318bb09d29e999d39f3d23d11864f;hpb=fb67d7201c231228f4a528f0c10d1cd17cddd6c5;p=pspp-builds.git diff --git a/src/math/ts/innovations.c b/src/math/ts/innovations.c index 3921ea1c..00d366ce 100644 --- a/src/math/ts/innovations.c +++ b/src/math/ts/innovations.c @@ -38,8 +38,8 @@ #include static void -get_mean_variance (const gsl_matrix *data, - struct innovations_estimate **est) +get_mean (const gsl_matrix *data, + struct innovations_estimate **est) { size_t n; @@ -49,9 +49,8 @@ get_mean_variance (const gsl_matrix *data, for (n = 0; n < data->size2; n++) { - est[n]->n_obs = 2.0; + est[n]->n_obs = 0.0; est[n]->mean = 0.0; - est[n]->variance = 0.0; } for (i = 0; i < data->size1; i++) { @@ -60,20 +59,36 @@ get_mean_variance (const gsl_matrix *data, tmp = gsl_matrix_get (data, i, n); if (!gsl_isnan (tmp)) { + est[n]->n_obs += 1.0; 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 < data->size2; n++) +} +static void +update_cov (struct innovations_estimate **est, gsl_vector_const_view x, + gsl_vector_const_view y, size_t lag) +{ + size_t j; + double xj; + double yj; + + for (j = 0; j < x.vector.size; j++) { - /* Maximum likelihood estimate of the variance. */ - est[n]->variance /= est[n]->n_obs; + xj = gsl_vector_get (&x.vector, j); + yj = gsl_vector_get (&y.vector, j); + if (!gsl_isnan (xj)) + { + if (!gsl_isnan (yj)) + { + xj -= est[j]->mean; + yj -= est[j]->mean; + *(est[j]->cov + lag) += xj * yj; + } + } } } - static int get_covariance (const gsl_matrix *data, struct innovations_estimate **est, size_t max_lag) @@ -81,54 +96,55 @@ get_covariance (const gsl_matrix *data, size_t lag; size_t j; size_t i; - double x; - double y; int rc = 1; assert (data != NULL); assert (est != NULL); - + + for (j = 0; j < data->size2; j++) + { + for (lag = 0; lag <= max_lag; lag++) + { + *(est[j]->cov + lag) = 0.0; + } + } + /* + The rows are in the outer loop because a gsl_matrix is stored in + row-major order. + */ for (i = 0; i < data->size1; i++) { - for (j = 0; j < data->size2; j++) + for (lag = 0; lag <= max_lag && lag < data->size1 - i; lag++) { - x = gsl_matrix_get (data, i, j); - - if (!gsl_isnan (x)) - { - x -= est[j]->mean; - for (lag = 1; lag <= max_lag && lag < (data->size1 - i); lag++) - { - y = gsl_matrix_get (data, i + lag, j); - if (!gsl_isnan (y)) - { - y -= est[j]->mean; - *(est[j]->cov + lag - 1) += y * x; - est[j]->n_obs += 1.0; - } - } - } + update_cov (est, gsl_matrix_const_row (data, i), + gsl_matrix_const_row (data, i + lag), lag); } } - for (lag = 1; lag <= max_lag; lag++) + for (j = 0; j < data->size2; j++) { - for (j = 0; j < data->size2; j++) + for (lag = 0; lag <= max_lag; lag++) { - *(est[j]->cov + lag) /= (est[j]->n_obs - lag); + *(est[j]->cov + lag) /= est[j]->n_obs; } } + return rc; } + static double -innovations_convolve (double **theta, struct innovations_estimate *est, - int i, int j) +innovations_convolve (double *x, double *y, struct innovations_estimate *est, + int i) { int k; double result = 0.0; - for (k = 0; k < j; k++) + assert (x != NULL && y != NULL); + assert (est != NULL); + assert (est->scale != NULL); + assert (i > 0); + for (k = 0; k < i; k++) { - result += theta[i-1][i-k-1] * theta[j][j-k-1] * est->scale[k]; + result += x[k] * y[k] * est->scale[i-k-1]; } return result; } @@ -140,14 +156,16 @@ innovations_update_scale (struct innovations_estimate *est, double *theta, size_t j; size_t k; - - result = est->variance; - for (j = 0; j < i; j++) + if (i < (size_t) est->max_lag) { - k = i - j - 1; - result -= theta[k] * theta[k] * est->scale[j]; + result = est->cov[0]; + for (j = 0; j < i; j++) + { + k = i - j - 1; + result -= theta[k] * theta[k] * est->scale[j]; + } + est->scale[i] = result; } - est->scale[i] = result; } static void init_theta (double **theta, size_t max_lag) @@ -170,16 +188,16 @@ innovations_update_coeff (double **theta, struct innovations_estimate *est, size_t i; size_t j; size_t k; - double v; for (i = 0; i < max_lag; i++) { - for (j = 0; j <= i; j++) + theta[i][i] = est->cov[i+1] / est->scale[0]; + for (j = 1; j <= i; j++) { k = i - j; - theta[i][k] = (est->cov[k] - - innovations_convolve (theta, est, i, j)) - / est->scale[k]; + 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); } @@ -195,7 +213,7 @@ get_coef (const gsl_matrix *data, theta = xnmalloc (max_lag, sizeof (*theta)); for (i = 0; i < max_lag; i++) { - theta[i] = xnmalloc (i + 1, sizeof (**(theta + i))); + theta[i] = xnmalloc (max_lag, sizeof (**(theta + i))); } for (n = 0; n < data->size2; n++) @@ -236,35 +254,51 @@ get_coef (const gsl_matrix *data, free (theta); } +static void +innovations_struct_init (struct innovations_estimate *est, + const struct design_matrix *dm, + size_t lag) +{ + size_t j; + + est->mean = 0.0; + /* COV[0] stores the lag 0 covariance (i.e., the variance), COV[1] + holds the lag-1 covariance, etc. + */ + est->cov = xnmalloc (lag + 1, sizeof (*est->cov)); + est->scale = xnmalloc (lag + 1, sizeof (*est->scale)); + est->coeff = xnmalloc (lag, sizeof (*est->coeff)); /* No intercept. */ + + /* + The loop below is an unusual use of PSPP_COEFF_INIT(). In a + typical model, one column of a DESIGN_MATRIX has one + coefficient. But in a time-series model, one column has many + coefficients. + */ + for (j = 0; j < lag; j++) + { + pspp_coeff_init (est->coeff + j, dm); + } + est->max_lag = (double) lag; +} + struct innovations_estimate ** -pspp_innovations (const gsl_matrix *data, size_t lag) +pspp_innovations (const struct design_matrix *dm, size_t lag) { struct innovations_estimate **est; size_t i; - size_t j; - est = xnmalloc (data->size2, sizeof *est); - for (i = 0; i < data->size2; i++) + est = xnmalloc (dm->m->size2, sizeof *est); + for (i = 0; i < dm->m->size2; i++) { est[i] = xmalloc (sizeof *est[i]); /* est[i]->variable = vars[i]; */ - est[i]->mean = 0.0; - est[i]->variance = 0.0; - /* COV does not the variance (i.e., the lag 0 covariance). So COV[0] - holds the lag 1 covariance, COV[i] holds the lag i+1 covariance. */ - est[i]->cov = xnmalloc (lag, sizeof (*est[i]->cov)); - est[i]->scale = xnmalloc (lag, sizeof (*est[i]->scale)); - est[i]->coeff = xnmalloc (lag, sizeof (*est[i]->coeff)); - est[i]->max_lag = (double) lag; - for (j = 0; j < lag; j++) - { - est[i]->coeff[j] = xmalloc (sizeof (*(est[i]->coeff + j))); - } + innovations_struct_init (est[i], dm, lag); } - get_mean_variance (data, est); - get_covariance (data, est, lag); - get_coef (data, est, lag); + get_mean (dm->m, est); + get_covariance (dm->m, est, lag); + get_coef (dm->m, est, lag); return est; } @@ -275,12 +309,13 @@ pspp_innovations_free_one (struct innovations_estimate *est) size_t i; assert (est != NULL); - free (est->cov); - free (est->scale); for (i = 0; i < (size_t) est->max_lag; i++) { pspp_coeff_free (est->coeff[i]); } + free (est->scale); + free (est->cov); + free (est); } void pspp_innovations_free (struct innovations_estimate **est, size_t n)