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);
    }
}