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