+ Reply to Thread
Page 4 of 4 FirstFirst 1 2 3 4
Results 46 to 52 of 52

Thread: Bayesian Statistics Question

  1. #46
    Points: 884, Level: 15
    Level completed: 84%, Points required for next Level: 16

    Posts
    79
    Thanks
    2
    Thanked 0 Times in 0 Posts

    Re: Bayesian Statistics Question




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

  2. #47
    Points: 884, Level: 15
    Level completed: 84%, Points required for next Level: 16

    Posts
    79
    Thanks
    2
    Thanked 0 Times in 0 Posts

    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
    Last edited by MegaMan; 12-12-2011 at 03:29 AM.

  3. #48
    TS Contributor
    Points: 22,410, Level: 93
    Level completed: 6%, Points required for next Level: 940

    Posts
    3,020
    Thanks
    12
    Thanked 565 Times in 537 Posts

    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 M is very easy to find.

    M = \sup_{x > 0} \frac {x^4e^{-3x}(1 - e^{-2x})^7} {\frac {1} {\Gamma(5)}3^5x^4e^{-3x}} = \sup_{x > 0} \frac {4!} {243} (1 - e^{-2x})^7 = \frac {1} {10.125}

    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. #49
    Points: 487, Level: 9
    Level completed: 74%, Points required for next Level: 13

    Posts
    4
    Thanks
    0
    Thanked 1 Time in 1 Post

    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. #50
    Devorador de queso
    Points: 95,940, Level: 100
    Level completed: 0%, Points required for next Level: 0
    Awards:
    Posting AwardCommunity AwardDiscussion EnderFrequent Poster
    Dason's Avatar
    Location
    Tampa, FL
    Posts
    12,937
    Thanks
    307
    Thanked 2,630 Times in 2,246 Posts

    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. #51
    TS Contributor
    Points: 22,410, Level: 93
    Level completed: 6%, Points required for next Level: 940

    Posts
    3,020
    Thanks
    12
    Thanked 565 Times in 537 Posts

    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. #52
    Points: 884, Level: 15
    Level completed: 84%, Points required for next Level: 16

    Posts
    79
    Thanks
    2
    Thanked 0 Times in 0 Posts

    Re: Bayesian Statistics Question


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

+ Reply to Thread
Page 4 of 4 FirstFirst 1 2 3 4

           




Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts






Advertise on Talk Stats