Categories &

Functions List

Function Reference: gpfit

statistics: paramhat = gpfit (x)
statistics: [paramhat, paramci] = gpfit (x)
statistics: [paramhat, paramci] = gpfit (x, alpha)
statistics: [paramhat, paramci] = gpfit (x, alpha, options)
statistics: [paramhat, paramci] = gpfit (x, alpha, options, freq)

Estimate parameters and confidence intervals for the generalized Pareto distribution.

paramhat = gpfit (x) returns the maximum likelihood estimates of the parameters of the generalized Pareto distribution given the data in x. paramhat(1) is the shape parameter, k, and paramhat(2) is the scale parameter, sigma.

gpfit does not estimate the location parameter theta and assumes it to be zero, so x must not contain negative values. To fit data with a known nonzero theta, subtract it from x before calling gpfit; the estimates of k and sigma are unchanged by the shift.

[paramhat, paramci] = gpfit (x) returns the 95% confidence intervals for the estimated parameters k and sigma as a 2-by-2 matrix whose first row holds the lower bounds and whose second row holds the upper bounds.

[…] = gpfit (x, alpha) also returns the 100 * (1 - alpha) percent confidence intervals for the parameter estimates. By default, the optional argument alpha is 0.05 corresponding to 95% confidence intervals. Pass in [] for alpha to use the default values.

[…] = gpfit (x, alpha, options) specifies control parameters for the iterative algorithm used to compute ML estimates with the fminsearch function. options is a structure with the following fields and their default values:

  • options.Display = "off"
  • options.MaxFunEvals = 400
  • options.MaxIter = 200
  • options.TolX = 1e-6

[…] = gpfit (x, alpha, options, freq) accepts a vector of the same size as x giving the number of times each element of x was observed. This fourth argument is an Octave extension; MATLAB’s gpfit takes three inputs at most.

When the shape parameter falls below -1 the likelihood is unbounded: the density at the upper endpoint of the support diverges as that endpoint closes onto the largest observation, so no maximum likelihood estimate exists and whatever is returned is an arbitrary point on that ridge. gpfit warns in this case and returns NaN confidence intervals. The estimate it does return always keeps every observation strictly inside the fitted support, since the likelihood is infinite outside it. This is a deliberate deviation: MATLAB has been measured returning parameters for such data under which the largest observation has zero density and its own gplike returns Inf.

When k = 0 and theta = 0, the Generalized Pareto is equivalent to the exponential distribution. When k > 0 and theta = k / k the Generalized Pareto is equivalent to the Pareto distribution. The mean of the Generalized Pareto is not finite when k >= 1 and the variance is not finite when k >= 1/2. When k >= 0, the Generalized Pareto has positive density for x > theta, or, when theta < 0, for 0 <= (x - theta) / sigma <= -1 / k.

Further information about the generalized Pareto distribution can be found at https://en.wikipedia.org/wiki/Generalized_Pareto_distribution

See also: gpcdf, gpinv, gppdf, gprnd, gplike, gpstat

Source Code: gpfit

Sample 2 populations from different generalized Pareto distributions Assume location parameter θ is known

 theta = 0;
 rand ('seed', 5);    # for reproducibility
 r1 = gprnd (1, 2, theta, 20000, 1);
 rand ('seed', 2);    # for reproducibility
 r2 = gprnd (3, 1, theta, 20000, 1);
 r = [r1, r2];

Plot them normalized and fix their colors

 hist (r, [0.1:0.2:100], 5);
 h = findobj (gca, 'Type', 'patch');
 set (h(1), 'facecolor', 'r');
 set (h(2), 'facecolor', 'c');
 ylim ([0, 1]);
 xlim ([0, 5]);
 hold on

Estimate their α and β parameters

 k_sigmaA = gpfit (r(:,1));
 k_sigmaB = gpfit (r(:,2));

Plot their estimated PDFs

 x = [0.01, 0.1:0.2:18];
 y = gppdf (x, k_sigmaA(1), k_sigmaA(2), theta);
 plot (x, y, '-pc');
 y = gppdf (x, k_sigmaB(1), k_sigmaB(2), theta);
 plot (x, y, '-sr');
 hold off
 legend ({'Normalized HIST of sample 1 with k=1 and σ=2', ...
          'Normalized HIST of sample 2 with k=2 and σ=2', ...
          sprintf("PDF for sample 1 with estimated k=%0.2f and σ=%0.2f", ...
                  k_sigmaA(1), k_sigmaA(2)), ...
          sprintf("PDF for sample 3 with estimated k=%0.2f and σ=%0.2f", ...
                  k_sigmaB(1), k_sigmaB(2))})
 title ('Two population samples from different generalized Pareto distributions')
 text (2, 0.7, 'Known location parameter θ = 0')
 hold off
plotted figure