2012/02/01

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

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

まえがき

今回試してみたのは、Y=C0+Cn×X^Cnの多項式回帰です。
n値を変えることで以下の様に変化します。
n=2の時
1m
n=3の時
3m

n=4の時
4m

n=10の時
10m

サンプル

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

DLLコード

※前回のコードをそのまま残し、修正しました。
#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       対象データ数(2の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);
}
//----------------------------------------------------------------------------

まとめ

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