頭と尻尾はくれてやる!

パソコンおやじのiPhoneアプリ・サイト作成・運営日記


n次の近似曲線を求める(2)

n次の近似曲線を求める(1)
↑この続きだよ。
今回は最小二乗法から得られる連立方程式の各係数を求めて行列を作るからね。

まず、ガウスジョルダン法で計算する行列を配列一つに入れておくんだけど、どういう順番で入れるかを書いおくよ。
例えば二次関数で近似する場合、最小二乗法から成り立つ式は次のようなものだよね(前回のリンク先にあったのを写しただけだけど)。

数式1

このご時世に手書き文字なんてなかなかキュートでしょ。新しかったiPadで書いたんだよ。どうでもいい情報だね。
これは書き直すと次のような形になるよね。あったでしょ、行列で連立方程式を記述するのって。
コードではこの一番左と右辺の行列を順番に一つの配列に入れておくんだ。赤い数字で配列に入れる時の順番を書いてるけど、ゼロから書いてるから順番というより配列の添字だね。

数式2

不思議な感じがするだろうからもう一度言うけど、右辺の行列も左辺のに続けて入れていくからね。


それじゃコードを書こう。まずは係数を算出するところ。
前回紹介した最小二乗法を説明するページには二次関数、三次関数で近似する場合の式があるけど、ここからだけでもn次関数で近似する場合の行列はどういう風になるかは想像がつくよね。
落ち着いて考えたらどうということはないので、俺が書いたコードを理解しようとするのが億劫なら自分で書いた方が速いと思うよ。
ともかく俺が書いたのはこんな感じ。まずは順番にxの何乗とかのΣを求めるよ。
{
    //左辺側は 0 から degree x 2次までの乗数の総和が必要
    double *xpow = calloc(degree*2+1, sizeof(double));
    
    //右辺に関してはdegree+1のデータが必要
    double *xpowy = calloc(degree+1, sizeof(double));
    
    
    for(int ii=0;ii<[dotArray count];ii++) {
        
        CGPoint point = [[dotArray objectAtIndex:ii] CGPointValue];
        
        //左辺側の必要な要素を計算
        for (int jj=0; jj<= degree*2; jj++) {
            xpow[jj]+=pow(point.x, jj);//jj乗を足していく
        }
        
        //右辺側の計算
        for (int jj=0; jj<=degree;jj++) {
            xpowy[jj] += point.y * pow(point.x, jj);//y x x^jj
        }
    }
}
ここでdegreeはint型で、degree次関数で近似するってことだよ。degreeが3なら三次関数で近似するってことね。
それから元々の近似させたい点がNSArray *dotArrayにCGPointで入っている。通常NSArrayにCGPoint型データはそのままでは入れられないからNSValueのvalueWithCGPoint:メソッドを使って入れてる。使う時は上のようにCGPointValueで取り出せるからね。


次に上で計算した値を行列に順番に入れていくよ。上の方にある図のような順番で入れるからね。もう一度言うけど、左辺のも右辺のも一つの配列に入れるよ、え?何回も言い過ぎ?
{
    double matrixArray[(degree+1)*(degree+2)];
    
    for (int jj=0;jj<=degree;jj++) {
        for (int ii=0;ii<=degree;ii++) {
            matrixArray[jj*(degree+1)+ii] = xpow[jj+ii];
        }
    }
    //右辺側も一緒に入れる
    for (int ii=0;ii<=degree;ii++) {
        matrixArray[(degree+1)*(degree+1)+ii] = xpowy[ii];
    }
}
これで配列matrixArrayに必要なデータは入ったね。
この配列をガウスジョルダン法で解いてくれる関数へ渡すんだけど、長くなったので今日はここまで!



<< n次の近似曲線を求める(3)   TopPage  n次の近似曲線を求める(1) >>

コメント

勘違いだったら申し訳ありませんが、
このプログラムだとmatrixArray[0]が1に鳴ってしまいませんか?

2015.03.14      編集


管理者にだけ表示を許可する
 

トラックバック

トラックバックURL
http://ringsbell.blog117.fc2.com/tb.php/701-3f3e3bbe




Copyright ©頭と尻尾はくれてやる!. Powered by FC2 Blog. Template by eriraha.

FC2Ad