Regression in Julia

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

######################################
## hipótese verdadeira não-linear 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))

######################################
## hipótese verdadeira não-linear 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))

######################################
## hipótese verdadeira linear
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.point, Theme(default_color=color("red"))),
    layer(x=x, y=y, Geom.point),
    layer(x=x, y=ysample, Geom.point, Theme(default_color=color("black"))))
end

######################################
# para hipotese linear
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()

######################################
# para hipotese não-linear 1
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()