# Motivation
library(hdrm)
Golub1999 <- read_data(Golub1999)
summ <- summary(lm(Golub1999$X~Golub1999$y))
tstat <- sapply(summ, function(s) s$coef[2,3])
p <- sapply(summ, function(s) s$coef[2,4])
z <- qnorm(p/2) * sign(-tstat)
h <- length(p)
q <- p.adjust(p, method='BH')
sum(q < 0.01)
range(p[q < 0.01])
range(abs(z[q < 0.01]))

# Illustrating the concept of kernel density estimation
x <- c(1, 3, 5, 6, 6.2, 6.6, 7, 8, 10)
n <- length(x)
xx <- seq(-2, 10, len=101)
hist(x, border="white", main="", freq=FALSE, las=1)
rug(x, col="red", lwd=3)
for (i in 1:n) lines(xx, dnorm(xx, x[i], 0.75)/10, col="#008DFFB2", lwd=3)
d <- density(x, bw=0.75)
lines(d, lwd=3)

# Density and choice of bandwidth
z6 <- z[z > -6 & z < 6]
brk <- seq(-6, 6, length = 99)
# First plot
par(mfrow=c(1,3))
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
# Second plot
par(mfrow=c(1,3))
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
density(z6, bw=.05) |> lines(col='slateblue', lwd=3)
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
# Third plot
par(mfrow=c(1,3))
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
density(z6, bw=.05) |> lines(col='slateblue', lwd=3)
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
density(z6, bw=1.5) |> lines(col='slateblue', lwd=3)
# Fourth plot
par(mfrow=c(1,3))
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
density(z6, bw=.05) |> lines(col='slateblue', lwd=3)
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
density(z6) |> lines(col='slateblue', lwd=3)
h <- hist(z6, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, xlab='z')
density(z6, bw=1.5) |> lines(col='slateblue', lwd=3)

# Local FDR illustration
f <- function(z, pi0=1, fig=4, ...) {
  # Calculation
  dens <- density(z, bw="nrd")
  f <- approxfun(dens$x, dens$y)
  f0 <- function(z) pi0*dnorm(z)
  lfdr <- pmin(f0(z)/f(z), 1)

  # Plot
  h <- hist(z, breaks=seq(min(z), max(z), length = 99), plot=FALSE)
  zz <- seq(min(z), max(z), len=299)
  ylim=c(0, max(c(h$density, dens$y, f0(0))))
  if (fig > 0) plot(h, main="", border=FALSE, col="lightgray", freq=FALSE, las=1, ylim=ylim)
  if (fig > 1) lines(dens, col=breheny::pal(2)[1], lwd=2)
  if (fig > 2) lines(zz, f0(zz), col=breheny::pal(2)[2], lwd=2)
  if (fig > 3) {
    fdr_zz <- f0(h$mids)/f(h$mids)
    y <- pmax(h$density * (1 - fdr_zz), 0)
    for (k in 1:length(h$mids)) lines(rep(h$mids[k],2), c(0, y[k]), lwd = 2, col = "red")
  }
  return(invisible(lfdr))
}
f(z, fig=1)
f(z, fig=2)
f(z, fig=3)
lfdr1 <- f(z, fig=4)

# Estimating pi0
lfdr2 <- f(z, pi0=.53)

par(mar=c(4,4,2,0.5))
f1 <- approxfun(z, lfdr1)
f2 <- approxfun(z, lfdr2)
zz <- seq(-5, 5, len=99)
plot(zz, f1(zz), type='l', lwd=3, bty='n', las=1, ylim=c(0, 1), col=hdrm::pal(2)[1],
     xlab='z', ylab='Local FDR')
lines(zz, f2(zz), lwd=3, col=hdrm::pal(2)[2])
hdrm::toplegend(legend=c(expression(hat(pi)[0]*'='*1), expression(hat(pi)[0]*'='*0.53)),
          lwd=3, col=hdrm::pal(2))
plot(zz, f1(zz), type='l', lwd=3, bty='n', las=1, ylim=c(0, 1), col=hdrm::pal(2)[1],
     xlab='z', ylab='Local FDR')
lines(zz, f2(zz), lwd=3, col=hdrm::pal(2)[2])
hdrm::toplegend(legend=c(expression(hat(pi)[0]*'='*1), expression(hat(pi)[0]*'='*0.53)),
          lwd=3, col=hdrm::pal(2))
lines(zz, rep(0.1, length(zz)), lty=2)

uniroot(function(x) f1(x)-0.1, lower=2, upper=4)$root
uniroot(function(x) f2(x)-0.1, lower=2, upper=4)$root
sum(lfdr1 < .1)
sum(lfdr2 < .1)