I’ve implemented a (slightly adapted) julia version of the Generalized Pareto estimation method described in this paper: DOI 10.1198/TECH.2009.08017

Am I missing any obvious performance improvements? Is it correct that this is not a good use case for StaticArrays, as the size of the small arrays is not known at compile time?

```
using Statistics
function fit_gpd(x::Vector{T}) where T <: Real
x = sort(x)
x_star = x[floor(Int, length(x) / 4.0 + 0.5)]
n = length(x)
m::Int = 30 + floor(sqrt(n))
θ = [1 / x[n] + (1 - sqrt(m / (j - 0.5))) / (3 * x_star) for j in 1:m]
lx(a, x) = let
a = -a
k = [mean(log1p.(b .* x)) for b in a]
@. log(a / k) - k - 1
end
l_θ = n * lx(θ, x)
w_θ = 1.0 ./ [sum(exp.(l_θ .- l_θ[j])) for j in 1:m]
θ_hat = sum(θ .* w_θ)
ξ = mean(log1p.(-θ_hat .* x))
σ = -ξ ./ θ_hat
ξ = ξ * n / (n + 10.0) + 10.0 * 0.5 / (n + 10.0)
return σ, ξ
end
s = rand(10000)
fit_gpd(s)
```

Thank you for your help