GSLで多項式回帰を試してみた■■■

今回は、前回の続きとして、多項式回帰を試してみました。

まえがき

今回試してみたのは、Ｙ＝Ｃ０＋Ｃｎ×Ｘ＾Ｃｎの多項式回帰です。
ｎ値を変えることで以下の様に変化します。
ｎ＝２の時

ｎ＝３の時


ｎ＝４の時


ｎ＝１０の時

サンプル

サンプルコードは、ここです。
※上記ｎ値は、変数名Ｐｗです。

ＤＬＬコード

※前回のコードをそのまま残し、修正しました。

#pragma once
#define WIN32_LEAN_AND_MEAN 
#include <Windows.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>//数学用
#include ".\include\gsl\gsl_sort.h"//配列並び替え用
#include ".\include\gsl\gsl_wavelet.h"//ウェーブレット用
#include ".\include\gsl\gsl_multifit.h"//重回帰

BOOL APIENTRY DllMain(HANDLE hModule,DWORD ul_reason_for_call,LPVOID lpReserved)
{
//----
switch(ul_reason_for_call)
{
 case DLL_PROCESS_ATTACH:
 case DLL_THREAD_ATTACH:
 case DLL_THREAD_DETACH:
 case DLL_PROCESS_DETACH:
 break;
}
//----
return(TRUE);
}

//+-------------------------------------------------------------------+
//|重回帰（Multifit_liner)Y=c0+c1X+c2X^2+c3X^3...cpwX^pw
//|input[]   入力用バッファ
//|pw        変数量
//|out[pw+1]    出力用バッファ[c0,c1,c2,c3,....cpw.chisq]
//|Ns       対象データ数
//+-------------------------------------------------------------------+
__declspec(dllexport) void __stdcall Multifit_liner(double *input,int pw,double *out,int Ns){
 int i,ipw;
 double    ix,chisq;
 gsl_matrix *X, *cov;
 gsl_vector *y, *w, *c;
 
 //配列（ベクトル配列）の確保
 X = gsl_matrix_alloc(Ns, pw);
 y = gsl_vector_alloc(Ns);
 c = gsl_vector_alloc(pw);
 cov = gsl_matrix_alloc(pw, pw);

 //配列（ベクトル配列）入力
 for(i=0;i<Ns;i++){
  ix = i*0.001;
  gsl_matrix_set (X, i, 0, 1.0);
  for(ipw=1;ipw<pw;ipw++)gsl_matrix_set (X, i, ipw, pow(ix,ipw)); 
  gsl_vector_set (y, i, input[i]);
 }
  
 //処理
 gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (Ns, pw);
 gsl_multifit_linear(X,y,c,cov,&chisq,work);
 gsl_multifit_linear_free (work);
 
 //抽出
 for(i=0;i<pw;i++)out[i]=gsl_vector_get(c,(i));
 out[pw]=chisq;

 //メモリ解放
 gsl_matrix_free(X);
 gsl_matrix_free(cov);
 gsl_vector_free(y);
 gsl_vector_free(c);

}


//+-------------------------------------------------------------------+
//|ウェーブレット変換（DWT）
//|input[]  入力用バッファ
//|out[]    出力用バッファ
//|Ns       対象データ数（２のNs乗)
//|D        ウェーブレット係数（2:Harr 4~daubechies
//|fs       フィルタ
//+-------------------------------------------------------------------+
__declspec(dllexport) void __stdcall DWT(double *input,double *out,int Ns,int D,int fs){
 int i;
 int siz = (int)pow(2.0,Ns);
 if(D%2==1)D++;
 if(D>20)D=20;
 
 //配列のコピー
 memcpy(out,input,sizeof(double)*siz); 
 
 //ウェーブレットの初期設定(ウェーブレット名,D値）
 gsl_wavelet *w;
 if(D==2){
  w = gsl_wavelet_alloc(gsl_wavelet_haar,D);
 }else{
  w = gsl_wavelet_alloc(gsl_wavelet_daubechies, D);
   }
  
 //作業領域の確保(数量）
 gsl_wavelet_workspace *work;
 work = gsl_wavelet_workspace_alloc(siz);
 size_t *p = (size_t *)malloc(siz * sizeof(size_t));
 double *abs = (double *)malloc(siz * sizeof (double));
 if (abs == NULL || p== NULL){printf("メモリを確保できません"); exit(0);}
 
 //変換(初期設定,データ,進み幅,作業エリア）
 gsl_wavelet_transform_forward(w, out, 1, siz,work);
    
 //フィルタ作業
 for(i = 0; i < siz; i++) abs[i] = fabs(out[i]);
 gsl_sort_index(p, abs, 1, siz);
 for(i = 0; (i + fs) < siz; i++) out[p[i]] = 0;
    
 //逆変換
 gsl_wavelet_transform_inverse(w, out, 1, siz, work);

    //メモリ解放
 free(abs);
 free(p);
 gsl_wavelet_workspace_free(work);
 gsl_wavelet_free(w);
}
//----------------------------------------------------------------------------

まとめ

適用範囲を広くとりすぎたり、ｎ値を大きくしすぎるとバグが発生する可能性がありますので注意してください。