参数估计-Weibull分布-两参数估计迭代算法
常用于为失效时间数据建模。例如,一个制造商希望计算某个部件在一年、两年或更多年后失效的概率。此分布广泛地应用于工程、医学研究、金融和气候学。
Weibull 分布由形状、尺度和阈值等参数描述。阈值参数为零的情况称为 2 参数 Weibull 分布。只为非负变量定义此分布。
取决于参数的值,Weibull 分布可以具有各种形状。
这种分布的主要优点之一在于它可以具有其他类型分布的特征,从而在拟合不同类型的数据时极其灵活。一般在可靠性分析中使用
常见数学统计算法包内包含各种分布的pdf,cdf,参数估计却很少提供,但是项目中必须要用,所以实现了一个经过优化的迭代算法(C#版本)
(其中有使用Gamma函数,正态分布等,比较常见,此处代码不提供了)
1 public sealed class WeibullDistribution 2 { 3 /// 形状参数 4 private double _alpha; 5 /// 尺度参数 6 private double _beta; 7 /// 正交化分布(方便计算) 8 private double _norm; 9 10 /// <summary> 11 /// 创建一个分布 12 /// </summary> 13 /// <param name="shape"></param> 14 /// <param name="scale"></param> 15 public WeibullDistribution(double shape, double scale) 16 { 17 if (shape <= 0) 18 throw new ArgumentOutOfRangeException( 19 "Shape parameter must be positive"); 20 if (scale <= 0) 21 throw new ArgumentOutOfRangeException( 22 "Scale parameter must be positive"); 23 DefineParameters(shape, scale); 24 } 25 public double ln(double x) { return Math.Log(x, Math.E); } 26 27 public double SigmaLnXi(IList<double> doubles) 28 { 29 double sum = 0; 30 foreach (var item in doubles) 31 { 32 sum += ln(item); 33 } 34 return sum; 35 } 36 37 public double SigmaPowXi(IList<double> doubles, double beta0) 38 { 39 double sum = 0; 40 foreach (var item in doubles) 41 { 42 sum += Math.Pow(item, beta0); 43 } 44 return sum; 45 46 } 47 48 public double SigmaPowXi2(IList<double> doubles, double beta0) 49 { 50 double sum = 0; 51 foreach (var item in doubles) 52 { 53 sum += Math.Pow(item, beta0) * ln(item); 54 } 55 return sum; 56 57 } 58 /// <summary> 59 /// 使用迭代计算数值解进行威布尔参数估计 60 /// </summary> 61 /// <param name="datas"></param> 62 public WeibullDistribution(IList<double> datas) 63 { 64 //参数估计 65 NumericalVariable n = new NumericalVariable(datas); 66 double xbar = n.Mean; 67 double sd = n.StandardDeviation; 68 double E = 0.001; 69 70 double b0 = 1.2 * xbar / sd; 71 double b = b0; 72 double Beta = int.MaxValue; 73 //迭代计算beta 74 while (Math.Abs(Beta - b) >= E) 75 { 76 Beta = 1.0 / ((SigmaPowXi2(datas, b) / SigmaPowXi(datas, b)) - (1.0 / datas.Count * SigmaLnXi(datas))); 77 b = (Beta + b) / 2; 78 } 79 // 80 //计算Alpha 81 double Alpha = Math.Pow(1.0 / datas.Count * SigmaPowXi(datas, Beta), 1.0 / Beta); 82 DefineParameters(Beta, Alpha); 83 } 84 85 86 87 88 public double Average 89 { 90 get { return Fn.Gamma(1 / _alpha) * _beta / _alpha; } 91 92 set 93 { 94 throw new InvalidOperationException( 95 "Can not set average on Weibull distribution"); 96 } 97 } 98 99 100 public void DefineParameters(double shape, double scale) 101 { 102 _alpha = shape; 103 _beta = scale; 104 _norm = _alpha / Math.Pow(_beta, _alpha); 105 } 106 107 108 public double DistributionValue(double x) 109 { 110 return 1.0 - Math.Exp(-Math.Pow(x / _beta, _alpha)); 111 } 112 113 114 public string Name 115 { 116 get { return "Weibull distribution"; } 117 } 118 119 public double[] Parameters 120 { 121 get { return new double[] { _alpha, _beta }; } 122 set { DefineParameters(value[0], value[1]); } 123 } 124 125 126 public double InverseDistributionValue(double x) 127 { 128 return Math.Pow(-Math.Log(1 - x), 1.0 / _alpha) * _beta; 129 } 130 131 132 public override string ToString() 133 { 134 return string.Format("Weibull distribution ({0:####0.00000},{1:####0.00000})", _alpha, _beta); 135 } 136 137 138 public double Value(double x) 139 { 140 return _norm * Math.Pow(x, _alpha - 1) * Math.Exp(-Math.Pow(x / _beta, _alpha)); 141 } 142 143 144 public double Variance 145 { 146 get 147 { 148 double s = Fn.Gamma(1 / _alpha); 149 return _beta * _beta * (2 * Fn.Gamma(2 / _alpha) 150 - s * s / _alpha) / _alpha; 151 } 152 } 153 154 155 }
http://kwok.io/