[ncl-talk] GEV assistance needed
Michael Notaro
mnotaro at wisc.edu
Tue Nov 23 13:07:32 MST 2021
As far as I can tell, the NCL-computed GEV parameters based on MLE do not match prior studies' results.
I say that with the huge caveat that I am not a GEV expert.
For example, based on the annual max values for 1950-2014 in Table 1 of this paper:
https://www.researchgate.net/publication/288991817_Predicted_Return_Level_of_Annual_Maximum_Earthquake_in_Northern_California/link/56881a6208aebccc4e155539/download
the location, scale, and shape parameters are 5.62, 0.54, and -0.098, respectively, in Table 2.
When I perform this task below with those data values,
yearmax=(/5.3,5.8,7.2,6.1,6.3,5.8,5.8,5.8,5.5,5.5,5.7,6.4,5.5,5.5,5.1,6.0,5.0, \
5.3,6.0,6.2,5.6,5.8,6.0,5.3,7.2,5.5,6.0,5.4,4.9,6.3,6.2,5.6,5.6,5.4, \
5.5,5.0,6.0,5.6,5.9,5.7,5.4,5.2,5.2,5.0,5.2,5.7,6.3,5.0,5.8,5.8,6.9, \
5.9,5.5,6.7,6.6,5.6,6.6,6.8,6.0,6.6,6.3,6.6,7.7,6.2,6.8,6.3,5.6,5.6, \
5.5,7.3,6.1,6.6,5.9,6.5,6.0,7.0,5.0,5.5,5.4,4.9,6.5,4.7,5.6,5.7,6.8/)
param=extval_mlegev(yearmax,0,False)
print(param)
I get the following, which does not match the results from the study at all.
(0) 3.7994 location
(1) 2.233068 scale
(2) 0.4374145 shape
Is there anyone with experience with GEV and NCL that can offer some insights?
Thanks,
Michael
Michael Notaro
Associate Director
Nelson Institute Center for Climatic Research
University of Wisconsin-Madison
Phone: (608) 261-1503
Email: mnotaro at wisc.edu
________________________________
From: Dennis Shea <shea at ucar.edu>
Sent: Wednesday, November 17, 2021 9:58 AM
To: Michael Notaro <mnotaro at wisc.edu>
Cc: ncl-talk at ucar.edu <ncl-talk at ucar.edu>
Subject: Re: [ncl-talk] GEV assistance needed
Hi Michael,
Not sure I can help on this. I am not very knowledgeable on MLE or GEV
As noted in the documentation, NCL invokes a August 1994 version of:
Algorithm 215: Maximum-Likelihood Estimation of the Parameters of the Generalized Extreme-Value Distribution
J. R. M. Hosking
Journal of the Royal Statistical Society. Series C (Applied Statistics)
Vol. 34, No. 3 (1985), pp. 301-310
URL: http://www.jstor.org/stable/2347483
Code: http://ftp.uni-bayreuth.de/math/statlib/apstat/215
My recollection is that we tested this on a few 'book' examples and compared the output to results generated by R's gev.fit function.
====
extval_mlegev<https://www.ncl.ucar.edu/Document/Functions/Built-in/extval_mlegev.shtml>: Example 1 illustrates a comparison. I believe you have used R. Perhaps trying your data using R
will provide another (NCL, MATLAB) viewpoint.
====
Again, I am not very knowledgeable on either MLE or GEV so I have no idea what the -1 does within the algorithm.
Regards
D
On Tue, Nov 16, 2021 at 3:24 PM Michael Notaro via ncl-talk <ncl-talk at mailman.ucar.edu<mailto:ncl-talk at mailman.ucar.edu>> wrote:
In the example below, I have 20 values of precipitation in inches.
Each value is the highest daily precipitation for one of 20 years.
I want to estimate the 1-day/100-year GEV precipitation level.
For some reason, my older script included multiplying the shape
parameter by -1. That gave a GEV value of 65757", definitely wrong.
If I get rid of the *-1 part, the GEV value becomes more reasonable,
1.73", but too low compared to the 20 annual max values. So I have
two questions.
1. Is there something wrong with the script below?
2. Is there ever a need to multiply any of the parameters by -1? I
3. did it once before but have no recollection of why I did it.
Thanks, Michael
yearmax=(/1.83,1.69,3.24,1.99,1.69,1.79,2.78,1.70,4.85,2.35,3.37,1.72,3.96,2.35,2.29,3.33,3.97,3.53,3.11,1.69/)
param=extval_mlegev(yearmax,0,False)
param_mle_location=param(0)
param_mle_scale=param(1)
param_mle_shape=param(2)*-1.
coeff=-1.*log(1.-1./100.)
gev=param_mle_scale/param_mle_shape*[coeff^(-1.*param_mle_shape)-1.]+param_mle_location
print("param_mle_location="+param_mle_location+" param_mle_scale="+param_mle_scale+" param_mle_shape="+param_mle_shape+" coeff="+coeff+" gev= "+gev)
param_mle_location=1.65715 param_mle_scale=0.216018 param_mle_shape=2.98224 coeff=0.0100503 gev= 65757.5
earmax=(/1.83,1.69,3.24,1.99,1.69,1.79,2.78,1.70,4.85,2.35,3.37,1.72,3.96,2.35,2.29,3.33,3.97,3.53,3.11,1.69/)
param=extval_mlegev(yearmax,0,False)
param_mle_location=param(0)
param_mle_scale=param(1)
param_mle_shape=param(2)*-1.
coeff=-1.*log(1.-1./100.)
gev=param_mle_scale/param_mle_shape*[coeff^(-1.*param_mle_shape)-1.]+param_mle_location
print("param_mle_location="+param_mle_location+" param_mle_scale="+param_mle_scale+" param_mle_shape="+param_mle_shape+" coeff="+coeff+" gev= "+gev)
param_mle_location=1.65715 param_mle_scale=0.216018 param_mle_shape=-2.98224 coeff=0.0100503 gev= 1.72958
Michael Notaro
Associate Director
Nelson Institute Center for Climatic Research
University of Wisconsin-Madison
Phone: (608) 261-1503
Email: mnotaro at wisc.edu<mailto:mnotaro at wisc.edu>
_______________________________________________
ncl-talk mailing list
ncl-talk at mailman.ucar.edu<mailto:ncl-talk at mailman.ucar.edu>
List instructions, subscriber options, unsubscribe:
https://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20211123/5fc43e25/attachment.html>
More information about the ncl-talk
mailing list