Simulation of Exp-Abs-xy
Target distribution:
f(x,y)exp(xyayx).f(x, y)\propto \exp(-\vert x\vert -\vert y\vert - a\vert y-x\vert)\,.
Gibbs:
    f(xy)exp(xayx)f(x\mid y)\propto \exp(-\vert x\vert - a\vert y-x\vert)
    f(yx)exp(yayx)f(y\mid x)\propto \exp(-\vert y\vert - a\vert y-x\vert)

Accept-Reject

f(xy)exp(xayx)exp(axy)f(x\mid y)\propto \exp(-\vert x\vert - a\vert y-x\vert)\le \exp(-a\vert x - y\vert)
1
using Distributions
Copied!
1
a = 10
2
f(x, y) = exp( -abs(x) - a * abs(y-x) )
3
g(x, y) = exp(-a * abs(x-y))
4
5
function gibbs_ar(N::Int)
6
x = ones(N+1)
7
y = ones(N+1)
8
for i = 2:(N+1)
9
xstar = rand(Laplace(y[i-1], 1/a))
10
if rand() < f(xstar, y[i-1]) / g(xstar, y[i-1])
11
x[i] = xstar
12
else
13
x[i] = x[i-1]
14
end
15
16
ystar = rand(Laplace(x[i], 1/a))
17
if rand() < f(x[i], ystar) / g(x[i], ystar)
18
y[i] = ystar
19
else
20
y[i] = y[i-1]
21
end
22
end
23
return x, y
24
end
Copied!
1
gibbs_ar (generic function with 1 method)
Copied!
1
x, y = gibbs_ar(1000)
Copied!
1
([1.0, 1.0, 1.0, 1.0, 1.0, 1.07195, 0.930769, 0.713139, 0.546008, 0.666174 … 0.0213817, 0.0968143, -0.0422362, 0.0506649, 0.315365, 0.0955934, -0.13926, -0.099665, -0.611045, -0.611045], [1.0, 0.958085, 0.958085, 0.958085, 0.958085, 0.958085, 0.906274, 0.560885, 0.61187, 0.734409 … 0.140089, 0.0351813, -0.0497683, 0.116464, 0.116464, -0.0539619, -0.0539619, -0.478698, -0.549537, -0.549537])
Copied!
1
using Plots
Copied!
1
plot(x, y, seriestype = :scatter, legend = false)
Copied!
svg

Inverse CDF

Refer to the handwritten supplementary for detailed derivation.
1
using StatsBase
Copied!
1
function one_gibbs(y)
2
p1 = exp(-a*y) * ( exp((a-1)*y) - 1 ) / (a - 1) * (y >= 0)
3
p2 = exp(-a*y) * exp((a+1)*min(0, y)) / (a + 1)
4
p3 = exp(a*y) * exp(-(1+a)*max(0, y)) / (a + 1)
5
p4 = exp(a*y) * ( exp((1-a)*y) - 1) / (a - 1) * (y < 0)
6
flag = sample(1:4, pweights([p1, p2, p3, p4]))
7
u = rand()
8
if flag == 1
9
x = log( 1 + u * ( exp((a-1)*y) - 1 ) ) / (a - 1)
10
elseif flag == 2
11
x = log(u) / (a + 1) + min(0, y)
12
elseif flag == 3
13
x = -log(u) / (a + 1) + max(0, y)
14
else
15
x = log( 1 - u * ( 1 - exp((1-a)*y) ) ) / (1 - a)
16
end
17
return x
18
end
19
20
function gibbs_inv(N::Int)
21
x = ones(N+1)
22
y = ones(N+1)
23
for i = 2:(N+1)
24
x[i] = one_gibbs(y[i-1])
25
y[i] = one_gibbs(x[i])
26
end
27
return x, y
28
end
Copied!
1
gibbs_inv (generic function with 1 method)
Copied!
1
one_gibbs(3)
Copied!
1
3.2469523161594838
Copied!
1
x, y = gibbs_inv(1000)
Copied!
1
([1.0, 1.01954, 1.25488, 1.31049, 1.25068, 1.32997, 1.38793, 1.24884, 1.16866, 1.22989 … 0.508065, 0.397666, 0.56615, 0.735377, 0.556286, 0.502009, 0.125919, 0.233455, 0.522778, 0.312981], [1.0, 1.24181, 1.271, 1.18791, 1.27124, 1.34111, 1.19227, 1.09944, 1.1339, 1.19598 … 0.54177, 0.36308, 0.677268, 0.514871, 0.640218, 0.316574, 0.256377, 0.382262, 0.270825, 0.0834128])
Copied!
1
plot(x, y, seriestype = :scatter, legend = false)
Copied!
svg

Slice Sampling

f(x,y)exp(xyaxy)=0Iu1exp(x)du10Iu2exp(y)du20Iu3exp(axy)du3f(x, y)\propto \exp(-\vert x\vert - \vert y\vert - a\vert x- y\vert) = \int_0^\infty \mathbb{I}_{u_1\le\exp(-\vert x\vert)}du_1\int_0^\infty \mathbb{I}_{u_2\le\exp(-\vert y\vert)}du_2\int_0^\infty \mathbb{I}_{u_3\le\exp(-a\vert x-y\vert)}du_3
Then
xy,u1,u2,u3f(xy,u1,u2,u3)Iu1exp{x}Iu3exp{axy}yx,u1,u2,u3f(yx,u1,u2,u3)Iu2exp{y}Iu3exp{axy}u1xU(0,exp{x})u3x,yU(0,exp{axy})u2yU(0,exp{x})\begin{align*} x|y,u_1,u_2,u_3&\sim f(x|y,u_1,u_2,u_3)\propto\mathbb{I}_{u_1\le\exp\{-|x|\}}\mathbb{I}_{u_3\le\exp\{-a|x-y|\}}\\ y|x,u_1,u_2,u_3&\sim f(y|x,u_1,u_2,u_3)\propto\mathbb{I}_{u_2\le\exp\{-|y|\}}\mathbb{I}_{u_3\le\exp\{-a|x-y|\}}\\ u_1|x&\sim\mathcal U(0,\exp\{-|x|\})\\ u_3|x,y&\sim\mathcal U(0,\exp\{-a|x-y|\})\\ u_2|y&\sim\mathcal U(0,\exp\{-|x|\})\\ \end{align*}
It is straightforward to use the following Julia program to simulate.
1
using Distributions
2
a = 10
3
N = 1000
4
x = ones(N)
5
y = ones(N)
6
for t = 2:N
7
u1 = rand() * exp(-abs(x[t-1]))
8
u2 = rand() * exp(-abs(y[t-1]))
9
u3 = rand() * exp(-a*abs(x[t-1]-y[t-1]))
10
x[t] = rand(Uniform(log(u1), -log(u1)))
11
while a*abs(x[t]-y[t-1]) > -log(u3)
12
x[t] = rand(Uniform(log(u1), -log(u1)))
13
end
14
y[t] = rand(Uniform(log(u2), -log(u2)))
15
while a*abs(x[t]-y[t]) > -log(u3)
16
y[t] = rand(Uniform(log(u2), -log(u2)))
17
end
18
end
19
using Plots
20
plot(x, y, seriestype=:scatter, legend=false)
Copied!

Reference

Last modified 2yr ago