Class: Statsample::TimeSeries::Arima::KF::LogLikelihood

Inherits:
Object
  • Object
show all
Defined in:
lib/bio-statsample-timeseries/arima/likelihood.rb

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(params, timeseries, p, q) ⇒ LogLikelihood

Returns a new instance of LogLikelihood.



20
21
22
23
24
25
26
 
# File 'lib/bio-statsample-timeseries/arima/likelihood.rb', line 20

def initialize(params, timeseries, p, q)
  @params = params
  @timeseries = timeseries
  @p = p
  @q = q
  ll
end

Instance Attribute Details

#aic ⇒ Object (readonly)

aic Gives AIC(Akaike Information Criterion) www.scss.tcd.ie/Rozenn.Dahyot/ST7005/13AICBIC.pdf



18
19
20
 
# File 'lib/bio-statsample-timeseries/arima/likelihood.rb', line 18

def aic
  @aic
end

#log_likelihood ⇒ Object (readonly)

log_likelihood Gives log likelihood value of an ARMA(p, q) process on given parameters



9
10
11
 
# File 'lib/bio-statsample-timeseries/arima/likelihood.rb', line 9

def log_likelihood
  @log_likelihood
end

#sigma ⇒ Object (readonly)

sigma Gives sigma value of an ARMA(p,q) process on given parameters



13
14
15
 
# File 'lib/bio-statsample-timeseries/arima/likelihood.rb', line 13

def sigma
  @sigma
end

Instance Method Details

#ll ⇒ Object

Log likelihood function. iteratively minimized by simplex algorithm via KalmanFilter.ks Not meant to be used directly. Will make it private later.



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
 
# File 'lib/bio-statsample-timeseries/arima/likelihood.rb', line 31

def ll
  params, timeseries = @params, @timeseries
  p, q = @p, @q

  phi = []
  theta = []
  phi = params[0...p] if p > 0
  theta = params[(p)...(p + q)] if q > 0

  [phi, theta].each do |v|
    if v.size>0 and v.map(&:abs).inject(:+) > 1
      return
    end
  end

  m = [p, q].max
  h = Matrix.column_vector(Array.new(m,0))
  m.times do |i|
    h[i,0] = phi[i] if i< p
    h[i,0] = h[i,0] + theta[i] if i < q
  end

  t = Matrix.zero(m)
  #set_column is available in utility.rb
  t = t.set_column(0, phi)
  if(m > 1)
    t[0...(m-1), 1...m] = Matrix.I(m-1)
    #chances of extra constant 0 values as unbalanced column, so:
    t = Matrix.columns(t.column_vectors)
  end

  g = Matrix[[1]]
  a_t = Matrix.column_vector(Array.new(m,0))
  n = timeseries.size
  z = Matrix.row_vector(Array.new(m,0))
  z[0,0] = 1
  p_t = Matrix.I(m)
  v_t, f_t = Array.new(n,0), Array.new(n, 0)

  n.times do |i|
    v_t[i] = (z * a_t).map { |x| timeseries[i] - x }[0,0]

    f_t[i] = (z * p_t * (z.transpose)).map { |x| x + 1 }[0,0]

    k_t = ((t * p_t * z.transpose) + h).map { |x| x / f_t[i] }

    a_t = (t * a_t) + (k_t * v_t[i])
    l_t = t - k_t * z
    j_t = h - k_t

    p_t = (t * p_t * (l_t.transpose)) + (h * (j_t.transpose))
  end

  pot = v_t.map(&:square).zip(f_t).map { |x,y| x / y}.inject(:+)
  sigma_2 = pot.to_f / n.to_f

  f_t_log_sum = f_t.map { |x| Math.log(x) }.inject(:+)
  @log_likelihood = -0.5 * (n*Math.log(2*Math::PI) + n*Math.log(sigma_2) + f_t_log_sum + n)
  
  @sigma = sigma_2
  @aic = -(2 * @log_likelihood - 2*(p+q+1))
  #puts ("ll = #{-ll}")
  return @log_likelihood
end

#to_s ⇒ Object 



96
97
98
99
 
# File 'lib/bio-statsample-timeseries/arima/likelihood.rb', line 96

def to_s
  sprintf("LogLikelihood(p = %d, q = %d) on params: [%s]",
          @p, @q, @params.join(', '))
end