graph/neural_network_linear_regression.R at main · patrickrd/graph · GitHub

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
###Caron's model
## f = function(x,y){1-exp(-2*x * y)}
## n=5000 ##because that way the readers have been the quality of the manifold approximation
## Z = rgamma(n,1,1)
## P = matrix(nrow=n,ncol=n)
## for (i in 1:n){
##         P[i,] = f(Z[i],Z)
## }
## d=100
## E = eigen(P)
## keep=sort(order(abs(E$values),decreasing=TRUE)[1:d])
## X = E$vectors[,keep] %*% diag(sqrt(abs(E$values[keep])))
## ##plot(X[,1:2], cex=.2)
## A=matrix(runif(n*n)<= P,nrow=n)
## ##A=matrix(rnorm(n*n,sd=0.05), nrow=n)+P
## A = sym(A)
## res=getXQ(A, X, sum(E$values[keep]>0), sum(E$values[keep]<0))
## XhatQ = as.matrix(res$Xhat %*% t(res$Q))
## save(Z,res,file="embedding_neural_network_data.Rout")


alpha = 2; beta = 5; std = 1
y = alpha + beta*Z + rnorm(n,sd=std)

ntrain = 100

library(tensorflow)
##install_tensorflow()
##tf$constant("Hellow Tensorflow")
##install.packages("keras")
library(keras)

model <- keras_model_sequential()
model %>%
    layer_dense(units = 256, activation = 'relu', input_shape = c(d)) %>%
    layer_dropout(rate = 0.4) %>%
    layer_dense(units = 128, activation = 'relu') %>%
    layer_dropout(rate = 0.3) %>%
    layer_dense(units = 1)


# Configure a model for mean-squared error regression.
model %>% compile(
  optimizer = 'adam',
  loss = 'mse',           # mean squared error
  metrics = list('mae')   # mean absolute error
)

model %>% summary()

history <- model %>% fit(
  XhatQ[1:ntrain,], y[1:ntrain],
  epochs = 30, batch_size = 128,
  validation_split = 0.2
  )


##points(Z, yhat, col="green")

##(model %>% evaluate(XhatQ, y, verbose = 0))

df = as.data.frame(cbind(XhatQ[1:ntrain,],rep(1,ntrain)))
modellm <- lm (y[1:ntrain] ~ ., df)

library(glmnet)
##fit.lasso <- glmnet(cbind(XhatQ,rep(1,n)), y, family="gaussian", lambda=1)
fit.lasso <- cv.glmnet(cbind(XhatQ[1:ntrain,],rep(1,ntrain)), y[1:ntrain], family="gaussian")


##yhatlassomin = predict(fit.lasso, cbind(XhatQ,rep(1,n)), lambda=fit.lasso$lambda.min)
##predict(fit.lasso, cbind(XhatQ,rep(1,n)))
## library(glmnet)
## fit.lasso <- cv.glmnet(cbind(XhatQ[1:ntrain],rep(1,ntrain)), y[1:ntrain], family="gaussian")
## fit.lasso <- cv.glmnet(cbind(XhatQ[1:n],rep(1,n)), y[1:n], family="gaussian")
## yhatlasso = predict(fit.lasso, cbind(XhatQ,rep(1,n)), lambda=fit.lasso$lambda.1se)
## ##yhatlassomin = predict(fit.lasso, cbind(XhatQ,rep(1,n)), lambda=fit.lasso$lambda.min)

## plot(Z,yhatlasso)
## plot(Z,yhatlassomin)
## mean((yhatlassomin - y)^2)

library(randomForest)
fit.rf = randomForest(cbind(XhatQ[1:ntrain,],rep(1,ntrain)), y[1:ntrain])
yhatrf = predict(fit.rf, cbind(XhatQ,rep(1,n)))
plot(Z, yhatrf, col="red")
mean((yhatrf - y)^2)


###the comparison
test = (ntrain+1):n

yhatnn=model %>% predict(XhatQ[test,])
yhatnna=model %>% predict(XhatQ)

dft = as.data.frame(cbind(XhatQ[test,],rep(1,length(test))))
dfta = as.data.frame(cbind(XhatQ,rep(1,n)))
yhatlm = predict(modellm, dft)
yhatlma = predict(modellm, dfta)

yhatlasso = predict(fit.lasso, cbind(XhatQ[test,],rep(1,length(test))), lambda=fit.lasso$lambda.1se)
yhatlassoa = predict(fit.lasso, cbind(XhatQ,rep(1,n)), lambda=fit.lasso$lambda.1se)

yhatrf = predict(fit.rf, cbind(XhatQ[test,],rep(1,length(test))))
yhatrfa = predict(fit.rf, cbind(XhatQ,rep(1,n)))

pdf("embedding_regression2.pdf", width=6, height=5)
par(mar=c(4.5, 4.5, 0.5, 0.5))
plot(Z, y, cex=.1, pch=16, xlab="Z", ylab="Y")
o = order(Z)
lines(Z[o], yhatnna[o], cex=.1, col="red")
lines(Z[o], yhatlma[o], cex=.1, col="purple")
lines(Z[o], yhatlassoa[o], cex=.1, col="orange")
lines(Z[o], yhatrfa[o], cex=.1, col="green")
lines(sort(Z),alpha + beta*sort(Z))
legend("bottomright",legend=c(paste0("Ideal regression (", round(mean((y[test]-(alpha+beta*Z[test]))^2), digit=2),")"), paste0("Random forest (",round(mean((yhatrf-y[test])^2), digit=2),")"), paste0("Neural network (",round(mean((yhatnn-y[test])^2), digit=2),")"), paste0("Lasso (",round(mean((yhatlasso-y[test])^2), digit=2),")"), paste0("Least squares (",round(mean((yhatlm-y[test])^2), digit=2),")")), col=c("black","green", "red","orange","purple"), lty=1)
dev.off()


####


##A full dimensional regression
beta_large = rnorm(d);
yfull = alpha + -X %*% beta_large + rnorm(n,sd=std)

history <- model %>% fit(
  XhatQ[1:ntrain,], yfull[1:ntrain],
  epochs = 30, batch_size = 128,
  validation_split = 0.2
  )

yhatnn=model %>% predict(XhatQ[test,])
mean((yhatnn-yfull[test])^2)

fit.rf = randomForest(cbind(XhatQ[1:ntrain,],rep(1,ntrain)), yfull[1:ntrain])
yhatrf = predict(fit.rf, cbind(XhatQ[test,],rep(1,length(test))))
mean((yhatrf-yfull[test])^2)