Skip to content

Instantly share code, notes, and snippets.

@IshitaTakeshi
Last active April 8, 2018 09:39
Show Gist options
  • Save IshitaTakeshi/f2a4a9b92af3ed83c2db5384821d6b66 to your computer and use it in GitHub Desktop.
Save IshitaTakeshi/f2a4a9b92af3ed83c2db5384821d6b66 to your computer and use it in GitHub Desktop.
type Lasso <: Model
w::AA{Float64, 1}
b::Float64
λ::Float64
η::Float64
Lasso(;λ = 1.0, η = 0.01) = new(zeros(Float64, 0), 0., λ, η)
end
proxₗ₁(λ, w) = sign.(w) .* max.(abs.(w) - λ, 0)
function update!{T <: Real}(model::Lasso, X::AA{Float64, 2}, y::AA{T, 1})
for i in 1:size(X, 2)
Δw, Δb = ∇L(model, view(X, :, i), y[i])
model.w = proxₗ₁(model.λ, model.w - model.η * Δw)
model.b = proxₗ₁(model.λ, model.b - model.η * Δb)
end
return model
end
function init_weights!(model::Lasso, X::AA{Float64, 2})
model.w = zeros(size(X, 1))
model.b = 0
return model
end
function fit!{T <: Real}(model::Lasso, X::AA{Float64, 2}, y::AA{T, 1};
n_iter = 200)
assert(size(X, 2) == size(y, 1))
init_weights!(model, X)
for i in 1:n_iter
model = update!(model, X, y)
end
return model
end
function ∇L{T <: Real}(model::Lasso, x::AA{Float64, 1}, y::T)
yₚ = predict(model, x)
Δ = yₚ - y
return x * Δ, Δ
end
function ∇L{T <: Real}(model::Lasso, X::AA{Float64, 2}, y::AA{T, 1})
n_features, n_samples = size(X)
Δw = zeros(n_features)
Δb = 0
for i in 1:n_samples
Δwᵢ, Δbᵢ = ∇L(model, view(X, :, i), y[i])
Δw += Δwᵢ
Δb += Δbᵢ
end
Δw = Δw / n_samples
Δb = Δb / n_samples
return Δw, Δb
end
predict(model::Lasso, x::AA{Float64, 1}) = dot(model.w, x) + model.b
predict(model::Lasso, X::AA{Float64, 2}) = vec(model.w' * X + model.b)
# data is available at https://archive.ics.uci.edu/ml/datasets/Wine+Quality
using CSV
srand(1234)
const AA = AbstractArray
include("model.jl")
include("utils.jl")
include("lasso.jl")
include("visualization.jl")
include("datasets/wine.jl")
include("datasets/make_regression.jl")
function search_param(X::AA{Float64, 2}, y::AA{Float64, 1}, λ::Float64)
min_error = Inf
argmin_η = 0
argmin_model = nothing
for η in [1e-6, 5e-6, 1e-5, 5e-5, 1e-4]
model = Lasso(λ = λ, η = η)
model = fit!(model, X, y, n_iter = 200)
E = mean_squared_error(y, predict(model, X))
if E < min_error
min_error = E
argmin_η = η
argmin_model = model
end
end
return min_error, argmin_η, argmin_model
end
X, y = make_regression(400, [4.0, 0.0], 0.8; noise = 3.0)
for λ in [0.0, 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4]
min_error, argmin_η, model = search_param(X, y, λ)
println("")
println("λ = $λ η = $argmin_η")
println("error = $min_error")
println(model.w)
println(model.b)
plot_regression_line(model, X, y)
clf()
plot_regression_surface(model, X, y)
clf()
end
abstract type Model end
using Base.Test
import Base.error
const AA = AbstractArray
include("utils.jl")
function test_mean_std_normalization()
X = Float64[
3 4 -1 3 0;
9 1 6 -8 1;
3 9 0 -4 -4;
5 -1 9 7 -1;
]
X = mean_std_normalization(X)
for x in mean(X, 1)
@test abs(x) ≈ 0.0 atol=1e-10
end
for x in std(X, 1)
@test abs(x) ≈ 1.0 atol=1e-10
end
end
function test_error()
yₜ = [1, 2, 1, 3, 4]
yₚ = [1, 3, 1, 1, 3]
@test error(yₜ, yₚ) == 1.2
end
test_mean_std_normalization()
test_error()
function meshgrid(xs1::AA{Float64, 1}, xs2::AA{Float64, 1})
return ([x₁ for x₁=xs1, x₂=xs2], [x₂ for x₁=xs1, x₂=xs2])
end
function train_test_split(X, y; train_ratio = 0.8)
assert(size(X, 2) == size(y, 1))
N = size(X, 2)
n_train = Int(round(N * train_ratio))
indices = shuffle(1:N)
train_indices = indices[1:n_train]
test_indices = indices[n_train+1:end]
X_train, X_test = X[:, train_indices], X[:, test_indices]
y_train, y_test = y[train_indices], y[test_indices]
return X_train, X_test, y_train, y_test
end
function mean_std_normalization(X::AA{Float64, 2})
X = X .- mean(X, 1)
X = X ./ std(X, 1)
return X
end
function mean_squared_error{T, U <: Real}(yₜ::AA{T, 1}, yₚ::AA{U, 1})
assert(length(yₜ) == length(yₚ))
N = length(yₜ)
dot(yₜ - yₚ, yₜ - yₚ) / N
end
using PyPlot
function plot_regression_line(model::Model,
X::AA{Float64, 2}, y::AA{Float64, 1})
N = 201
n_features = size(X, 1)
for i in 1:n_features
subplot(n_features, 1, i)
u = view(X, i, :)
scatter(u, y)
# draw the regression line
Xₚ = zeros(n_features, N)
v = linspace(minimum(u), maximum(u), N)
Xₚ[i, :] = v
plot(v, predict(model, Xₚ))
end
show()
end
function plot_regression_surface(model::Model, X::AA{Float64, 2}, y::AA{Float64, 1})
assert(size(X, 1) == 2)
x1, x2 = view(X, 1, :), view(X, 2, :)
scatter3D(x1, x2, y)
N = 201
L₁, H₁ = minimum(x1), maximum(x1)
L₂, H₂ = minimum(x2), maximum(x2)
xs1 = linspace(L₁, H₁, N)
xs2 = linspace(L₂, H₂, N)
X₁, X₂ = meshgrid(xs1, xs2)
Y = [predict(model, [x₁, x₂]) for x₁=xs1, x₂=xs2]
plot_surface(X₁, X₂, Y)
show()
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment