-
Notifications
You must be signed in to change notification settings - Fork 5
/
demo_gendirmnreg.m
81 lines (71 loc) · 2.03 KB
/
demo_gendirmnreg.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
%% Generalized Dirichlet-Multinomial regression and sparse regression
% A demo of gen. Dirichlet-Multinomial regression and sparse regression
%% Generate generalized Dirichlet-Multinomial random vectors from covariates
clear;
% reset random seed
s = RandStream('mt19937ar','Seed',1);
RandStream.setGlobalStream(s);
% sample size
n = 500;
% # covariates
p = 15;
% # bins
d = 5;
% design matrix
X = randn(n,p);
% true regression coefficients
A = zeros(p,d-1);
B = zeros(p,d-1);
nzidx = [1 3 5];
A(nzidx,:) = 0.5.*ones(length(nzidx),d-1);
B(nzidx,:) = 0.5.*ones(length(nzidx),d-1);
alpha = exp(X*A);
beta = exp(X*B);
batchsize = 25+unidrnd(25,n,1);
Y = gendirmnrnd(batchsize,alpha, beta);
%% Fit generalized Dirichlet-Multinomial regression
tic;
[Bhat1,Bhat2,stats_gdm] = gendirmnreg(X,Y);
toc;
display(Bhat1);
display(Bhat2);
display(stats_gdm);
display(stats_gdm.se);
display(stats_gdm.wald_pvalue);
%% Fit generalized Dirichlet-Multinomial sparse regression - lasso/group/nuclear penalty
penalty = {'sweep','group','nuclear'};
ngridpt = 10;
dist = 'gendirmn';
for i = 1:length(penalty)
pen = penalty{i};
[~, stats] = mglm_sparsereg(X,Y,inf,'penalty',pen,'dist',dist);
maxlambda = stats.maxlambda;
lambdas = exp(linspace(log(maxlambda),log(maxlambda/100),ngridpt));
BICs = zeros(1,ngridpt);
tic;
for j=1:ngridpt
if j==1
B0 = zeros(p,2*(d-1));
else
B0 = B_hat;
end
[B_hat, stats] = mglm_sparsereg(X,Y,lambdas(j),'penalty',pen, ...
'dist',dist,'B0',B0);
BICs(j) = stats.BIC;
end
toc;
% True signal versus estimated signal
[bestbic,bestidx] = min(BICs);
lambdas(bestidx)
B_best = mglm_sparsereg(X,Y,lambdas(bestidx),'penalty',pen,'dist',dist);
figure;
subplot(1,3,1);
semilogx(lambdas,BICs);
ylabel('BIC');
xlabel('\lambda');
xlim([min(lambdas) max(lambdas)]);
subplot(1,3,2);
imshow(mat2gray(-[A,B])); title('True B');
subplot(1,3,3);
imshow(mat2gray(-B_best)); title([pen ' estimate']);
end