参数估计Weibull分布两参数估计迭代算法

常用于为失效时间数据建模。例如，一个制造商希望计算某个部件在一年、两年或更多年后失效的概率。此分布广泛地应用于工程、医学研究、金融和气候学。 
Weibull 分布由形状、尺度和阈值等参数描述。阈值参数为零的情况称为 2 参数 Weibull 分布。只为非负变量定义此分布。 
取决于参数的值，Weibull 分布可以具有各种形状。

这种分布的主要优点之一在于它可以具有其他类型分布的特征，从而在拟合不同类型的数据时极其灵活。一般在可靠性分析中使用

常见数学统计算法包内包含各种分布的pdf，cdf，参数估计却很少提供，但是项目中必须要用，所以实现了一个经过优化的迭代算法（C#版本）

（其中有使用Gamma函数，正态分布等，比较常见，此处代码不提供了）

public sealed class WeibullDistribution
    {
        /// 形状参数
        private double _alpha;
        /// 尺度参数
        private double _beta;
        /// 正交化分布（方便计算）
        private double _norm;

        /// <summary>
        /// 创建一个分布
        /// </summary>
        /// <param name="shape"></param>
        /// <param name="scale"></param>
        public WeibullDistribution(double shape, double scale)
        {
            if (shape <= 0)
                throw new ArgumentOutOfRangeException(
                                        "Shape parameter must be positive");
            if (scale <= 0)
                throw new ArgumentOutOfRangeException(
                                        "Scale parameter must be positive");
            DefineParameters(shape, scale);
        }
        public double ln(double x) { return Math.Log(x, Math.E); }

        public double SigmaLnXi(IList<double> doubles)
        {
            double sum = 0;
            foreach (var item in doubles)
            {
                sum += ln(item);
            }
            return sum;
        }

        public double SigmaPowXi(IList<double> doubles, double beta0)
        {
            double sum = 0;
            foreach (var item in doubles)
            {
                sum += Math.Pow(item, beta0);
            }
            return sum;

        }

        public double SigmaPowXi2(IList<double> doubles, double beta0)
        {
            double sum = 0;
            foreach (var item in doubles)
            {
                sum += Math.Pow(item, beta0) * ln(item);
            }
            return sum;

        }
        /// <summary>
        /// 使用迭代计算数值解进行威布尔参数估计
        /// </summary>
        /// <param name="datas"></param>
        public WeibullDistribution(IList<double> datas)
        {
            //参数估计
            NumericalVariable n = new NumericalVariable(datas);
            double xbar = n.Mean;
            double sd = n.StandardDeviation;
            double E = 0.001;

            double b0 = 1.2 * xbar / sd;
            double b = b0;
            double Beta = int.MaxValue;
            //迭代计算beta
            while (Math.Abs(Beta - b) >= E)
            {
                Beta = 1.0 / ((SigmaPowXi2(datas, b) / SigmaPowXi(datas, b)) - (1.0 / datas.Count * SigmaLnXi(datas)));
                b = (Beta + b) / 2;
            }
            //
           //计算Alpha
            double Alpha = Math.Pow(1.0 / datas.Count * SigmaPowXi(datas, Beta), 1.0 / Beta);
            DefineParameters(Beta, Alpha);
        }



     
        public double Average
        {
            get { return Fn.Gamma(1 / _alpha) * _beta / _alpha; }

            set
            {
                throw new InvalidOperationException(
                                        "Can not set average on Weibull distribution");
            }
        }

       
        public void DefineParameters(double shape, double scale)
        {
            _alpha = shape;
            _beta = scale;
            _norm = _alpha / Math.Pow(_beta, _alpha);
        }

      
        public double DistributionValue(double x)
        {
            return 1.0 - Math.Exp(-Math.Pow(x / _beta, _alpha));
        }

       
        public string Name
        {
            get { return "Weibull distribution"; }
        }

        public double[] Parameters
        {
            get { return new double[] { _alpha, _beta }; }
            set { DefineParameters(value[0], value[1]); }
        }

       
        public double InverseDistributionValue(double x)
        {
            return Math.Pow(-Math.Log(1 - x), 1.0 / _alpha) * _beta;
        }

      
        public override string ToString()
        {
            return string.Format("Weibull distribution ({0:####0.00000},{1:####0.00000})", _alpha, _beta);
        }

       
        public double Value(double x)
        {
            return _norm * Math.Pow(x, _alpha - 1) * Math.Exp(-Math.Pow(x / _beta, _alpha));
        }

       
        public double Variance
        {
            get
            {
                double s = Fn.Gamma(1 / _alpha);
                return _beta * _beta * (2 * Fn.Gamma(2 / _alpha)
                                                        - s * s / _alpha) / _alpha;
            }
        }

      
    }