{"id":10177,"date":"2018-02-18T20:44:18","date_gmt":"2018-02-18T19:44:18","guid":{"rendered":"http:\/\/djalil.chafai.net\/blog\/?p=10177"},"modified":"2018-06-28T15:53:08","modified_gmt":"2018-06-28T13:53:08","slug":"gumbel-fit-with-julia","status":"publish","type":"post","link":"https:\/\/djalil.chafai.net\/blog\/2018\/02\/18\/gumbel-fit-with-julia\/","title":{"rendered":"Gumbel fit with Julia"},"content":{"rendered":"

 <\/p>\n

\"Gumbel<\/a>
Gumbel fit example<\/figcaption><\/figure>\n

The\u00a0Gumbel distribution<\/a> appears in the modelling of extreme values phenomena and the famous\u00a0Fisher\u2013Tippett\u2013Gnedenko theorem<\/a>. Its cumulative distribution function is given by<\/p>\n

$$x\\in\\mathbb{R}\\mapsto\\exp\\left(-\\exp\\left(-\\frac{x-\\mu}{\\sigma}\\right)\\right)$$<\/p>\n

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

$$\\mu+\\gamma\\sigma\\quad\\text{and}\\quad\\frac{\\pi^2\\sigma^2}{6},$$<\/p>\n

where $\\gamma=0.577\\ldots$ is\u00a0Euler's constant<\/a>. The maximum likelihood estimators<\/a> $\\hat{\\mu}$ and $\\hat{\\sigma}$ satisfy<\/p>\n

$$
\n\\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}}}.
\n$$<\/p>\n

Note that the second formula is implicit. Unfortunately, the current version of the Julia<\/a> package\u00a0Distributions.jl<\/a>\u00a0does 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<\/a> implements a Gumbel fit in the fitdistrplus<\/a> package with a notable contribution from my colleague Christophe Dutang<\/a>.<\/p>\n

# Pkg.add(\"Roots\") # https:\/\/github.com\/JuliaMath\/Roots.jl\r\nusing Roots\r\n\"\"\" \r\nComputes approximately the Maximum Likelihood Estimator (MLE) of \r\nthe position and scale parameters of a Gumbel distribution with\r\ncumulative distribution function x-> exp(-exp(-(x-mu)\/sigma)).\r\nHas mean mu+eulergamma*sigma and variance sigma^2*pi^2\/6.\r\nThe MLE for mu is a function of data and the MLE for sigma.\r\nThe MLE for sigma is the root of a nonlinear function of data,\r\ncomputed numerically with function fzero() from package Roots.\r\n\r\nReference: page 24 of book ISBN 1-86094-224-5.\r\nSamuel Kotz and Saralees Nadarajah\r\nExtreme value distributions.Theory and applications. \r\nImperial College Press, London, 2000. viii+187 pp. \r\nhttps:\/\/ams.org\/mathscinet-getitem?mr=1892574 \r\n\r\n``no distribution should be stated without an explanation of how the\r\nparameters are estimated even at the risk that the methods used will \r\nnot stand up to the present rigorous requirements of mathematically \r\nminded statisticians''. E. J. Gumbel (1958).\r\n\"\"\"\r\nfunction gumbel_fit(data)\r\n  f(x) = x - mean(data) + sum(data.*exp(-data\/x)) \/ sum(exp(-data\/x)) \r\n  sigma = fzero(f,sqrt(6)*std(data)\/pi) # moment estimator as initial value\r\n  mu = - sigma * log(sum(exp(-data\/sigma))\/length(data))\r\n  return mu , sigma\r\nend #function\r\n<\/pre>\n
# Pkg.add(\"Plots\") # https:\/\/github.com\/JuliaPlots\/Plots.jl\r\n# Pkg.add(\"PyPlot\") # https:\/\/github.com\/JuliaPy\/PyPlot.jl\r\n# http:\/\/docs.juliaplots.org\/\r\nusing Plots\r\nfunction gumbel_fit_plots(data,label,filename) \r\n  n = length(data)\r\n  data = reshape(data,n,1)\r\n  mu , sigma = gumbel_fit(data)\r\n  X = collect(linspace(minimum(data),maximum(data),1000))\r\n  tmp = (X-mu)\/sigma\r\n  Y = exp(-tmp-exp(-tmp))\/sigma\r\n  pyplot()\r\n  histogram!(data,\r\n             nbins = round(Int,sqrt(n)),\r\n             normed = true,\r\n             label = @sprintf(\"%s histogram n=%i\",label,n),\r\n             color = :white)\r\n  plot!(X,Y,\r\n        title = @sprintf(\"Gumbel density mu=%2.2f sigma=%2.2f\",mu,sigma),\r\n        titlefont = font(\"Times\", 10),\r\n        label = \"Gumbel density\",\r\n        lw = 3,\r\n        linestyle = :solid, \r\n        linecolor = :darkblue,\r\n        grid = false,\r\n        border = false)\r\n        gui()\r\n  #savefig(string(filename,\".svg\"))\r\n  savefig(string(filename,\".png\"))\r\nend #function\r\n<\/pre>\n
function gumbel_fit_example()\r\n  gumbel_fit_plots(-log(-log(rand(1,1000))),\"Sample\",\"gumbelfitexample\") \r\nend #function<\/pre>\n","protected":false},"excerpt":{"rendered":"
  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}$…<\/p>\n
Continue readingGumbel fit with Julia<\/span><\/a><\/div>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"iawp_total_views":282},"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/10177"}],"collection":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/comments?post=10177"}],"version-history":[{"count":37,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/10177\/revisions"}],"predecessor-version":[{"id":10413,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/10177\/revisions\/10413"}],"wp:attachment":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/media?parent=10177"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/categories?post=10177"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/tags?post=10177"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}