Julia による Lévy flight のプロット

levyflight.jl

#=
Lévy flight
Copyright (c) 2021 Daichi Furukawa, github.com/Loliver1224
License: MIT
=#

using SpecialFunctions  # gamma()
using Plots

N = 200  # max iter
Dim = 2  # dimension

struct LevyFlight
    β
    σ
    LevyFlight(; β=1.5) = new(
        β,
        (gamma(1 + β) * sin(π * β / 2) / (gamma((1 + β) / 2) * β * 2 ^ ((β-1)/2))) ^ (1/β)
    )
end

function Levy(Dim, β, σ)
    u = randn(1, Dim) * σ
    v = randn(1, Dim)

    u ./ abs.(v) .^ (1/β)
end

function step(flight::LevyFlight, Dim)
    θ = 2π * rand()
    x, y = [cos(θ) sin(θ)] .* Levy(Dim, flight.β, flight.σ)

    append!(xs, xs[end] + x)
    append!(ys, ys[end] + y)
end

# --- Main ---
xs = [randn()]
ys = [randn()]

flight = LevyFlight()
for i=1:N
    step(flight, Dim)
end

plot(
    xs, ys,
    xlabel="x",
    ylabel="y",
    title="Lévy flight",
    legend=false
)