1. ## Re: Bayesian Statistics Question

The thing is, I tried to use gamma density as proposal and couldn't find the M

2. ## Re: Bayesian Statistics Question

Hey, do you know why this code does not seem to work and seems to be running infinitely and I have to stop the process.

Code:
v<-runif(1)
n<-50 # sample size

y<-numeric(n)
accepted.samples<-0
while (accepted.samples<n) {
x<-rgamma(n,5,3)
f.at.x<-x^4*exp(-3*x)*(1-exp(-2*x))^7
g.at.x<-dgamma(x, 5, 3)
accept<-(f.at.x) / (g.at.x)
if (v<accept) {
accepted.samples<-accepted.samples+1
y[accepted.samples]<-x
}
}
It's something to do with the M * missing

3. ## Re: Bayesian Statistics Question

You cannot use 1 uniform random variable to decide whether all proposal sample are accepted or not. Each sample require 1 independent uniform random variable.

Also the optimal is very easy to find.

So according to Dason's code, just need to change

Code:
 id.accept <- (u <= 10.125*r)
then will improve the acceptance rate by almost 10 times from 0.06 to 0.65.

4. ## Re: Bayesian Statistics Question

I used another proposal function which can be integrated analytically and inverted. Therefore I only need uniform sampling in my code. Gamma sampling also uses rejection sampling.

Code:
tmp <- function(x){
x^4*exp(-3*x)*(1-exp(-2*x))^7
}

k <- 1/integrate(tmp, 0, Inf)$value desired.density <- function(x){k*tmp(x)} n <- 100000 proposals <- tan(pi*runif(n)) + 1.6 u <- runif(n) r <- tmp(proposals)/0.045*(1 + (proposals - 1.6)^2) id.accept <- (u <= r) samp <- proposals[id.accept] accept.rate <- mean(id.accept) print(accept.rate) plot(density(samp)) curve(desired.density, add = T, col = "red") 5. ## Re: Bayesian Statistics Question This probably wouldn't be good for your assignment but there is also brute force sampling. The basic idea is that you grid up the support of your distribution. Since the distribution of interest has infinite support we need to find where essentially all of the mass is contained. Then you sample points from your grid with probability proportional to the density of interest. Code: tmp <- function(x){ x^4*exp(-3*x)*(1-exp(-2*x))^7 } integrate(tmp, 0, Inf)$value
# [1] 0.06483012
integrate(tmp, 0, Inf)\$value -> k
f <- function(x){1/k*tmp(x)}
integrate(f, 0, 5)
# 0.9986952 with absolute error < 1e-05
integrate(f, 0, 10)
# 1 with absolute error < 1e-05
grid <- seq(0, 10, .00001)
p <- tmp(grid)
?sample
# starting httpd help server ... done
samp <- sample(grid, 10000, replace = TRUE, prob = p)

6. ## Re: Bayesian Statistics Question

To fanky: you try to use the Cauchy as the proposal, which will generate negative random variates. Now the target density evaluated at negative points is negative, which make the code u <= r also reject those invalid case as the ratio r < 0; but in general it is not true and need to be more careful.

7. ## Re: Bayesian Statistics Question

Thanks to everybody who helped. You don't know how much this was appreciated!