威布尔分布代码实现

public class WeibullDistribution : ContinuousDistribution
    {
        private class MMSystem
        {
            private double mean;

            private double variance;

            public MMSystem(NumericalVariable variable)
            {
                this.mean = variable.Mean;
                this.variance = variable.Variance;
            }

            private Vector Evaluate(Vector input, Vector output)
            {
                double num;
                num = input[0];
                double num2;
                num2 = input[1];
                double num3;
                num3 = GammaFunctions.Gamma(1.0 + 1.0 / num2);
                output[0] = this.mean - num * num3;
                output[1] = this.variance - num * num * (GammaFunctions.Gamma(1.0 + 2.0 / num2) - num3 * num3);
                return output;
            }

            public void Solve(ref double a, ref double b)
            {
                DoglegSystemSolver doglegSystemSolver;
                doglegSystemSolver = new DoglegSystemSolver();
                doglegSystemSolver.TargetFunction = Evaluate;
                doglegSystemSolver.InitialGuess = new GeneralVector(a, b);
                Vector vector;
                vector = doglegSystemSolver.Solve();
                a = vector[0];
                b = vector[1];
            }
        }

        private class MLSystem
        {
            private NumericalVariable variable;

            private NumericalVariable log;

            public MLSystem(NumericalVariable variable)
            {
                this.variable = variable;
                this.log = variable.Transforms.Log();
            }

            private Vector Evaluate(Vector input, Vector output)
            {
                double num;
                num = input[0];
                double num2;
                num2 = input[1];
                int length;
                length = this.variable.Length;
                double num3;
                num3 = 0.0;
                double num4;
                num4 = 0.0;
                double num5;
                num5 = 0.0;
                double num6;
                num6 = Math.Log(num);
                for (int i = 0; i < length; i++)
                {
                    double num7;
                    num7 = Math.Pow(this.variable[i] / num, num2);
                    double num8;
                    num8 = this.log[i] - num6;
                    num3 += num8;
                    num5 += num8 * num7;
                    num4 += num7;
                }
                output[0] = (double)length + num2 * num3 - num2 * num5;
                output[1] = (double)length - num4;
                return output;
            }

            public void Solve(ref double a, ref double b)
            {
                DoglegSystemSolver doglegSystemSolver;
                doglegSystemSolver = new DoglegSystemSolver();
                doglegSystemSolver.TargetFunction = Evaluate;
                doglegSystemSolver.InitialGuess = new GeneralVector(a, b);
                doglegSystemSolver.Solve();
                a = doglegSystemSolver.CurrentPoint[0];
                b = doglegSystemSolver.CurrentPoint[1];
            }
        }

        private double scale;

        private double shape;

        private double location;

        public override double Mean => this.location + this.scale * GammaFunctions.Gamma(1.0 + 1.0 / this.shape);

        public override double Variance => this.scale * this.scale * (GammaFunctions.Gamma(1.0 + 2.0 / this.shape) - Math.Pow(GammaFunctions.Gamma(1.0 + 1.0 / this.shape), 2.0));

        public override double Skewness
        {
            get
            {
                double num;
                num = GammaFunctions.Gamma(1.0 + 1.0 / this.shape);
                double num2;
                num2 = GammaFunctions.Gamma(1.0 + 2.0 / this.shape);
                double num3;
                num3 = num2 - num * num;
                return (num * (2.0 * num * num - 3.0 * num2) + GammaFunctions.Gamma(1.0 + 3.0 / this.shape)) / (num3 * Math.Sqrt(num3));
            }
        }

        public override double Kurtosis
        {
            get
            {
                double num;
                num = GammaFunctions.Gamma(1.0 + 1.0 / this.shape);
                double num2;
                num2 = GammaFunctions.Gamma(1.0 + 2.0 / this.shape);
                double num3;
                num3 = num2 - num * num;
                return (6.0 * num * num * (2.0 * num2 - num * num) - 3.0 * num2 * num2 - 4.0 * num * GammaFunctions.Gamma(1.0 + 3.0 / this.shape) + GammaFunctions.Gamma(1.0 + 4.0 / this.shape)) / (num3 * num3);
            }
        }

        public double LocationParameter => this.location;

        public double ScaleParameter => this.scale;

        public double ShapeParameter => this.shape;

        public WeibullDistribution(double shape, double scale, double location)
        {
            if (shape <= 0.0)
            {
                ThrowException.ArgumentOutOfRange("shape");
            }
            if (scale <= 0.0)
            {
                ThrowException.ArgumentOutOfRange("scale");
            }
            this.scale = scale;
            this.shape = shape;
            this.location = location;
        }

        public WeibullDistribution(double shape, double scale)
            : this(shape, scale, 0.0)
        {
        }

        public WeibullDistribution(double shape)
            : this(shape, 1.0, 0.0)
        {
        }

        public WeibullDistribution(NumericalVariable variable, EstimationMethod method)
        {
            if (variable == null)
            {
                ThrowException.ArgumentNull("variable");
            }
            double a;
            a = variable.Mean;
            double b;
            b = 1.0;
            switch (method)
            {
            case EstimationMethod.Default:
            case EstimationMethod.MatchingMoments:
                new MMSystem(variable).Solve(ref a, ref b);
                break;
            case EstimationMethod.MaximumLikelihood:
                new MLSystem(variable).Solve(ref a, ref b);
                break;
            default:
                ThrowException.ArgumentOutOfRange("method");
                break;
            }
            this.scale = a;
            this.shape = b;
        }

        public WeibullDistribution(NumericalVariable variable)
            : this(variable, EstimationMethod.MatchingMoments)
        {
        }

        public override double ProbabilityDensityFunction(double x)
        {
            if (x < this.location)
            {
                return 0.0;
            }
            if (x == this.location)
            {
                if (this.shape == 1.0)
                {
                    return 1.0 / this.scale;
                }
                if (this.shape < 1.0)
                {
                    return double.PositiveInfinity;
                }
                return 0.0;
            }
            if (this.shape == 1.0)
            {
                return Math.Exp((0.0 - (x - this.location)) / this.scale) / this.scale;
            }
            return this.shape / this.scale * Math.Exp(0.0 - Math.Pow((x - this.location) / this.scale, this.shape) + (this.shape - 1.0) * Math.Log((x - this.location) / this.scale));
        }

        public override double DistributionFunction(double x)
        {
            if (x <= this.location)
            {
                return 0.0;
            }
            return 1.0 - Math.Exp(0.0 - Math.Pow((x - this.location) / this.scale, this.shape));
        }

        public override double MomentFunction(int order, double x)
        {
            if (order < 0)
            {
                ThrowException.ArgumentOutOfRange("order");
            }
            if (order == 0)
            {
                return this.DistributionFunction(x) - this.DistributionFunction(0.0);
            }
            if (x <= this.location)
            {
                return 0.0;
            }
            if (this.shape == 1.0)
            {
                return Math.Exp((0.0 - this.location) / this.scale) * ElementaryFunctions.Pow(0.0 - this.scale, order) * GammaFunctions.IncompleteGamma(1 + order, (0.0 - x) / this.scale);
            }
            if (this.location == 0.0)
            {
                return ElementaryFunctions.Pow(this.scale, order) * (GammaFunctions.Gamma(1.0 + (double)order / this.shape) - GammaFunctions.IncompleteGamma(1.0 + (double)order / this.shape, Math.Pow(x / this.scale, this.shape)));
            }
            return base.MomentFunction(order, x);
        }

        public override double InverseDistributionFunction(double probability)
        {
            if (probability < 0.0 || probability > 1.0)
            {
                ThrowException.ArgumentOutOfRange("probability");
            }
            if (probability == 0.0)
            {
                return this.location;
            }
            if (probability == 1.0)
            {
                return double.PositiveInfinity;
            }
            return this.location + this.scale * Math.Pow(0.0 - Math.Log(1.0 - probability), 1.0 / this.shape);
        }

        public override double GetRandomVariate(System.Random random)
        {
            if (random == null)
            {
                ThrowException.ArgumentNull("random");
            }
            return this.location + this.scale * Math.Pow(0.0 - Math.Log(random.NextDouble()), 1.0 / this.shape);
        }
    }

 这样很棒！

 如果对您有帮忙，非常感谢您支持一下创造者的付出！

 感谢支持技术分享，请扫码点赞支持：

技术合作交流qq:2401315930