zoukankan      html  css  js  c++  java
  • C# ——计算线性关系kb值、R平方,类似于excel的趋势线线性关系功能

    这些功能Excel上都有,原理一模一样,现在需要C#的实现代码;
    各函数的线性拟合,相关系数、截距为0(即强制过原点)等等

    源码:https://download.csdn.net/download/a1062484747/15769519

    相关系数R²的公式引用:http://blog.csdn.net/huwei2003/article/details/18553775(验证过)
    1.一次线性、二次曲线、指数、对数、幂等函数拟合及相关系数R²的代码实现(指数函数拟合的相关系数R²和Excel有出入);
    2.一次线性的截距为0(即强制过原点)的代码实现;
    3.代码三次乃至多项以上的函数拟合有问题,不会改,望有大神补充修改一下;
    4.有没有大神补充一下二次曲线、指数这2个函数拟合时截距为0(即强制过原点)的拟合代码或者数学公式。

    代码实现在Excel验证过!

    代码如下:

    C#
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
     
    namespace 高斯消元法
    {
        class Program
        {
            static void Main(string[] args)
            {
               /* double[,] xArray = new double[,]
                {
                     
                        { 2.000000 ,-1.000000 , 3.000000,  1.000000},
                        { 4.000000 , 2.000000 , 5.000000,  4.000000},
                        { 1.000000 , 2.000000 , 0.000000 , 7.000000}
                };*/
     
                System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
                double[] y = new double[] { 29152.3, 47025.3, 86852.3, 132450.6, 200302.3, 284688.1, 396988.3 };
                double[] x = new double[] { 1.24, 2.37, 5.12, 8.12, 12.19, 17.97, 24.99 };
     
               // double[,] xArray;
                double[] ratio;
                double[] yy = new double[y.Length];
     
                Console.WriteLine("一次拟合:");
                sw.Start();
                ratio = FittingFunct.Linear(y, x);
                sw.Stop();
     
                foreach (double num in ratio)
                {
                    Console.WriteLine(num);
                }
                for (int i = 0; i < x.Length; i++)
                {
                    yy[i] = ratio[0] + ratio[1] * x[i];
                }
                Console.WriteLine("R²=: " + FittingFunct.Pearson(y, yy) + "
    ");
                //Console.WriteLine("一次拟合计算时间:");
                //Console.WriteLine(sw.ElapsedMilliseconds);
     
                Console.WriteLine("一次拟合(截距为0,即强制过原点):");
                sw.Start();
                ratio = FittingFunct.LinearInterceptZero(y, x);
                sw.Stop();
     
                foreach (double num in ratio)
                {
                    Console.WriteLine(num);
                }
                for (int i = 0; i < x.Length; i++)
                {
                    yy[i] = ratio[0] * x[i];
                }
                Console.WriteLine("R²=: " + FittingFunct.Pearson(y, yy) + "
    ");
                //Console.WriteLine("一次拟合计算时间:");
                //Console.WriteLine(sw.ElapsedMilliseconds);
     
                Console.WriteLine("二次拟合:");
                sw.Start();
                ratio = FittingFunct.TowTimesCurve(y, x);
                sw.Stop();
     
                foreach (double num in ratio)
                {
                    Console.WriteLine(num);
                }
                for (int i = 0; i < x.Length; i++)
                {
                    yy[i] = ratio[0] + ratio[1] * x[i] + ratio[2] * x[i] * x[i];
                }
                Console.WriteLine("R²=: " + FittingFunct.Pearson(y, yy) + "
    ");
                //Console.WriteLine("二次拟合计算时间:");
                //Console.WriteLine(sw.ElapsedMilliseconds);
     
                Console.WriteLine("对数拟合计算时间:");
                sw.Start();
                ratio = FittingFunct.LOGEST(y, x);
                sw.Stop();
     
                foreach (double num in ratio)
                {
                    Console.WriteLine(num);
                }
                for (int i = 0; i < x.Length; i++)
                {
                    yy[i] = ratio[1]*Math.Log10(x[i]) + ratio[0];
                }
                Console.WriteLine("R²=: " + FittingFunct.Pearson(y, yy) + "
    ");
                //Console.WriteLine("对数拟合计算时间:");
                //Console.WriteLine(sw.ElapsedMilliseconds);
     
                Console.WriteLine("幂级数拟合:");
                sw.Start();
                ratio=FittingFunct.PowEST(y,x);
                sw.Stop();
     
                 foreach (double num in ratio)
                {
                    Console.WriteLine(num);
                }
                 for (int i = 0; i < x.Length; i++)
                 {
                     yy[i] = ratio[0] * Math.Pow(x[i], ratio[1]);
                 }
                 Console.WriteLine("R²=: " + FittingFunct.Pearson(y, yy) + "
    ");
                 //Console.WriteLine("幂级数拟合计算时间:");
                 //Console.WriteLine(sw.ElapsedMilliseconds);
     
                Console.WriteLine("指数函数拟合:");
                sw.Start();
                ratio = FittingFunct.IndexEST(y, x);
                sw.Stop();
                foreach (double num in ratio)
                {
                    Console.WriteLine(num);
                }
                for (int i = 0; i < x.Length; i++)
                {
                    yy[i] = ratio[0] * Math.Exp(x[i] * ratio[1]);
                }
                Console.WriteLine("R²=: " + FittingFunct.Pearson(y, yy));
                //Console.WriteLine("指数函数拟合计算时间:");
                //Console.WriteLine(sw.ElapsedMilliseconds);
     
                Console.ReadKey();
            
            }
        }
    }
    C#
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
     
    namespace 高斯消元法
    {
        class FittingFunct
        {
              #region 多项式拟合函数,输出系数是y=a0+a1*x+a2*x*x+.........,按a0,a1,a2输出
            static public double[] Polyfit(double[] y, double[] x, int order)
            {
                  double[,] guass = Get_Array(y, x, order);
                
                  double[] ratio = Cal_Guass(guass, order + 1);
                 
                  return ratio;
            }
            #endregion
     
              #region 一次拟合函数,y=a0+a1*x,输出次序是a0,a1
            static public double[] Linear(double[] y,double[] x)
            {
                double[] ratio = Polyfit(y, x, 1);
                return ratio;
            }
             #endregion
     
            #region 一次拟合函数,截距为0,y=a0x,输出次序是a0
            static public double[] LinearInterceptZero(double[] y, double[] x)
            {
                double divisor = 0; //除数
                double dividend = 0; //被除数
                for (int i = 0; i < x.Length;i++ )
                {
                    divisor += x[i] * x[i];
                    dividend += x[i] * y[i];
                }
                if (divisor == 0)
                {
                    throw (new Exception("除数不为0!"));
                    return null;
                }
                return new double[] { dividend / divisor };
     
            }
            #endregion
     
            #region 二次拟合函数,y=a0+a1*x+a2x²,输出次序是a0,a1,a2
            static public double[] TowTimesCurve(double[] y, double[] x)
            {
                double[] ratio = Polyfit(y, x, 2);
                return ratio;
            }
            #endregion
     
              #region 对数拟合函数,.y= c*(ln x)+b,输出为b,c
            static public double[] LOGEST(double[] y, double[] x)
            {
                double[] lnX = new double[x.Length];
     
                for (int i = 0; i < x.Length; i++)
                {
                    if (x[i] == 0 || x[i] < 0)
                    {
                        throw (new Exception("正对非正数取对数!"));
                        return null;
                    }
                    lnX[i] = Math.Log(x[i]);
                }
     
                return Linear(y, lnX);
            }
            #endregion
     
              #region 幂函数拟合模型, y=c*x^b,输出为c,b
            static public double[] PowEST(double[] y, double[] x)
            {
                double[] lnX = new double[x.Length];
                double[] lnY = new double[y.Length];
                double[] dlinestRet;
     
                for (int i = 0; i < x.Length; i++)
                {
                    lnX[i] = Math.Log(x[i]);
                    lnY[i] = Math.Log(y[i]);
                }
     
                dlinestRet = Linear(lnY, lnX);
     
               dlinestRet[0] = Math.Exp(dlinestRet[0]);
                
                return dlinestRet;
            }
             #endregion
     
             #region 指数函数拟合函数模型,公式为 y=c*m^x;输出为 c,m
             static  public double[] IndexEST(double[] y, double[] x)
            {
                double[] lnY = new double[y.Length];
                double[] ratio;
                for (int i = 0; i < y.Length; i++)
                {
                    lnY[i] = Math.Log(y[i]);
                 }
     
                ratio = Linear(lnY, x);
                for (int i = 0; i < ratio.Length; i++)
                {
                    if (i == 0)
                    {
                        ratio[i] = Math.Exp(ratio[i]);
                    }
                 }
                return ratio;
            }
             #endregion
     
             #region 相关系数R²部分
             public static double Pearson(IEnumerable<double> dataA, IEnumerable<double> dataB)
             {
                 int n = 0;
                 double r = 0.0;
     
                 double meanA = 0;
                 double meanB = 0;
                 double varA = 0;
                 double varB = 0;
                 int ii = 0;
                 using (IEnumerator<double> ieA = dataA.GetEnumerator())
                 using (IEnumerator<double> ieB = dataB.GetEnumerator())
                 {
                     while (ieA.MoveNext())
                     {
                         if (!ieB.MoveNext())
                         {
                             //throw new ArgumentOutOfRangeException("dataB", Resources.ArgumentArraysSameLength);
                         }
                         ii++;
                         //Console.WriteLine("FF00::  " + ii + " --  " + meanA + " -- " + meanB + " -- " + varA + "  ---  " + varB);
                         double currentA = ieA.Current;
                         double currentB = ieB.Current;
     
                         double deltaA = currentA - meanA;
                         double scaleDeltaA = deltaA / ++n;
     
                         double deltaB = currentB - meanB;
                         double scaleDeltaB = deltaB / n;
     
                         meanA += scaleDeltaA;
                         meanB += scaleDeltaB;
     
                         varA += scaleDeltaA * deltaA * (n - 1);
                         varB += scaleDeltaB * deltaB * (n - 1);
                         r += (deltaA * deltaB * (n - 1)) / n;
                         //Console.WriteLine("FF00::  " + ii + " --  " + meanA + " -- " + meanB + " -- " + varA + "  ---  " + varB);
                     }
     
                     if (ieB.MoveNext())
                     {
                         //throw new ArgumentOutOfRangeException("dataA", Resources.ArgumentArraysSameLength);
                     }
                 }
                 return (r / Math.Sqrt(varA * varB)) * (r / Math.Sqrt(varA * varB));
             }
             #endregion 
     
    #region 最小二乘法部分
     
              #region 计算增广矩阵
            static  private double[] Cal_Guass(double[,] guass,int count)
            {
                double temp;
                double[] x_value;
     
                for (int j = 0; j < count; j++)
                {
                    int k = j;
                    double min = guass[j,j];
     
                    for (int i = j; i < count; i++)
                    {
                        if (Math.Abs(guass[i, j]) < min)
                        {
                            min = guass[i, j];
                            k = i;
                        }
                    }
     
                    if (k != j)
                    {
                        for (int x = j; x <= count; x++)
                        {
                            temp = guass[k, x];
                            guass[k, x] = guass[j, x];
                            guass[j, x] = temp;
                        }
                    }
     
                    for (int m = j + 1; m < count; m++)
                    {
                        double div = guass[m, j] / guass[j, j];
                        for (int n = j; n <= count; n++)
                        {
                            guass[m, n] = guass[m, n] - guass[j, n] * div;
                        }
                    }
     
                   /* System.Console.WriteLine("初等行变换:");
                    for (int i = 0; i < count; i++)
                    {
                        for (int m = 0; m < count + 1; m++)
                        {
                            System.Console.Write("{0,10:F6}", guass[i, m]);
                        }
                        Console.WriteLine();
                    }*/
                }
                x_value= Get_Value(guass, count);
     
                return x_value;
     
                /*if (x_value == null)
                    Console.WriteLine("方程组无解或多解!");
                else
                {
                    foreach (double x in x_value)
                    {
                        Console.WriteLine("{0:F6}", x);
                    }
                }*/
            }
     
            #endregion
     
              #region 回带计算X值
            static private double[] Get_Value(double[,] guass,int count)
            {
                double[] x = new double[count];
                double[,] X_Array = new double[count, count];
                int rank = guass.Rank;//秩是从0开始的
     
                for (int i = 0; i < count; i++)
                    for (int j = 0; j < count; j++)
                        X_Array[i, j] = guass[i, j];
     
                if (X_Array.Rank < guass.Rank)//表示无解
                {
                    return null;
                }
     
                if (X_Array.Rank < count-1)//表示有多解
                {
                    return null;
                }
                //回带计算x值
                x[count - 1] = guass[count - 1, count] / guass[count-1, count-1];
                for (int i = count - 2; i >= 0; i--)
                {
                    double temp=0;
                    for (int j = i; j < count; j++)
                    {
                        temp += x[j] * guass[i, j];
                    }
                    x[i] = (guass[i, count] - temp) / guass[i, i];
                }
     
                    return x;
            }
              #endregion
     
              #region  得到数据的法矩阵,输出为发矩阵的增广矩阵
            static private double[,] Get_Array(double[] y, double[] x, int n)
            {
                double[,] result = new double[n + 1, n + 2];
     
                if (y.Length != x.Length)
                {
                    throw (new Exception("两个输入数组长度不一!"));
                    //return null;
                }
     
                for (int i = 0; i <= n; i++)
                {
                    for (int j = 0; j <= n; j++)
                    {
                        result[i, j] = Cal_sum(x, i + j);
                    }
                    result[i, n + 1] = Cal_multi(y, x, i);
                }
     
                return result;
            }
     
         #endregion
     
              #region 累加的计算
            static private double Cal_sum(double[] input,int order)
            {
                double result=0;
                int length = input.Length;          
     
                for (int i = 0; i < length; i++)
                {
                    result += Math.Pow(input[i], order);
                }
     
               return result;
            }
            #endregion
     
              #region 计算∑(x^j)*y
            static private double Cal_multi(double[] y,double[] x,int order)
            {
                double result = 0;
     
                int length = x.Length;
     
                for (int i = 0; i < length; i++)
                {
                    result += Math.Pow(x[i], order) * y[i];
                }
     
                return result;
            }
             #endregion
     
    #endregion
        }
    }

     转http://www.skcircle.com/?id=362

  • 相关阅读:
    js自动小轮播
    js字符串
    工资
    可是姑娘,你为什么要编程呢?
    程序猿媳妇儿注意事项
    js勾选时显示相应内容
    js点击显示隐藏
    js选项卡
    js数组
    js旋转V字俄罗斯方块
  • 原文地址:https://www.cnblogs.com/eve612/p/14528500.html
Copyright © 2011-2022 走看看