开启辅助访问

amobbs.com 阿莫电子技术论坛

登录 注册 找回密码
搜索
amobbs.com 阿莫电子技术论坛»论坛首页 › 单片机 › DSP › 算法——纯C语言最小二乘法曲线拟合
bottom↓
返回列表 发新帖
回复: 93

算法——纯C语言最小二乘法曲线拟合

  [复制链接]
zouzhichao

出10入23汤圆
发表于 2013-5-20 16:45:26 | 显示全部楼层 |阅读模式
纯手写，刚写完，还没来得及写注释，已通过Matlab的polyfit验证（阶数高或者数据量太大会有double数据溢出的危险，低阶的都吻合），时间有点紧，程序注释，数学推导等时间充裕一点再贴上来
// 最小二乘法拟合.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"
//
#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))
//
/***********************************************************************************
***********************************************************************************/
static int GetXY(const char* FileName, double* X, double* Y, int* Amount)
{
        FILE* File = fopen(FileName, "r");
        if (!File)
                return -1;
        for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)
                if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))
                        break;
        fclose(File);
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int PrintPara(double* Para, int SizeSrc)
{
        int i, j;
        for (i = 0; i < SizeSrc; i++)
        {
                for (j = 0; j <= SizeSrc; j++)
                        printf("%10.6lf ", ParaBuffer(Para, i, j));
                printf("\r\n");
        }
        printf("\r\n");
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int ParalimitRow(double* Para, int SizeSrc, int Row)
{
        int i;
        double Max, Min, Temp;
        for (Max = abs(ParaBuffer(Para, Row, 0)), Min = Max, i = SizeSrc; i; i--)
        {
                Temp = abs(ParaBuffer(Para, Row, i));
                if (Max < Temp)
                        Max = Temp;
                if (Min > Temp)
                        Min = Temp;
        }
        Max = (Max + Min) * 0.000005;
        for (i = SizeSrc; i >= 0; i--)
                ParaBuffer(Para, Row, i) /= Max;
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int Paralimit(double* Para, int SizeSrc)
{
        int i;
        for (i = 0; i < SizeSrc; i++)
                if (ParalimitRow(Para, SizeSrc, i))
                        return -1;
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int ParaPreDealA(double* Para, int SizeSrc, int Size)
{
        int i, j;
        for (Size -= 1, i = 0; i < Size; i++)
        {
                for (j = 0; j < Size; j++)
                        ParaBuffer(Para, i, j) = ParaBuffer(Para, i, j) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, j) * ParaBuffer(Para, i, Size);
                ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, SizeSrc) * ParaBuffer(Para, i, Size);
                ParaBuffer(Para, i, Size) = 0;
                ParalimitRow(Para, SizeSrc, i);
        }
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int ParaDealA(double* Para, int SizeSrc)
{
        int i;
        for (i = SizeSrc; i; i--)
                if (ParaPreDealA(Para, SizeSrc, i))
                        return -1;
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int ParaPreDealB(double* Para, int SizeSrc, int OffSet)
{
        int i, j;
        for (i = OffSet + 1; i < SizeSrc; i++)
        {
                for (j = OffSet + 1; j <= i; j++)
                        ParaBuffer(Para, i, j) *= ParaBuffer(Para, OffSet, OffSet);
                ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, OffSet, OffSet) - ParaBuffer(Para, i, OffSet) * ParaBuffer(Para, OffSet, SizeSrc);
                ParaBuffer(Para, i, OffSet) = 0;
                ParalimitRow(Para, SizeSrc, i);
        }
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int ParaDealB(double* Para, int SizeSrc)
{
        int i;
        for (i = 0; i < SizeSrc; i++)
                if (ParaPreDealB(Para, SizeSrc, i))
                        return -1;
        for (i = 0; i < SizeSrc; i++)
        {
                if (ParaBuffer(Para, i, i))
                {
                        ParaBuffer(Para, i, SizeSrc) /= ParaBuffer(Para, i, i);
                        ParaBuffer(Para, i, i) = 1.0;
                }
        }
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int ParaDeal(double* Para, int SizeSrc)
{
        PrintPara(Para, SizeSrc);
        Paralimit(Para, SizeSrc);
        PrintPara(Para, SizeSrc);
        if (ParaDealA(Para, SizeSrc))
                return -1;
        PrintPara(Para, SizeSrc);
        if (ParaDealB(Para, SizeSrc))
                return -1;
        return 0;
}
/***********************************************************************************
***********************************************************************************/
static int GetParaBuffer(double* Para, const double* X, const double* Y, int Amount, int SizeSrc)
{
        int i, j;
        for (i = 0; i < SizeSrc; i++)
                for (ParaBuffer(Para, 0, i) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, 0, i) += pow(*(X + j), 2 * (SizeSrc - 1) - i);
        for (i = 1; i < SizeSrc; i++)
                for (ParaBuffer(Para, i, SizeSrc - 1) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, i, SizeSrc - 1) += pow(*(X + j), SizeSrc - 1 - i);
        for (i = 0; i < SizeSrc; i++)
                for (ParaBuffer(Para, i, SizeSrc) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, i, SizeSrc) += (*(Y + j)) * pow(*(X + j), SizeSrc - 1 - i);
        for (i = 1; i < SizeSrc; i++)
                for (j = 0; j < SizeSrc - 1; j++)
                        ParaBuffer(Para, i, j) = ParaBuffer(Para, i - 1, j + 1);
        return 0;
}
/***********************************************************************************
***********************************************************************************/
int Cal(const double* BufferX, const double* BufferY, int Amount, int SizeSrc, double* ParaResK)
{
        double* ParaK = (double*)malloc(SizeSrc * (SizeSrc + 1) * sizeof(double));
        GetParaBuffer(ParaK, BufferX, BufferY, Amount, SizeSrc);
        ParaDeal(ParaK, SizeSrc);
        for (Amount = 0; Amount < SizeSrc; Amount++, ParaResK++)
                *ParaResK = ParaBuffer(ParaK, Amount, SizeSrc);
        free(ParaK);
        return 0;
}
/***********************************************************************************
***********************************************************************************/
int main(int argc, char* argv[])
{
        int Amount;
        double BufferX[1024], BufferY[1024], ParaK[6];
        if (GetXY((const char*)"test.txt", (double*)BufferX, (double*)BufferY, &Amount))
                return 0;
        Cal((const double*)BufferX, (const double*)BufferY, Amount, sizeof(ParaK) / sizeof(double), (double*)ParaK);
        for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)
                printf("ParaK[%d] = %lf!\r\n", Amount, ParaK[Amount]);
        return 0;
}

测试test.txt里的内容
      0.995119      -7.620000
      2.001185      -2.460000
      2.999068     108.760000
      4.001035     625.020000
      4.999859    2170.500000
      6.004461    5814.580000
      6.999335   13191.840000
      7.999433   26622.060000
      9.002257   49230.220000
     10.003888   85066.500000
     11.004076  139226.280000
     12.001602  217970.140000
     13.003390  328843.860000
     14.001623  480798.420000
     15.003034  684310.000000
     16.002561  951499.980000
     17.003010 1296254.940000
     18.003897 1734346.660000
     19.002563 2283552.120000
     20.003530 2963773.500000
matlab拟合m代码：
load test.txt;
x = test(:, 1)';
y = test(:, 2)';
para = polyfit(x, y, 5);

阿莫论坛20周年了！感谢大家的支持与爱护！！

知道什么是神吗？其实神本来也是人，只不过神做了人做不到的事情 所以才成了神。 （头文字D, 杜汶泽）
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 16:54:51 | 显示全部楼层
VC工程文件

本帖子中包含更多资源

您需要 登录 才可以下载或查看，没有帐号？注册

x
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 16:58:12 | 显示全部楼层
二阶的很好用（二阶的话还可以优化很多，为了适应普遍性牺牲了一些优化），高阶的还需要注意很多地方（主要是限制幅度），程序尚在改进中...
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 17:01:10 | 显示全部楼层
从数学算法上理论上可以实现很高的阶数拟合，但是受到数据溢出的限制，高阶比较麻烦，可能该成大数运算会好一些
回复

   举报本楼层   

mhw

出0入0汤圆
发表于 2013-5-20 17:06:43 | 显示全部楼层
不错……想起以前做数值处理时，参考的那本经典的C算法，通篇变量都是a,b,c,d……aaa,bbb……看得令人崩溃
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 17:08:54 | 显示全部楼层
mhw 发表于 2013-5-20 17:06
不错……想起以前做数值处理时，参考的那本经典的C算法，通篇变量都是a,b,c,d……aaa,bbb……看得令人崩溃 ...

刚写完，写了一整天了，有点累了
回复

   举报本楼层   

nnnkey

出0入0汤圆
发表于 2013-5-20 17:20:34 来自手机 | 显示全部楼层
楼主，这个主要用途是什么=-O不懂
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 17:24:34 | 显示全部楼层
nnnkey 发表于 2013-5-20 17:20
楼主，这个主要用途是什么=-O不懂

目前没有用武之地，看到有个老帖子（http://www.amobbs.com/thread-5044969-1-1.html），突发兴致写的，刚开始写了二阶的，然后升级到高阶
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 17:30:15 | 显示全部楼层
nnnkey 发表于 2013-5-20 17:20
楼主，这个主要用途是什么=-O不懂

之前用PCI1714采集数据时用到过拟合，但是阶数很高，用Matlab混合编程做的，自己写的这个也就一个原理性的东西，拟合低阶少数据的场合还行，高阶暂时不考虑它
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 17:31:54 | 显示全部楼层
还有一个老帖子http://www.amobbs.com/thread-753268-1-1.html
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 22:35:27 | 显示全部楼层
写得很辛苦，自己顶一下，别沉了，为它添加简单的注释，详细的注释不好添加，准备把数学推导整理出来再添加详细注释，
简单思路如下：
1，采用目标函数对多项式系数求偏导，得到最优值条件，组成一个方程组；
2，方程组的解法采用行列式变换（两次变换：普通行列式——三角行列式——对角行列式——求解），行列式的求解算法上优化过一次了，目前还没有更好的思路再优化运算方法，限幅和精度准备再修改修改
目前存在的问题：
1，代码还是太粗糙
2，数学原理可行，但是计算机运算有幅度溢出和精度问题，这方面欠考虑，导致高阶大数据可能拟合失败
下面附简单注释版代码：
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"

//把二维的数组与一维数组的转换，也可以直接用二维数组，只是我的习惯是不用二维数组
#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))
//

/***********************************************************************************
从txt文件里读取double型的X，Y数据
txt文件里的存储格式为
X1  Y1
X2  Y2
X3  Y3
X4  Y4
X5  Y5
X6  Y6
X7  Y7
X8  Y8
函数返回X，Y，以及数据的数目（以组为单位）
***********************************************************************************/
static int GetXY(const char* FileName, double* X, double* Y, int* Amount)
{
        FILE* File = fopen(FileName, "r");
        if (!File)
                return -1;
        for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)
                if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))
                        break;
        fclose(File);
        return 0;
}

/***********************************************************************************
打印系数矩阵，只用于调试，不具备运算功能
对于一个N阶拟合，它的系数矩阵大小是（N + 1）行（N + 2）列
double* Para：系数矩阵存储地址
int SizeSrc：系数矩阵大小（SizeSrc）行（SizeSrc + 1）列
***********************************************************************************/
static int PrintPara(double* Para, int SizeSrc)
{
        int i, j;
        for (i = 0; i < SizeSrc; i++)
        {
                for (j = 0; j <= SizeSrc; j++)
                        printf("%10.6lf ", ParaBuffer(Para, i, j));
                printf("\r\n");
        }
        printf("\r\n");
        return 0;
}

/***********************************************************************************
系数矩阵的限幅处理，防止它溢出，目前这个函数很不完善，并不能很好地解决这个问题
原理：矩阵解行列式，同一行乘以一个系数，行列式的解不变
当然，相对溢出问题，还有一个精度问题，也是同样的思路，现在对于这两块的处理很不完善，有待优化
以行为单位处理
***********************************************************************************/
static int ParalimitRow(double* Para, int SizeSrc, int Row)
{
        int i;
        double Max, Min, Temp;
        for (Max = abs(ParaBuffer(Para, Row, 0)), Min = Max, i = SizeSrc; i; i--)
        {
                Temp = abs(ParaBuffer(Para, Row, i));
                if (Max < Temp)
                        Max = Temp;
                if (Min > Temp)
                        Min = Temp;
        }
        Max = (Max + Min) * 0.000005;
        for (i = SizeSrc; i >= 0; i--)
                ParaBuffer(Para, Row, i) /= Max;
        return 0;
}

/***********************************************************************************
同上，以矩阵为单位处理
***********************************************************************************/
static int Paralimit(double* Para, int SizeSrc)
{
        int i;
        for (i = 0; i < SizeSrc; i++)
                if (ParalimitRow(Para, SizeSrc, i))
                        return -1;
        return 0;
}

/***********************************************************************************
系数矩阵行列式变换
***********************************************************************************/
static int ParaPreDealA(double* Para, int SizeSrc, int Size)
{
        int i, j;
        for (Size -= 1, i = 0; i < Size; i++)
        {
                for (j = 0; j < Size; j++)
                        ParaBuffer(Para, i, j) = ParaBuffer(Para, i, j) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, j) * ParaBuffer(Para, i, Size);
                ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, SizeSrc) * ParaBuffer(Para, i, Size);
                ParaBuffer(Para, i, Size) = 0;
                ParalimitRow(Para, SizeSrc, i);
        }
        return 0;
}

/***********************************************************************************
系数矩阵行列式变换，与ParaPreDealA配合
完成第一次变换，变换成三角矩阵
***********************************************************************************/
static int ParaDealA(double* Para, int SizeSrc)
{
        int i;
        for (i = SizeSrc; i; i--)
                if (ParaPreDealA(Para, SizeSrc, i))
                        return -1;
        return 0;
}

/***********************************************************************************
系数矩阵行列式变换
***********************************************************************************/
static int ParaPreDealB(double* Para, int SizeSrc, int OffSet)
{
        int i, j;
        for (i = OffSet + 1; i < SizeSrc; i++)
        {
                for (j = OffSet + 1; j <= i; j++)
                        ParaBuffer(Para, i, j) *= ParaBuffer(Para, OffSet, OffSet);
                ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, OffSet, OffSet) - ParaBuffer(Para, i, OffSet) * ParaBuffer(Para, OffSet, SizeSrc);
                ParaBuffer(Para, i, OffSet) = 0;
                ParalimitRow(Para, SizeSrc, i);
        }
        return 0;
}

/***********************************************************************************
系数矩阵行列式变换，与ParaPreDealB配合
完成第一次变换，变换成对角矩阵，变换完毕
***********************************************************************************/
static int ParaDealB(double* Para, int SizeSrc)
{
        int i;
        for (i = 0; i < SizeSrc; i++)
                if (ParaPreDealB(Para, SizeSrc, i))
                        return -1;
        for (i = 0; i < SizeSrc; i++)
        {
                if (ParaBuffer(Para, i, i))
                {
                        ParaBuffer(Para, i, SizeSrc) /= ParaBuffer(Para, i, i);
                        ParaBuffer(Para, i, i) = 1.0;
                }
        }
        return 0;
}

/***********************************************************************************
系数矩阵变换
***********************************************************************************/
static int ParaDeal(double* Para, int SizeSrc)
{
        PrintPara(Para, SizeSrc);
        Paralimit(Para, SizeSrc);
        PrintPara(Para, SizeSrc);
        if (ParaDealA(Para, SizeSrc))
                return -1;
        PrintPara(Para, SizeSrc);
        if (ParaDealB(Para, SizeSrc))
                return -1;
        return 0;
}

/***********************************************************************************
最小二乘法的第一步就是从XY数据里面获取系数矩阵
double* Para：系数矩阵地址
const double* X：X数据地址
const double* Y：Y数据地址
int Amount：XY数据组数
int SizeSrc：系数矩阵大小（SizeSrc）行（SizeSrc + 1）列
***********************************************************************************/
static int GetParaBuffer(double* Para, const double* X, const double* Y, int Amount, int SizeSrc)
{
        int i, j;
        for (i = 0; i < SizeSrc; i++)
                for (ParaBuffer(Para, 0, i) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, 0, i) += pow(*(X + j), 2 * (SizeSrc - 1) - i);
        for (i = 1; i < SizeSrc; i++)
                for (ParaBuffer(Para, i, SizeSrc - 1) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, i, SizeSrc - 1) += pow(*(X + j), SizeSrc - 1 - i);
        for (i = 0; i < SizeSrc; i++)
                for (ParaBuffer(Para, i, SizeSrc) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, i, SizeSrc) += (*(Y + j)) * pow(*(X + j), SizeSrc - 1 - i);
        for (i = 1; i < SizeSrc; i++)
                for (j = 0; j < SizeSrc - 1; j++)
                        ParaBuffer(Para, i, j) = ParaBuffer(Para, i - 1, j + 1);
        return 0;
}

/***********************************************************************************
整个计算过程
***********************************************************************************/
int Cal(const double* BufferX, const double* BufferY, int Amount, int SizeSrc, double* ParaResK)
{
        double* ParaK = (double*)malloc(SizeSrc * (SizeSrc + 1) * sizeof(double));
        GetParaBuffer(ParaK, BufferX, BufferY, Amount, SizeSrc);
        ParaDeal(ParaK, SizeSrc);
        for (Amount = 0; Amount < SizeSrc; Amount++, ParaResK++)
                *ParaResK = ParaBuffer(ParaK, Amount, SizeSrc);
        free(ParaK);
        return 0;
}

/***********************************************************************************
测试main函数
数据组数：20
阶数：5
***********************************************************************************/
int main(int argc, char* argv[])
{
        //数据组数
        int Amount;
        //XY缓存，系数矩阵缓存
        double BufferX[1024], BufferY[1024], ParaK[6];
        if (GetXY((const char*)"test.txt", (double*)BufferX, (double*)BufferY, &Amount))
                return 0;
        //运算
        Cal((const double*)BufferX, (const double*)BufferY, Amount, sizeof(ParaK) / sizeof(double), (double*)ParaK);
        //输出系数
        for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)
                printf("ParaK[%d] = %lf!\r\n", Amount, ParaK[Amount]);
        //屏幕暂留
        system("pause");
        return 0;
}
回复

   举报本楼层   

wcm_e

出0入0汤圆
发表于 2013-5-20 22:46:05 | 显示全部楼层
学习，谢谢
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-20 23:00:35 | 显示全部楼层
wcm_e 发表于 2013-5-20 22:46
学习，谢谢

呵呵，共同学习
回复

   举报本楼层   

bigallium

出0入0汤圆
发表于 2013-5-21 00:07:39 | 显示全部楼层
这个要顶
话说我到现在还搞不懂最小二乘法的原理...
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:11:25 | 显示全部楼层
bigallium 发表于 2013-5-21 00:07
这个要顶
话说我到现在还搞不懂最小二乘法的原理...

今天有点晚了，数学推导就不详细发了，不过在之前一个老帖子里提到过二阶的（高阶的处理方法类似）http://www.amobbs.com/thread-5044969-1-1.html，介绍了一丁点的数学原理，很不详细，更详细的正在整理，office用得不熟，很吃力，希望明天能贴出来
回复

   举报本楼层   

rifjft

出0入0汤圆
发表于 2013-5-21 00:13:22 | 显示全部楼层
三寸金…… 很特别
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:42:16 | 显示全部楼层
C:\Documents and Settings\Administrator\桌面\新建文件夹\1.jpg
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:46:17 | 显示全部楼层
不擅长Office，手写了一个简单数学推导（以二阶为例），手机拍的，不会发图，图片在pic.rar里，共三张

本帖子中包含更多资源

您需要 登录 才可以下载或查看，没有帐号？注册

x
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:48:01 | 显示全部楼层
发错了，漏了一张，重发

本帖子中包含更多资源

您需要 登录 才可以下载或查看，没有帐号？注册

x
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:51:28 | 显示全部楼层
系数矩阵很容易看规律来，高阶的也符合这个规律，然后写一个通用的行列式的解法，问题就OK了一半
行列变换的顺序按照图中给出的顺序可以很方便地设计出运算量比较小的高阶行列式变换矩阵算法（类似迭代），上文中的几个Deal函数就是这个思路
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:52:06 | 显示全部楼层
数学上理解不难，时基实施会有很大问题，主要是数据精度，数据的溢出问题
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:53:48 | 显示全部楼层
程序中设计有类似增益调节的函数（Paralimit)，可以部分的解决一些问题，但是目前没有很完美的思路，所以高阶，多数据，大数据慎用此程序
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:55:23 | 显示全部楼层
程序的调试思路也比较简单，一步步跟踪系数矩阵的变换，每次变换后输出系数矩阵，看看问题在哪，就知道有哪些地方不完善了
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:56:36 | 显示全部楼层
矩阵变换思路：
普通矩阵——下三角矩阵——对角矩阵，之所以这么设计，考虑到高阶采用嵌套和迭代
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 00:56:54 | 显示全部楼层
宿舍没电了，呜呜
回复

   举报本楼层   

DashMBO

出0入0汤圆
发表于 2013-5-21 08:20:23 | 显示全部楼层
waiting for 下文。
回复

   举报本楼层   

xjf616

出0入0汤圆
发表于 2013-5-21 08:23:05 | 显示全部楼层
mark，支持lz
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 08:32:41 | 显示全部楼层
xjf616 发表于 2013-5-21 08:23
mark，支持lz

谢谢支持
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 08:35:13 | 显示全部楼层
转载一个别人写的，转载地址：http://wenwen.soso.com/z/q190332592.htm，刚在搜搜看到的，还没细读，验证，不过貌似我们用的原理是一样的，先转了再度吧，博采众长

#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <process.h>
#define N 5//N个点
#define T 3 //T次拟合
#define W 1//权函数
#define PRECISION 0.00001
float pow_n(float a,int n)
{
int i;
if(n==0)
return(1);
float res=a;
for(i=1;i<n;i++)
{
res*=a;
}
return(res);
}
void mutiple(float a[][N],float b[][T+1],float c[][T+1])
{
float res=0;
int i,j,k;
for(i=0;i<T+1;i++)
for(j=0;j<T+1;j++)
{
res=0;
for(k=0;k<N;k++)
{
res+=a[k]*b[k][j];
c[j]=res;
}
}
}
void matrix_trans(float a[][T+1],float b[][N])
{
int i,j;
for(i=0;i<N;i++)
{
for(j=0;j<T+1;j++)
{
b[j]=a[j];
}
}
}
void init(float x_y[][2],int n)
{
int i;
printf("请输入%d个已知点：\n",N);
for(i=0;i<n;i++)
{
printf("(x%d y%d):",i,i);
scanf("%f %f",&x_y[0],&x_y[1]);
}
}
void get_A(float matrix_A[][T+1],float x_y[][2],int n)
{
int i,j;
for(i=0;i<N;i++)
{
for(j=0;j<T+1;j++)
{
matrix_A[j]=W*pow_n(x_y[0],j);
}
}
}
void print_array(float array[][T+1],int n)
{
int i,j;
for(i=0;i<n;i++)
{
for(j=0;j<T+1;j++)
{
printf("%-g",array[j]);
}
printf("\n");
}
}
void convert(float argu[][T+2],int n)
{
int i,j,k,p,t;
float rate,temp;
for(i=1;i<n;i++)
{
for(j=i;j<n;j++)
{
if(argu[i-1][i-1]==0)
{
for(p=i;p<n;p++)
{
if(argu[p][i-1]!=0)
break;
}
if(p==n)
{
printf("方程组无解!\n");
exit(0);
}
for(t=0;t<n+1;t++)
{
temp=argu[i-1][t];
argu[i-1][t]=argu[p][t];
argu[p][t]=temp;
}
}
rate=argu[j][i-1]/argu[i-1][i-1];
for(k=i-1;k<n+1;k++)
{
argu[j][k]-=argu[i-1][k]*rate;
if(fabs(argu[j][k])<=PRECISION)
argu[j][k]=0;
}
}
}
}
void compute(float argu[][T+2],int n,float root[])
{
int i,j;
float temp;
for(i=n-1;i>=0;i--)
{
temp=argu[n];
for(j=n-1;j>i;j--)
{
temp-=argu[j]*root[j];
}
root=temp/argu;
}
}
void get_y(float trans_A[][N],float x_y[][2],float y[],int n)
{
int i,j;
float temp;
for(i=0;i<n;i++)
{
temp=0;
for(j=0;j<N;j++)
{
temp+=trans_A[j]*x_y[j][1];
}
y=temp;
}
}
void cons_formula(float coef_A[][T+1],float y[],float coef_form[][T+2])
{
int i,j;
for(i=0;i<T+1;i++)
{
for(j=0;j<T+2;j++)
{
if(j==T+1)
coef_form[j]=y;
else
coef_form[j]=coef_A[j];
}
}
}
void print_root(float a[],int n)
{
int i,j;
printf("%d个点的%d次拟合的多项式系数为:\n",N,T);
for(i=0;i<n;i++)
{
printf("a[%d]=%g,",i+1,a);
}
printf("\n");
printf("拟合曲线方程为:\ny(x)=%g",a[0]);
for(i=1;i<n;i++)
{
printf(" + %g",a);
for(j=0;j<i;j++)
{
printf("*X");
}
}
printf("\n");
}
void process()
{
float x_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef_formu[T+1][T+2],y[T+1],a[T+1];
init(x_y,N);
get_A(matrix_A,x_y,N);
printf("矩阵A为：\n");
print_array(matrix_A,N);
matrix_trans(matrix_A,trans_A);
mutiple(trans_A,matrix_A,coef_A);
printf("法矩阵为：\n");
print_array(coef_A,T+1);
get_y(trans_A,x_y,y,T+1);
cons_formula(coef_A,y,coef_formu);
convert(coef_formu,T+1);
compute(coef_formu,T+1,a);
print_root(a,T+1);
}
void main()
{
process();
}
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 08:39:59 | 显示全部楼层
DashMBO 发表于 2013-5-21 08:20
waiting for 下文。

二阶的推导思路在那个pic.rar里面，高阶的你找那个二阶拟合的系数矩阵，会发现规律，高阶拟合的数学证明最好用矩阵运算去理解，二阶的参数不多（3个）直接强行运算就好了
矩阵的行变换算法顺序，按照图中的顺序就好，运算次数不是很多
数据溢出，数据精度这些目前还没思路
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 08:59:10 | 显示全部楼层
二阶的系数矩阵方程组：
(目标函数y = A * x^2 + b * x^1 + C *x^0)
注意！
(X^N表示的是X1^N + X2^N + X3^N + X4^N + …… + Xn^N)
(Y * X^N表示的是Y1 * X1^N + Y1 * X2^N + Y3 * X3^N + Y4 * X4^N + …… + Yn * Xn^N)
X^4 * A + X^3 * B + X^2 * C = Y * X ^2
X^3 * A + X^2 * B + X^1 * C = Y * X ^1
X^2 * A + X^1 * B + X^0 * C = Y * X ^0

三阶的系数矩阵方程组：
注意！
(X^N表示的是X1^N + X2^N + X3^N + X4^N + …… + Xn^N)
(Y * X^N表示的是Y1 * X1^N + Y1 * X2^N + Y3 * X3^N + Y4 * X4^N + …… + Yn * Xn^N)
(目标函数y = A * x^3 + b * x^2 + C *x^1 + D * x^0)
X^6 * A + X^5 * B + X^4 * C + X^3 * D = Y * X ^3
X^5 * A + X^4 * B + X^3 * C + X^2 * D = Y * X ^2
X^4 * A + X^3 * B + X^2 * C + X^1 * D = Y * X ^1
X^3 * A + X^2 * B + X^1 * C + X^0 * D = Y * X ^0

更高阶的规律也跟上面一样，所以第一步就是获取系数矩阵，需要算的值有X^0， X^1， …… X^（2N）， Y * X^0, Y * X^1,……Y* X^N，N表示阶数
在本程序中用如下代码实现：
/***********************************************************************************
最小二乘法的第一步就是从XY数据里面获取系数矩阵
double* Para：系数矩阵地址
const double* X：X数据地址
const double* Y：Y数据地址
int Amount：XY数据组数
int SizeSrc：系数矩阵大小（SizeSrc）行（SizeSrc + 1）列
***********************************************************************************/
static int GetParaBuffer(double* Para, const double* X, const double* Y, int Amount, int SizeSrc)
{
        int i, j;
//先把第0行的搞定
        for (i = 0; i < SizeSrc; i++)
                for (ParaBuffer(Para, 0, i) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, 0, i) += pow(*(X + j), 2 * (SizeSrc - 1) - i);
//把第（SizeSrc - 1）列的搞定
        for (i = 1; i < SizeSrc; i++)
                for (ParaBuffer(Para, i, SizeSrc - 1) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, i, SizeSrc - 1) += pow(*(X + j), SizeSrc - 1 - i);
//把第（SizeSrc）列的搞定
        for (i = 0; i < SizeSrc; i++)
                for (ParaBuffer(Para, i, SizeSrc) = 0, j = 0; j < Amount; j++)
                        ParaBuffer(Para, i, SizeSrc) += (*(Y + j)) * pow(*(X + j), SizeSrc - 1 - i);
//剩下的全是拷贝，不存在复杂计算了，省了很多事
        for (i = 1; i < SizeSrc; i++)
                for (j = 0; j < SizeSrc - 1; j++)
                        ParaBuffer(Para, i, j) = ParaBuffer(Para, i - 1, j + 1);
        return 0;
}
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:01:07 | 显示全部楼层
推导残差方程——对残差方程偏导——求出系数矩阵——解系数矩阵——OK了
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:03:32 | 显示全部楼层
残差方程的推导以二阶为例（有耐心的话三阶也能，更高阶的推导换思路），再找出矩阵规律（有点投机取巧了，要是有耐心，还可以证明一下，不过貌似不简单）
回复

   举报本楼层   

loongsuns

出0入0汤圆
发表于 2013-5-21 09:05:26 | 显示全部楼层
这个偶尔会用到，mark一下，经常都是matlab算好，如果可以直接自行解算那就有意思了，多谢分享
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:18:29 | 显示全部楼层
解系数矩阵的算法：行变换
1，矩阵变换有个特点，尽可能快地变换出更多的0，而且0的位置最好有规律，因为可以在写循环的时候直接PASS掉，省掉很多浮点运算，快速性变好
2，高阶运算的计算步骤很多，最好是实现迭代或者递推，这样写代码轻松很多
3，矩阵系数中可能本来就含有0，因此慎用除法，当然用乘法也有一个严重的问题，很容易溢出（这一条现在很不严谨，可能有疏漏谬误）
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:25:12 | 显示全部楼层
以三阶为例，行变换顺序：
K00   K01   K02   K03   K04
K10   K11   K12   K13   K14
K20   K21   K22   K23   K24
K30   K31   K32   K33   K34
利用第3行消去第0,1,2行的第3列
K00   K01   K02   000   K04
K10   K11   K12   000   K14
K20   K21   K22   000   K24
K30   K31   K32   K33   K34
利用第2行消去第0,1行的第2列
K00   K01   000   000   K04
K10   K11   000   000   K14
K20   K21   K22   000   K24
K30   K31   K32   K33   K34
利用第1行消去第0行的第1列
K00   000   000   000   K04
K10   K11   000   000   K14
K20   K21   K22   000   K24
K30   K31   K32   K33   K34
上面的三个步骤用的函数，运用了递推的算法
/***********************************************************************************
系数矩阵行列式变换
***********************************************************************************/
static int ParaPreDealA(double* Para, int SizeSrc, int Size)
{
        int i, j;
        for (Size -= 1, i = 0; i < Size; i++)
        {
                for (j = 0; j < Size; j++)
                        ParaBuffer(Para, i, j) = ParaBuffer(Para, i, j) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, j) * ParaBuffer(Para, i, Size);
                ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, SizeSrc) * ParaBuffer(Para, i, Size);
                ParaBuffer(Para, i, Size) = 0;
                ParalimitRow(Para, SizeSrc, i);
        }
        return 0;
}

/***********************************************************************************
系数矩阵行列式变换，与ParaPreDealA配合
完成第一次变换，变换成三角矩阵
***********************************************************************************/
static int ParaDealA(double* Para, int SizeSrc)
{
        int i;
        for (i = SizeSrc; i; i--)
                if (ParaPreDealA(Para, SizeSrc, i))
                        return -1;
        return 0;
}
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:27:59 | 显示全部楼层
K00   000   000   000   K04
K10   K11   000   000   K14
K20   K21   K22   000   K24
K30   K31   K32   K33   K34
用第0行消去第1,2,3行
K00   000   000   000   K04
000   K11   000   000   K14
000   K21   K22   000   K24
000   K31   K32   K33   K34
用第1行消去第2,3行
K00   000   000   000   K04
000   K11   000   000   K14
000   000   K22   000   K24
000   000   K32   K33   K34
用第2行消去第3行
K00   000   000   000   K04
000   K11   000   000   K14
000   000   K22   000   K24
000   000   000   K33   K34
这个操作用的函数如下，也运用了递推
/***********************************************************************************
系数矩阵行列式变换
***********************************************************************************/
static int ParaPreDealB(double* Para, int SizeSrc, int OffSet)
{
        int i, j;
        for (i = OffSet + 1; i < SizeSrc; i++)
        {
                for (j = OffSet + 1; j <= i; j++)
                        ParaBuffer(Para, i, j) *= ParaBuffer(Para, OffSet, OffSet);
                ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, OffSet, OffSet) - ParaBuffer(Para, i, OffSet) * ParaBuffer(Para, OffSet, SizeSrc);
                ParaBuffer(Para, i, OffSet) = 0;
                ParalimitRow(Para, SizeSrc, i);
        }
        return 0;
}

/***********************************************************************************
系数矩阵行列式变换，与ParaPreDealB配合
完成第一次变换，变换成对角矩阵，变换完毕
***********************************************************************************/
static int ParaDealB(double* Para, int SizeSrc)
{
        int i;
        for (i = 0; i < SizeSrc; i++)
                if (ParaPreDealB(Para, SizeSrc, i))
                        return -1;
        for (i = 0; i < SizeSrc; i++)
        {
                if (ParaBuffer(Para, i, i))
                {
                        ParaBuffer(Para, i, SizeSrc) /= ParaBuffer(Para, i, i);
                        ParaBuffer(Para, i, i) = 1.0;
                }
        }
        return 0;
}
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:28:26 | 显示全部楼层
至此，这个三阶拟合的推导已经结束
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:34:52 | 显示全部楼层
loongsuns 发表于 2013-5-21 09:05
这个偶尔会用到，mark一下，经常都是matlab算好，如果可以直接自行解算那就有意思了，多谢分享 ...

工程上我也是用Matlab跟VC混合编程实现，自己写主要是怕溢出，精度这些，可靠性有点担忧，Matlab考虑的肯定比这个要全面的多，至少是好几个数量级，用着放心
写这个只是突发奇想，在思路上还是有值得借鉴的地方
回复

   举报本楼层   

winterw

出0入0汤圆
发表于 2013-5-21 09:40:58 | 显示全部楼层
感谢楼主。
mark一下回头看看
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 09:44:47 | 显示全部楼层
winterw 发表于 2013-5-21 09:40
感谢楼主。
mark一下回头看看

共同进步
回复

   举报本楼层   

yywin

出0入0汤圆
发表于 2013-5-21 13:25:57 | 显示全部楼层
LZ能不能推荐一本关于最小二乘原理的书
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 13:31:04 | 显示全部楼层
yywin 发表于 2013-5-21 13:25
LZ能不能推荐一本关于最小二乘原理的书

这个好像没有很特别的，大学教程《线性代数》和《高等数学》估计还不错，在百度文库上看到有几个PPT讲了一下最小二乘法的一些数学知识（全部用的矩阵运算，线性代数和高等数学忘了的话理解不容易），我没有细看
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 13:32:31 | 显示全部楼层
yywin 发表于 2013-5-21 13:25
LZ能不能推荐一本关于最小二乘原理的书

这东西主要是一个数学原理知识，还有一个就是计算机的具体实现，这两个都比较重要
回复

   举报本楼层   

15085362

出0入0汤圆
发表于 2013-5-21 13:35:53 | 显示全部楼层
正需要这个，多谢分享
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 13:38:40 | 显示全部楼层
15085362 发表于 2013-5-21 13:35
正需要这个，多谢分享

不客气，目前还不够完善，注意上面提到的一些限制，有些地方慎用，不保证任何情况100%拟合正确，只是提供一个思路
回复

   举报本楼层   

无级电工

出0入0汤圆
发表于 2013-5-21 13:52:07 | 显示全部楼层
赞一个：有功底。
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2013-5-21 13:53:48 | 显示全部楼层
无级电工 发表于 2013-5-21 13:52
赞一个：有功底。

惭愧难当，大一的时候《线性代数》零分来着
回复

   举报本楼层   

993430676

出0入0汤圆
发表于 2013-5-23 10:06:45 | 显示全部楼层
最小二乘法是数值分析里面的离散情形下的最佳平方逼近问题。
回复

   举报本楼层   

电子小生

出0入0汤圆
发表于 2014-4-11 19:09:17 | 显示全部楼层
mark                                          
回复

   举报本楼层   

zdhhpc

出0入0汤圆
发表于 2014-7-16 21:38:10 | 显示全部楼层
群主真是有心啊！！
回复

   举报本楼层   

emobbs

出0入0汤圆
发表于 2014-9-4 17:52:34 | 显示全部楼层
曲线拟合
回复

   举报本楼层   

bli19

出0入0汤圆
发表于 2014-9-12 03:19:57 来自手机 | 显示全部楼层
擦，我当时弄直线提取都差点吐了，单精度溢出，拟合出来了的倒是能看，但用LCD实时显示是相当地慢动作
回复

   举报本楼层   

nrwhc

出0入0汤圆
发表于 2014-9-23 00:40:57 | 显示全部楼层
最小二乘QR分解
回复

   举报本楼层   

lgy026

出0入0汤圆
发表于 2014-9-24 16:21:01 | 显示全部楼层
好厉害，顶一个
回复

   举报本楼层   

manyman

出0入0汤圆
发表于 2015-3-2 23:08:05 | 显示全部楼层
楼主，小白问一个问题哈，这种曲线拟合是不是可以用来做预测
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2015-3-2 23:14:16 | 显示全部楼层
manyman 发表于 2015-3-2 23:08
楼主，小白问一个问题哈，这种曲线拟合是不是可以用来做预测

可以的，不过有很多限制
回复

   举报本楼层   

wszyjsw2

出0入0汤圆
发表于 2015-3-2 23:22:56 | 显示全部楼层
mark  最小二乘法曲线拟合
回复

   举报本楼层   

migrant

出0入0汤圆
发表于 2015-3-3 09:34:46 | 显示全部楼层
多谢楼主贡献
回复

   举报本楼层   

xf331785508

出0入0汤圆
发表于 2015-3-3 10:25:44 | 显示全部楼层
zouzhichao 发表于 2013-5-20 22:35
写得很辛苦，自己顶一下，别沉了，为它添加简单的注释，详细的注释不好添加，准备把数学推导整理出来再添加 ...
  1. //把二维的数组与一维数组的转换，也可以直接用二维数组，只是我的习惯是不用二维数组
  2. #define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))
复制代码


第一句的宏，跳出个没有的参数SizeSrc ，不如把参数加进去。
我没理解错的话，SizeSrc应该是二维数组的第一维总长度吧。

  1. #define  ARRAY_CONVERT(Buffer, row, col, n)  (*((Buffer) + (row)*((n)+1) + (col)))
复制代码
回复

   举报本楼层   

manyman

出0入0汤圆
发表于 2015-3-3 22:13:03 | 显示全部楼层
zouzhichao 发表于 2015-3-2 23:14
可以的，不过有很多限制

有没有 关于预测方面的资料啊
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2015-3-3 22:17:56 来自手机 | 显示全部楼层
xf331785508 发表于 2015-3-3 10:25
第一句的宏，跳出个没有的参数SizeSrc ，不如把参数加进去。
我没理解错的话，SizeSrc应该是二维数组的 ...

是的，你说的没错
回复

   举报本楼层   

kisscontrol

出0入0汤圆
发表于 2015-3-4 14:06:34 | 显示全部楼层
楼主很用心，很细致。   
回复

   举报本楼层   

zzzhpzhp

出0入0汤圆
发表于 2015-5-2 08:20:31 | 显示全部楼层
谢谢分享
回复

   举报本楼层   

baoqi

出0入0汤圆
发表于 2015-5-2 08:37:13 | 显示全部楼层
谢谢分享
回复

   举报本楼层   

悟空间66

出0入0汤圆
发表于 2015-5-2 09:40:38 | 显示全部楼层
收藏，技术型好贴
回复

   举报本楼层   

symcu

出0入0汤圆
发表于 2015-6-12 11:29:19 | 显示全部楼层
正是我要用到的 太感谢了。
回复

   举报本楼层   

ycwjl728

出0入24汤圆
发表于 2015-6-12 14:40:37 | 显示全部楼层
Mark!学习一下
回复

   举报本楼层   

jsszdfdn

出0入0汤圆
发表于 2015-8-25 13:43:05 | 显示全部楼层
zouzhichao 发表于 2015-3-3 22:17
是的，你说的没错

你好 楼主 请教一个问题 如果用最小二乘法 根据图形上的一些点 求得一个函数 y=p(x)， 然后我将这些点以其中某点为原心 逆时针 或 顺时针 方向转10度 此时我再用 最小二乘法去求这些点的 拟和函数 会和原先的函数一样么？
                              这个问题主要是想知道 角度对拟合函数有何影响
                                                                                                                        谢谢。
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2015-8-25 20:46:31 | 显示全部楼层
jsszdfdn 发表于 2015-8-25 13:43
你好 楼主 请教一个问题 如果用最小二乘法 根据图形上的一些点 求得一个函数 y=p(x)， 然后我将这些点以 ...

肯定有影响啊！！！
回复

   举报本楼层   

jsszdfdn

出0入0汤圆
发表于 2015-8-26 10:41:54 | 显示全部楼层
zouzhichao 发表于 2015-8-25 20:46
肯定有影响啊！！！

你好 楼主  ， 这样的话 原图中的一系列点的之间的位置关系虽没变化（比如 在16*16的框中 有个“2” 这个2上的每个点之间都有位置关系）但如果整体位置相对于框有变化后 经过最小二乘法所得的曲线公式 就肯定不同可以这么理解么？ 谢谢。
回复

   举报本楼层   

jsszdfdn

出0入0汤圆
发表于 2015-8-26 10:44:24 | 显示全部楼层
最小二乘法如果对一个圆 或则 "8" 这样的封闭图形能求出拟合公式么？
回复

   举报本楼层   

jsszdfdn

出0入0汤圆
发表于 2015-8-26 10:47:06 | 显示全部楼层
当位置变化后（拟合前的内容不变，比如 16*16中的“2”，平移或旋转一点角度）这时 拟合之前和拟合之后的公式差别会在哪里？
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2015-8-26 19:54:07 | 显示全部楼层
jsszdfdn 发表于 2015-8-26 10:41
你好 楼主  ， 这样的话 原图中的一系列点的之间的位置关系虽没变化（比如 在16*16的框中 有个“2” 这个 ...

这里的最小二乘法是多项式拟合，目标函数是多项式，不是任意图形
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2015-8-26 19:55:12 | 显示全部楼层
jsszdfdn 发表于 2015-8-26 10:44
最小二乘法如果对一个圆 或则 "8" 这样的封闭图形能求出拟合公式么？

这个不太好实现
回复

   举报本楼层   

jsszdfdn

出0入0汤圆
发表于 2015-8-27 09:13:23 | 显示全部楼层
zouzhichao 发表于 2015-8-26 19:55
这个不太好实现

谢谢 楼主的回复，对最小二乘法的使用范有点了解了。
回复

   举报本楼层   

tkloveyang

出0入0汤圆
发表于 2015-9-27 17:10:31 | 显示全部楼层
zouzhichao 发表于 2013-5-20 17:31
还有一个老帖子http://www.amobbs.com/thread-753268-1-1.html

看过你之前的帖子，不知道我理解的对不对，对N个数据M阶拟合计算量是N组数的合并化简运算+解(M+1)元方程组吗？
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2015-9-27 17:34:32 来自手机 | 显示全部楼层
tkloveyang 发表于 2015-9-27 17:10
看过你之前的帖子，不知道我理解的对不对，对N个数据M阶拟合计算量是N组数的合并化简运算+解(M+1)元方程 ...

是的，就是这样的
回复

   举报本楼层   

eliterxzgxu

出0入0汤圆
发表于 2015-9-27 23:15:56 | 显示全部楼层
Mark!学习一下
回复

   举报本楼层   

MRRIGHT

出0入0汤圆
发表于 2015-10-11 09:28:59 来自手机 | 显示全部楼层
要学习一下
回复

   举报本楼层   

bi大痣

出0入0汤圆
发表于 2016-10-14 11:40:52 | 显示全部楼层
mark一下，楼主加油！
回复

   举报本楼层   

captain_steve

出0入0汤圆
发表于 2016-10-14 13:02:42 | 显示全部楼层
mark一下，加油加油
回复

   举报本楼层   

philix

出0入0汤圆
发表于 2017-6-25 08:58:15 | 显示全部楼层
高手，谢谢分享
回复

   举报本楼层   

ddcour

出10入95汤圆
发表于 2017-11-12 12:42:51 | 显示全部楼层
像高手致敬！多谢分享！
回复

   举报本楼层   

wind2100

出0入0汤圆
发表于 2018-5-11 07:58:46 | 显示全部楼层
顶, 看起来简单, 写和调起来有花上几天时间,是不行的.
回复

   举报本楼层   

zhongjp85

出0入0汤圆
发表于 2018-5-11 14:59:07 | 显示全部楼层
火钳刘明 多谢LZ
回复

   举报本楼层   

taojie

出0入0汤圆
发表于 2018-6-19 09:18:21 | 显示全部楼层
这个的理论知识挺深的
回复

   举报本楼层   

zouzhichao

出10入23汤圆
 楼主| 发表于 2018-6-19 09:19:28 来自手机 | 显示全部楼层
taojie 发表于 2018-6-19 09:18
这个的理论知识挺深的

这个还深啊？
回复

   举报本楼层   

ywlzh

出0入0汤圆
发表于 2018-6-19 09:33:51 | 显示全部楼层
http://www.gnu.org/software/gsl/manual/html_node/index_old.html

一看时间 老帖了

回复

   举报本楼层   

bg6agf

出0入0汤圆
发表于 2018-6-19 10:24:20 | 显示全部楼层
真老帖也！
回复

   举报本楼层   

zf12862177

出0入0汤圆
发表于 2018-6-19 11:03:52 | 显示全部楼层
谢谢分享
回复

   举报本楼层   

wang0529

出0入0汤圆
发表于 2018-7-15 12:55:40 | 显示全部楼层
多谢楼主贡献
回复

   举报本楼层   

eliterxzgxu

出0入0汤圆
发表于 2018-7-15 13:18:13 来自手机 | 显示全部楼层
感谢楼主分享
回复

   举报本楼层   

jacktau

出0入0汤圆
发表于 2018-7-15 14:23:01 | 显示全部楼层
非常感谢，学习一下~
回复

   举报本楼层   

返回列表 发新帖
回帖提示： 反政府言论将被立即封锁ID 。在按“提交”前，请自问一下：我这样表达会给举报吗，会给自己惹麻烦吗？ 另外：尽量不要使用Mark、顶等没有意义的回复。不得大量使用大字体和彩色字。【本论坛不允许直接上传手机拍摄图片，浪费大家下载带宽和论坛服务器空间，请压缩后（图片小于1兆）才上传。压缩方法可以在微信里面发给自己（不要勾选“原图），然后下载，就能得到压缩后的图片。注意：要连续压缩2次才能满足要求！！】。另外，手机版只能上传图片，要上传附件需要切换到电脑版（不需要使用电脑，手机上切换到电脑版就行，页面底部）。
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 注册

本版积分规则

手机版|Archiver|amobbs.com 阿莫电子技术论坛 ( 粤ICP备2022115958号, 版权所有：东莞阿莫电子贸易商行 创办于2004年 (公安交互式论坛备案：44190002001997 ) )

GMT+8, 2024-7-23 11:27

© Since 2004 www.amobbs.com, 原www.ourdev.cn, 原www.ouravr.com

快速回复 返回顶部 返回列表