-
-
Notifications
You must be signed in to change notification settings - Fork 368
/
autocovariance.hpp
135 lines (120 loc) · 4.46 KB
/
autocovariance.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#ifndef STAN_ANALYZE_MCMC_AUTOCOVARIANCE_HPP
#define STAN_ANALYZE_MCMC_AUTOCOVARIANCE_HPP
#include <stan/math/prim.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/variance.hpp>
#include <unsupported/Eigen/FFT>
#include <complex>
#include <vector>
namespace stan {
namespace analyze {
/**
* Write autocorrelation estimates for every lag for the specified
* input sequence into the specified result using the specified FFT
* engine. Normalizes lag-k autocorrelation estimators by N instead
* of (N - k), yielding biased but more stable estimators as
* discussed in Geyer (1992); see
* https://projecteuclid.org/euclid.ss/1177011137. The return vector
* will be resized to the same length as the input sequence with
* lags given by array index.
*
* <p>The implementation involves a fast Fourier transform,
* followed by a normalization, followed by an inverse transform.
*
* <p>An FFT engine can be created for reuse for type double with:
*
* <pre>
* Eigen::FFT<double> fft;
* </pre>
*
* @tparam T Scalar type.
* @param y Input sequence.
* @param ac Autocorrelations.
* @param fft FFT engine instance.
*/
template <typename T, typename DerivedA, typename DerivedB>
void autocorrelation(const Eigen::MatrixBase<DerivedA>& y,
Eigen::MatrixBase<DerivedB>& ac, Eigen::FFT<T>& fft) {
size_t N = y.size();
size_t M = math::internal::fft_next_good_size(N);
size_t Mt2 = 2 * M;
// centered_signal = y-mean(y) followed by N zeros
Eigen::Matrix<T, Eigen::Dynamic, 1> centered_signal(Mt2);
centered_signal.setZero();
centered_signal.head(N) = y.array() - y.mean();
Eigen::Matrix<std::complex<T>, Eigen::Dynamic, 1> freqvec(Mt2);
fft.fwd(freqvec, centered_signal);
// cwiseAbs2 == norm
freqvec = freqvec.cwiseAbs2();
Eigen::Matrix<std::complex<T>, Eigen::Dynamic, 1> ac_tmp(Mt2);
fft.inv(ac_tmp, freqvec);
// use "biased" estimate as recommended by Geyer (1992)
ac = ac_tmp.head(N).real().array() / (N * N * 2);
ac /= ac(0);
}
/**
* Write autocovariance estimates for every lag for the specified
* input sequence into the specified result using the specified FFT
* engine. Normalizes lag-k autocovariance estimators by N instead
* of (N - k), yielding biased but more stable estimators as
* discussed in Geyer (1992); see
* https://projecteuclid.org/euclid.ss/1177011137. The return vector
* will be resized to the same length as the input sequence with
* lags given by array index.
*
* <p>The implementation involves a fast Fourier transform,
* followed by a normalization, followed by an inverse transform.
*
* <p>This method is just a light wrapper around the three-argument
* autocovariance function
*
* @tparam T Scalar type.
* @param y Input sequence.
* @param acov Autocovariances.
*/
template <typename T, typename DerivedA, typename DerivedB>
void autocovariance(const Eigen::MatrixBase<DerivedA>& y,
Eigen::MatrixBase<DerivedB>& acov) {
Eigen::FFT<T> fft;
autocorrelation(y, acov, fft);
using boost::accumulators::accumulator_set;
using boost::accumulators::stats;
using boost::accumulators::tag::variance;
accumulator_set<double, stats<variance>> acc;
for (int n = 0; n < y.size(); ++n) {
acc(y(n));
}
acov = acov.array() * boost::accumulators::variance(acc);
}
/**
* Write autocovariance estimates for every lag for the specified
* input sequence into the specified result using the specified FFT
* engine. Normalizes lag-k autocovariance estimators by N instead
* of (N - k), yielding biased but more stable estimators as
* discussed in Geyer (1992); see
* https://projecteuclid.org/euclid.ss/1177011137. The return vector
* will be resized to the same length as the input sequence with
* lags given by array index.
*
* <p>The implementation involves a fast Fourier transform,
* followed by a normalization, followed by an inverse transform.
*
* <p>This method is just a light wrapper around the three-argument
* autocovariance function
*
* @tparam T Scalar type.
* @param y Input sequence.
* @param acov Autocovariances.
*/
template <typename T>
void autocovariance(const std::vector<T>& y, std::vector<T>& acov) {
size_t N = y.size();
acov.resize(N);
const Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, 1>> y_map(&y[0], N);
Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1>> acov_map(&acov[0], N);
autocovariance<T>(y_map, acov_map);
}
} // namespace analyze
} // namespace stan
#endif