-
Notifications
You must be signed in to change notification settings - Fork 11
/
demo_lsq.m
209 lines (189 loc) · 6.91 KB
/
demo_lsq.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
%% Sparse Linear Regression
% A demonstration of sparse linear regression using SparseReg toolbox.
% Sparsity is in the general sense: variable selection, total variation
% regularization, polynomial trend filtering, and others. Various penalties
% are implemented: elestic net (enet), power family (bridge regression),
% log penalty, SCAD, and MCP.
%% Sparse linear regression (n>p)
% Simulate a sample data set (n=500, p=100)
clear;
s = RandStream('mt19937ar','Seed',1);
RandStream.setGlobalStream(s);
n = 500;
p = 100;
X = randn(n,p); % generate a random design matrix
X = bsxfun(@rdivide, X, sqrt(sum(X.^2,1))); % normalize predictors
X = [ones(size(X,1),1) X]; % add intercept
b = zeros(p+1,1); % true signal:
b(2:6) = 3; % first 5 predictors are 3
b(7:11) = -3; % next 5 predictors are -3
y = X*b+randn(n,1); % response vector
%%
% Sparse regression at a fixed tuning parameter value
penalty = 'enet'; % set penalty function to lasso
penparam = 1;
penidx = ... % leave intercept unpenalized
[false; true(size(X,2)-1,1)];
lambdastart = ... % find the maximum tuning parameter to start
max(lsq_maxlambda(sum(X(:,penidx).^2),-y'*X(:,penidx),penalty,penparam));
display(lambdastart);
lambda = 0.9*lambdastart; % tuning parameter value
betahat = ... % lasso regression
lsq_sparsereg(X,y,lambda,'penalty',penalty,'penparam',penparam,'penidx',penidx);
figure; % plot penalized estimate
bar(0:length(betahat)-1,betahat);
xlabel('j');
ylabel('\beta_j');
xlim([-1,length(betahat)]);
title([penalty '(' num2str(penparam) '), \lambda=' num2str(lambda,2)]);
lambda = 0.5*lambdastart; % try a smaller tuning parameter value
betahat = ... % sparse regression
lsq_sparsereg(X,y,lambda,'penalty',penalty,'penparam',penparam,'penidx',penidx);
figure; % plot penalized estimate
bar(0:length(betahat)-1,betahat);
xlabel('j');
ylabel('\beta_j');
xlim([-1,length(betahat)]);
title([penalty '(' num2str(penparam) '), \lambda=' num2str(lambda,2)]);
%%
% Solution path for lasso
penalty = 'enet'; % set penalty function
penparam = 1;
penidx = [false; true(size(X,2)-1,1)]; % leave intercept unpenalized
tic;
[rho_path,beta_path] = ... % compute solution path
lsq_sparsepath(X,y,'penalty',penalty,'penparam',penparam,'penidx',penidx);
timing = toc;
figure;
plot(rho_path,beta_path);
xlabel('\rho');
ylabel('\beta(\rho)');
xlim([min(rho_path),max(rho_path)]);
title([penalty '(' num2str(penparam) '), ' num2str(timing,2) ' sec']);
%%
% Solution path for power (0.5)
penalty = 'power'; % set penalty function to power
penparam = 0.5;
tic;
[rho_path,beta_path] = ...
lsq_sparsepath(X,y,'penalty',penalty,'penparam',penparam,'penidx',penidx);
timing = toc;
figure;
plot(rho_path,beta_path);
xlabel('\rho');
ylabel('\beta(\rho)');
xlim([min(rho_path),max(rho_path)]);
title([penalty '(' num2str(penparam) '), ' num2str(timing,2) ' sec']);
%%
% Compare solution paths from different penalties
penalty = {'enet' 'enet' 'enet' 'power' 'power' 'log' 'log' 'log' 'scad'};
penparam = [1 1.5 2 0.5 1 0 1 5 3.7];
penidx = [false; true(size(X,2)-1,1)]; % leave intercept unpenalized
figure
for i=1:length(penalty)
tic;
[rho_path,beta_path] = lsq_sparsepath(X,y,...
'penalty',penalty{i},'penparam',penparam(i),'penidx',penidx);
timing = toc;
subplot(3,3,i);
plot(rho_path,beta_path);
if (i==8)
xlabel('\rho');
end
if (i==4)
ylabel('\beta(\rho)');
end
xlim([min(rho_path),max(rho_path)]);
title([penalty{i} '(' num2str(penparam(i)) '), ' num2str(timing,1) 's']);
end
%% Fused linear regression
% Fused lasso (fusing the first 10 predictors)
D = zeros(9,size(X,2)); % regularization matrix for fusing first 10 preds
D(10:10:90) = 1;
D(19:10:99) = -1;
display(D(1:9,1:11));
penalty = 'enet'; % set penalty function to lasso
penparam = 1;
tic;
[rho_path, beta_path] = lsq_regpath(X,y,D,'penalty',penalty,'penparam',penparam);
timing = toc;
figure;
plot(rho_path,beta_path(2:11,:));
xlabel('\rho');
ylabel('\beta(\rho)');
xlim([min(rho_path),max(rho_path)]);
title([penalty '(' num2str(penparam) '), ' num2str(timing,2) ' sec']);
%%
% Same fusion problem, but with power, log, MCP, and SCAD penalty
penalty = {'enet' 'power' 'log' 'mcp'};
penparam = [1.5 0.5 1 1];
for i=1:length(penalty)
tic;
[rho_path, beta_path] = lsq_regpath(X,y,D,'penalty',penalty{i},...
'penparam',penparam(i));
timing = toc;
subplot(2,2,i);
plot(rho_path,beta_path(2:11,:));
xlim([min(rho_path),max(rho_path)]);
title([penalty{i} '(' num2str(penparam(i)) '), ' num2str(timing,1) 's']);
end
%% Sparse linear regression (n<p)
% Simulate another sample data set (n=100, p=1000)
clear;
n = 100;
p = 1000;
X = randn(n,p); % generate a random design matrix
X = bsxfun(@rdivide, X, sqrt(sum(X.^2,1))); % normalize predictors
X = [ones(size(X,1),1) X]; % add intercept
b = zeros(p+1,1); % true signal
b(2:6) = 3; % first 5 predictors are 3
b(7:11) = -3; % next 5 predictors are -3
y = X*b+randn(n,1); % response vector
%%
% Solution path for lasso
maxpreds = 51; % run solution path until 50 predictors are in
penalty = 'enet'; % set penalty function
penparam = 1;
penidx = [false; true(size(X,2)-1,1)]; % leave intercept unpenalized
tic;
[rho_path,beta_path,eb_path] = lsq_sparsepath(X,y,'penidx',penidx, ...
'maxpreds',maxpreds,'penalty',penalty,'penparam',penparam);
timing = toc;
[~,ebidx] = min(eb_path);
figure;
plot(rho_path,eb_path);
xlabel('\rho');
ylabel('EBC');
xlim([min(rho_path),max(rho_path)]);
title([penalty '(' num2str(penparam) '), ' num2str(timing,2) ' sec']);
line([rho_path(ebidx), rho_path(ebidx)], ylim);
figure;
plot(rho_path,beta_path);
xlabel('\rho');
ylabel('\beta(\rho)');
xlim([min(rho_path),max(rho_path)]);
title([penalty '(' num2str(penparam) '), ' num2str(timing,2) ' sec']);
line([rho_path(ebidx), rho_path(ebidx)], ylim);
%%
% Solution path for power (0.5)
penalty = 'power'; % set penalty function to power
penparam = 0.5;
tic;
[rho_path,beta_path,eb_path] = lsq_sparsepath(X,y,'penidx',penidx, ...
'maxpreds',maxpreds,'penalty',penalty,'penparam',penparam);
timing = toc;
[~,ebidx] = min(eb_path);
figure;
plot(rho_path,eb_path);
xlabel('\rho');
ylabel('EBC');
xlim([min(rho_path),max(rho_path)]);
title([penalty '(' num2str(penparam) '), ' num2str(timing,2) ' sec']);
line([rho_path(ebidx), rho_path(ebidx)], ylim);
figure;
plot(rho_path,beta_path);
xlabel('\rho');
ylabel('\beta(\rho)');
xlim([min(rho_path),max(rho_path)]);
title([penalty '(' num2str(penparam) '), ' num2str(timing,2) ' sec']);
line([rho_path(ebidx), rho_path(ebidx)], ylim);