pomp
Inference for partially observed Markov processes
•All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
probe_acf.c File Reference
#include "internal.h"
#include <stdio.h>
Include dependency graph for probe_acf.c:

Go to the source code of this file.

Functions

static void pomp_acf_compute (double *acf, double *x, int n, int nvars, int *lags, int nlag)
 
static void pomp_ccf_compute (double *ccf, double *x, double *y, int n, int *lags, int nlag)
 
SEXP probe_acf (SEXP x, SEXP lags, SEXP corr)
 
SEXP probe_ccf (SEXP x, SEXP y, SEXP lags, SEXP corr)
 

Function Documentation

◆ pomp_acf_compute()

static void pomp_acf_compute ( double *  acf,
double *  x,
int  n,
int  nvars,
int *  lags,
int  nlag 
)
static

Definition at line 11 of file probe_acf.c.

11 {
12 double *p, *p0, *p1, *p2;
13 long double xx;
14 int i, j, k, lag, ct;
15
16 // first center each row
17 for (j = 0, p = x; j < nvars; j++, p++) {
18 for (k = 0, p0 = p, xx = 0, ct = 0; k < n; p0 += nvars, k++) {
19 if (R_FINITE(*p0)) {
20 xx += *p0;
21 ct++;
22 }
23 }
24 if (ct < 1) err("series %d has no data",j+1);
25 xx /= ct; // mean of x[j,]
26 for (k = 0, p0 = p; k < n; p0 += nvars, k++)
27 if (R_FINITE(*p0)) *p0 -= xx;
28 }
29
30 // compute covariances
31 for (j = 0, p0 = x, p = acf; j < nvars; j++, p0++) { // loop over series
32 for (i = 0; i < nlag; i++, p++) { // loop over lags
33 lag = lags[i]; // i-th lag
34 for (k = 0, ct = 0, xx = 0, p1 = p0, p2 = p0+lag*nvars; k < n-lag; k++, p1 += nvars, p2 += nvars)
35 if (R_FINITE(*p1) && R_FINITE(*p2)) {
36 xx += (*p1)*(*p2);
37 ct++;
38 }
39 *p = (ct > 0) ? xx/ct : R_NaReal;
40 }
41 }
42
43}
err
#define err(...)
Definition pomp.h:21
nvars
int nvars
Definition trajectory.c:131
Here is the caller graph for this function:

◆ pomp_ccf_compute()

static void pomp_ccf_compute ( double *  ccf,
double *  x,
double *  y,
int  n,
int *  lags,
int  nlag 
)
static

Definition at line 47 of file probe_acf.c.

47 {
48 double *p, *p1, *p2;
49 long double xx;
50 int j, k, lag, ct;
51
52 // first center x
53 for (k = 0, xx = 0, ct = 0, p = x; k < n; k++, p++) {
54 if (R_FINITE(*p)) {
55 xx += *p;
56 ct++;
57 }
58 }
59 if (ct < 1) err("series 1 has no data");
60 xx /= ct; // mean of x[j]
61 for (k = 0, p = x; k < n; k++, p++) {
62 if (R_FINITE(*p)) *p -= xx;
63 }
64
65 // now center y
66 for (k = 0, xx = 0, ct = 0, p = y; k < n; k++, p++) {
67 if (R_FINITE(*p)) {
68 xx += *p;
69 ct++;
70 }
71 }
72 if (ct < 1) err("series 2 has no data");
73 xx /= ct; // mean of y[j]
74 for (k = 0, p = y; k < n; k++, p++) {
75 if (R_FINITE(*p)) *p -= xx;
76 }
77
78 // compute covariances
79 for (j = 0, p = ccf; j < nlag; j++, p++) { // loop over lags
80 lag = lags[j];
81 if (lag < 0) {
82 for (k = 0, xx = 0, ct = 0, p1 = x-lag, p2 = y; k < n+lag; k++, p1++, p2++)
83 if (R_FINITE(*p1) && R_FINITE(*p1)) {
84 xx += (*p1)*(*p2);
85 ct++;
86 }
87 *p = (ct > 0) ? xx/ct : R_NaReal;
88 } else {
89 for (k = 0, xx = 0, ct = 0, p1 = x, p2 = y+lag; k < n-lag; k++, p1++, p2++)
90 if (R_FINITE(*p1) && R_FINITE(*p1)) {
91 xx += (*p1)*(*p2);
92 ct++;
93 }
94 *p = (ct > 0) ? xx/ct : R_NaReal;
95 }
96 }
97
98}
Here is the caller graph for this function:

◆ probe_acf()

SEXP probe_acf ( SEXP  x,
SEXP  lags,
SEXP  corr 
)

Definition at line 100 of file probe_acf.c.

100 {
101 SEXP ans, ans_names;
102 SEXP X;
103 int nlag, correlation, nvars, n;
104 int j, k, l;
105 double *p, *p1, *cov;
106 int *lag;
107 char tmp[BUFSIZ];
108
109 nlag = LENGTH(lags); // number of lags
110 PROTECT(lags = AS_INTEGER(lags));
111 lag = INTEGER(lags);
112 correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
113
114 nvars = INTEGER(GET_DIM(x))[0]; // number of variables
115 n = INTEGER(GET_DIM(x))[1]; // number of observations
116
117 PROTECT(X = duplicate(AS_NUMERIC(x)));
118
119 PROTECT(ans = NEW_NUMERIC(nlag*nvars));
120
121 pomp_acf_compute(REAL(ans),REAL(X),n,nvars,lag,nlag);
122
123 if (correlation) {
124 l = 0;
125 cov = (double *) R_alloc(nvars,sizeof(double));
126 pomp_acf_compute(cov,REAL(X),n,nvars,&l,1); // compute lag-0 covariance
127 for (j = 0, p = REAL(ans), p1 = cov; j < nvars; j++, p1++)
128 for (k = 0; k < nlag; k++, p++)
129 *p /= *p1;
130 }
131
132 PROTECT(ans_names = NEW_STRING(nlag*nvars));
133 for (j = 0, l = 0; j < nvars; j++) {
134 for (k = 0; k < nlag; k++, l++) {
135 snprintf(tmp,BUFSIZ,"acf[%d]",lag[k]);
136 SET_STRING_ELT(ans_names,l,mkChar(tmp));
137 }
138 }
139 SET_NAMES(ans,ans_names);
140
141 UNPROTECT(4);
142 return ans;
143}
X
#define X
Definition gompertz.c:14
pomp_acf_compute
static void pomp_acf_compute(double *acf, double *x, int n, int nvars, int *lags, int nlag)
Definition probe_acf.c:11
cov
SEXP cov
Definition trajectory.c:129
Here is the call graph for this function:

◆ probe_ccf()

SEXP probe_ccf ( SEXP  x,
SEXP  y,
SEXP  lags,
SEXP  corr 
)

Definition at line 145 of file probe_acf.c.

145 {
146 SEXP ans, ans_names;
147 SEXP X, Y;
148 double cov[2], xx, *p;
149 int nlag, n, correlation;
150 int k;
151 char tmp[BUFSIZ];
152
153 nlag = LENGTH(lags);
154 PROTECT(lags = AS_INTEGER(lags));
155 correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
156
157 n = LENGTH(x); // n = # of observations
158 if (n != LENGTH(y)) err("'x' and 'y' must have equal lengths"); // #nocov
159
160 PROTECT(X = duplicate(AS_NUMERIC(x)));
161 PROTECT(Y = duplicate(AS_NUMERIC(y)));
162
163 PROTECT(ans = NEW_NUMERIC(nlag));
164
165 pomp_ccf_compute(REAL(ans),REAL(X),REAL(Y),n,INTEGER(lags),nlag);
166
167 if (correlation) {
168 k = 0;
169 pomp_acf_compute(&cov[0],REAL(X),n,1,&k,1); // compute lag-0 covariance of x
170 pomp_acf_compute(&cov[1],REAL(Y),n,1,&k,1); // compute lag-0 covariance of y
171 xx = sqrt(cov[0]*cov[1]);
172 for (k = 0, p = REAL(ans); k < nlag; k++, p++) *p /= xx; // convert to correlation
173 }
174
175 PROTECT(ans_names = NEW_STRING(nlag));
176 for (k = 0; k < nlag; k++) {
177 snprintf(tmp,BUFSIZ,"ccf[%d]",INTEGER(lags)[k]);
178 SET_STRING_ELT(ans_names,k,mkChar(tmp));
179 }
180 SET_NAMES(ans,ans_names);
181
182 UNPROTECT(5);
183 return ans;
184}
Y
#define Y
Definition gompertz.c:13
pomp_ccf_compute
static void pomp_ccf_compute(double *ccf, double *x, double *y, int n, int *lags, int nlag)
Definition probe_acf.c:47
Here is the call graph for this function: