;+ ; NAME: ; gev_pdf ; ; PURPOSE: ; This function calculates the probability density function of the ; Generalised Extreme Value distribution for a given set of parameters. ; ; CATEGORY: ; Statistics ; ; CALLING SEQUENCE: ; gev_pdf, values, location=location, scale=scale, shape=shape ; ; INPUTS: ; VALUES: A required floating vector of length N_VALUES specifying the ; values at which to calculate the probability density function. ; LOCATION, LOG, SCALE, SHAPE ; ; KEYWORD PARAMETERS: ; LOCATION: A required float scalar or N_VALUES vector specifying the ; location parameter of the GEV distribution. If an N_VALUES vector, ; then each of the N_VALUES elements is matched with the corresponding ; element in the VALUES input. ; LOG: If set, then the log likelihood (the logarithm of the probability ; density function) is returned. ; SCALE: A required float scalar or N_VALUES vector specifying the ; scale parameter of the GEV distribution. If an N_VALUES vector, ; then each of the N_VALUES elements is matched with the corresponding ; element in the VALUES input. ; SHAPE: A required float scalar or N_VALUES vector specifying the ; shape parameter of the GEV distribution. If an N_VALUES vector, ; then each of the N_VALUES elements is matched with the corresponding ; element in the VALUES input. ; ; OUTPUTS: ; RESULT: A floating vector of length N_VALUES containing the probability ; density function for each values specified in VALUES. ; ; USES: ; - ; ; PROCEDURE: ; This function calculates the cumulative density function of the ; Generalised Extreme Value distribution analytically. ; ; EXAMPLE: ; ; Plot the cumulative distribution function for a given set of parameters ; values = findgen( 501.) / 100. ; pdfgev = gev_pdf( values, location=2.7, scale=0.30, shape=-0.11 ) ; plot, values, pdfgev ; ; LICENSE: ; This code was written as part of the Extreme Weather Event Real-Time ; Attribution Machine (EWERAM) project supported by the New Zealand Ministry ; of Business, Innovation and Employment. The code is free for use under ; the terms of the Creative Commons License v2.0 ; (http://creativecommons.org/licenses/by-nc-sa/2.0/). ; ; MODIFICATION HISTORY: ; Written by: Daithi A. Stone (dastone@runbox.com), 2019-11-06 ;- ;*********************************************************************** FUNCTION GEV_PDF, $ VALUES, $ LOCATION=location, SCALE=scale, SHAPE=shape, $ LOG=log_opt ;*********************************************************************** ; Constants and checks ; The number of values for which to calculate the CDF n_values = n_elements( values ) if n_values eq 0 then stop ; Check the GEV parameters n_location = n_elements( location ) if max( n_location eq [ 1, n_values ] ) ne 1 then stop n_scale = n_elements( scale ) if max( n_scale eq [ 1, n_values ] ) ne 1 then stop if min( scale ) lt 0. then stop n_shape = n_elements( shape ) if max( n_shape eq [ 1, n_values ] ) ne 1 then stop ; The maximum number of parameters n_param_max = max( [ n_location, n_scale, n_shape ] ) ; The option to return the log likelihood log_opt = keyword_set( log_opt ) ;*********************************************************************** ; Calculate the CDF ; Copy parameters to temporary variables, making all consistent in length if n_param_max eq 1 then begin temp_location = location[0] temp_scale = scale[0] temp_shape = shape[0] endif else begin if n_location eq 1 then begin temp_location = location + fltarr( n_param_max ) endif else begin temp_location = location endelse if n_scale eq 1 then begin temp_scale = scale + fltarr( n_param_max ) endif else begin temp_scale = scale endelse if n_shape eq 1 then begin temp_shape = shape + fltarr( n_param_max ) endif else begin temp_shape = shape endelse endelse ; Initialise the output scalar/vector result = 0. * values ; The case when the shape parameter is zero id = where( temp_shape eq 0., n_id ) if n_id gt 0 then begin ; Calculate the standard PDF if log_opt eq 0 then begin if n_param_max eq 1 then begin result = 1. / temp_scale $ * exp( ( values - temp_location ) / temp_scale $ - exp( ( values - temp_location ) / temp_scale ) ) endif else begin result[id] = 1. / temp_scale[id] $ * exp( ( values[id] - temp_location[id] ) / temp_scale[id] $ - exp( ( values[id] - temp_location[id] ) / temp_scale[id] ) ) endelse ; Calculate the log likelihood endif else begin if n_param_max eq 1 then begin result = -alog( temp_scale ) $ + ( values - temp_location ) / temp_scale $ - exp( ( values - temp_location ) / temp_scale ) endif else begin result[id] = -alog( temp_scale[id] ) $ + ( values[id] - temp_location[id] ) / temp_scale[id] $ - exp( ( values[id] - temp_location[id] ) / temp_scale[id] ) endelse endelse endif ; The case when the shape parameter is not zero id = where( temp_shape ne 0., n_id ) if n_id gt 0 then begin ; If we have legal parameter combinations if n_param_max eq 1 then begin temp = 1. + temp_shape * ( values - temp_location ) / temp_scale endif else begin temp = 1. $ + temp_shape[id] * ( values[id] - temp_location[id] ) / temp_scale[id] endelse id_good = where( temp gt 0, n_id_good ) if n_id_good gt 0 then begin ; Calculate standard PDF if log_opt eq 0 then begin if n_param_max eq 1 then begin result[id_good] $ = 1. / temp_scale * ( temp[id_good] ^ ( -1. / temp_shape - 1 ) ) $ * ( exp( -( temp[id_good] ^ ( -1. / temp_shape ) ) ) ) endif else begin result[id_good] $ = 1. / temp_scale[id[id_good]] $ * ( temp[id_good] ^ ( -1. / temp_shape[id[id_good]] - 1. ) ) $ * ( exp( -( temp[id_good] ^ ( -1. / temp_shape[id[id_good]] ) ) ) ) endelse ; Calculate the log likelihood endif else begin if n_param_max eq 1 then begin result[id_good] = -alog( temp_scale ) $ - ( 1. + 1. / temp_shape ) * alog( temp[id_good] ) $ - temp[id_good] ^ ( -1. / temp_shape ) endif else begin result[id[id_good]] = -alog( temp_scale[id[id_good]] ) $ - ( 1. + 1. / temp_shape[id[id_good]] ) * alog( temp[id_good] ) $ - temp[id_good] ^ ( -1. / temp_shape[id[id_good]] ) endelse endelse endif endif ;*********************************************************************** ; The end return, result END