C#最小二乘法进行曲线拟合及相关系数
两个类：
类1：
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) + "\r\n");
//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) + "\r\n");
//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) + "\r\n");
//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) + "\r\n");
//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) + "\r\n");
//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();
}
}
}
类2：
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
}
}
C#最小二乘法进行曲线拟合及相关系数相关推荐
- 最小二乘法多项式曲线拟合原理与实现--转
原文地址:http://blog.csdn.net/jairuschan/article/details/7517773/ 概念 最小二乘法多项式曲线拟合,根据给定的m个点,并不要求这条曲线精确地经过 ...
- 最小二乘法多项式曲线拟合数学原理及其C++实现
目录 0 前言 1 最小二乘法概述 2 最小二乘法求解多项式曲线系数向量的数学推导 2.1 代数法 2.2 矩阵法 3 代码实现 4 总结 参考 0 前言 自动驾驶开发中经常涉及到多项式曲线拟合,本文 ...
- 最小二乘法多项式曲线拟合原理与实现
概念 最小二乘法多项式曲线拟合,根据给定的m个点,并不要求这条曲线精确地经过这些点,而是曲线y=f(x)的近似曲线y= φ(x). 原理 [原理部分由个人根据互联网上的资料进行总结,希望对大家能有用] ...
- python最小二乘法拟合三维曲线_python_numpy最小二乘法的曲线拟合
在了解了最小二乘法的基本原理之后python_numpy实用的最小二乘法理解,就可以用最小二乘法做曲线拟合了 1.直线拟合 直线拟合 已知图中拟合数据的坐标,对图中的拟合数据进行直线拟合. 依旧使用最 ...
- 最小二乘法多项式曲线拟合原理与实现（错误地方已经修改底层补充自己写的java实现）
目录(?) [-] 概念 原理 运行前提 代码 运行效果 概念 最小二乘法多项式曲线拟合,根据给定的m个点,并不要求这条曲线精确地经过这些点,而是曲线y=f(x)的近似曲线y= φ(x). 原理 [原 ...
- c++ 三次多项式拟合_最小二乘法多项式曲线拟合数学原理及其C++实现
本文使用 Zhihu On VSCode 创作并发布 0 前言 自动驾驶开发中经常涉及到多项式曲线拟合,本文详细描述了使用最小二乘法进行多项式曲线拟合的数学原理,通过样本集构造范德蒙德矩阵,将一元 N ...
- 最小二乘法多项式曲线拟合及其python实现
最小二乘法多项式曲线拟合及其python实现 多项式曲线拟合问题描述 最小二乘法 针对overfitting,加入正则项 python实现 运行结果 多项式曲线拟合问题描述 问题描述:给定一些数据点, ...
- 最小二乘法进行曲线拟合
工作需求,这里记录一下数值插值和数值分析方面的算法,希望和大家一起进步. 曲线拟合的最小二乘定义 求一条曲线,使数据点均在离此曲线的上方或下方不远处,所求的曲线称为拟合曲线, 它既能反映数据的总体分布 ...
- 最小二乘法实现曲线拟合
说明,本文章的源代码来着于网络,本人已在实际项目中反复使用过,证明没问题. 1.简介 已知曲线上的n个点,可以使用某条曲线去拟合,使得整体上所有的点都逼近曲线,可以使用不同的角度去判断整体逼近,最小二 ...
最新文章
- 【转】【iOS知识学习】_视图控制对象生命周期-init、viewDidLoad、viewWillAppear、viewDidAppear、viewWillDisappear等的区别及用途...
- 设置textview背景色为透明
- Subversion之路－－实现精细的目录访问权限控制（v1.0 更新于2006.12.05）(二)
- tomcat5下jsp出现getOutputStream() has already been called for this response异常的原因和解决方法...
- 前端之同源策略 Jsonp 与 CORS
- Master-Detail（主表明细），确认可以出货的SQL指令 -- Not Exists
- python 调用mysql_Python调用Mysql
- sql删除语句_推荐强大开源的数据库SQL语句审核平台，再也不用担心删除跑路了！...
- 光源时间_D65光源对色灯箱的操作步骤及作业标准
- Flink从入门到入土
- erlang的简单模拟半包的产生
- AD排错最佳实践—— 技巧与捷径-笔记
- 正则表达式 两个符号的字段_Tableau正则提取字段部分内容
- [原创]独立模式安装Hive
- 易语言解析html实例,易语言完整示例（单设备）
- 【QT】QCustomPlot图表控件
- 用什么系统搭建nas服务器,自己家里搭建NAS服务器有什么好方案
- 弱水三千只取一瓢，Forcepoint的变与不变
- Shell脚本之shift用法
- mysql查询top10_各个数据库中TOP10记录的查询方法