-
Notifications
You must be signed in to change notification settings - Fork 0
/
Accept_Reject.R
105 lines (85 loc) · 3.18 KB
/
Accept_Reject.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#--------------------------------
# Acceptance-Rejection sampling #
#--------------------------------
#----------
# Example 1
#----------
# We use Accept-Reject to estimate the mean and variance of a beta distribution
n = 10000 # number of draws
f = function(x) {dbeta(x,6.3,2.7)} # target distribution
c = 2.67 # maximum (found using for ex:
optimize(f, interval=c(0, 1), maximum=TRUE)
g = dunif # proposal distribution
# Accept-Reject algorithm
theta=c()
for (i in 1:n){
x <- runif(1, 0, 1)
y <- runif(1, 0, c)
acceptance_rate <- f(x) / (c*g(x)) # acceptance probability
if (y <= acceptance_rate){
theta[i] <- x
}
else {next}
}
# Plotting the results along with the true distribution
hist(theta, breaks=100, freq=FALSE, main="Random draws from beta(6.3, 2.7)")
lines(x, dbeta(x, 6.3, 2.7))
# and let's see if we get the same results
cat("Estimated mean is equal to", mean(theta, na.rm=TRUE))
quantile(theta, probs=c(0.05,0.95), na.rm=TRUE)
cat("Standard deviation is equal to", sd(theta, na.rm=TRUE))
# check
mean(rbeta(10000, 6.3, 2.7))
sd(rbeta(10000, 6.3, 2.7))
# same results as when using the inbuilt function rbeta()
points = runif(n)
uniforms = runif(n)
accept = uniforms < (f(points)/(c*g(points)))
# plotting
curve(f, lwd=2, main="Acceptance - Rejection method to sample random variates from a Beta(2.7, 6.3)")
points(points, c*uniforms, pch=ifelse(accept,4,4), col=ifelse(accept,"darkred","gray60"), cex=0.5)
legend("bottomleft", c("f","accepted","rejected"),
lwd=c(2,NA,NA), col=c("black","darkred","gray60"),
pch=c(NA,4,4), bg="white")
#----------
# Example 2
#----------
set.seed(2023)
# we want to learn about the mean survival time of butterflies using bayesian analysis
n = 10000 # number of draws
# density of exponential distribution with parameter theta = 1
g <- function(x) return( exp(-x) ) # proposal, easy to sample from
# posterior distribution of the data
# 13 = 6 + 3 + 2 + 2 : number of butterflies alive
# 10 = 4 + 3 + 1 + 0 + 2 : number of butterflies dead
# 1/x is the Jeffrey's prior
# exp(-13*x) * (1-exp(-x))^10 ) is the likelihood of the data
f <- function(x) return( 1/x *exp(-13*x) * (1-exp(-x))^10 ) # target (posterior), hard to sample from
optimize(f, interval=c(0, 100), maximum=TRUE)
c <- 5 * 10^(-7)
# Accept-Reject algorithm
theta=c()
for (i in 1:n){
x <- runif(1, 0, 1)
y <- runif(1, 0, c)
acceptance_rate <- f(x) / (c*g(x)) # acceptance probability
if (y <= acceptance_rate){
theta[i] <- x
}
else {next}
}
# and let's see if we get the same results
cat("Estimated EST is equal to", 1/mean(theta, na.rm=TRUE))
# Estimated EST is equal to 1.946919
points = rexp(n)
uniforms = runif(n, min = 0, max = c*g(points))
accepted = uniforms < f(points)
# plotting
curve(f, lwd=2, main="Acceptance - Rejection method to sample random variates from the target distribution")
points(points, uniforms, pch=ifelse(accepted,4,4), col=ifelse(accepted,"darkred","gray60"), cex=0.5)
legend("bottomleft", c("f","accepted","rejected"),
lwd=c(2,NA,NA), col=c("black","darkred","gray60"),
pch=c(NA,4,4), bg="white")
#----
# end
#----