In [1]:
function setrand()
s = rand(1:1000)
println("s = ", s)
srand(s)
end
function loaddata(n)
d = 2 # dimension of u space
s = 1.4 # mean separation factor
fracs = [0.5, 0.5] # fraction in the different classes
K = length(fracs) # number of classes
# square roots of covariances
covs = Any[]
for i=1:K
A = randn(2,2)
D,V = eig(A*A')
D2 = [randn()+0.7 0; 0 randn()+0.7]
push!(covs, V*D2*V')
end
# assign categories according to fractions
cfracs = cumsum(fracs)
cats = zeros(Int64,n)
c = 1
for i=1:n
if i/n > cfracs[c]
c += 1
end
cats[i] = c
end
# choose the means of each cluster
means = zeros(d,K)
for i=1:K
means[:,i] = s*randn(d)
end
# assign U and V
U = zeros(n,d)
V = zeros(n)
for i=1:n
c = cats[i]
V[i] = c
x = (means[:,c] + covs[c]*randn(d))'
U[i,:] = x
end
Unew, Vnew = normalizedata(U,V)
return Unew, Vnew, K
end
function normalizedata(U,V)
n = size(U,1)
xmin = minimum(U[:,1]) - 0.1
xmax = maximum(U[:,1]) + 0.1
ymin = minimum(U[:,2]) - 0.1
ymax = maximum(U[:,2]) + 0.1
m = 0.5*[ xmax+xmin ymax+ymin]
Unew = U - repmat(m, n, 1)
Unew[:,1] = Unew[:,1] * 8/(xmax-xmin)
Unew[:,2] = Unew[:,2] * 8/(ymax-ymin)
return Unew, V
end
Out[1]:
In [2]:
using PyPlot, Convex, ECOS
srand(323)
#srand(999)
U,V,K = loaddata(100)
n = length(V)
X = [ones(n) U]
y = V
y[V.==2] = -1
U_p = U[y.==1,:]
U_n = U[y.==-1,:]
## least squares classifier
theta = Variable(3);
obj = sumsquares(1-y.*(X*theta)) + 0.1*sumsquares(theta[2:end])
#obj = sumsquares(X*theta-y) + 0.1*sumsquares(theta[2:end])
problem = minimize(obj)
solve!(problem, ECOSSolver())
theta_opt = evaluate(theta)
# separating plane
s1 = -5:0.1:5
s2 = - ( theta_opt[1] + s1*theta_opt[2] ) / theta_opt[3]
s2p = - ( 1 + theta_opt[1] + s1*theta_opt[2] ) / theta_opt[3]
s2n = - ( -1 + theta_opt[1] + s1*theta_opt[2] ) / theta_opt[3]
s2_ls = s2
figure()
plot(U_p[:,1], U_p[:,2], "o", alpha=0.5)
plot(U_n[:,1], U_n[:,2], "v", alpha=0.5)
plot(s1,s2)
plot(s1,s2p, ":")
plot(s1,s2n, ":")
grid("on")
axis("square")
xlim(-5, 5)
ylim(-5, 5)
title("Least squares classifier")
## logistic regression
theta = Variable(3);
obj = logisticloss(-y.*(X*theta)) + 0.1*sumsquares(theta[2:end])
problem = minimize(obj)
solve!(problem, ECOSSolver())
theta_opt = evaluate(theta)
# separating plane
s1 = -5:0.1:5
s2 = - ( theta_opt[1] + s1*theta_opt[2] ) / theta_opt[3]
s2_lr = s2
figure()
plot(U_p[:,1], U_p[:,2], "o", alpha=0.5)
plot(U_n[:,1], U_n[:,2], "v", alpha=0.5)
plot(s1,s2)
grid("on")
axis("square")
xlim(-5, 5)
ylim(-5, 5)
title("Logistic regression")
## support vector machine
theta = Variable(3);
obj = sum(max(0,1-y.*(X*theta))) + 0.1*sumsquares(theta[2:end])
problem = minimize(obj)
solve!(problem, ECOSSolver())
theta_opt = evaluate(theta)
# separating plane
s1 = -5:0.1:5
s2 = - ( theta_opt[1] + s1*theta_opt[2] ) / theta_opt[3]
s2p = - ( 1 + theta_opt[1] + s1*theta_opt[2] ) / theta_opt[3]
s2n = - ( -1 + theta_opt[1] + s1*theta_opt[2] ) / theta_opt[3]
s2_svm = s2
figure()
plot(U_p[:,1], U_p[:,2], "o", alpha=0.5, label="Positive data")
plot(U_n[:,1], U_n[:,2], "v", alpha=0.5, label="Negative data")
plot(s1,s2, label=L"$\theta^Tx=0$")
plot(s1,s2p, ":", label=L"$\theta^Tx=-1$")
plot(s1,s2n, ":", label=L"$\theta^Tx=1$")
grid("on")
axis("square")
legend()
xlim(-5, 5)
ylim(-5, 5)
title("Support vector machine")
figure()
plot(U_p[:,1], U_p[:,2], "o", alpha=0.5, label="Positive data")
plot(U_n[:,1], U_n[:,2], "v", alpha=0.5, label="Negative data")
plot(s1,s2_ls, label="Least squares classifier")
plot(s1,s2_lr, label="Logistic regression")
plot(s1,s2_svm, label="Support vector machine")
grid("on")
axis("square")
legend()
xlim(-5, 5)
ylim(-5, 5)
Out[2]: