############################################
# STAT 740 R examples
# Tim Hanson
# U. South Carolina
# September 18, 2017
############################################

############################################
# compute E(x^4) for x ~ N(1,4)
# Note: error bands are pointwise!
############################################

m=1000
g=Vectorize(function(x){x^4}) # function g(x)
x=g(rnorm(m,1,2)) # X ~ N(1,4)
estint=cumsum(x)/(1:m) 
esterr=sqrt(cumsum((x-estint)^2))/(1:m)
plot(estint,xlab="Mean and error range",type="l",lwd=2,
ylim=mean(x)+20*c(-esterr[m],esterr[m]),ylab="") 
lines(estint+1.96*esterr,col="gold",lwd=2) 
lines(estint-1.96*esterr,col="gold",lwd=2) 
abline(a=73,b=0,col="blue",lwd=2)

############################################
# compute 90th percentile of x^4 for x ~ N(1,4)
# & 95% equal-tailed probability interval (PI)
# & 95% highest-density PI
############################################

m=100000
x=rnorm(m,1,2)
hist(g(x),freq=F,main="Histogram of g(X)",xlim=c(0,1000),breaks=100)
quantile(g(x),0.9)

cat("equal-tailed 95% PI=(",quantile(g(x),0.025),",",quantile(g(x),0.975),") \n")

intervals=20
a=seq(0,0.05,length=intervals)
m=matrix(intervals*2,intervals,2)
for(i in 1:intervals){m[i,]=quantile(g(x),c(a[i],a[i]+0.95))}
a.min=which.min(m[,2]-m[,1])
cat("highest density 95% PI=(",quantile(g(x),a[a.min]),",",
 quantile(g(x),0.95+a[a.min]),") \n")

############################################
# different quantiles for bumpy density
############################################

f=Vectorize(function(x){exp(-0.5*x^2)*(sin(6*x)^2+3*cos(x)^2*sin(4*x)^2+1)})
x=seq(-4,4,length=1000)
plot(x,f(x),type="l")
coi=integrate(f,-5,5)
coi

M=max(f(x)/dnorm(x,0,1))
 
sim=function(n){
 y=1:n
 for(i in 1:n){
  u=-1; x=0
  repeat{
   u=runif(1); x=rnorm(1)
   if(u<=f(x)/(M*dnorm(x,0,1))) break
  } # end repeat
 y[i]=x
 } # i=1,...,n
 y
}

d=sim(100000)
quantile(d,0.9) # weighted version
# all identical for p*m that is an integer
sort(d)[round(0.9*100000)]
sort(d)[floor(0.9*100000)]
sort(d)[ceiling(0.9*100000)]