Merge commit '4966036b86efe39c7a96db67b375f2d7467fa6c1'

cdrnet · cdrnet · commit 320c1cc1d1e7 · 2015-08-09T18:54:53.000+02:00
diff --git a/src/Numerics/Distributions/Weibull.cs b/src/Numerics/Distributions/Weibull.cs
@@ -420,6 +420,49 @@ public static double CDF(double shape, double scale, double x)
             return -SpecialFunctions.ExponentialMinusOne(-Math.Pow(x, shape)*Math.Pow(scale, -shape));
         }
 
+        /// <summary>
+        /// Implemented according to: Parameter estimation of the Weibull probability distribution, 1994, Hongzhu Qiao, Chris P. Tsokos 
+        /// </summary>
+        /// <param name="samples"></param>
+        /// <param name="randomSource"></param>
+        /// <returns>Returns a Weibull distribution.</returns>
+        public static Weibull Estimate(IEnumerable<double> samples, System.Random randomSource = null)
+        {
+            var samp = samples as double[] ?? samples.ToArray();
+            double n = samp.Count(), s1 = 0, s2 = 0, s3 = 0, previousC = Int32.MinValue, QofC = 0;
+
+            if (n <= 1) throw new Exception("Observations not sufficient");
+
+            // Start values
+            double c = 10; double b = 0;
+
+            while (Math.Abs(c - previousC) >= 0.0001)
+            {
+                s1 = s2 = s3 = 0;
+                foreach (double x in samp)
+                {
+                    if (x > 0)
+                    {
+                        s1 += Math.Log(x);
+                        s2 += Math.Pow(x, c);
+                        s3 += Math.Pow(x, c) * Math.Log(x);
+                    }
+                }
+                QofC = n * s2 / (n * s3 - s1 * s2);
+
+                previousC = c;
+                c = (c + QofC) / 2;
+            }
+
+            foreach (double x in samp)
+                if (x > 0)
+                    b += Math.Pow(x, c);
+
+            b = Math.Pow(b / n, 1 / c);
+
+            return new Weibull(c, b, randomSource);
+        }
+
         /// <summary>
         /// Generates a sample from the Weibull distribution.
         /// </summary>
