From fbbc4d9655835c0d6fdf58f231e59d7007a99407 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 23 Apr 2018 18:32:49 -0400 Subject: lots of changes --- lib/yule_walker.c | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 lib/yule_walker.c (limited to 'lib/yule_walker.c') diff --git a/lib/yule_walker.c b/lib/yule_walker.c new file mode 100644 index 0000000..77b3e96 --- /dev/null +++ b/lib/yule_walker.c @@ -0,0 +1,41 @@ + +#include "yule_walker.h" + +double yule_walker(const autocorr_t *a) { + gsl_vector *rhos = gsl_vector_alloc(a->W); + gsl_matrix *mat = gsl_matrix_alloc(a->W, a->W); + gsl_vector *phis = gsl_vector_alloc(a->W); + + for (count_t i = 0; i < a->W; i++) { + gsl_vector_set(rhos, i, rho(a, i)); + } + + for (count_t i = 0; i < a->W; i++) { + gsl_matrix_set(mat, i, i, 1.0); + + for (count_t j = 0; j < i; j++) { + gsl_matrix_set(mat, i, j, gsl_vector_get(rhos, i - 1 - j)); + } + + for (count_t j = 0; j < a->W - 1 - i; j++) { + gsl_matrix_set(mat, i, i + 1 + j, gsl_vector_get(rhos, j)); + } + } + + int out = gsl_linalg_cholesky_solve(mat, rhos, phis); + + double rhophi = 0; + double onephi = 0; + + for (count_t i = 0; i < a->W; i++) { + rhophi += gsl_vector_get(rhos, i) * gsl_vector_get(phis, i); + onephi += gsl_vector_get(phis, i); + } + + gsl_vector_free(rhos); + gsl_matrix_free(mat); + gsl_vector_free(phis); + + return (1 - rhophi) / pow(1 - onephi, 2); +} + -- cgit v1.2.3-70-g09d2