求一对数曲线拟合算法或现成的免费的库

[varlj]

做一个东西，需要曲线拟合算法，这东西涉及到的数学知识对我而言有点难，网上找到的也多是对于直线的拟合，短时间估计搞不出来，又比较急
在网上找了不少地方，都没有合适的算法，麻烦大家帮一下忙
最终的曲线公式应该是
Y = A*LN(X)+B
需要求出A和B

如有多组数据Y（分别对应X从1到23）,大家可以在EXCEL里面用散点图求得公式
9999992.348
9999992.35
9999992.354
9999992.359
9999992.361
9999992.365
9999992.366
9999992.37
9999992.371
9999992.374
9999992.376
9999992.377
9999992.379
9999992.38
9999992.382
9999992.384
9999992.386
9999992.387
9999992.389
9999992.39
9999992.39
9999992.392
9999992.392

得到公式应该是y = 0.016Ln(x) + 1E7
即A=0.016 B=10000000

我的要求是：
求一算法，可以返回式中的A和B或者完整的公式
或者是网上的现成的，免费的库（不要破解的）


谢谢大家

[CCCCCCCCCCCCCCC]

那你还是直接用matlab或mathematica吧

[wuyq11]

C#与Matlab混合编程

[varlj]

matlab和mathematica都是收费的吧？有没有免费版本（我指的是没有日期限制的）

[varlj]

顶一个……网上找了半天，几乎全部是收费的

[gomoku]

其实你可以用直线拟合来做：

设X' = Ln(X)，则问题转化为
Y = AX' + B
求A和B是标准的直线拟合问题。

注：该方法求出来的AB，可以使得拟合点(Y,Ln(X))与实际点的总平方差最小，
但不见得能使拟合(Y,X)与实际(Y,X)的总平方差最小。

[varlj]

[Quote=引用 4 楼 gomoku 的回复:]
其实你可以用直线拟合来做：

设X' = Ln(X)，则问题转化为
Y = AX' + B
求A和B是标准的直线拟合问题。

注：该方法求出来的AB，可以使得拟合点(Y,Ln(X))与实际点的总平方差最小，
但不见得能使拟合(Y,X)与实际(Y,X)的总平方差最小。

[/Quote]

谢谢，因为这个函数是实物的数学模型，我试了一下，如果用直接来拟合，和对数拟合得到的结果相差有点大

[netstray]

ding

[varlj]

谢谢gomoku
真的是感慨万分，虽然一直知道数学及算法对程序的重要性，但是一直以来做的业务相关比较多，对于算法及数学难免有些不上心了
这拟合算法的说明网上其实也挺多，原理也介绍的挺清楚，但是以我的基础，根本没有信心去实现这个算法

再次谢谢gomoku及楼上各位

[varlj]

测试了几组数据，和EXCEL和其他专门的拟合工具得到的结果，在有效位上（小数后5位）完全一致

[gomoku]

昨天晚上写了一个“高斯-牛顿法”用于非线性拟合。

可惜中文wiki没有高斯-牛顿的解释，英文版的链接在这：Gauss–Newton algorithm

大概的说这是一种迭代方法，先随便估计一个初始参数，算出误差。对参数们进行偏导，估计它们对误差的影响，然后对估计参数在收敛误差的最大方向进行修正。

对那些多重峰值的情况，它可能会被不好的初始值拖入陷阱，只得到局部最优。但在大部分情况下，它还是能又快又好地得到最优解。（注：Levenberg–Marquardt法稳定性会更好一些）。

贴有三个文件，Test.cs，GaussNewton.cs，MyMatrix.cs。使用方法参考Test.cs，比如其中的PolynomialTest()用来解 y = Ax*x + Bx + C的方程，Main()用来解 y = Aln(x) + B的方程。

// Test.cs
//
using System;

class Test
{
    static void Main()
    {
        double[] Y = 
        {
            9999992.348,
            9999992.35,
            9999992.354,
            9999992.359,
            9999992.361,
            9999992.365,
            9999992.366,
            9999992.37,
            9999992.371,
            9999992.374,
            9999992.376,
            9999992.377,
            9999992.379,
            9999992.38,
            9999992.382,
            9999992.384,
            9999992.386,
            9999992.387,
            9999992.389,
            9999992.39,
            9999992.39,
            9999992.392,
            9999992.392
        };
        double[] X = new double[Y.Length];
        for (int i = 0; i < X.Length; i++) X[i] = i + 1;

        // f(x) = A ln(x) + B
        GaussNewton.F f = delegate(double[] coefficients, double x)
        {
            return coefficients[0] * Math.Log(x) + coefficients[1];
        };
        GaussNewton gaussNewton = new GaussNewton(2);
        gaussNewton.Initialize(Y, X, f);
        double[] answer = gaussNewton.Coefficients;     //answer[0] = 0.016 answer[1]=9999992.3386
    }

    static void PolynomialTest()
    {
        // y = -x*x + 2*x + 3
        double[] X = { 1, 2, 3,   8 };
        double[] Y = { 4, 3, 0, -45 };

        // f(x) = A*x*x + B*x + C
        GaussNewton.F f = delegate(double[] coefficients, double x)
        {
            return coefficients[0] * x * x + coefficients[1] * x + coefficients[2];
        };

        GaussNewton gaussNewton = new GaussNewton(3);
        gaussNewton.Initialize(Y, X, f);
        double[] answer = gaussNewton.Coefficients;     //A=-1 B=2 C=3
    }
}

[gomoku]

// GaussNewton.cs
//
using System;

public class GaussNewton
{
    public delegate double F(double[] coefficients, double x);

    /// <summary> Construct a GaussNewton solver </summary>
    public GaussNewton(int coefficientCount)
    {
        this.coefficientCount = coefficientCount;
        this.coefficients = new double[coefficientCount];
    }

    /// <summary> Initialize the solver without a guess. Y = f(X) </summary>
    public void Initialize(double[] Y, double[] X, F f)
    {
        Initialize(Y, X, f, null);
    }

    /// <summary> Initialize the solver: Y = f(X) with a guessed approximation</summary>
    public void Initialize(double[] Y, double[] X, F f, double[] coefficientGuess)
    {
        if (X == null || Y == null || f == null) throw new ArgumentNullException();
        if (X.Length != Y.Length) throw new ArgumentException("Y and X not in pairs");

        this.xs = X;
        this.ys = Y;
        this.function = f;
        this.solved = false;

        if (coefficientGuess != null && coefficientGuess.Length == this.coefficientCount)
        {
            Array.Copy(coefficientGuess, this.coefficients, this.coefficientCount);
        }
    }

    /// <summary> Get the solved coefficents </summary>
    /// <exception cref="InvalidOperationException">Throw when not answer can be found</exception>
    public double[] Coefficients
    {
        get { if (!this.solved) { Solve(); } return this.coefficients; }
    }

    /// <summary> Get the residual error (SSE) </summary>
    public double SumOfSquredError
    {
        get { return this.sse; }
    }

    /// <summary> Iteratively solve the coefficients using Gauss-Newton method </summary>
    /// <remarks> http://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm </remarks>
    public double[] Solve()
    {
        if( this.ys == null) throw new InvalidOperationException("Not yet initialized");

        double lastSSE = double.MaxValue;
        double[] errors = new double[ this.ys.Length ];

        // let C0 be the current coefficient guess, C1 be the better answer we are after
        // let E0 be the error using current guess
        // we have:
        // JacT * Jac * (C1 - C0) = JacT * E0
        //
        MyMatrix jacobian = Jacobian();
        MyMatrix jacobianT = jacobian.Transpose();
        MyMatrix product = jacobianT * jacobian;
        SquareMatrix inverse = SquareMatrix.FromMatrix(product).Inverse();

        for (int iteration = 0; iteration < GaussNewton.MaxIteration; iteration++)
        {
            this.sse = 0;

            for (int i = 0; i < this.ys.Length; i++)
            {
                double y = function(this.coefficients, this.xs[i]);
                errors[i] = this.ys[i] - y;
                sse += errors[i] * errors[i];
            }

            if (lastSSE - sse < GaussNewton.ConvergeThreshold)
            {
                this.solved = true;
                return this.coefficients;
            }

            double[] shift = inverse * (jacobianT * errors);

            for (int i = 0; i < this.coefficientCount; i++)
            {
                this.coefficients[i] += shift[i];
            }

            lastSSE = sse;
        }
        throw new InvalidOperationException("No answer can be found");
    }

    /// <summary> Calculate a Jacobian matrix. </summary>
    /// <remarks> http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant </remarks>
    private MyMatrix Jacobian()
    {
        double[][] p = new double[this.coefficientCount][];
        for (int i = 0; i < p.Length; i++)
        {
            p[i] = new double[this.coefficientCount];
            p[i][i] = 1;
        }

        MyMatrix jacobian = new MyMatrix(this.ys.Length, this.coefficientCount);
        for (int i = 0; i < this.ys.Length; i++)
        {
            for (int j = 0; j < this.coefficientCount; j++)
            {
                jacobian[i, j] = function(p[j], this.xs[i]);
            }
        }
        return jacobian;
    }

    private bool solved;
    private int coefficientCount;
    private double[] coefficients;
    private double[] xs;
    private double[] ys;
    private double sse;
    private F function;

    public static readonly int MaxIteration = 16;
    public static readonly double Epsilon = 1e-10;
    public static readonly double ConvergeThreshold = 1e-8;
}

[tzs2304]

up

[gomoku]

// MyMatrix.cs
//
using System;

public class MyMatrix
{
    private int cols, rows;
    protected double[,] data;

    /// <summary>Construct a matrix</summary>
    public MyMatrix(int rows, int cols)
    {
        this.rows = rows;
        this.cols = cols;
        this.data = new double[rows, cols];
    }

    /// <summary>Construct a matrix by copying the data</summary>
    public MyMatrix(double[,] data) : this( data.GetLength(0), data.GetLength(1) )
    {
        System.Buffer.BlockCopy(data, 0, this.data, 0, sizeof(double) * this.rows * this.cols);
    }

    /// <summary> Gets the number of rows in the matrix</summary>
    public int Rows { get { return this.rows; } }

    /// <summary> Gets the number of columns in the matrix</summary>
    public int Cols { get { return this.cols; } }

    /// <summary> Gets the raw data of the matrix. Should be treated as readonly</summary>
    public double[,] RawData { get { return this.data; } }

    /// <summary> Gets the sets the matrix element</summary>
    public double this[int row, int col]
    {
        get { return this.data[row, col]; }
        set { this.data[row, col] = value; }
    }

    /// <summary> Multiply by a vector. The width of the matrix should equal the height of the vector</summary>
    public double[] Multiply(double[] x)
    {
        if (x == null || x.Length != this.cols) throw new ArgumentException("Invalid vector");
        double[] result = new double[this.rows];

        for (int i = 0; i < result.Length; i++)
        {
            double sum = 0;
            for (int j = 0; j < this.cols; j++)
            {
                sum += this.data[i, j] * x[j];
            }
            result[i] = sum;
        }
        return result;
    }

    /// <summary> Multiply by another matrix. The width of left matrix should equal the height of right matrix </summary>
    public MyMatrix Multiply(MyMatrix m)
    {
        if (m == null || m.rows != this.cols) throw new ArgumentException("Invalid matrix to multiply");
        MyMatrix result = new MyMatrix(this.rows, m.cols);
        int inner = this.cols;

        for (int row = 0; row < result.rows; row++)
        {
            for (int col = 0; col < result.cols; col++)
            {
                double sum = 0;
                for (int i = 0; i < inner; i++) sum += this[row, i] * m[i, col];
                result[row, col] = sum;
            }
        }
        return result;
    }

    /// <summary> Create a transposed matrix </summary>
    public MyMatrix Transpose()
    {
        MyMatrix result = new MyMatrix(this.cols, this.rows);
        for (int row = 0; row < this.rows; row++)
        {
            for (int col = 0; col < this.cols; col++)
            {
                result[col, row] = this.data[row, col];
            }
        }
        return result;
    }

    public static double[] operator *(MyMatrix m, double[] vector)
    {
        return m.Multiply(vector);
    }

    public static MyMatrix operator *(MyMatrix m1, MyMatrix m2)
    {
        return m1.Multiply(m2);
    }
}

public class SquareMatrix : MyMatrix
{
    int rank;

    /// <summary>Construct an identity matrix</summary>
    public SquareMatrix(int rank) : base(rank, rank)
    {
        this.rank = rank;
        for (int i = 0; i < rank; i++) this.data[i, i] = 1;
    }

    /// <summary>Construct a square matrix by copying the data</summary>
    public SquareMatrix(double[,] data) : base(data)
    {
        if ( data.GetLength(0) != data.GetLength(1))
        {
            throw new ArgumentException("Only square matrix supported");
        }
        this.rank = data.GetLength(0);
    }

    /// <summary>Use Laplace expansion to calculate the determinant of a square matrix</summary>
    public double Determinant()
    {
        if (this.rank == 1)
        {
            return data[0, 0];
        }
        if (this.rank == 2)
        {
            return data[0, 0] * data[1, 1] - data[1, 0] * data[0, 1];
        }

        double det = 0;
        for (int i = 0; i < this.rank; i++)
        {
            det += this.data[i, 0] * this.Cofactor(i, 0).Determinant() * (i % 2 == 0 ? 1 : -1);
        }
        return det;
    }

    /// <summary>Use Cramer's rule to recursively find the inverse of a square matrix</summary>
    public SquareMatrix Inverse()
    {
        double det = this.Determinant();
        if (Math.Abs(det) < 1e-8) throw new InvalidOperationException("Cannot inverse a singular matrix");

        SquareMatrix result = new SquareMatrix(this.rank);
        for (int i = 0; i < this.rank; i++)
        {
            for (int j = 0; j < this.rank; j++)
            {
                result.data[j, i] = this.Cofactor(i, j).Determinant() * ((i + j) % 2 == 0 ? 1 : -1) / det;
            }
        }
        return result;
    }

    /// <summary>Get a Cofactor matrix </summary>
    public SquareMatrix Cofactor(int i, int j)
    {
        if (this.rank < 2) throw new InvalidOperationException("Rank is less than two");

        SquareMatrix result = new SquareMatrix(this.rank - 1);
        double[,] buf = result.data;
        for (int y = 0; y < this.rank; y++)
        {
            for (int x = 0; x < this.rank; x++)
            {
                if (y != i && x != j)
                {
                    buf[y - (y > i ? 1 : 0), x - (x > j ? 1 : 0)] = this.data[y, x];
                }
            }
        }
        return result;
    }

    /// <summary> Create a square matrix from a general matrix</summary>
    public static SquareMatrix FromMatrix(MyMatrix m)
    {
        if (m == null || m.Cols != m.Rows) throw new ArgumentException("Not a valid square matrix");
        return new SquareMatrix(m.RawData);
    }
}

[lvdazhuang]

顶了，还有什么其他的算法没有？

[超级大笨狼]

已经揭贴了，但是不妨碍我传播技巧。

用matlab+C#做这个事情比较明智。

我现在用这个办法做多元线性回归，就一行代码regress(Y,X)
就可以求出Y=CC+AX!+BX2+CX3....中的系数

曲线拟合也有现成的例子，只要找本书大概看一下，什么牛顿法之类的都可以秒杀。

matlab+C#百度视频就有，再搜一下网页，很容易做到。

dll加载的时候慢，需要提前加载，计算的时候很快。
比自己写C#代码速度要快。
因为我们可能只是用C#实现了，未必去优化它。

[超级大笨狼]

关于拟合和插值，我觉得还是利用现成的库比较好，因为matlab有多种现成的拟合函数，你可以逐一实验，快速挑选最适合的。还可以画出图来。而且，正确性和速度性都有保证。

用C#，相当于自己发明轮子，会给你带来多少个不眠之夜啊！~~

[超级大笨狼]

赞一下，gomoku的技术和知识！~~

[SimpleGIS]

gomoku写得确实不错，我也和EXCEL比对了，确实一样
不过我用它来做类似y = 231.3e0.049x的拟合时，就报错了，麻烦gomoku再看一下呢，难道你的类就运行对数拟合和二元一次拟合？

[IT老人]

好帖子，正好用到，谢谢