GeneralizedParetoDistribution.cs

[FREENQ]

/*

• Fredrik Enqvist
• [email protected]
• 2016-02-17
• Not implemented (lack of time/interest/understanding): Variance(), Entropy(), Support(), Fit(), LogProbabilityDensityFunction(), Generate()
• Implemeted and tested against Extreme Optimization with identical results: Mean(), Median(), ProbabilityDensityFunction(), DistributionFunction()
  */

namespace Accord.Statistics.Distributions.Univariate
{
using Accord.Math;
using Accord.Math.Optimization;
using Accord.Statistics.Distributions;
using Accord.Statistics.Distributions.Fitting;
using AForge;
using System;
using System.ComponentModel;

public class GeneralizedParetoDistribution : UnivariateContinuousDistribution,
    IFormattable, ISampleableDistribution<double>
{

    double location;
    double scale;
    double shape;    


    /// <summary>
    ///   Creates new Pareto distribution.
    /// </summary>
    /// 
    public GeneralizedParetoDistribution()
        : this(1, 1, 1)
    {
    }

    /// <summary>
    /// Creates a Generalized Pareto Distribution.
    /// </summary>
    /// <param name="location">The x location.</param>
    /// <param name="scale">The scale parameter. Must be > 0.</param>
    /// <param name="shape">The shape of the distribution.</param>
    public GeneralizedParetoDistribution(double location, [Positive] double scale, [Positive] double shape)
    {
        if (scale <= 0)
        {
            throw new ArgumentOutOfRangeException("scale",
                "Scale must be positive.");
        }

        init(location, scale, shape);
    }


    private void init(double location, double scale, double shape)
    {
        this.location = location;
        this.scale = scale;
        this.shape = shape;
    }

    /// <summary>
    /// Get scale parameter.
    /// </summary>
    public double Scale
    {
        get { return scale; }
    }

    /// <summary>
    /// Get shape parameter.
    /// </summary>
    public double Shape
    {
        get { return shape; }
    }


    public override double Variance
    {
        get { return 999; }
    }


    public override double Entropy
    {
        get { return 999; }
    }


    public override DoubleRange Support
    {
        get { return new DoubleRange(scale, Double.PositiveInfinity); }
    }

    /// <summary>
    /// Mode. (Always == location???).
    /// </summary>
    public override double Mode
    {
        get { return location; }
    }


    public override double Mean
    {
        get { return location + (scale / (1 - shape)); }
    }


    public override double Median
    {
        get { return location + (scale * (Math.Pow(2, shape) - 1) / shape); }
    }

    /// <summary>
    /// The cumulative distribution function (CDF) evaluated at point x.
    /// </summary>
    /// <param name="x">Point x.</param>
    /// <returns>The CDF at point.</returns>
    public override double DistributionFunction(double x)
    {
        // PDF components
        double m = (x - location) / scale;
        double k = 1 + (shape * m);
        double l = -1 / shape;

        // domain logic
        bool shapePos = shape > 0;
        bool shapeNeg = shape < 0;
        bool shapeZero = shape == 0; // special case 

        bool xA = x >= location;
        bool xB = x <= (location - (scale / shape));

        bool A = shapePos && xA;
        bool B = shapeNeg && xA && xB;
        bool C = shapeZero && xA; // special case

        // CDF function
        if (A || B)
            return 1 - Math.Pow(k, l);
        if (C)
            return 1 - Math.Exp(-1 * m);

        return 0;
    }


    /// <summary>
    /// The probability distribution function (PDF) evaluated at point x.
    /// </summary>
    /// <param name="x">Point x.</param>
    /// <returns>PDF at point x.</returns>
    public override double ProbabilityDensityFunction(double x)
    {
        // PDF components
        double m = (x - location) / scale;
        double k = 1 + (shape * m);
        double l = -1 * ((1 / shape) + 1);

        // domain logic
        bool shapePos = shape > 0; 
        bool shapeNeg = shape < 0;
        bool shapeZero = shape == 0; // special case 

        bool xA = x >= location;
        bool xB = x <= (location - (scale / shape));

        bool A = shapePos && xA;
        bool B = shapeNeg && xA && xB;
        bool C = shapeZero && xA; // special case

        // PDF function
        if (A || B)
            return Math.Pow(k, l) / scale;
        if (C)
            return Math.Exp(-1 * m) / scale;

        return 0;
    }


    // Log PDF
    public override double LogProbabilityDensityFunction(double x)
    {
        if (x >= scale)
            return 999;
        return 0;
    }


    /// <summary>
    ///   Fits the underlying distribution to a given set of observations.
    /// </summary>
    /// <param name="observations">
    ///   The array of observations to fit the model against. The array
    ///   elements can be either of type double (for univariate data) or
    ///   type double[] (for multivariate data).
    /// </param>
    /// <param name="weights">
    ///   The weight vector containing the weight for each of the samples.</param>
    /// <param name="options">
    ///   Optional arguments which may be used during fitting, such
    ///   as regularization constants and additional parameters.</param>
    public override void Fit(double[] observations, double[] weights, IFittingOptions options)
    {
        if (options != null)
            throw new ArgumentException("This method does not accept fitting options.");

        double xm = observations.Min();

        double lnx = Math.Log(xm);

        double alpha;

        if (weights == null)
        {
            double den = 0;
            for (int i = 0; i < observations.Length; i++)
                den += Math.Log(observations[i]) - lnx;
            alpha = observations.Length / den;
        }
        else
        {
            double den = 0;
            for (int i = 0; i < observations.Length; i++)
                den += (Math.Log(observations[i]) - lnx) * weights[i];
            alpha = weights.Sum() / den;
        }

        //init(xm, alpha);
    }


    public override object Clone()
    {
        return new ParetoDistribution(scale, shape);
    }


    public override string ToString(string format, IFormatProvider formatProvider)
    {
        return String.Format(formatProvider, "Pareto(x; xm = {0}, α = {1})",
            scale.ToString(format, formatProvider),
            shape.ToString(format, formatProvider));
    }

    /// <summary>
    ///   Generates a random vector of observations from the current distribution.
    /// </summary> 
    /// <param name="samples">The number of samples to generate.</param>
    /// <returns>A random vector of observations drawn from this distribution.</returns>
    public override double[] Generate(int samples)
    {
        double[] U = UniformContinuousDistribution.Standard.Generate(samples);

        for (int i = 0; i < U.Length; i++)
            U[i] = scale / Math.Pow(U[i], 1.0 / shape);

        return U;
    }

    /// <summary>
    ///   Generates a random observation from the current distribution.
    /// </summary>
    /// <returns>A random observations drawn from this distribution.</returns> 
    public override double Generate()
    {
        double U = UniformContinuousDistribution.Standard.Generate();

        return scale / Math.Pow(U, 1.0 / shape);
    }
}

}

[cesarsouza]

Thanks a lot for the contribution!! I will fill in some of the missing methods and commit it to the repository.

Regards,
Cesar

[FREENQ]

Happy to contribute! Great work Cesar!

Besides the 999 usd Extreme Optimization framework, Accord.NET is the only one providing the Gen Pareto distribution class.

/FE

---

From: César Souza [email protected]
Sent: Sunday, August 21, 2016 5:04 PM
To: accord-net/framework
Cc: Fredrik Enqvist; Author
Subject: Re: [accord-net/framework] GeneralizedParetoDistribution.cs (#201)

Thanks a lot for the contribution!! I will fill in some of the missing methods and commit it to the repository. Thansk again!

Regards,
Cesar

You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHubhttps://github.com//issues/201#issuecomment-241262766, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AQf9sAzTOIeemQQ0O-xeKN0YvYDjiMQuks5qiGjugaJpZM4HccJf.