class Statsample::GLM::MLE::Normal

Normal Distribution MLE estimation. Usage: TODO : Document this properly

mle=Statsample::MLE::Normal.new
mle.newton_raphson(x,y)
beta=mle.coefficients
likelihood=mle.likelihood(x,y,beta)
iterations=mle.iterations

Protected Instance Methods

_log_likelihood(x,y,b) click to toggle source

Total MLE for given X, Y and B matrices (overridden over the Base version)

# File lib/statsample-glm/glm/mle/normal.rb, line 21
def _log_likelihood(x,y,b)
  n      = x.row_size.to_f
  sigma2 = b[b.row_size-1,0]
  betas  = Matrix.columns([b.column(0). to_a[0...b.row_size-1]])
  e      = y - (x * betas)
  last   = (1 / (2*sigma2)) * e.t * e
  (-(n / 2.0) * Math::log(2*Math::PI))-((n / 2.0)*Math::log(sigma2)) - last[0,0]
end
first_derivative(x,y,p) click to toggle source

First derivative for Normal Model. p should be [k+1,1], because the last parameter is sigma^2

# File lib/statsample-glm/glm/mle/normal.rb, line 31
def first_derivative(x,y,p)
  raise "x.rows != y.rows" if x.row_size != y.row_size
  raise "x.columns + 1 != p.rows" if x.column_size + 1 != p.row_size
  
  n = x.row_size
  k = x.column_size
  b = Array.new(k)

  k.times{|i| b[i]=[p[i,0]]}
  beta   = Matrix.rows(b)
  sigma2 = p[k,0]
  sigma4 = sigma2 * sigma2
  e = y-(x * (beta))
  xte = x.transpose*(e)
  ete = e.transpose*(e)
  #rows of the Jacobian
  rows = Array.new(k+1)
  k.times{|i| rows[i] = [xte[i,0] / sigma2]}
  rows[k] = [ete[0,0] / (2*sigma4) - n / (2*sigma2)]
  Matrix.rows(rows, true)
end
measurement(data_set, coefficients) click to toggle source
# File lib/statsample-glm/glm/mle/normal.rb, line 17
def measurement data_set, coefficients
  (data_set * coefficients[0..-2,0]).map { |xb| xb }
end
second_derivative(x,y,p) click to toggle source

second derivative for normal model p should be [k+1,1], because the last parameter is sigma^2

# File lib/statsample-glm/glm/mle/normal.rb, line 55
def second_derivative(x,y,p)
  raise "x.rows != y.rows"        if x.row_size != y.row_size
  raise "x.columns + 1 != p.rows" if x.column_size + 1 != p.row_size
  #n = x.row_size
  k = x.column_size
  b = Array.new(k)
  k.times{|i| b[i] = [p[i,0]]}
  beta = Matrix.rows(b)
  sigma2 = p[k,0]
  sigma4 = sigma2*sigma2
  sigma6 = sigma2*sigma2*sigma2
  e = y-(x*(beta))
  xtx = x.transpose*(x)
  xte = x.transpose*(e)
  ete = e.transpose*(e)
  #rows of the Hessian
  rows = Array.new(k+1)

  k.times do |i|
    row = Array.new(k+1)
    k.times do |j|
        row[j] = -xtx[i,j] / sigma2
    end
    row[k] = -xte[i,0] / sigma4
    rows[i] = row
  end

  last_row = Array.new(k+1)
  k.times do |i|
    last_row[i] = -xte[i,0] / sigma4
  end

  last_row[k] = 2*sigma4 - ete[0,0] / sigma6
  rows[k] = last_row
  Matrix.rows(rows, true)
end