diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-04-23 18:32:49 -0400 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-04-23 18:32:49 -0400 |
commit | fbbc4d9655835c0d6fdf58f231e59d7007a99407 (patch) | |
tree | cf534240022ef4f3f39f3e865dbfa6032666f04e /lib/yule_walker.c | |
parent | c9f253b3987ee3e5ae88a99ddd2757801026efb2 (diff) | |
download | c++-fbbc4d9655835c0d6fdf58f231e59d7007a99407.tar.gz c++-fbbc4d9655835c0d6fdf58f231e59d7007a99407.tar.bz2 c++-fbbc4d9655835c0d6fdf58f231e59d7007a99407.zip |
lots of changes
Diffstat (limited to 'lib/yule_walker.c')
-rw-r--r-- | lib/yule_walker.c | 41 |
1 files changed, 41 insertions, 0 deletions
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); +} + |