##################################################
### SIR, Sampling Importance Resampling, simple example

### f is N(0,1), g is t_5


## draw x from g
n = 1000
set.seed(99)
tdf = 5
x = rt(n,df=tdf)

## look at weights function on a grid a values
xv = seq(from=-4,to=4,length.out=1000)
dnx = dnorm(xv)
dtx = dt(xv,df=tdf)
par(mfrow=c(2,2))

plot(xv,dnx,type="l")
title(main="f=N(0,1)",cex.main=1.5)

plot(xv,dtx,type="l")
title(main="g=t5",cex.main=1.5)

wv = dnx/dtx
plot(xv,wv,type="l")
title(main="not normalized w= f/g = nor/t",cex.main=1.5)

wsv = wv/sum(wv)
plot(xv,wsv,type="l")
title(main="normalized w",cex.main=1.5)

#weights on draws

dnx = dnorm(x)
dtx = dt(x,df=tdf)
wv = dnx/dtx
wsv = wv/sum(wv)

## plot weights on draws
par(mfrow=c(1,1))
plot(wsv)

## plot normalized weight at each x from g
plot(range(x),range(wsv),type="n")
for(i in 1:n) {
   lines(c(x[i],x[i]),c(0,wsv[i]))
}

## resample
m = 5000
xrs = sample(x,m,replace=TRUE,prob=wsv)

## plot the resampled draws, should look iid standard normal
par(mfrow=c(1,2))
hist(xrs,prob=TRUE)
qqnorm(xrs)
abline(0,1,col="red")