scipy genextreme.fit gives: "runtimewarning: divide by zero"

The problem is that your code has some calculations in which it is trying to divide by zero or NaN.
You can either fix that by diving into the individual calculations and handling the zeroes and NaNs such that your algorithm skips them or you could do the following:

import numpy as np
from scipy.stats
import genextreme as gev
with np.errstate(divide = 'ignore', invalid = 'ignore'):
   t = np.array([3.8482, 3.6435, 3.0417, 4.0329, 3.2967, 3.3535, 3.6179, 3.3042, 3.6164, 3.5855, 2.7932, 2.8833, 2.6513, 2.7794, 3.2649, 3.2613, 3.1736, 3.1131, 3.3896, 4.2891])
a = gev.fit(t)

Suggestion : 2

I'm trying to do a GEV-fit using the anycodings_python genextreme package in SciPy. Although I can anycodings_python get it to estimate the parameters, I get a anycodings_python warning that it is dividing by zero.,These values are close to what I expected, anycodings_python but I don't understand the warning. ,I'm new to Python and wanted to do a quick anycodings_python first test with my existing installation. ,The problem is that your code has some anycodings_scipy calculations in which it is trying to anycodings_scipy divide by zero or NaN. You can either anycodings_scipy fix that by diving into the individual anycodings_scipy calculations and handling the zeroes and anycodings_scipy NaNs such that your algorithm skips them anycodings_scipy or you could do the following:

The data given in the code is roughly Gumbel anycodings_python distributed (c=0), since I have checked it anycodings_python on a Gumbel plot and made a Gumbel fit using anycodings_python linear regression.

import numpy as np
from scipy.stats
import genextreme as gev
t = np.array([3.8482, 3.6435, 3.0417, 4.0329, 3.2967, 3.3535, 3.6179, 3.3042, 3.6164, 3.5855, 2.7932, 2.8833, 2.6513, 2.7794, 3.2649, 3.2613, 3.1736, 3.1131, 3.3896, 4.2891])
a = gev.fit(t)

Based on my linear regression I expected a anycodings_python shape parameter close to c=0, a location anycodings_python parameter close to 3.15 and a scale anycodings_python parameter close to 0.39. The actual result anycodings_python was:

However, a value is apparently assigned to anycodings_python a, since:

>>> a
   (0.16458924337692377, 3.1800328240261857, 0.37894174199431357)

The problem is that your code has some anycodings_scipy calculations in which it is trying to anycodings_scipy divide by zero or NaN. You can either anycodings_scipy fix that by diving into the individual anycodings_scipy calculations and handling the zeroes and anycodings_scipy NaNs such that your algorithm skips them anycodings_scipy or you could do the following:

import numpy as np
from scipy.stats
import genextreme as gev
with np.errstate(divide = 'ignore', invalid = 'ignore'):
   t = np.array([3.8482, 3.6435, 3.0417, 4.0329, 3.2967, 3.3535, 3.6179, 3.3042, 3.6164, 3.5855, 2.7932, 2.8833, 2.6513, 2.7794, 3.2649, 3.2613, 3.1736, 3.1131, 3.3896, 4.2891])
a = gev.fit(t)

Suggestion : 3

Posted on Mon 10 September 2018 in analysis

import numpy as np
from scipy.stats
import genextreme as gev
from scipy.stats
import norm
from scipy.integrate
import trapz
import matplotlib.pyplot as plt
np.random.seed(523)

def get_sinusoidal_signal(n, amplitude, period, phase, mu0 = 4.0, coeff_var = 0.1):
   time = np.linspace(0, n, n)
n_predictors = len(amplitude)
assert len(period) == n_predictors
assert len(phase) == n_predictors
predictors = np.array([np.sin(2 * np.pi * time / period[i] + phase[i]) for i in range(n_predictors)]).T
mu = np.repeat(mu0, time.size)
for i in range(n_predictors):
   mu += amplitude[i] * predictors[: , i]
sigma = coeff_var * mu
sigma[np.where(sigma < 0.1)] = 0.1
observed = np.exp(np.random.normal(loc = mu, scale = sigma))
return time, observed

n_predictors = 3
time, observed = get_sinusoidal_signal(
   n = 100000,
   period = [9, 13, 24],
   amplitude = np.random.normal(loc = 0, scale = 0.5, size = n_predictors),
   phase = np.random.uniform(low = 0, high = 2 * np.pi, size = n_predictors),
)

muhat, sigmahat = norm.fit(np.log(observed))
shapehat, lochat, scalehat = gev.fit(observed)
sflo = np.linspace(0, 600, 1000) # we will make predictions here
norm_prob = norm.pdf(np.log(sflo), loc = muhat, scale = sigmahat)
norm_prob /= trapz(x = sflo, y = norm_prob)
gev_prob = gev.pdf(sflo, shapehat, scale = scalehat, loc = lochat)

/usr/local / miniconda3 / envs / jupyter - blog / lib / python3 .7 / site - packages / ipykernel_launcher.py: 4: RuntimeWarning: divide by zero encountered in log
after removing the cwd from sys.path.

fig, axes = plt.subplots(nrows = 1, ncols = 2, figsize = (9, 4.5), sharey = True, gridspec_kw = {
   'width_ratios': [2, 1]
})
ax = axes[0]
ax.plot(time[0: 1000], observed[0: 1000])
ax.grid()
ax.set_title('Time Series of Observed Streamflow')
ax.set_xlabel('Time [year]')
ax.set_ylabel('Observed Streamflow')

ax = axes[1]
ax.plot(norm_prob, sflo, label = 'Naive Log-Normal')
ax.plot(gev_prob, sflo, label = 'GEV')
ax.hist(observed, orientation = "horizontal", density = True, bins = 30, alpha = 0.5, label = 'Observed')
ax.grid()
ax.set_xlabel('Probability Distribution Function')
ax.set_title('Histogram and GEV Fit')
ax.legend()

plt.show()