####################################################
#Importance sampling: Gamma

#version I - know normalizing constants Cp, Cq
#Suppose we have samples from Gamma(a,b) (say)
a<-7
b<-2

#simulate n from the ground up
n<-100000
U<-matrix(-log(runif(a*n)),a,n)
Y<-apply(U,2,'sum')/b

#now suppose need mean at
alpha<-7.8
beta<-3.1
exp_Y1<-mean(Y*Y^(alpha-a)*exp(-(beta-b)*Y)*(gamma(a)/gamma(alpha))*(beta^alpha/b^a))

exp_Y1
alpha/beta

#################################################################
#see the estimator break down as alpha moves away from a

m<-15
alpha<-seq(from=1,to=15,length.out=m)
beta<-3.1

n<-100000
a<-7
b<-2
U<-matrix(-log(runif(a*n)),a,n)
Y<-apply(U,2,'sum')/b

ex<-rep(NA,m)
#calulate the mean of f(X)=1 in Gamma(alpha,beta) using Gamma(a,b)
for (i in 1:m) {
	ex[i]<-mean(Y^(alpha[i]-a)*exp(-(beta-b)*Y)
			*(gamma(a)/gamma(alpha[i]))*(beta^alpha[i]/b^a))
}
plot(alpha,ex)
lines(alpha,alpha^0)

#####################################################
#version II dont know normalizations

num_Y<-mean(Y*Y^(alpha-a)*exp(-(beta-b)*Y))
den_Y<-mean(Y^(alpha-a)*exp(-(beta-b)*Y))
exp_Y2<-num_Y/den_Y

exp_Y2
alpha/beta

den_Y
(gamma(alpha)/gamma(a))*(b^a/beta^alpha)