欧美成人精品视频一区二区三区,亚洲日韩a￠￥在线视频,欧美特大特白屁股xxxx

標(biāo)題: 最小二乘法曲線擬合 C語(yǔ)言實(shí)現(xiàn) [打印本頁(yè)]

作者: jialinx 時(shí)間: 2015-9-15 14:58
標(biāo)題: 最小二乘法曲線擬合 C語(yǔ)言實(shí)現(xiàn)
簡(jiǎn)單思路如下：
1，采用目標(biāo)函數(shù)對(duì)多項(xiàng)式系數(shù)求偏導(dǎo)，得到最優(yōu)值條件，組成一個(gè)方程組；
2，方程組的解法采用行列式變換（兩次變換：普通行列式——三角行列式——對(duì)角行列式——求解），行列式的求解算法上優(yōu)化過(guò)一次了，目前還沒(méi)有更好的思路再優(yōu)化運(yùn)算方法，限幅和精度準(zhǔn)備再修改修改
目前存在的問(wèn)題：
1，代碼還是太粗糙
2，數(shù)學(xué)原理可行，但是計(jì)算機(jī)運(yùn)算有幅度溢出和精度問(wèn)題，這方面欠考慮，導(dǎo)致高階大數(shù)據(jù)可能擬合失敗
下面附簡(jiǎn)單注釋版代碼：
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"

//把二維的數(shù)組與一維數(shù)組的轉(zhuǎn)換，也可以直接用二維數(shù)組，只是我的習(xí)慣是不用二維數(shù)組
#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))
//

/***********************************************************************************
從txt文件里讀取double型的X，Y數(shù)據(jù)
txt文件里的存儲(chǔ)格式為
X1  Y1
X2  Y2
X3  Y3
X4  Y4
X5  Y5
X6  Y6
X7  Y7
X8  Y8
函數(shù)返回X，Y，以及數(shù)據(jù)的數(shù)目（以組為單位）
***********************************************************************************/
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;
}

/***********************************************************************************
打印系數(shù)矩陣，只用于調(diào)試，不具備運(yùn)算功能
對(duì)于一個(gè)N階擬合，它的系數(shù)矩陣大小是（N + 1）行（N + 2）列
double* Para：系數(shù)矩陣存儲(chǔ)地址
int SizeSrc：系數(shù)矩陣大小（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;
}

/***********************************************************************************
系數(shù)矩陣的限幅處理，防止它溢出，目前這個(gè)函數(shù)很不完善，并不能很好地解決這個(gè)問(wèn)題
原理：矩陣解行列式，同一行乘以一個(gè)系數(shù)，行列式的解不變
當(dāng)然，相對(duì)溢出問(wèn)題，還有一個(gè)精度問(wèn)題，也是同樣的思路，現(xiàn)在對(duì)于這兩塊的處理很不完善，有待優(yōu)化
以行為單位處理
***********************************************************************************/
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;
}

/***********************************************************************************
系數(shù)矩陣行列式變換
***********************************************************************************/
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;
}

/***********************************************************************************
系數(shù)矩陣行列式變換，與ParaPreDealA配合
完成第一次變換，變換成三角矩陣
***********************************************************************************/
static int ParaDealA(double* Para, int SizeSrc)
{
      int i;
      for (i = SizeSrc; i; i--)
            if (ParaPreDealA(Para, SizeSrc, i))
                     return -1;
      return 0;
}

/***********************************************************************************
系數(shù)矩陣行列式變換
***********************************************************************************/
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;
}

/***********************************************************************************
系數(shù)矩陣行列式變換，與ParaPreDealB配合
完成第一次變換，變換成對(duì)角矩陣，變換完畢
***********************************************************************************/
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;
}

/***********************************************************************************
系數(shù)矩陣變換
***********************************************************************************/
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數(shù)據(jù)里面獲取系數(shù)矩陣
double* Para：系數(shù)矩陣地址
const double* X：X數(shù)據(jù)地址
const double* Y：Y數(shù)據(jù)地址
int Amount：XY數(shù)據(jù)組數(shù)
int SizeSrc：系數(shù)矩陣大�。⊿izeSrc）行（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;
}

/***********************************************************************************
整個(gè)計(jì)算過(guò)程
***********************************************************************************/
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;
}

/***********************************************************************************
測(cè)試main函數(shù)
數(shù)據(jù)組數(shù)：20
階數(shù)：5
***********************************************************************************/
int main(int argc, char* argv[])
{
      //數(shù)據(jù)組數(shù)
      int Amount;
      //XY緩存，系數(shù)矩陣緩存
      double BufferX[1024], BufferY[1024], ParaK[6];
      if (GetXY((const char*)"test.txt", (double*)BufferX, (double*)BufferY, &Amount))
            return 0;
      //運(yùn)算
      Cal((const double*)BufferX, (const double*)BufferY, Amount, sizeof(ParaK) / sizeof(double), (double*)ParaK);
      //輸出系數(shù)
      for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)
            printf("ParaK[%d] = %lf!\r\n", Amount, ParaK[Amount]);
      //屏幕暫留
      system("pause");
      return 0;
}

作者: weiyida 時(shí)間: 2021-1-15 10:14
static int GetXY(const char* FileName, double* X, double* Y, int* Amount)
{
   FILE* File = fopen(FileName, "r");//這里的FILE在哪里有定義�。�
      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;
}

作者: rundstedt 時(shí)間: 2021-1-23 17:50
最后擬合出來(lái)的是個(gè)啥？

作者: 點(diǎn)擊小子 時(shí)間: 2022-7-1 20:49
其實(shí)有更好的算法的，在圖形處理領(lǐng)域，如果按你這個(gè)算法，cpu會(huì)崩潰的。計(jì)算殘差可以參考：strchr的這篇文章 "Calculating standard deviation in one pass"

作者: Hephaestus 時(shí)間: 2022-7-1 22:21

rundstedt 發(fā)表于 2021-1-23 17:50
最后擬合出來(lái)的是個(gè)啥？

顯然是個(gè)n階多項(xiàng)式。

歡迎光臨 (http://www.torrancerestoration.com/bbs/)