1 /*
2 * R : A Computer Language for Statistical Data Analysis
3
4 * Copyright (C) 1999-2016 The R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/.
19 */
20
21 #include <R.h>
22
23 static void
burg(int n,double * x,int pmax,double * coefs,double * var1,double * var2)24 burg(int n, double*x, int pmax, double *coefs, double *var1, double *var2)
25 {
26 double d, phii, *u, *v, *u0, sum;
27
28 u = (double *) R_alloc(n, sizeof(double));
29 v = (double *) R_alloc(n, sizeof(double));
30 u0 = (double *) R_alloc(n, sizeof(double));
31
32 for(int i = 0; i < pmax*pmax; i++) coefs[i] = 0.0;
33 sum = 0.0;
34 for(int t = 0; t < n; t++) {
35 u[t] = v[t] = x[n - 1 - t];
36 sum += x[t] * x[t];
37 }
38 var1[0] = var2[0] = sum/n;
39 for(int p = 1; p <= pmax; p++) { /* do AR(p) */
40 sum = 0.0;
41 d = 0;
42 for(int t = p; t < n; t++) {
43 sum += v[t]*u[t-1];
44 d += v[t]*v[t] + u[t-1]*u[t-1];
45 }
46 phii = 2*sum/d;
47 coefs[pmax*(p-1) + (p-1)] = phii;
48 if(p > 1)
49 for(int j = 1; j < p; j++)
50 coefs[p-1 + pmax*(j-1)] =
51 coefs[p-2 + pmax*(j-1)] - phii* coefs[p-2 + pmax*(p-j-1)];
52 /* update u and v */
53 for(int t = 0; t < n; t++)
54 u0[t] = u[t];
55 for(int t = p; t < n; t++) {
56 u[t] = u0[t-1] - phii * v[t];
57 v[t] = v[t] - phii * u0[t-1];
58 }
59 var1[p] = var1[p-1] * (1 - phii * phii);
60 d = 0.0;
61 for(int t = p; t < n; t++) d += v[t]*v[t] + u[t]*u[t];
62 var2[p] = d/(2.0*(n-p));
63 }
64 }
65
66 #include <Rinternals.h>
67
Burg(SEXP x,SEXP order)68 SEXP Burg(SEXP x, SEXP order)
69 {
70 x = PROTECT(coerceVector(x, REALSXP));
71 int n = LENGTH(x), pmax = asInteger(order);
72 SEXP coefs = PROTECT(allocVector(REALSXP, pmax * pmax)),
73 var1 = PROTECT(allocVector(REALSXP, pmax + 1)),
74 var2 = PROTECT(allocVector(REALSXP, pmax + 1));
75 burg(n, REAL(x), pmax, REAL(coefs), REAL(var1), REAL(var2));
76 SEXP ans = PROTECT(allocVector(VECSXP, 3));
77 SET_VECTOR_ELT(ans, 0, coefs);
78 SET_VECTOR_ELT(ans, 1, var1);
79 SET_VECTOR_ELT(ans, 2, var2);
80 UNPROTECT(5);
81 return ans;
82 }
83