Class: Statsample::TimeSeries::Pacf::Pacf
- Inherits:
-
Object
- Object
- Statsample::TimeSeries::Pacf::Pacf
- Defined in:
- lib/bio-statsample-timeseries/timeseries/pacf.rb
Class Method Summary collapse
- .diag(mat) ⇒ Object
-
.levinson_durbin(series, nlags = 10, is_acovf = false) ⇒ Object
returns: -sigma_v: estimate of the error variance -arcoefs: AR coefficients -pacf: pacf function -sigma: some function.
- .pacf_yw(timeseries, max_lags, method = 'yw') ⇒ Object
- .solve_matrix(matrix, out_vector) ⇒ Object
- .toeplitz(arr) ⇒ Object
-
.yule_walker(ts, k = 1, method = 'yw') ⇒ Object
- returns: -rho
- autoregressive coefficients -sigma
-
sigma parameter.
Class Method Details
.diag(mat) ⇒ Object
63 64 65 66 67 |
# File 'lib/bio-statsample-timeseries/timeseries/pacf.rb', line 63 def self.diag(mat) #returns array of diagonal elements of a matrix. #will later abstract it to matrix.rb in Statsample return mat.each_with_index(:diagonal).map { |x, r, c| x } end |
.levinson_durbin(series, nlags = 10, is_acovf = false) ⇒ Object
returns: -sigma_v: estimate of the error variance -arcoefs: AR coefficients -pacf: pacf function -sigma: some function
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 |
# File 'lib/bio-statsample-timeseries/timeseries/pacf.rb', line 28 def self.levinson_durbin(series, nlags = 10, is_acovf = false) if is_acovf series = series.map(&:to_f) else #nlags = order(k) of AR in this case series = series.acvf.map(&:to_f)[0..nlags] end #phi = Array.new((nlags+1), 0.0) { Array.new(nlags+1, 0.0) } order = nlags phi = Matrix.zero(nlags + 1) sig = Array.new(nlags+1) #setting initial point for recursion: phi[1,1] = series[1]/series[0] #phi[1][1] = series[1]/series[0] sig[1] = series[0] - phi[1, 1] * series[1] 2.upto(order).each do |k| phi[k, k] = (series[k] - (Statsample::Vector.new(phi[1...k, k-1]) * series[1...k].reverse.to_ts).sum) / sig[k-1] #some serious refinement needed in above for matrix manipulation. Will do today 1.upto(k-1).each do |j| phi[j, k] = phi[j, k-1] - phi[k, k] * phi[k-j, k-1] end sig[k] = sig[k-1] * (1-phi[k, k] ** 2) end sigma_v = sig[-1] arcoefs_delta = phi.column(phi.column_size - 1) arcoefs = arcoefs_delta[1..arcoefs_delta.size] pacf = diag(phi) pacf[0] = 1.0 return [sigma_v, arcoefs, pacf, sig, phi] end |
.pacf_yw(timeseries, max_lags, method = 'yw') ⇒ Object
6 7 8 9 10 11 12 13 14 |
# File 'lib/bio-statsample-timeseries/timeseries/pacf.rb', line 6 def self.pacf_yw(timeseries, max_lags, method = 'yw') #partial autocorrelation by yule walker equations. #Inspiration: StatsModels pacf = [1.0] (1..max_lags).map do |i| pacf << yule_walker(timeseries, i, method)[0][-1] end pacf end |
.solve_matrix(matrix, out_vector) ⇒ Object
143 144 145 146 147 148 149 150 151 152 153 154 |
# File 'lib/bio-statsample-timeseries/timeseries/pacf.rb', line 143 def self.solve_matrix(matrix, out_vector) solution_vector = Array.new(out_vector.size, 0) matrix = matrix.to_a k = 0 matrix.each do |row| row.each_with_index do |element, i| solution_vector[k] += element * 1.0 * out_vector[i] end k += 1 end solution_vector end |
.toeplitz(arr) ⇒ Object
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 |
# File 'lib/bio-statsample-timeseries/timeseries/pacf.rb', line 113 def self.toeplitz(arr) #Generates Toeplitz matrix - #http://en.wikipedia.org/wiki/Toeplitz_matrix #Toeplitz matrix are equal when they are stored in row & #column major #=> arr = [0, 1, 2, 3] #=> result: #[[0, 1, 2, 3], # [1, 0, 1, 2], # [2, 1, 0, 1], # [3, 2, 1, 0]] eplitz_matrix = Array.new(arr.size) { Array.new(arr.size) } 0.upto(arr.size - 1) do |i| j = 0 index = i while i >= 0 do eplitz_matrix[index][j] = arr[i] j += 1 i -= 1 end i = index + 1; k = 1 while i < arr.size do eplitz_matrix[index][j] = arr[k] i += 1; j += 1; k += 1 end end eplitz_matrix end |
.yule_walker(ts, k = 1, method = 'yw') ⇒ Object
returns:
- -rho
-
autoregressive coefficients
- -sigma
-
sigma parameter
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 |
# File 'lib/bio-statsample-timeseries/timeseries/pacf.rb', line 83 def self.yule_walker(ts, k = 1, method='yw') ts = ts - ts.mean n = ts.size if method.downcase.eql? 'yw' #unbiased => denominator = (n - k) denom =->(k) { n - k } else #mle #denominator => (n) denom =->(k) { n } end r = Array.new(k + 1) { 0.0 } r[0] = ts.map { |x| x**2 }.inject(:+).to_f / denom.call(0).to_f 1.upto(k) do |l| r[l] = (ts[0...-l].zip(ts[l...ts.size])).map do |x| x.inject(:*) end.inject(:+).to_f / denom.call(l).to_f end r_R = toeplitz(r[0...-1]) mat = Matrix.columns(r_R).inverse() phi = solve_matrix(mat, r[1..r.size]) phi_vector = Statsample::Vector.new(phi, :scale) r_vector = Statsample::Vector.new(r[1..r.size], :scale) sigma = r[0] - (r_vector * phi_vector).sum return [phi, sigma] end |