Press "Enter" to skip to content

Month: February 2018

Gumbel fit with Julia

 

Gumbel fit example
Gumbel fit example

The Gumbel distribution appears in the modelling of extreme values phenomena and the famous Fisher–Tippett–Gnedenko theorem. Its cumulative distribution function is given by

$$x\in\mathbb{R}\mapsto\exp\left(-\exp\left(-\frac{x-\mu}{\sigma}\right)\right)$$

where $\mu\in\mathbb{R}$ and $\sigma\in(0,\infty)$ are location and scale parameters. Its mean and variance are

$$\mu+\gamma\sigma\quad\text{and}\quad\frac{\pi^2\sigma^2}{6},$$

where $\gamma=0.577\ldots$ is Euler’s constant. The maximum likelihood estimators $\hat{\mu}$ and $\hat{\sigma}$ satisfy

$$
\hat{\mu} =-\hat{\sigma}\log\left(\frac{1}{n}\sum_{i=1}^{n}\mathrm{e}^{-x_i/\hat{\sigma}}\right)\quad\text{and}\quad\hat{\sigma} = \frac{\sum_{i=1}^nx_i}{n}-\frac{\sum_{i=1}^{n} x_i\mathrm{e}^{-x_i/\hat{\sigma}}}{ \sum_{i=1}^{n} \mathrm{e}^{-x_i/\hat{\sigma}}}.
$$

Note that the second formula is implicit. Unfortunately, the current version of the Julia package Distributions.jl does not implement yet the Gumbel distro in fit()/fit_mle(). Here is a home brewed program to fit a Gumbel distribution on data using maximum likelihood and moments. The produced graphic is above. Note that the GNU R software implements a Gumbel fit in the fitdistrplus package with a notable contribution from my colleague Christophe Dutang.

# Pkg.add("Roots") # https://github.com/JuliaMath/Roots.jl
using Roots
""" 
Computes approximately the Maximum Likelihood Estimator (MLE) of 
the position and scale parameters of a Gumbel distribution with
cumulative distribution function x-> exp(-exp(-(x-mu)/sigma)).
Has mean mu+eulergamma*sigma and variance sigma^2*pi^2/6.
The MLE for mu is a function of data and the MLE for sigma.
The MLE for sigma is the root of a nonlinear function of data,
computed numerically with function fzero() from package Roots.

Reference: page 24 of book ISBN 1-86094-224-5.
Samuel Kotz and Saralees Nadarajah
Extreme value distributions.Theory and applications. 
Imperial College Press, London, 2000. viii+187 pp. 
https://ams.org/mathscinet-getitem?mr=1892574 

``no distribution should be stated without an explanation of how the
parameters are estimated even at the risk that the methods used will 
not stand up to the present rigorous requirements of mathematically 
minded statisticians''. E. J. Gumbel (1958).
"""
function gumbel_fit(data)
  f(x) = x - mean(data) + sum(data.*exp(-data/x)) / sum(exp(-data/x)) 
  sigma = fzero(f,sqrt(6)*std(data)/pi) # moment estimator as initial value
  mu = - sigma * log(sum(exp(-data/sigma))/length(data))
  return mu , sigma
end #function
# Pkg.add("Plots") # https://github.com/JuliaPlots/Plots.jl
# Pkg.add("PyPlot") # https://github.com/JuliaPy/PyPlot.jl
# http://docs.juliaplots.org/
using Plots
function gumbel_fit_plots(data,label,filename) 
  n = length(data)
  data = reshape(data,n,1)
  mu , sigma = gumbel_fit(data)
  X = collect(linspace(minimum(data),maximum(data),1000))
  tmp = (X-mu)/sigma
  Y = exp(-tmp-exp(-tmp))/sigma
  pyplot()
  histogram!(data,
             nbins = round(Int,sqrt(n)),
             normed = true,
             label = @sprintf("%s histogram n=%i",label,n),
             color = :white)
  plot!(X,Y,
        title = @sprintf("Gumbel density mu=%2.2f sigma=%2.2f",mu,sigma),
        titlefont = font("Times", 10),
        label = "Gumbel density",
        lw = 3,
        linestyle = :solid, 
        linecolor = :darkblue,
        grid = false,
        border = false)
        gui()
  #savefig(string(filename,".svg"))
  savefig(string(filename,".png"))
end #function
function gumbel_fit_example()
  gumbel_fit_plots(-log(-log(rand(1,1000))),"Sample","gumbelfitexample") 
end #function
1 Comment

Playing a bit with Julia

The Julia Language
This short post is devoted to a couple of Julia programs.

Gaussian Unitary Ensemble. It is the Boltzmann-Gibbs measure with density proportional to
$$
(x_1,\ldots,x_N)\in\mathbb{R}^N\mapsto\mathrm{e}^{-N\sum_{i=1}^N x_i^2}\prod_{i<j}(x_i-x_j)^2.
$$
It is the distribution of the eigenvalues of a random $N\times N$ Hermitian matrix distributed according to the Gaussian probability measure with density proportional to
$$
H\mapsto\mathrm{e}^{-N\mathrm{Tr}(H^2)}.
$$
It is a famous exactly solvable model of mathematical physics. There is a nice formula for the mean empirical spectral distribution $\mathbb{E}\mu_N$ where $\mu_N=\frac{1}{N}\sum_{i=1}^N\delta_{x_i}$, namely
$$
x\in\mathrm{R}\mapsto\frac{\mathrm{e}^{-\frac{N}{2}x^2}}{\sqrt{2\pi N}}\sum_{\ell=0}^{N-1}P_\ell^2(\sqrt{N}x)
$$
where ${(P_\ell)}_{\ell\geq0}$ are the Hermite polynomials which are the orthonormal polynomials for the standard Gaussian distribution $\mathcal{N}(0,1)$. The computation of the Laplace transform and a subtle reasonning, see this paper, reveal that it converges as $N\to\infty$ towards the the Wigner semicircle distribution with density with respect to the Lebesgue measure
$$
x\in\mathbb{R}\mapsto\frac{\sqrt{4-x^2}}{2\pi}\mathbf{1}_{x\in[-2,2]}.
$$
Here is a nice plot followed by the Julia code used to produce it.

function normalized_hermite_polynomials_numeric(n,x)
    if size(x,2) != 1
        error("Second argument must be column vector.")
    elseif n == 0
        return ones(length(x),1)
    elseif n == 1
        return x
    else
        P = [ ones(length(x),1) x ]
        for k = 1:n-1
            P = [ P  x.*P[:,k+1]/sqrt(k+1)-P[:,k]*sqrt(1/(1+1/k)) ]
        end
        return P # matrix with columns P_0(x),...,P_n(x)
    end
end #function

function gue_numeric(n,x)
    if size(x,2) != 1
        error("Second argument must be column vector.")
    else
        y = sqrt(n)*x
        p = normalized_hermite_polynomials_numeric(n-1,y).^2
        return exp(-y.^2/2) .* sum(p,2) / sqrt(2*pi*n)
    end
end #function

# Pkg.add("Plots") # http://docs.juliaplots.org/latest/tutorial/
using Plots
pyplot()
N = 8
X = collect(-3:.01:3)
Y = gue_numeric(N,X)
function semicircle(x)
    if (abs(x)>2) return 0 else return sqrt(4-x^2)/(2*pi) end
end # function
Y = [Y [semicircle(x) for x in X]]
plot(X,Y,
     title = @sprintf("GUE mean Empirical Spectral Distribution beta=2 N=%i",N),
     titlefont = font("Times", 10),
     label = ["GUE mean ESD" "SemiCircle"],
     lw = 2,
     linestyle = [:dash :solid],     
     linecolor = [:darkblue :darkred],
     aspect_ratio = 2*pi,
     grid = false,
     border = false)
gui()
#savefig("gue.svg")
savefig("gue.png")

Here is, just for fun, a way to produce symbolically Hermite polynomials.

# Pkg.add("Polynomials") # http://juliamath.github.io/Polynomials.jl/latest/
using Polynomials

function unnormalized_hermite_polynomial_symbolic(n) 
    if n == 0
        return Poly([1])
    elseif n == 1
        return Poly([0,1])
    else
        P = [Poly([1]) Poly([0,1])]
        for k = 1:n-1
            P = [ P Poly([0,1])*P[k+1]-k*P[k] ]
        end
        return P # better than just P[n+1]
    end
end #function

Complex Ginibre Ensemble. It is the Boltzmann-Gibbs measure with density proportional to
$$
(x_1,\ldots,x_N)\in\mathbb{C}^N\mapsto\mathrm{e}^{-N\sum_{i=1}^N|x_i|^2}\prod_{i<j}|x_i-x_j|^2.
$$
It is the distribution of the eigenvalues of a random $N\times N$ complex matrix distributed according to the Gaussian probability measure with density proportional to
$$
M\mapsto\mathrm{e}^{-N\mathrm{Tr}(MM^*)}
$$
where $M^*=\overline{M}^\top$. Yet another exactly solvable model of mathematical physics, with a nice formula for the mean empirical spectral distribution $\mathbb{E}\mu_N$ where $\mu_N=\frac{1}{N}\sum_{i=1}^N\delta_{x_i}$, given by
$$
z\in\mathbb{C}\mapsto\frac{\mathrm{e}^{-N|z|^2}}{\pi}\sum_{\ell=0}^{N-1}\frac{|\sqrt{N}z|^{2\ell}}{\ell!}
$$
which is the analogue of the one given above for the Gaussian Unitary Ensemble. Using induction and integration by parts, it turns out that this density can be rewritten as
$$
z\in\mathbb{C}\mapsto\int_{N|z|^2}^\infty\frac{u^{N-1}\mathrm{e}^{-u}}{\pi(N-1)!}\mathrm{d}u
=\frac{\Gamma(N,N|z|^2)}{\pi}
$$
where $\Gamma$ is the normalized incomplete Gamma function and where we used the identity
$$
\mathrm{e}^{-r}\sum_{\ell=0}^n\frac{r^\ell}{\ell!}=\frac{1}{n!}\int_r^\infty u^n\mathrm{e}^{-u}\mathrm{d}u.
$$
Note that the function $t\mapsto 1-\Gamma(N,t)$ is the cumulative distribution function of the Gamma distribution with shape parameter $N$ and scale parameter $1$. As a consequence, denoting $X_1,\ldots,X_N$ i.i.d. exponential random variables of unit mean, the density can be written as
$$
z\in\mathbb{C}\mapsto\frac{1}{\pi}\mathbb{P}\left(\frac{X_1+\cdots+X_N}{N}\geq|z|^2\right).
$$
Now the law of large numbers implies that this density converges pointwise as $N\to\infty$ towards the density of the uniform distribution on the unit disk of the complex plane
$$
z\in\mathbb{C}\mapsto\frac{\mathbf{1}_{|z|\leq1}}{\pi}.
$$
One can also use i.i.d. Poisson random variables $Y_1,\ldots,Y_N$ of mean $|z|^2$, giving
$$
z\in\mathbb{C}\mapsto\frac{1}{\pi}\mathbb{P}\left(\frac{Y_1+\cdots+Y_N}{N}<1\right),
$$
see this former post. Here are the plots with respect to $|z|$, followed by the Julia code.

# Pkg.add("Distributions") # https://juliastats.github.io/Distributions.jl
# Pkg.add("Plots") # http://docs.juliaplots.org/latest/tutorial/
using Plots
using Distributions 
N = 8
X = collect(0:.01:3)
Y = ccdf(Gamma(N,1),N*X.^2)/pi
function unit(x)
    if (x>1) return 0 else return 1/pi end
end # function
Y = [Y [unit(x) for x in X]]
pyplot()
plot(X,Y,
     title = @sprintf("Ginibre mean Empirical Spectral Distribution beta=2 N=%i",N),
     titlefont = font("Times", 10),
     label = ["Ginibre mean ESD" "Equilibrium"],
     lw = 2,
     linestyle = [:dash :solid],     
     linecolor = [:darkblue :darkred],
     grid = false,
     border = false)
gui()
#savefig("ginibre.svg")
savefig("ginibre.png")

Further reading. The Introducing Julia wikibook on Wikibooks.

Leave a Comment
Syntax · Style · Tracking & Privacy.