rm(list=ls())

n=260
x_pars=c(0,1);
y_pars=c(1,.3,7);

#observed data follows  gamma-shifted distribution
y=rgamma(n,y_pars[1],y_pars[2])-y_pars[3]
sy=sort(y)
#quantiles come from the standard normal distribution
p=pnorm(sy,x_pars[1],x_pars[2])

############################
# Linear Interpolation
############################
w=abs(diff(sy)) #interval width
dp=abs(diff(p)) #differenced quantiles
fx=dp/w         #estimated hight of pdf
windows()
par(mfrow=c(2,1))
plot(sy[-n],fx,typ='l',main="Linear Interpolation",
      ylab="density",xlab="y",xlim=c(-3,3),ylim=c(0,.5)) #estimated pdf
points(sy,dnorm(sy,x_pars[1],x_pars[2]),type='l',
       col=2,lty=2) #actual pdf
points(sy,dgamma(sy+y_pars[3],y_pars[1],y_pars[2]),
       type='l',col=4,lty=4) #pdf of observations
legend("topright",legend=c("Actual PDF","Estimated PDF","Obs. PDF"),
       lty=c(2,1,4),col=c(2,1,4))

############################
# Monte Carlo Smoothing
############################
precision=100 #if this is to high the KDE may have difficulty + longer run-time
m=n*precision


u=sort(runif(m))

ind=which(u>p[1] & u<=p[n]) #throw these out for now (assume  0 density on either side of the observations)
m2=length(ind)
u=u[ind];

#simulate from linear interp. cdf using probability integral transform
z=numeric(m2)
st=1
for(i in 1:m2){
  for(j in 2:n){
    if(u[i]<=p[j]){
      st=j;
      interp=(u[i]-p[j-1])/dp[j-1]
      z[i]=sy[j-1]+w[j-1]*interp
      break
    }
  }
}


plot(density(z),type='l',main="Linear Interpolation + MC Smoothing",
     ylab="density",xlab="y",xlim=c(-3,3),ylim=c(0,.5)) #estimated pdf
points(sy,dnorm(sy,x_pars[1],x_pars[2]),type='l',
       col=2,lty=2) #actual pdf
points(sy,dgamma(sy+y_pars[3],y_pars[1],y_pars[2]),
       type='l',col=4,lty=4) #pdf of observations
legend("topright",legend=c("Actual PDF","Estimated PDF","Obs. PDF"),
       lty=c(2,1,4),col=c(2,1,4))