Importance sampling

################################
#IMPORTANCE SAMPLING
#
# int{ h(x)*p(x) }dx = int{ h(x)*(p(x)/g(x))*g(x) }dx = int{ h(x)*w(x)*g(x) }dx
#
################################

#Importance SAMPLING INTEGRAL
w <- function(x) dunif(x,0.01,1)/dbeta(x,0.7,1)
f <- function(x) x^(-0.5)
X <- rbeta(1000,0.7,1)
Y <- w(X)*f(X)
c(mean(Y),var(Y))  #TRUE IS 1.8

## [1] 1.810358 0.359825

# http://www.statistik.tuwien.ac.at/public/filz/students/seminar/ss11/kuerner_vortrag.pdf

ch=function(la){
integrate(function(x){x^(la-1)*exp(-x)},0,Inf)$val
}

plot(
lgamma(seq(.01,10,le=100)),
log(apply(as.matrix(seq(.01,10,le=100)),1,ch)),
xlab="log(integrate(f))",
ylab="log(gamma(lambda))",pch=19,cex=.6
)

cac=rcauchy(10)+350

lik=function(the){
    u=dcauchy(cac[1]-the)
    for (i in 2:10)
    u=u*dcauchy(cac[i]-the)
    return(u)
}

integrate(lik,-Inf,Inf)

## 4.826682e-44 with absolute error < 9.5e-44

integrate(lik,200,400)

## 1.953913e-08 with absolute error < 3.6e-08

##

 h=function(x){(cos(50*x)+sin(20*x))^2}
 par(mar=c(2,2,2,1),mfrow=c(2,1))
 curve(h,xlab="Function",ylab="",lwd=2)
 integrate(h,0,1)

## 0.9652009 with absolute error < 1.9e-10

x=h(runif(10^4))
 estint=cumsum(x)/(1:10^4)
 esterr=sqrt(cumsum((x-estint)^2))/(1:10^4)
 plot(estint,xlab="Mean and error range",type="l",lwd=2,
 ylim=mean(x)+20*c(-esterr[10^4],esterr[10^4]),ylab="")
 lines(estint+2*esterr,col="gold",lwd=2)
 lines(estint-2*esterr,col="gold",lwd=2)

#http://www.stats.ox.ac.uk/~winkel/ASim11.pdf

p<-0.5
n<-1000000
X<-rep(NA,n)
X[1]<-0
for (j in 2:n) {
X[j]<-X[j-1]+2*(runif(1)<=p)-1
}
time<-1:n
plot(time,X,type="l")