In [2]:
import numpy as np
import pandas as pd
from scipy import stats

import pyreto
from pyreto.tests import utilities

Pareto distributions


In [3]:
p = lambda alpha, xmin: stats.pareto(alpha-1, loc=0, scale=xmin)

In [4]:
prng = np.random.RandomState(42)
alpha, xmin = 3.0, 1.0
pareto_data = pd.Series(pyreto.distributions.Pareto._rvs(prng, 2500, xmin, alpha), name='samples')

In [5]:
# check that I get same estimate for alpha given reported xmin...
result1 = pyreto.distributions.Pareto.fit(pareto_data, xmin=xmin)

In [6]:
result1.params['alpha']


Out[6]:
2.9998119111215154

In [7]:
result1.xmin


Out[7]:
1.0

In [8]:
f = p(result1.params['alpha'], 1.0)

In [9]:
f.pdf(5.0)


Out[9]:
0.016003339044528487

In [10]:
pyreto.distributions.Pareto._pdf(5.0, 1.0, result1.params['alpha'])


Out[10]:
0.016003339044528487

In [12]:
f = stats.pareto(2.0, 1.0)

In [13]:
stats.pareto.

In [28]:
pyreto.distributions.Pareto._pdf??

In [22]:
# check that I get same estimate for alpha given reported xmin...
result2 = pyreto.distributions.Pareto.fit(pareto_data, xmin=None, quantile=0.9, method='brute')

In [23]:
result2.params['alpha']


Out[23]:
1.9909069323487545

In [24]:
result2.xmin


Out[24]:
1.1331781782153696

In [25]:
# check that I get same estimate for alpha given reported xmin...
result3 = pyreto.distributions.Pareto.fit(pareto_data, xmin=None, quantile=0.9, method='bounded')

In [26]:
result3.params['alpha']


Out[26]:
1.9924312888829019

In [27]:
result3.xmin


Out[27]:
1.3693896782872266

Cauchy distribution


In [7]:
prng = np.random.RandomState(1234)
cauchy_rvs = stats.cauchy.rvs(size=2500, random_state=prng)
cauchy_data = pd.Series(cauchy_rvs, name='samples')

In [8]:
# check that I get same estimate for alpha given reported xmin...
result = pyreto.distributions.Pareto.fit(cauchy_data, xmin=None, quantile=0.99, method='brute')

In [12]:
# will raise assertion error if test fails!
desired_alpha = 2.0
utilities.test_scaling_exponent_estimation(desired_alpha, result)

Levy-Stable distributions


In [2]:
prng = np.random.RandomState(4321)
desired_alpha, beta = 1.5, 1.0 
stable_rvs = stats.levy_stable.rvs(desired_alpha - 1, beta, size=1000, random_state=prng)
stable_data = pd.Series(stable_rvs, name='samples')

In [3]:
result = pyreto.distributions.Pareto.fit(stable_data, xmin=None, quantile=0.99, method='bounded')

In [4]:
# will raise assertion error if test fails!
utilities.test_scaling_exponent_estimation(desired_alpha, result)

In [84]:
# goodness of fit test should fail to reject Pareto tails!
test_kwargs = {'seed': None, 'result': result, 'data': stable_data,
               'xmin': None, 'quantile': 0.99, 'method': 'bounded'}
pvalue, _ = pyreto.distributions.Pareto.test_goodness_of_fit(**test_kwargs)
assert pvalue > 0.1, "Goodness of fit test Type I error!"

In [ ]: