Skip to content

Commit 45f3d82

Browse files
author
mauriciogtec
committed
#1 changes to README for 1.0, and added elastic net example
1 parent 34bd1cd commit 45f3d82

File tree

1 file changed

+41
-12
lines changed

1 file changed

+41
-12
lines changed

README.md

Lines changed: 41 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -36,9 +36,10 @@ Let's verify the result
3636

3737

3838
```julia
39-
x = -linspace(-10.0, 10.0, 100)
40-
envelop = eval_envelop.(sampler.envelop, x)
41-
target = f.(x)
39+
# Plot the results and compare to target distribution
40+
x = range(-10.0, 10.0, length=100)
41+
envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
42+
target = [f(xi) for xi in x]
4243

4344
histogram(sim, normalize = true, label = "Histogram")
4445
plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])
@@ -56,12 +57,14 @@ plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])
5657
f(x) = β^α * x^-1) * exp(-β*x) / gamma(α)
5758
support = (0.0, Inf)
5859

60+
# Build the sampler and simulate 10,000 samples
5961
sampler = RejectionSampler(f, support)
6062
@time sim = run_sampler!(sampler, 10000)
6163

62-
x = linspace(0.0, 5.0, 100)
63-
envelop = eval_envelop.(sampler.envelop, x)
64-
target = f.(x)
64+
# Plot the results and compare to target distribution
65+
x = range(0.0, 5.0, length=100)
66+
envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
67+
target = [f(xi) for xi in x]
6568

6669
histogram(sim, normalize = true, label = "Histogram")
6770
plot!(x, [target envelop], width = 2, label = ["Gamma(α, β)" "Envelop"])
@@ -71,8 +74,6 @@ plot!(x, [target envelop], width = 2, label = ["Gamma(α, β)" "Envelop"])
7174

7275

7376

74-
75-
7677
![](img/example2.png)
7778

7879
## Truncated distributions and unknown normalisation constant
@@ -85,12 +86,14 @@ We don't to provide an exact density--it will sample up to proportionality--and
8586
f(x) = β^α * x^-1) * exp(-β*x) / gamma(α)
8687
support = (1.0, 3.5)
8788

89+
# Build the sampler and simulate 10,000 samples
8890
sampler = RejectionSampler(f, support)
8991
@time sim = run_sampler!(sampler, 10000)
9092

91-
x = linspace(0.01, 8.0, 100)
92-
envelop = eval_envelop.(sampler.envelop, x)
93-
target = f.(x)
93+
# Plot the results and compare to target distribution
94+
x = range(0.01, 8.0, length=100)
95+
envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
96+
target = [f(xi) for xi in x]
9497

9598
histogram(sim, normalize = true, label = "histogram")
9699
plot!(x, [target envelop], width = 2, label = ["target density" "envelop"])
@@ -99,7 +102,33 @@ plot!(x, [target envelop], width = 2, label = ["target density" "envelop"])
99102
0.007766 seconds (181.82 k allocations: 3.024 MiB)
100103

101104

105+
![](img/example3.png)
106+
107+
## Elastic Net distribution
108+
109+
The following example arises from elastic net regression and smoothing problems. In these cases, the integration constants are not available analytically.
110+
111+
```julia
112+
# Define function to be sampled
113+
function f(x, μ, λ1, λ2)
114+
δ = x - μ
115+
nl = λ1 * abs(δ) + λ2 * δ^2
116+
return exp(-nl)
117+
end
118+
support = (-Inf, Inf)
119+
120+
# Build the sampler and simulate 10,000 samples
121+
μ, λ1, λ2 = 0.0, 2.0, 1.0
122+
sampler = RejectionSampler(x -> f(x, μ, λ1, λ2), support, max_segments = 5)
123+
@time sim = run_sampler!(sampler, 10000);
102124

125+
# Plot the results and compare to target distribution
126+
x = range(-2.3, 2.3, length=100)
127+
envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
128+
target = [f(xi, μ, λ1, λ2) for xi in x]
103129

130+
histogram(sim, normalize = true, label = "histogram")
131+
plot!(x, [target envelop], width = 2, label = ["target density" "envelop"])
132+
```
104133

105-
![](img/example3.png)
134+
![](img/example4.png)

0 commit comments

Comments
 (0)