灰色预测模型 c# 算法实现

public class GrayModel    {        private double a0, a1, a2;        private int size;        private double error;        public GrayModel()        {        }        public void build(double[] x0)        {            size = x0.Length;            double[] x1 = new double[size];            x1[0] = x0[0];            for (int i = 1; i < size; i++)            {                x1[i] = x0[i] + x1[i - 1];            }            double[,] b = new double[size - 1, 2];            double[,] bt = new double[2, size - 1];            double[,] y = new double[size - 1, 1];            for (int i = 0; i < b.GetLength(0); i++)            {                b[i, 0] = -(x1[i] + x1[i + 1]) / 2;                b[i, 1] = 1;                bt[0, i] = b[i, 0];                bt[1, i] = 1;                y[i, 0] = x0[i + 1];            }            double[,] t = new double[2, 2];            multiply(bt, b, t);            t = inverse(t);            double[,] t1 = new double[2, size - 1];            multiply(t, bt, t1);            double[,] t2 = new double[2, 1];            multiply(t1, y, t2);            a0 = t2[0, 0];            double u = t2[1, 0];            a2 = u / a0;            a1 = x0[0] - a2;            a0 = -a0;            error = 0;            for (int i = 0; i < x0.Length; i++)            {                double d = (x0[i] - getX0(i));                error += d * d;            }            error /= x0.Length;        }        /// <summary>        /// 误差        /// </summary>        /// <returns></returns>        public double getError()        {            return error;        }        double getX1(int k)        {            return a1 * Math.Exp(a0 * k) + a2;        }        double getX0(int k)        {            // return a0 * a1 * Math.exp(a0 * k);              if (k == 0)                return a1 * Math.Exp(a0 * k) + a2;            else                return a1 * (Math.Exp(a0 * k) - Math.Exp(a0 * (k - 1)));        }       /// <summary>       /// 预测后续的值       /// </summary>       /// <param name="index"></param>       /// <returns></returns>        public double nextValue(int index)        {            if (index < 0)                throw new Exception("超出索引范围");            return getX0(size + index);        }       /// <summary>       /// 预测下一个值       /// </summary>       /// <returns></returns>        public double nextValue()        {            return nextValue(0);        }        static double[,] inverse(double[,] t)        {            double[,] a = new double[2, 2];            double det = t[0, 0] * t[1, 1] - t[0, 1] * t[1, 0];            a[0, 0] = t[1, 1] / det;            a[0, 1] = -t[1, 0] / det;            a[1, 0] = -t[0, 1] / det;            a[1, 1] = t[0, 0] / det;            return a;        }        static void multiply(double[,] left, double[,] right, double[,] dest)        {            int n1 = left.GetLength(0);            int m1 = left.GetLength(1);            int m2 = right.GetLength(1);            for (int k = 0; k < n1; k++)            {                for (int s = 0; s < m2; s++)                {                    dest[k, s] = 0;                    for (int i = 0; i < m1; i++)                    {                        dest[k, s] += left[k, i] * right[i, s];                    }                }            }        }    }

https://www.cnblogs.com/njcxwz/p/5198187.html

http://cdmd.cnki.com.cn/Article/CDMD-84502-2009157569.htm

static void Main(string[] args)
        {
           double[] a=new double[5]{2.874,3.278,3.337,3.390,3.679};
           GrayModel gm = new GrayModel();
           gm.build(a);
            Console.WriteLine(gm.nextValue());
            Console.Read();
        }

public class GrayModel
    {
        private double a0, a1, a2;
        private int size;
        private double error;

        public GrayModel()
        {
        }

        public void build(double[] x0)
        {
            size = x0.Length;
            double[] x1 = new double[size];
            x1[0] = x0[0];
            for (int i = 1; i < size; i++)
            {
                x1[i] = x0[i] + x1[i - 1];
            }
            double[,] b = new double[size - 1, 2];
            double[,] bt = new double[2, size - 1];
            double[,] y = new double[size - 1, 1];
            for (int i = 0; i < b.GetLength(0); i++)
            {
                b[i, 0] = -(x1[i] + x1[i + 1]) / 2;
                b[i, 1] = 1;
                bt[0, i] = b[i, 0];
                bt[1, i] = 1;
                y[i, 0] = x0[i + 1];
            }
            double[,] t = new double[2, 2];
            multiply(bt, b, t);
            t = inverse(t);
            double[,] t1 = new double[2, size - 1];
            multiply(t, bt, t1);
            double[,] t2 = new double[2, 1];
            multiply(t1, y, t2);
            a0 = t2[0, 0];
            double u = t2[1, 0];
            a2 = u / a0;
            a1 = x0[0] - a2;
            a0 = -a0;

            error = 0;
            for (int i = 0; i < x0.Length; i++)
            {
                double d = (x0[i] - getX0(i));
                error += d * d;
            }
            error /= x0.Length;
        }

        /// <summary>
        /// 误差
        /// </summary>
        /// <returns></returns>
        public double getError()
        {
            return error;
        }

        double getX1(int k)
        {
            return a1 * Math.Exp(a0 * k) + a2;
        }

        double getX0(int k)
        {
            // return a0 * a1 * Math.exp(a0 * k);  
            if (k == 0)
                return a1 * Math.Exp(a0 * k) + a2;
            else
                return a1 * (Math.Exp(a0 * k) - Math.Exp(a0 * (k - 1)));
        }

       /// <summary>
       /// 预测后续的值
       /// </summary>
       /// <param name="index"></param>
       /// <returns></returns>
        public double nextValue(int index)
        {
            if (index < 0)
                throw new Exception("超出索引范围");
            return getX0(size + index);
        }

       /// <summary>
       /// 预测下一个值
       /// </summary>
       /// <returns></returns>
        public double nextValue()
        {
            return nextValue(0);
        }

        static double[,] inverse(double[,] t)
        {
            double[,] a = new double[2, 2];
            double det = t[0, 0] * t[1, 1] - t[0, 1] * t[1, 0];
            a[0, 0] = t[1, 1] / det;
            a[0, 1] = -t[1, 0] / det;
            a[1, 0] = -t[0, 1] / det;
            a[1, 1] = t[0, 0] / det;
            return a;
        }

        static void multiply(double[,] left, double[,] right, double[,] dest)
        {
            int n1 = left.GetLength(0);
            int m1 = left.GetLength(1);
            int m2 = right.GetLength(1);
            for (int k = 0; k < n1; k++)
            {
                for (int s = 0; s < m2; s++)
                {
                    dest[k, s] = 0;
                    for (int i = 0; i < m1; i++)
                    {
                        dest[k, s] += left[k, i] * right[i, s];
                    }
                }
            }
        }
    }