Press "Enter" to skip to content

10 search results for "julia"

Parallel computing with Julia

The Julia Language

Julia is particularly attractive for its ability to run in parallel, not only on clusters, but also on the multiple cores of your laptop processor. Here is a minimalistic example using pmap().

# For parallel computing with julia -p NumberOfAdditionalCPU

@everywhere begin # this part will be available on all CPUs

function computation(runid_is_useless)
  result = rand()
  return result
end # function computation()

end # @everywhere

### Main part - runs only on main instance of Julia.

Nprocs = nprocs()
print("Number of processors ",Nprocs,".\n")

## Launching computations on Nprocs parallel processors
results = @time pmap(computation,1:Nprocs)

## Post-processing
for i in 1:Nprocs
  print(@sprintf("Result of CPU %i is %f\n",i,results[i]))
end # for i

On my laptop, my processor has 4 cores, and I run this code with julia -p 3 code.jl which gives

Number of processors 4.
0.548238 seconds (214.75 k allocations: 9.788 MB, 5.76% gc time)
Result of CPU 1 is 0.725449
Result of CPU 2 is 0.377958
Result of CPU 3 is 0.443606
Result of CPU 4 is 0.016641

The same code with -p 39 can also be run on the 40 CPUs cluster of my department! For more material on parallel computing with Julia, you may take a look at Julia documentation.

Leave a Comment

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
2 Comments

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

The Julia Language

The Julia Language

I played a bit this spring with Julia (the software, not the girl next door ;-). For Debian : apt-get install julia julia-doc. This was the source of great enthusiasm! Julia is a promising software for scientific computing. Yet another interpreter around Lapack? Yes and no. Some people believe that Julia is a revolution in scientific computing. For most of my daily needs, Julia is already almost preferable to Matlab/Octave/Scilab/R/C++/Python/… due to multiple bleeding edge technical features, including a high performance just-in-time compiler. Julia is flexible like Matlab, is fast like C, is free software like Octave/Scilab/R (MIT license), is text friendly like Python, and much more. Julia is also good for its competitors: it will help them to evolve. The reading of «Why we created Julia» by the creators of Julia is informative. An impressive summary of Julia features is given on julialang.org, which provides many resources, including videos featuring the great Alan Edelman (MIT group leader, co-creator of Julia).

Learning Julia is a pretty easy task. Here is one of my first Julia program:

Pkg.add("Winston") # Pkg.add("Winston") is needed once for all.
using Winston # for simple graphics à la Matlab
n = 1000;
(D,U) = eig((randn(n,n)+im*randn(n,n))/sqrt(2*n));
I = [-1:.01:1]; J = sqrt(1-I.^2);
hold(true)
plot(real(D),imag(D),"b.",I,J,"r--",I,-J,"r--")
title(@sprintf("Complex Ginibre Ensemble %dx%d",n,n))
file("circ.png")

(this is one of my standard tests: Ginibre ensemble!) which produces the following image: Ginibre with JuliaJulia is still young – born in 2012, version 0.2 – and needs time to gain maturity, notably for plotting support among other things. One can regret the lack of a graphical user interface, but one can use Emacs+ESS, or the wonderful IPython notebooks in web browser via the Julia package IJulia. For Debian: apt-get install ess ipython-notebook. Give it a try!

 

4 Comments

Can't find what you're looking for? Try refining your search:

Syntax · Style · .