Stan Math Library  2.15.0
reverse mode automatic differentiation
autocorrelation.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_AUTOCORRELATION_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_AUTOCORRELATION_HPP
3 
5 #include <unsupported/Eigen/FFT>
6 #include <complex>
7 #include <vector>
8 
9 namespace stan {
10  namespace math {
11 
12  namespace {
17  size_t fft_next_good_size(size_t N) {
18  if (N <= 2) return 2;
19  while (true) {
20  size_t m = N;
21  while ((m % 2) == 0) m /= 2;
22  while ((m % 3) == 0) m /= 3;
23  while ((m % 5) == 0) m /= 5;
24  if (m <= 1)
25  return N;
26  N++;
27  }
28  }
29  }
30 
51  template <typename T>
52  void autocorrelation(const std::vector<T>& y,
53  std::vector<T>& ac,
54  Eigen::FFT<T>& fft) {
55  using std::vector;
56  using std::complex;
57 
58  size_t N = y.size();
59  size_t M = fft_next_good_size(N);
60  size_t Mt2 = 2 * M;
61 
62  vector<complex<T> > freqvec;
63 
64  // centered_signal = y-mean(y) followed by N zeroes
65  vector<T> centered_signal(y);
66  centered_signal.insert(centered_signal.end(), Mt2-N, 0.0);
67  T mean = stan::math::mean(y);
68  for (size_t i = 0; i < N; i++)
69  centered_signal[i] -= mean;
70 
71  fft.fwd(freqvec, centered_signal);
72  for (size_t i = 0; i < Mt2; ++i)
73  freqvec[i] = complex<T>(norm(freqvec[i]), 0.0);
74 
75  fft.inv(ac, freqvec);
76  ac.resize(N);
77 
78  /*
79  vector<T> mask_correction_factors;
80  vector<T> mask;
81  mask.insert(mask.end(), N, 1.0);
82  mask.insert(mask.end(), N, 0.0);
83 
84  freqvec.resize(0);
85  fft.fwd(freqvec, mask);
86  for (size_t i = 0; i < Nt2; ++i)
87  freqvec[i] = complex<T>(norm(freqvec[i]), 0.0);
88 
89  fft.inv(mask_correction_factors, freqvec);
90 
91  for (size_t i = 0; i < N; ++i) {
92  ac[i] /= mask_correction_factors[i];
93  }
94  */
95  for (size_t i = 0; i < N; ++i) {
96  ac[i] /= (N - i);
97  }
98  T var = ac[0];
99  for (size_t i = 0; i < N; ++i)
100  ac[i] /= var;
101  }
102 
119  template <typename T>
120  void autocorrelation(const std::vector<T>& y,
121  std::vector<T>& ac) {
122  Eigen::FFT<T> fft;
123  return autocorrelation(y, ac, fft);
124  }
125 
126  }
127 }
128 #endif
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:30
boost::math::tools::promote_args< T >::type mean(const std::vector< T > &v)
Returns the sample mean (i.e., average) of the coefficients in the specified standard vector...
Definition: mean.hpp:23
void autocorrelation(const std::vector< T > &y, std::vector< T > &ac, Eigen::FFT< T > &fft)
Write autocorrelation estimates for every lag for the specified input sequence into the specified res...

     [ Stan Home Page ] © 2011–2016, Stan Development Team.