X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fts%2Finnovations.c;h=75b88d6cb4d9bfbd0f9a73b4b734bd0d7464cbbd;hb=bd9396a24e862f444e1d3199f3a9dfd0d65daceb;hp=131284459096b2056b79c86e92b3b290c09bfea4;hpb=4dc2ebcfd1a113b25f6997ff3b66fa52ac41158b;p=pspp-builds.git diff --git a/src/math/ts/innovations.c b/src/math/ts/innovations.c index 13128445..75b88d6c 100644 --- a/src/math/ts/innovations.c +++ b/src/math/ts/innovations.c @@ -97,26 +97,24 @@ get_covariance (const gsl_matrix *data, if (!gsl_isnan (x)) { x -= est[j]->mean; - for (lag = 1; lag <= max_lag && lag < data->size1 - max_lag; lag++) + 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) += y * x; - est[i]->n_obs += 1.0; + *(est[j]->cov + lag - 1) += y * x; + est[j]->n_obs += 1.0; } } } } } - for (lag = 0; lag <= max_lag && lag < data->size1 - max_lag; lag++) + for (j = 0; j < data->size2; j++) { - for (j = 0; j < data->size2; j++) - { - *(est[j]->cov + lag) /= (est[j]->n_obs - lag); - } + *(est[j]->cov + lag - 1) /= est[j]->n_obs; } + return rc; } static double @@ -126,9 +124,9 @@ innovations_convolve (double **theta, struct innovations_estimate *est, int k; double result = 0.0; - for (k = 0; k < i; k++) + for (k = 0; k < j; k++) { - result += theta[i-1][i-k-1] * theta[j-1][j-k-1] * est->scale[k]; + result += theta[i-1][i-k-1] * theta[j][j-k-1] * est->scale[k]; } return result; } @@ -140,54 +138,70 @@ innovations_update_scale (struct innovations_estimate *est, double *theta, size_t j; size_t k; + if (i < (size_t) est->max_lag) + { + result = est->variance; + for (j = 0; j < i; j++) + { + k = i - j - 1; + result -= theta[k] * theta[k] * est->scale[j]; + } + est->scale[i] = result; + } +} +static void +init_theta (double **theta, size_t max_lag) +{ + size_t i; + size_t j; - result = est->cov[0]; - for (j = 0; j < i; j++) + for (i = 0; i < max_lag; i++) { - k = i - j; - result -= theta[k] * theta[k] * est->scale[j]; + for (j = 0; j <= i; j++) + { + theta[i][j] = 0.0; + } } - est->scale[i] = result; } +static void +innovations_update_coeff (double **theta, struct innovations_estimate *est, + size_t max_lag) +{ + size_t i; + size_t j; + size_t k; + for (i = 0; i < max_lag; i++) + { + for (j = 0; j <= i; j++) + { + k = i - j; + theta[i][k] = (est->cov[k] - + innovations_convolve (theta, est, i, j)) + / est->scale[k]; + } + innovations_update_scale (est, theta[i], i + 1); + } +} static void get_coef (const gsl_matrix *data, struct innovations_estimate **est, size_t max_lag) { - size_t j; size_t i; - size_t k; size_t n; - double v; double **theta; 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++) { - for (i = 0; i < max_lag; i++) - { - for (j = 0; j < i; j++) - { - theta[i][j] = 0.0; - } - } + init_theta (theta, max_lag); innovations_update_scale (est[n], theta[0], 0); - for (i = 0; i < max_lag; i++) - { - v = est[n]->cov[i]; - for (j = 0; j < i; j++) - { - k = i - j; - theta[i-1][k-1] = est[n]->cov[k] - - innovations_convolve (theta, est[n], i, j); - } - innovations_update_scale (est[n], theta[i], i); - } + innovations_update_coeff (theta, est[n], max_lag); /* Copy the final row of coefficients into EST->COEFF.*/ for (i = 0; i < max_lag; i++) { @@ -213,6 +227,7 @@ get_coef (const gsl_matrix *data, pspp_coeff_set_estimate (est[n]->coeff[i], theta[max_lag - 1][i]); } } + for (i = 0; i < max_lag; i++) { free (theta[i]); @@ -220,32 +235,69 @@ get_coef (const gsl_matrix *data, free (theta); } +static void +innovations_struct_init (struct innovations_estimate *est, size_t lag) +{ + size_t j; + + est->mean = 0.0; + est->variance = 0.0; + est->cov = xnmalloc (lag, sizeof (*est->cov)); + est->scale = xnmalloc (lag + 1, sizeof (*est->scale)); + est->coeff = xnmalloc (lag, sizeof (*est->coeff)); + est->max_lag = (double) lag; + /* 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. */ + for (j = 0; j < lag; j++) + { + est->coeff[j] = xmalloc (sizeof (*(est->coeff[j]))); + } +} + 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); + est[i] = xmalloc (sizeof *est[i]); /* 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]->scale = xnmalloc (lag, sizeof (*est[i]->scale)); - 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))); - } + innovations_struct_init (est[i], lag); } - get_mean_variance (data, est); - get_covariance (data, est, lag); - get_coef (data, est, lag); + get_mean_variance (dm->m, est); + get_covariance (dm->m, est, lag); + get_coef (dm->m, est, lag); return est; } + +static void +pspp_innovations_free_one (struct innovations_estimate *est) +{ + size_t i; + + assert (est != NULL); + 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) +{ + size_t i; + + assert (est != NULL); + for (i = 0; i < n; i++) + { + pspp_innovations_free_one (est[i]); + } + free (est); +}