# This code will reproduce the analysis, takes some time
set.seed(406406406)
n <- 200
df <- tibble(
x = seq(.05, 1, length = n),
y = sin(1 / x) + rnorm(n, 0, .1) # Doppler function
)
testx <- matrix(seq(.05, 1, length.out = 1e3), ncol = 1)
library(neuralnet)
library(splines)
fstar <- sin(1 / testx)
spline_test_err <- function(k) {
fit <- lm(y ~ bs(x, df = k), data = df)
yhat <- predict(fit, newdata = tibble(x = testx))
mean((yhat - fstar)^2)
}
Ks <- 1:15 * 10
SplineErr <- map_dbl(Ks, ~ spline_test_err(.x))
Jgrid <- c(5, 10, 15)
NNerr <- double(length(Jgrid)^3)
NNplot <- character(length(Jgrid)^3)
sweep <- 0
for (J1 in Jgrid) {
for (J2 in Jgrid) {
for (J3 in Jgrid) {
sweep <- sweep + 1
NNplot[sweep] <- paste(J1, J2, J3, sep = " ")
nn_out <- neuralnet(y ~ x, df,
hidden = c(J1, J2, J3),
threshold = 0.01, rep = 3
)
nn_results <- sapply(1:3, function(x) {
compute(nn_out, testx, x)$net.result
})
# Run them through the neural network
Yhat <- rowMeans(nn_results)
NNerr[sweep] <- mean((Yhat - fstar)^2)
}
}
}
bestK <- Ks[which.min(SplineErr)]
bestspline <- predict(lm(y ~ bs(x, bestK), data = df), newdata = tibble(x = testx))
besthidden <- as.numeric(unlist(strsplit(NNplot[which.min(NNerr)], " ")))
nn_out <- neuralnet(y ~ x, df, hidden = besthidden, threshold = 0.01, rep = 3)
nn_results <- sapply(1:3, function(x) compute(nn_out, testdata, x)$net.result)
# Run them through the neural network
bestnn <- rowMeans(nn_results)
plotd <- data.frame(
x = testdata, spline = bestspline, nnet = bestnn, truth = fstar
)
save.image(file = "data/nnet-example.Rdata")