Giter Club home page Giter Club logo

hill_simu's Introduction

Bachelor thesis and simulations

This repository was for my bachelor thesis and simulations for the thesis. My thesis is related to the field of extreme value theory. In the kandi folder there are all the files regarding the thesis itself:

  • the LaTeX source code kandi.tex.
  • thesis as a pdf file kandi.pdf.
  • figures pareto.pdf and cauchy.pdf.

In the simu folder there are all the files concerning the simulations. Purpose of the simulations is to generate independent samples and calculate values for Hill estimator. Hill estimator estimates parameter called extreme value index, that describes the tail behaviour of a distribution. Simulated data is used for two figures which show finite sample behaviour of the Hill estimator. Simulations were executed with the Aalto University's Triton computing cluster.

Here are all the files used for simulations:

  • hillSim.py contains the simulation code. It is run in fashion python hillSim.py [arguments], where [arguments] specify different simulation scenarios.

  • gen_args.py is used for creating command line arguments of form python hillSim.py [arguments].

  • simujob.py submits the job to the Triton and starts the execution.

Following files are for plotting the figures in my thesis:

  • hillPlot.py contains code for plotting and is run in fashion python hillPlot.py [arguments].

  • plot.sh executes hillPlot.py with various arguments and produces the figures. Figures appear in the kandi folder.

Running the simulations

Run

python gen_args.py

This generates text files containing command line arguments in the folder args.

After this submit the job by running

sbatch simujob.sh

Results will appear in the folder data.

Huge thanks for my thesis advisor Matias Heikkilä for all the advice both for mathematics and coding. Simulations follow very closely the same format as Matias Heikkilä's code in repository https://github.com/mapehe/evt-ica-simu.

hill_simu's People

Contributors

perej1 avatar mapehe avatar

Watchers

 avatar

Forkers

mapehe

hill_simu's Issues

Runtime error

Moi, yritin ajaa python hillSim.py, ja scripti tuotti seuraavan tulosteen

hillSim.py:43: RuntimeWarning: divide by zero encountered in power
return np.power(estimate,-1)
Traceback (most recent call last):
File "hillSim.py", line 184, in
main()
File "hillSim.py", line 155, in main
plt.hist(N, 50)
File "/usr/local/lib/python2.7/dist-packages/matplotlib/pyplot.py", line 3132, in hist
stacked=stacked, normed=normed, data=data, **kwargs)
File "/usr/local/lib/python2.7/dist-packages/matplotlib/init.py", line 1855, in inner
return func(ax, *args, **kwargs)
File "/usr/local/lib/python2.7/dist-packages/matplotlib/axes/_axes.py", line 6530, in hist
m, bins = np.histogram(x[i], bins, weights=w[i], **hist_kwargs)
File "/usr/local/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 670, in histogram
'range parameter must be finite.')
ValueError: range parameter must be finite.

Plotit

Tein tuossa nyt ensimmäiset plotit. Ajattelin, että kuvasta voi tehdä kauniimman kun itse skripti on kunnossa. Itselleni ongelmaksi muodostui miten k(n) kannattaisi antaa parametrina kun se on funktio n:stä. Nyt se on vaan määritelty erikseen for-loopissa, mikä ei vissiin ole kovin järkevää. Toiseksi noita argumentteja on aika paljon. Niitä pitäisi osata ehkä vähentää jollain tavalla. Mutta tosiaan seuraavanlaisella komennolla sain koodin ajettua: "python hillPlot.py -n1 100 -n2 30000 -s 100 -f1 data -r 500 -p pplot -g 0.33 -d pareto -f2 convergence.png". Tuon p parametrin ideana on määrittää mitä plotataan siinä vaiheessa kun vaihtoehtoja on enemmän kuin yksi.

Toinen versio: Datan generointi csv tiedostoon

Tein toisen version noiden korjausehdotusten perusteella. Nyt en plottaile vielä mitään vaan generoin vain estimaatteja parametrien n (otoskoko) ja r (estimaattien määrä) mukaan. Parametrit syötetetään vasta ajettaessa komentorivilla tyyliin "python3 hillSim.py -n 1000 -r 500" argparse kirjaston avulla. Tuloksena syntyy csv tiedosto, jossa on hill estimaatorin arvoja pareto ja cauchy jakaumille. Skripti on vielä varsin alkeellinen mutta ajattelin nyt lähteä koodaamaan vähän rauhallisemmin pala kerrallaan (toisin kuin ensimmäisessä versiossa).

Mean of empty slice

Kun koitan ajaa python hillSim.py -n 2000 -r 100, saan seuraavan tulosteen

/usr/local/lib/python2.7/dist-packages/numpy/core/fromnumeric.py:2957: RuntimeWarning: Mean of empty slice.
out=out, **kwargs)
/usr/local/lib/python2.7/dist-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars
ret = ret.dtype.type(ret / rcount)

Luulen että tuossa

return np.mean(sample[1:])

lasketaan jostain syystä keskiarvoa tyhjän joukon yli. Tsekkaahan vielä että skripti toimii sulla. Jos toimii sulla mutta ei mulla niin lähdetään yhdessä metsästämään bugia.

Muutamia korjausehdotuksia

Laitan tähän joitain parannusehdotuksia.

Ominaisuudet

Tällä hetkellä jako näyttäisi siltä että koodi

  • Tekee "hill-plotin"
  • Simuloi stokastista suppenemista
  • Simuloi jakaumasuppenemista
  • Plottaa tuloksia

Ehdotan että lähtisit mielummin liikkeelle tällaisella jaolla:

  • Simuloi hill-estimaattorin arvoja
  • Plottaa simuloidut arvot

Tuosta ensimmäisestä pointista voisi sitten varsin helposti saada nuo molemmat suppenemismoodit tarkasteltua. Hill-plot sisältyisi mahdollisesti tuohon jälkimmäiseen. Lähde siitä että toteuta skripti mikä vain tuottaa Hill estimaattorin arvoja näiden

hill_simu/hillSim.py

Lines 53 to 55 in 64f3caf

b = 3
n = 10000
r= 5000

vakioiden määrämällä tavalla ja sylkee tuotetut arvot vaikka vaan csv-tiedostoon. Katsotaan tuota plottailua sitten kun se generointijuttu toimii solidisti.

Modulaarisuus (ei tärkeä tässä vaiheessa)

Tällä hetkellä skriptiin on rakennettu sisään kaikenlaista toiminnallisuutta mikä kyllä liittyy läheisesti siihen mikä meitä kiinnostaa. Tämä pitäisi nyt kuitenkin organisoida siten, että skriptin eri vaiheita voisi käyttää toisistaan riippumatta.

Pari esimerkkiä: Tässä vissiin tehdään "hill plot"

hill_simu/hillSim.py

Lines 60 to 71 in 64f3caf

#Creating everithing necessary for the hill plot.
est = np.array([])
K = np.arange(1,int(n/3))
# Estimates for different values for k.
for k in K:
est=np.append(est, hill(k, X))
# Line y=b represents the true value for tail index
B = b*np.ones(K.shape)

toisaalta tässä

hill_simu/hillSim.py

Lines 79 to 83 in 64f3caf

nest = np.array([])
nvec = np.arange(10, n+1)
for i in nvec:
nest = np.append(nest, hill(kf(i), pareto.rvs(b, size=i)))

tarkastellaan estimaattorin asymptoottisia ominaiisuuksia. Nämä ovat selkeästi erillisiä juttuja emmekä välttämättä halua aina tehdä molempia. Loppuvaiheessa haluamme että simulointia ja plottailua pystyy tekemään erikseen, mikä onnistuu varsin helposti argparse-kirjastolla.

Funktio hill

hill_simu/hillSim.py

Lines 28 to 43 in 64f3caf

def hill(k, data):
# Data is sorted into descending order.
sorted_data = np.sort(data)[::-1]
# k largest elements.
kdata = sorted_data[:k]
# Elementwise natural logarithm.
logdata = np.log(kdata)
# The smallest element is subtracted from every element of this list.
sumdata = logdata - logdata[-1]
# Sum of elements and normalize with 1/k.
# This is estimate for extreme value index.
estimate = sum(sumdata)/k
# estimate for tail index
return np.power(estimate,-1)

Voi olla että tässä on jopa virhe (ainakin tässä #1 näyttäisi siltä että yritetään jakaa nollalla). Joka tapauksessa tuossa lienee ihan oikea idea, mutta sen voisi ilmaista simppelimmin, voit kopsata suoraan tämän pätkän mitä käytin jossain muualla

def hill(sample, k):
    sample = sorted(sorted(np.array(sample), reverse=True)[:k])
    sample = np.log(sample)-np.log(sample[0])
    return np.mean(sample[1:])

Tämä nyt tuottaa alfan sijaan gamman, mutta voidaan ajatella että ollaan kiinnostuneita ääriarvoindeksistä eikä häntäindeksistä.

Sisäänrakennetut vakiot

hill_simu/hillSim.py

Lines 53 to 55 in 64f3caf

b = 3
n = 10000
r= 5000

Näistä sisäänrakennetuista vakioista pitäisi päästä kokonaan eroon. Haluttaisiin, että python hillSim.py muotoisen komennon sijaan simukierros ajettaisiin komennolla python hillSim.py -b 3 -n 10000 -r 5000, jotta noita vakioita olisi sitten joustavaa säätää halutulla tavalla koskematta itse skriptiin. Tähän soveltuu argparse-kirjasto. Copypasteet tuolta vaan jonkin esimerkin ja muokkaat siitä sellaisen mitä tarvitset tähän tilanteeseen.

Tyyliseikkoja

Kommenteista

hill_simu/hillSim.py

Lines 15 to 18 in 64f3caf

# Own definition for k, which is a parameter in Hill estimator
# kf(n) is such a function that k-> inf and k/n -> 0
# k is now quite large because we want to see some bias
# for small k variance is great but bias would be nonexistent.

# Tämä on ihan hyvä syntaksi python kommenteille. Tehdään me tällaiset block-commentit kuitenkin näin:

"""
Own definition for k, which is a parameter in Hill estimator
kf(n) is such a function that k-> inf and k/n -> 0
k is now quite large because we want to see some bias
for small k variance is great but bias would be nonexistent.
"""

PEP8

Hyvä tapa pitää huolta koodin luettavuudesta on noudattaa PEP8 standardia. Järkevin tapa pitää huolta siitä että koodi on standarin mukaista on ajaa flake8 hillSim.py ja tehdä muutoksia niin kauan ettei komento tulosta mitään. On myös olemassa online-työkaluja tähän.

Korjausehdotuksia tekstiin

Hyvin olit saanut tekstiä aikaiseksi. Laitoin alle joitain juttuja, mitkä voisit katsoa kuntoon ennen kuin tapaatte Pauliinan kanssa.

Asiavirheet ja täsmennykset

  • Tämä lause kaipaa täsmennystä.

    hill_simu/kandi/kandi.tex

    Lines 454 to 467 in 85b028d

    \begin{theorem}
    There exists real constants $a_n>0$ and $b_n \in \mathbb{R}$ such that
    \begin{equation}
    \lim_{n\to\infty} F^n(a_nx + b_n) = G_{\gamma}(x) =
    \begin{cases}
    exp(-(1 + \gamma x)^{-\frac{1}{\gamma}}), \gamma \neq 0 \\
    exp(-e^{-x}), \gamma = 0,
    \end{cases}
    \label{mdaEq}
    \end{equation}
    for all x with $1+\gamma x > 0$ where $\gamma \in \mathbb{R}$.
    \end{theorem}

    Huomaathan, että suppenemisen F(a_n x + b_n) -> G(x) sijaan tulisi todeta, että F(a_n x + b_n) -> G(ax + b) (nykyisellään jakauman G sisältä puuttuu skaalauskertoimet a ja b). Vertaa tätä kirjaan ja muotoile tämä uudestaan siten, että korjaus sopii siihen.

  • Lähetin tähän

    If $Y_1, Y_2,...$ are i.d.d random variables from Pareto distribution with cdf $F_Y(y) = 1 - \frac{1}{y}$ then $U(Y_i) \overset{d}{=} X_i$, since

    liittyen sähköpostiin täsmennyksen josta puhuimme kun tavattiin viimeksi. Ehdotan, että erotat luettavuuden nimissä koko tämän väitteen omaksi lemmakseen, joka esitetään ennen tätä todistusta.

Tyyliseikkoja

  • Pitäsikö tänne

    \subsection*{Symbols}

    lisätä järjestyslukuihin viittaava merkintä X_{i, n}?

  • Poistathan numeroinnin yhtälöistä joihin ei viitata jälkeenpäin.

Pieniä kommentteja

Listaan alle joitain pieniä juttuja. Käythän tekstin läpi sillä silmällä etteivät nämä toistu muualla.

  • Tässä pitäisi olla \sup

    $x^*=sup\{x : F(x)<1\}$ & right endpoint of the distribution \\

  • Tämän voisi korvata täällä ehdotetulla tavalla.

    $1(p)=

  • Näissä pitäisi olla i.i.d. ja a.s. (yksi piste puuttuu)

    hill_simu/kandi/kandi.tex

    Lines 403 to 404 in 85b028d

    i.d.d & independent and identically distributed \\
    a.s & almost surely\\

  • Tässä pitäisi olla \max

    P(max(X_1, X_2, ... , X_n) \leq x) = P(X_1 \leq x, X_2 \leq x,..., X_n \leq x) = \\

  • Tällä rivillä pitäisi käyttää \left( ja \right) pelkkien sulkujen sijaan, jotta ne skaalautuvat oikein.

    \log(1 - \varepsilon) + (\gamma - \delta) \log(\frac{Y_{n-i,n}}{Y_{n-k,n}}) < \log(U(Y_{n-i,n})) - \log(U(Y_{n-k,n})) \\

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.