dist_density_fda

library(fda)

## Loading required package: splines

## Loading required package: Matrix

## 
## Attaching package: 'fda'

## The following object is masked from 'package:graphics':
## 
##     matplot

load("/home/boyazhang/repos/unifdist/code/space-fillingness_functions.RData")
l <- length(out)
msev <- rep(NA, l)
nv <- rep(NA, l)
dimv <- rep(NA, l)
design <- list()
for(i in 1:l){
  design[[i]] <- out[[i]][[1]]
  nv[i] <- out[[i]][[2]]
  dimv[i] <- out[[i]][[3]]
  msev[i] <- out[[i]][[4]]
}


# returns distance density y values for certain dim and n designs
f <- function(n, dim){
  mse1 <- msev[nv==n & dimv == dim]
  subind <- c(1:l)[nv==n & dimv == dim]
  yy <- matrix(rep(NA, 512*length(mse1)), ncol = 512)
  for(i in 1:length(subind)){
   x <- design[[subind[i]]]
   des <- density(dist(x), from = 0, to = sqrt(dim))
   yy[i,] <- des$y
  }
  return(list(yy = t(yy), mse = mse1))
}

out <- f(8,2)
#x <- out$yy
x <- out$yy-rowMeans(out$yy) ## centralize the data
mse <- log(out$mse) ## log transformation of response
data1=Data2fd(argvals=1:512, x)
fdparobj2=fdPar(data1,lambda=0.03)
smoothlist=smooth.basis(c(1:512), x, fdparobj2)

regressorl=vector("list",2)
regressorl[[1]]=rep(1,ncol(x))
regressorl[[2]]=smoothlist$fd

#beta list
conbasis = create.constant.basis(c(1,512))
betabasis = create.fourier.basis(c(1,512),11)

betalist = vector("list",2)
betalist[[1]] = conbasis
betalist[[2]] = betabasis

#----beta list with penalty---
Lcoef = c(0, (2*pi/512)^2, 0)
harmaccelLfd = vec2Lfd(Lcoef, c(0,512))
conbasis = create.constant.basis(c(1,512))
betabasis = create.fourier.basis(c(1, 512), 11)
lambda = 0.1
betafdPar = fdPar(betabasis, harmaccelLfd, lambda)
betalist[[1]] = conbasis
betalist[[2]] = betafdPar
#----cv for choosing lambda----
# loglam = seq(0,15,5)/10
# nlam = length(loglam)
# SSE.CV = matrix(0,nlam,1)
# for (ilam in 1:nlam) {
#   lambda = 10^loglam[ilam]
#   betalisti = betalist
#   betafdPar2 = betalisti[[2]]
#   betafdPar2$lambda = lambda
#   betalisti[[2]] = betafdPar2
#   fRegi = fRegress.CV(mse,regressorl,betalist)
#   SSE.CV[ilam] = fRegi$SSE.CV
# }
# plot(loglam,SSE.CV)

fRegressList1 = fRegress(mse,regressorl,betalist)

names(fRegressList1)

##  [1] "yfdPar"      "xfdlist"     "betalist"    "betaestlist" "yhatfdobj"  
##  [6] "Cmatinv"     "wt"          "df"          "OCV"         "gcv"

plot(fRegressList1$betaestlist[[2]]$fd,xlab="grid id on (0, sqrt(2))",
     ylab="Beta")

## [1] "done"

coef(fRegressList1$betaestlist[[1]])

##              [,1]
## [1,] -0.002081462

plot(mse,fRegressList1$yhatfdobj,xlab="true log-mse",
     ylab="estimated log-mse")
abline(coef = c(0,1))

#--------------to be modified---------------------------------------
#F test
yds <- mse
yhat = fRegressList1$yhatfdobj
# SSE1 = sum((yhat-yds)^2)
# SSE0 = sum((yds - mean(yds))^2)
# RSQ = (SSE0 - SSE1)/SSE0
# Fratio = ((SSE0-SSE1)/11)/(SSE1/(length(mse)-11-1))
# F.res = Fperm.fd(yds, regressorl, betalist)
#CI
resid = yds - yhat
SigmaE.= sum(resid^2)/(length(mse)-fRegressList1$df)
SigmaE = SigmaE.*diag(rep(1,length(mse)))
y2cMap = diag(rep(1,length(mse)))
stderrList = fRegress.stderr(fRegressList1, y2cMap,
                             SigmaE)
betafdPar = fRegressList1$betaestlist[[2]]
betafd = betafdPar$fd
betastderrList = stderrList$betastderrlist
betastderrfd = betastderrList[[2]]
plot(betafd, xlab="grid id on (0, sqrt(2))",
     ylab="Beta", main="Confidence interval of beta", lwd=2,ylim = c(-0.04,0.01))

## [1] "done"

lines(betafd+2*betastderrfd, lty=2, lwd=1, col = "red")
lines(betafd-2*betastderrfd, lty=2, lwd=1)
abline(h = 0)