-
Notifications
You must be signed in to change notification settings - Fork 0
/
GLMNET.Rmd
108 lines (88 loc) · 2.45 KB
/
GLMNET.Rmd
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
---
title: "SCADD"
author: "Bo Zhang"
date: "5/24/2019"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r}
library(picasso)
library(MASS)
library(Matrix)
library(glmnet)
ds = c(256,512,1024)
ns = c(100,200)
l2s = c()
FPs = c()
FNs = c()
for (i in 1:length(ns)) {
n = ns[i]
print(n)
for (j in 1:length(ds)) {
d = ds[j]
print(d)
true_theta = c(2,2,2,-1.5,-1.5,-1.5,2,2,2,2, rep(0, d-10))
l2 = 0
FP = 0
FN = 0
for (i in 1:200){
# generate data using rnorm
X = matrix(rnorm(n*d), nrow = n)
# generate labels with errors
Y = X %*% true_theta + rnorm(n, 0, 1.5)
fold = 5
result = glmnet(X, Y, lambda.min.ratio = 0.001, alpha = 1)
lambdas = result$lambda
mean_val_error = c()
for (j in 1:length(lambdas)){
lambda = lambdas[j]
errors = c()
for (i in 1:fold){
start = as.integer((i-1)/fold*n)
end = as.integer(i/fold*n)
X_vali = X[(start+1):end,]
Y_vali = Y[(start+1):end,]
#print(start)
#print(end)
if (start == 0){
X_train = X[(end+1):n,]
Y_train = Y[(end+1):n,]
}
else if (end == n){
X_train = X[1:start,]
Y_train = Y[1:start,]
}
else{
X_train = rbind(X[1:start,], X[(end+1):n,])
Y_train = c(Y[1:start,], Y[(end+1):n,])
}
#print(dim(X_train))
#print(length(Y_train))
est_beta = glmnet(X, Y, lambda = lambda, alpha = 1)$beta
est_Y = X_vali %*% est_beta
# mean sum of squared error
errors[i] = sum((est_Y - Y_vali) ** 2) / length(Y_vali)
}
# mean cross validation errors
mean_val_error[j] = mean(errors)
}
min_index = which.min(mean_val_error)
lambda = lambdas[min_index]
## Useing the lambda after cross validation to predict theta
result_optimal = glmnet(X, Y, lambda = lambda, alpha = 1)
pred_theta = result_optimal$beta
l2 = l2 + norm(as.matrix(pred_theta - true_theta), type = '2')
FP = FP + sum((true_theta == 0) & (pred_theta != 0))
FN = FN + sum((true_theta != 0) & (pred_theta == 0))
}
l2s = c(l2s, l2/200)
FPs = c(FPs, FP/200)
FNs = c(FNs, FN/200)
}
}
```
```{r}
glmnet.lasso = c(l2s, FPs, FNs)
```