From e64adc609d1cc0612851815d37f7d6759268a28c Mon Sep 17 00:00:00 2001
From: Jon Larborn <jon.larborn@outlook.com>
Date: Wed, 29 Jul 2015 14:10:45 +0200
Subject: [PATCH] Estimate extention

---
 src/Numerics/Distributions/Weibull.cs | 44 +++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/src/Numerics/Distributions/Weibull.cs b/src/Numerics/Distributions/Weibull.cs
index 896cb3cb..7acbb36f 100644
--- a/src/Numerics/Distributions/Weibull.cs
+++ b/src/Numerics/Distributions/Weibull.cs
@@ -420,6 +420,50 @@ namespace MathNet.Numerics.Distributions
             return -SpecialFunctions.ExponentialMinusOne(-Math.Pow(x, shape)*Math.Pow(scale, -shape));
         }
 
+        /// <summary>
+        /// Returns a Weibull distribution.
+        /// 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>
+        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>