# Simple simulation
**URL:** https://discourse.julialang.org/t/simple-simulation/135701
**Category:** Performance
**Created:** [February 17, 2026, 9:57am UTC](https://discourse.julialang.org/t/simple-simulation/135701 "2026-02-17T09:57:18Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@Lincoln\_Hannah](https://discourse.julialang.org/u/Lincoln_Hannah)
#### Post date: [February 17, 2026, 9:57am UTC](https://discourse.julialang.org/t/simple-simulation/135701/1 "2026-02-17T09:57:18Z")
Is there an easy way to speed this up?
```julia-auto
sims = 10_000
T = 1
σ = .1
Δ = .000001
σrΔ = σ√Δ
steps = (Int64 ∘ round)(T/Δ)
L = .9
H = 1.1
x = ones( sims )
alive = trues( sims )
for t ∈ 1:steps
x .+= randn(sims)σrΔ
@. alive &= L < x < H
end
mean(alive)
```
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [February 17, 2026, 10:15am UTC](https://discourse.julialang.org/t/simple-simulation/135701/2 "2026-02-17T10:15:25Z")
around 2x faster
```julia
julia> function main()
sims = 10_000
T = 1
σ = .1
Δ = .000001
σrΔ = σ√Δ
steps = (Int64 ∘ round)(T/Δ)
L = .9
H = 1.1
x = ones( sims )
alive = trues( sims )
cache = similar(x)
for t ∈ 1:steps
x .+= Random.randn!(cache) .* σrΔ
@. alive &= L < x < H
end
mean(alive)end
main (generic function with 1 method)
julia> @time main()
18.912667 seconds (9 allocations: 157.811 KiB)
0.3624
```
you could also try to switch random generator to check.
---
### Author:  [@JADekker](https://discourse.julialang.org/u/JADekker)
#### Post date: [February 17, 2026, 10:28am UTC](https://discourse.julialang.org/t/simple-simulation/135701/4 "2026-02-17T10:28:53Z")
How fast you can go depends on what you want to record: if you only care about `mean(alive)`, you can get very fast (if I didn’t make a mistake…).
`sim_loop` simply uses a loop instead of a vector of normals. You can then strip this further by using some problem-specific insights if you only care about `mean(alive)`: 1) you can stop path-wise as soon as `alive` is false, 2) you don’t need to record `alive` as a vector, 3) instead of multiplying by a Float each time, you can also divide L and H by a Float once, which I would expect to be faster (but isn’t?)
```julia-auto
using Random, BenchmarkTools
function sim_opt(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x = ones(sims)
alive = trues(sims)
rnd = Vector{Float64}(undef, sims)
for t ∈ 1:steps
randn!(rnd)
@. x += rnd * σrΔ
@. alive &= (L < x < H)
end
mean(alive)
end
function sim_loop(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x = ones(sims)
alive = trues(sims)
Threads.@threads for i in 1:sims
for _ in 1:steps
x[i] += randn() * σrΔ
alive[i] &= (L < x[i] < H)
end
end
mean(alive)
end
function sim_loop_only_alive(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x = ones(sims)
alive = trues(sims)
Threads.@threads for i in 1:sims
for _ in 1:steps
x[i] += randn() * σrΔ
if !(L < x[i] < H)
alive[i] = false
break
end
end
end
mean(alive)
end
function sim_loop_only_alive_no_vecs(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x₀ = 1.
n_alive = sims
Threads.@threads for i in 1:sims
x = x₀
for _ in 1:steps
x += randn() * σrΔ
if !(L < x < H)
n_alive -= 1
break
end
end
end
n_alive / sims
end
function sim_loop_only_alive_no_vecs_no_mul(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L_sc = 0.9 / σrΔ
H_sc = 1.1 / σrΔ
x₀ = 1. / σrΔ
n_alive = sims
Threads.@threads for i in 1:sims
x = x₀
for _ in 1:steps
x += randn()
if !(L_sc < x < H_sc)
n_alive -= 1
break
end
end
end
n_alive / sims
end
function tests(N)
@time sim_opt(N);
@time sim_loop(N);
@time sim_loop_only_alive(N);
@time sim_loop_only_alive_no_vecs(N);
@time sim_loop_only_alive_no_vecs_no_mul(N);
end
tests(1000)
```
gives
```julia-auto
3.063867 seconds (1.00 M allocations: 4.053 GiB, 6.50% gc time)
1.090280 seconds (16.19 k allocations: 1.043 MiB, 6.25% compilation time)
0.722454 seconds (14.61 k allocations: 995.672 KiB, 10.36% compilation time)
0.699596 seconds (13.65 k allocations: 907.141 KiB, 10.81% compilation time)
0.731058 seconds (13.49 k allocations: 896.547 KiB, 10.69% compilation time)
0.38
```
---
### Author:  [@JADekker](https://discourse.julialang.org/u/JADekker)
#### Post date: [February 17, 2026, 10:35am UTC](https://discourse.julialang.org/t/simple-simulation/135701/5 "2026-02-17T10:35:49Z")
Note that if you only care about `alive`, a loop-based version will get relatively faster when you add more and more `sims`: a vectorized approach requires all calculations to be done for all paths, also for the ones that are no longer `alive`.
---
### Author:  [@Lincoln\_Hannah](https://discourse.julialang.org/u/Lincoln_Hannah)
#### Post date: [February 17, 2026, 10:38am UTC](https://discourse.julialang.org/t/simple-simulation/135701/6 "2026-02-17T10:38:34Z")
much appreciated.
could @turbo be used? or does the sequential calculation make it ineffective
---
### Author:  [@technocrat](https://discourse.julialang.org/u/technocrat)
#### Post date: [February 17, 2026, 10:41am UTC](https://discourse.julialang.org/t/simple-simulation/135701/7 "2026-02-17T10:41:12Z")
When dealing with unreasonably slow code, it’s best to consider the algorithm. The example takes 24 seconds on my well provisioned Mac.
For a Brownian motion with constant drift=0, σ, and absorbing barriers at L and H starting from x₀=1 (the midpoint), there are closed-form survival probabilities via eigenfunction expansion. That would replace the entire simulation. Your Δ=1e-6 is very small, likely to avoid missing barrier crossings between discrete steps. Brownian bridge correction lets you use much larger Δ (e.g., 1e-3) while still accounting for intra-step barrier hits, reducing steps from 1M to 1K:
```julia-auto
function barrier_sim_bb(; sims=10_000, T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
σrΔ = σ * √Δ
steps = round(Int, T / Δ)
x = ones(sims)
alive = trues(sims)
for t ∈ 1:steps
@inbounds for i ∈ 1:sims
alive[i] || continue
x_old = x[i]
x[i] += randn() * σrΔ
x_new = x[i]
# direct exit
if !(L < x_new < H)
alive[i] = false
continue
end
# Brownian bridge: probability of crossing barrier between steps
p_low = exp(-2(x_old - L) * (x_new - L) / (σ^2 * Δ))
p_high = exp(-2(H - x_old) * (H - x_new) / (σ^2 * Δ))
if rand() < p_low + p_high
alive[i] = false
end
end
end
mean(alive)
end
@time barrier_sim_bb(sims=10_000, T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
```
quite a bit faster
```julia-auto
julia> @time barrier_sim_bb(sims=10_000, T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
0.093072 seconds (5 allocations: 81.422 KiB)
0.3762
```
---
### Author:  [@JADekker](https://discourse.julialang.org/u/JADekker)
#### Post date: [February 17, 2026, 10:42am UTC](https://discourse.julialang.org/t/simple-simulation/135701/8 "2026-02-17T10:42:51Z")
I was thinking the same but figured that OP might have more complicated dynamics in mind in the long run. If not, directly simulating hitting times of fixed barriers is also an approach: first simulate the first hitting time of either barrier; then simulate the state `x` at the first of {hitting time, T} using a conditional distribution. If you also need the terminal state and hitting occurred before time T, simulate x\_T using one normal random variate.
---
### Author:  [@Lincoln\_Hannah](https://discourse.julialang.org/u/Lincoln_Hannah)
#### Post date: [February 17, 2026, 10:45am UTC](https://discourse.julialang.org/t/simple-simulation/135701/9 "2026-02-17T10:45:56Z")
thank you.
I’ll combine this with JADekker’s approach.
---
### Author:  [@JADekker](https://discourse.julialang.org/u/JADekker)
#### Post date: [February 17, 2026, 11:00am UTC](https://discourse.julialang.org/t/simple-simulation/135701/10 "2026-02-17T11:00:07Z")
Combining, I get
```julia-auto
function barrier_sim_bb(sims=10_000; T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
σrΔ = σ * √Δ
steps = round(Int, T / Δ)
x = ones(sims)
alive = trues(sims)
for t ∈ 1:steps
@inbounds for i ∈ 1:sims
alive[i] || continue
x_old = x[i]
x[i] += randn() * σrΔ
x_new = x[i]
# direct exit
if !(L < x_new < H)
alive[i] = false
continue
end
# Brownian bridge: probability of crossing barrier between steps
p_low = exp(-2(x_old - L) * (x_new - L) / (σ^2 * Δ))
p_high = exp(-2(H - x_old) * (H - x_new) / (σ^2 * Δ))
if rand() < p_low + p_high
alive[i] = false
end
end
end
mean(alive)
end
function barrier_sim_bb_quick(sims=10_000; T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
σrΔ = σ * √Δ
steps = round(Int, T / Δ)
x = 1.
n_alive = sims
fr_inv = 2 / (σ^2 * Δ)
for i ∈ 1:sims
x_old = x
x_new = x_old
for t ∈ 1:steps
x_old = x_new
x_new += randn() * σrΔ
# direct exit
if !(L < x_new < H)
n_alive -= 1
break
end
# Brownian bridge: probability of crossing barrier between steps
p_low = exp(-(x_old - L) * (x_new - L) * fr_inv)
p_high = exp(-(H - x_old) * (H - x_new) * fr_inv)
if rand() < p_low + p_high
n_alive -= 1
break
end
end
end
n_alive / sims
end
```
gives with 1000 simulations
```julia-auto
0.010199 seconds (3 allocations: 8.219 KiB) # barrier_sim_bb
0.006918 seconds # barrier_sim_bb_quick
```
and 10.000 simulations
```julia-auto
0.126836 seconds (4 allocations: 79.531 KiB)
0.070636 seconds
```