##D. Kang, 2010-01
##Importance sampler


rm(list=ls())

fnlike <- function(ksam,t,z,sigi,sigs) {
	y <- exp(-1*ksam*t)
	sig <- sigi+ sigs*y

	arg1 <- (z-exp(-1*ksam*t))^2
	arg2 <- 2*sig*sig
	arg3 <- -1*arg1 / arg2

	denom <- sig*sqrt(2*3.14)
	like <- exp(arg3) / denom
	ln <- like
	ln
}


fnsamnorm <- function(kavg,ksig, n) {
	norsam <- ksig*rnorm(n, 0, 1)+kavg
	ksam <- norsam
	ksam
}



fnsamln <- function(klogavg, klogsig, n) {
	norsam <- klogsig*rnorm(n, 0, 1)+klogavg
	ksam <- exp(norsam)
	ksam
}



#set.seed(12345)

model <- 'y=exp(-Kel*t)'
variancemodel <- 'sig=sigi+sigs*y'
prior <- 'N'

nsam <- 1000		# posterior samples to generate
kavg <- 1			# mean of prm k, prior
ksig <- .3		      # std of prm k, prior
sigi <- .05		      # std prm. of obs. error model, intercept (additive)
sigs <- .15	            # std prm. of obs. error model, slope (proportional)
#z <- 0.905              # obs. for Kel <- 0.1, exp(-0.1)
#z <- .607			# obs. for Kel <- 0.5, exp(-0.5)
z <- .368			# obs. for Kel <- 1, exp(-1)
#z <- 0.223              # obs. for Kel <- 1.5, exp(-1.5)
#z <- .135			# obs. for Kel <- 2, exp(-2)
t <- 1		      # obs. time

par(mfrow=c(3,1))

#generate as many candidate samples as specified 
ksam <- fnsamnorm(kavg,ksig,nsam)
hist(ksam, freq=FALSE, prob=TRUE, main=paste('Prior ~ N(', kavg, ', ', ksig,'^2',')', sep=''), xlim=c(0,3))

like.ksam <- fnlike(ksam,t,z,sigi,sigs)
plot(ksam, like.ksam, xlim=c(0,3))
title('Likelihood')

weit <- like.ksam/sum(like.ksam)
#print(sum(weit))

print(post.mean <- sum(weit*ksam))
print(post.var <- sum(weit*ksam^2)-post.mean^2)
print(post.sig <- sqrt(post.var))

post.dist <- 1/(post.sig*sqrt(2*3.14))*exp(-(ksam-post.mean)^2/(post.sig*post.sig))
plot(ksam, post.dist, xlim=c(0,3))
title(paste('Posterior ~ N(', signif(post.mean,4), ', ', signif(post.sig,4),'^2',')', sep=''))