Linear Regression implementation with Gradient descent and linear and polynomial hypotheses

linear_regression.jl

#     - Linear Regression implementation with Gradient descent and linear and polynomial hypotheses (in Julia)
#     Copyright (C) 2015  Eric Aislan Antonelo

#     This program is free software: you can redistribute it and/or modify
#     it under the terms of the GNU General Public License as published by
#     the Free Software Foundation, either version 3 of the License, or
#     (at your option) any later version.

######################################
using Gadfly
using Distributions

######################################
## true nonlinear hypothesis 1
n = 3
x = 1:0.05:n
y =  sqrt(30x ) + 5sin(4x)

d = Normal(0,1.1)
ysample = y + rand(d,length(x))

plot(layer(x=x,y=y,Geom.line),
    layer(x=x,y=ysample, Geom.point))

######################################
## true nonlinear hypothesis 2
n = 9
x = 1:0.1:n
y =  0.1x + 0.4(x.^2) - 0.5(x.^3)

d = Normal(0,29.5)
ysample = y + rand(d,length(x))

plot(layer(x=x,y=y,Geom.line),
    layer(x=x,y=ysample, Geom.point))

######################################
## true linear hypothesis
n = 10
x = 1:0.1:n
d = Normal(0,8)
y = 5.4x
ysample = y + rand(d,length(x))
plot(x=x,y=ysample,Geom.point)

######################################
function cost(examples, theta)
    err = 0
    for (x,y) in examples
        err += (y - hypothesis(x, theta))^2
    end
    err/2length(examples)
end

function gradient_descent(examples, theta; learning_rate=0.05)
    err = 10000
#     learning_rate = 0.05
    max_episodes = 25
    i = 0
    while err > 0.5 && i < max_episodes
        println("Episode $i")
        theta = theta - learning_rate * gradient(examples, theta)
        println("theta: $theta")
        err = cost(examples, theta)
        println("cost: $err")
        i += 1
    end
    theta
end

function gradient(examples, theta)
    grad = zeros(size(theta))
    for (x,y) in examples
        #println(x," __ ",y)
        #println(hypothesis(x, theta)[1] -y)
        #println((hypothesis(x, theta) - y) * [1 x]')
        grad = grad + (hypothesis(x, theta) - y) * hypothesis_features(x)
    end
    grad/length(examples)
end

function hypothesis(x, theta)
    y =  hypothesis_features(x)' * theta
    y[1]
end

# Plota resultados
function plot_results()
    h_ = zeros(length(x))
    for i in 1:length(x)
        h_[i] = hypothesis(x[i], theta)
    end
    # hypothesis in red
    plot(layer(x=x, y=h_, Geom.line, Theme(default_color=color("red"))),
    layer(x=x, y=y, Geom.line),
    layer(x=x, y=ysample, Geom.point, Theme(default_color=color("black"))))
end


######################################
# linear hypothesis
function hypothesis_features(x)
    [1; x;]
end
theta_ini = ones(2) * 0.001
theta = gradient_descent(zip(x,ysample), theta_ini, learning_rate=0.01)
println("Theta $theta")
plot_results()

######################################
# # nonlinear hypothesis
# function hypothesis_features(x)
#     [1; x; sqrt(x); 1/(1+exp(-x)); x^2; x^3]
# end
# theta_ini = ones(6) * 0.00001
# theta = gradient_descent(zip(x,ysample), theta_ini, learning_rate=0.000001)
# println("Theta $theta")
# plot_results()