The main module for accessing objects is regreg.api
.
In [1]:
import numpy as np
import regreg.api as rr
from IPython.display import Math
We will compare our fitted values to R
s output from lars
.
In [2]:
%load_ext rmagic
In [3]:
%%R -o X,Y
library(lars)
data(diabetes)
X = diabetes$x
Y = diabetes$y
For a given $X, Y$, the squared error loss
In [15]:
loss = rr.squared_error(X, Y)
Math(loss.latexify())
Out[15]:
In [16]:
penalty_lagrange = rr.l1norm(10, lagrange=200)
Math(penalty_lagrange.latexify())
problem_lagrange = rr.simple_problem(loss, penalty_lagrange)
coef_lagrange = problem_lagrange.solve(tol=1.e-12)
print(coef_lagrange)
In [17]:
Math(problem_lagrange.latexify())
Out[17]:
In [18]:
implied_bound = np.fabs(coef_lagrange).sum()
print(implied_bound)
In [19]:
bound_constraint = rr.l1norm(10, bound=implied_bound)
problem_bound = rr.simple_problem(loss, bound_constraint)
Math(problem_bound.latexify())
Out[19]:
In [12]:
coef_bound = problem_bound.solve(tol=1.e-12)
print(coef_bound)
In [13]:
np.linalg.norm(coef_bound - coef_lagrange) / np.linalg.norm(coef_lagrange)
Out[13]:
In [14]:
Math(problem_bound.latexify())
Out[14]:
In [9]:
%%R -o L,B
lars_fit = lars(X, Y, type='lasso', normalize=FALSE, intercept=FALSE)
L = lars_fit$lambda
B = lars_fit$beta
In [10]:
solns = []
for lagrange in L:
penalty_lagrange.lagrange = lagrange
soln = problem_lagrange.solve(tol=1.e-12)
solns.append(soln.copy())
solns.append(np.dot(np.linalg.pinv(X), Y))
solns = np.array(solns)
In [11]:
print [np.linalg.norm(B[i] - solns[i]) / max(np.linalg.norm(B[i]), np.linalg.norm(solns[i]),1) for i in range(solns.shape[0])]
In [12]:
%%R -o B2
library(glmnet)
A=glmnet(X,Y, lambda=L/nrow(X),standardize=FALSE, intercept=FALSE)
B2 = as.matrix(A$beta)
print(L)
print(A$beta-B[1:12])
In [ ]:
print B2[-1]
print B[-1]
print solns[-2]
In [ ]:
print B2
In [ ]:
beta = np.array(soln['beta'].todense())
soln['lagrange'].shape
beta[1].shape
f1, ax1 = plt.subplots()
[ax1.plot(soln['lagrange'],beta[i]) for i in range(1,11)]
l1norm = np.fabs(beta[1:]).sum(0)
f2, ax2= plt.subplots()
[ax2.plot(l1norm,beta[i]) for i in range(1,11)];
In [ ]:
%%R
l = lars(X, Y)
plot(l)
#library(glmnet)
#g = glmnet(X, Y)
In [ ]: