Extracting confidence intervals on fit parameters 'c', 'loc', and 'scale'

[arfogg]

Hi,

Thanks for creating this super helpful python package, it's really easy to use and the documentation are really helpful. I use this package for space science research.

I'm interested in returning the confidence intervals or some other measure of error on the fit parameters 'c', 'loc', and 'scale'.

Currently I access the fit parameters by doing: EVA.model.fit_parameters . Having confidence intervals (or any other measure of error) on these parameters would be a great addition to the package. Apologies if this feature has already been implemented - if it has please could you guide me on how to extract the confidence intervals / errors?

Thanks,
Alexandra Fogg

[georgebv]

There are no confidence intervals available right away for the parameters, but you can calculate those yourself from the data.

For the Emcee model you can use the _trace attribute of shape number of free parameters x number of samples.

For the MLE model you can use the fit_parameter_cache attribute which is a list of samples of free parameters. Before it's available, however, you would need to call get_return_value with alpha argument and n_samples set to number you need (100 by default).

[georgebv]

pyextremes/src/pyextremes/models/model_base.py

Line 71 in a0b4baa

self._trace: typing.Optional[np.ndarray] = None

pyextremes/src/pyextremes/models/model_mle.py

Line 50 in a0b4baa

self.fit_parameter_cache: typing.List[typing.Tuple[float, ...]] = []

pyextremes/src/pyextremes/models/model_mle.py

Lines 69 to 98 in a0b4baa

	def get_return_value(
	self, exceedance_probability, alpha: typing.Optional[float] = None, **kwargs
	) -> tuple:
	"""
	Calculate return value and confidence interval bounds.
	Parameters
	----------
	exceedance_probability : array-like
	Exceedance probability or 1D array of exceedance probabilities.
	Each exceedance probability must be in the [0, 1) range.
	alpha : float, optional
	Width of confidence interval (0, 1).
	If None (default), return None
	for upper and lower confidence interval bounds.
	kwargs
	n_samples : int, optional
	Number of bootstrap samples used to estimate
	confidence interval bounds (default=100).
	Returns
	-------
	return_value : array-like
	Return values.
	ci_lower : array-like
	Lower confidence interval bounds.
	ci_upper : array-like
	Upper confidence interval bounds.
	"""

[arfogg]

Hello!

Thanks for your reply. I have been working on this using your suggestion above.

Here's an example of what I'm running to extract confidence intervals on the fit parameters:

Initiating the model:
eva_model=pyextremes.EVA(org_df.b)
block_size=pd.Timedelta(pd.to_timedelta("365.2425D"))
eva_model.get_extremes("BM", block_size=block_size )
eva_model.fit_model(model='Emcee')

Extract the fit parameters:
b_gevd_fit_params=pd.DataFrame(eva_model.model.fit_parameters, index=[0])
b_gevd_fit_params.rename(columns={'c':'shape_', 'loc':'location'}, inplace=True)

Plot the distribution/histogram of all the samplers for the shape_ fit parameters
fig,ax=plt.subplots()
ax.hist(eva_model.model.trace[:,:,0].flatten(), bins=100, color='grey')

Plot a vertical line to show where the determined shape_ lies
ax.axvline(x=b_gevd_fit_params.shape_[0], color='violet')

Some labels
ax.set_xlabel('shape')
ax.set_ylabel('occurrence')

Results in the following plot:

Now let's calculate the confidence interval
shape_ci_width=0.95*((np.std(eva_model.model.trace[:,:,0].flatten()))/(np.sqrt(eva_model.model.trace[:,:,0].flatten().size)))

lower_interval=np.nanmean(eva_model.model.trace[:,:,0].flatten())-shape_ci_width
upper_interval=np.nanmean(eva_model.model.trace[:,:,0].flatten())+shape_ci_width

This gives very small shape_ci_width=0.000538, leading to lower_interval=-0.292882 and upper_interval=-0.291805 for a b_gevd_fit_params.shape_[0]=-0.272371.

This doesn't seem realistic to me based on my experience with the confidence intervals in the rest of pyextremes. Should I be adding/subtracting shape_ci_width from b_gevd_fit_params.shape_[0]? Still shape_ci_width is very small.

Please could you let me know if there is something obvious I'm doing wrong here. I'm really grateful for your help with this!!

Thanks,
Alexandra

[georgebv]

You extracted data correctly, which is the key part. What you got wrong is:

• You got mean but what (I think) you are really interested it is mode (histogram peak).
• You assume that data is normally distributed, which is clearly not the case. You should get quantiles of your sample - 95% confidence interval lies between 0.025 and 0.975 quantiles

[arfogg]

Thanks for your quick reply!

I think I have estimated the confidence intervals correctly now - thank you for your help. I'm just putting the final few code lines here in case anyone else comes across the same problem:

Adding onto the code in my comment above:
Calculate the 0.025 and 0.975 quantiles
quantiles=np.quantile(eva_model.model.trace[:,:,0].flatten(), [0.025, 0.975])

Overplot these quantiles onto the distribution
ax.axvline(x=quantiles[0], color='violet', linestyle='dashed')
ax.axvline(x=quantiles[1], color='violet', linestyle='dashed')

Results in the following plot:

Thanks again for creating and maintaining this package - it's super helpful and well documented so I'm able to use for my novel space science research!

Thanks,
Alexandra