<< フィッティングで方程式を解いてみる(2)連立方程式 | top | 小林泰三「セピア色の凄惨」読了 >>

スポンサーサイト

一定期間更新がないため広告を表示しています

スポンサードリンク | - | | - | - |

フィッティングで方程式を解いてみる(3)デノイズ

前回、前々回、フィッティングで方程式を解いてみた。

例が簡単過ぎてあまり面白みを感じない方もいそうなので、もう少し凝ったことをしてみたい。

 

今回は、記事タイトルに偽りありであるが関数の最小化をしてみたい。

フィッティングが二乗和で表される関数の最小化をしていることを利用して方程式を解いていたのだが、今回は方程式ではなく、関数の最小化問題を解いてみたいと思う。

 

今回行うのは簡単なデノイズである。

いま、あるデータ{Xi}がある。本来、データは滑らかでなければならないとする。この仮定の下でノイズを含まないデータの推定を行ってみよう。

デノイズ後の値を{ai}とする。滑らかさの指標として、

¥[
¥sum_i(a_{i+1}-a_i)^2
¥]

を用いる。これが小さければ滑らかである。最も滑らかなのは全ての点で同一の値を持つ場合である。これではデノイズ出来たとはいい難い。{ai}は{Xi}にそれなりに近い必要がある。そこで

¥[
¥sum_i(a_i-X_i)^2 + ¥lambda¥sum_i(a_{i+1}-a_i)^2
¥]

を最小化する{ai}を求めることにする。λは定数で、これを大きくすると、滑らかであるが元のデータとの乖離が大きくなる。

この関数は二乗和となっているので、フィッティング機能を利用して解くことが出来る。

 

フィッティング対象のデータは次のようなものを用意する。データの次元をNとする。

 

0 X[0]

1 X[1]

・・・・

N-1 X[N-1]

N 0

・・・

2N-2 0

 

最初のN行は

¥[
¥sum_i(a_i-X_i)^2 
¥]

の部分である。後半は

¥[
¥lambda¥sum_i(a_{i+1}-a_i)^2
¥]

の部分である。

 

フィッティング関数は

¥[
¥sum_{i=0}^{N-1}¥delta_{xi}a_i + ¥sqrt{¥lambda}¥sum_{i=N}^{2N-2}¥delta_{xi}(a_{i-N+1}-a_{i-N})
¥]

と定義する。これで元の問題をフィッティング問題として定義できた。

 

それではやってみよう。今回は問題が少し複雑で流石にGnuplotに手打ちでは辛いので、次のコードでテキストを生成してGnuplotで読み込んだ。

 

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <vector>
#include <random>
#include <cmath>
double M_PI = 3.14159265358979;
using namespace std;

int main()
{

    FILE *fp = fopen("data.txt", "w");
    int N = 64;
    vector<double> val(N);
    std::mt19937 mt(0);
    std::normal_distribution<double> gauss(0, 0.1);
    for (int n = 0; n < N; n++)
    {
        val[n] = sin(2 * M_PI*n / (double)N) +gauss(mt);
        fprintf(fp, "%d %e¥n",n,val[n]);
    }
    for (int n = N; n < 2*N-1; n++)
    {
        fprintf(fp, "%d 0¥n",n);
    }
    fclose(fp);

    fp = fopen("macro.txt", "w");
    fprintf(fp,"d(i,j) = (i==j) ? 1 : 0¥n");
    fprintf(fp,"fit ");
    for (int n = 0; n < N; n++)
    {
        fprintf(fp,"d(x,%d)*a%d + ",n,n);
    }
    double lambda = 10;
    for (int n = N; n < 2*N-1; n++)
    {
        fprintf(fp,"%f*d(x,%d)*(a%d-a%d)", sqrt(lambda),n, n+1-N,n-N);
        if (n != 2 * N - 2) fprintf(fp," + ");
        else fprintf(fp," ");
    }
    fprintf(fp,"¥"data.txt¥" via ");
    for (int n = 0; n < N; n++)
    {

        fprintf(fp,"a%d ",n);
        if (n != N - 1) fprintf(fp,", ");

    }
    fclose(fp);

    return 0;
}

 

このコードは、オリジナルデータとしてsin関数を用意し、σ=0.1のガウスノイズを足し合わせてデータを生成している。

その後、Gnuplot用のコマンドを生成している。

具体的にはN=64, λ=1では次のコマンドを生成している。

d(i,j) = (i==j) ? 1 : 0
fit d(x,0)*a0 + d(x,1)*a1 + d(x,2)*a2 + d(x,3)*a3 + d(x,4)*a4 + d(x,5)*a5 + d(x,6)*a6 + d(x,7)*a7 + d(x,8)*a8 + d(x,9)*a9 + d(x,10)*a10 + d(x,11)*a11 + d(x,12)*a12 + d(x,13)*a13 + d(x,14)*a14 + d(x,15)*a15 + d(x,16)*a16 + d(x,17)*a17 + d(x,18)*a18 + d(x,19)*a19 + d(x,20)*a20 + d(x,21)*a21 + d(x,22)*a22 + d(x,23)*a23 + d(x,24)*a24 + d(x,25)*a25 + d(x,26)*a26 + d(x,27)*a27 + d(x,28)*a28 + d(x,29)*a29 + d(x,30)*a30 + d(x,31)*a31 + d(x,32)*a32 + d(x,33)*a33 + d(x,34)*a34 + d(x,35)*a35 + d(x,36)*a36 + d(x,37)*a37 + d(x,38)*a38 + d(x,39)*a39 + d(x,40)*a40 + d(x,41)*a41 + d(x,42)*a42 + d(x,43)*a43 + d(x,44)*a44 + d(x,45)*a45 + d(x,46)*a46 + d(x,47)*a47 + d(x,48)*a48 + d(x,49)*a49 + d(x,50)*a50 + d(x,51)*a51 + d(x,52)*a52 + d(x,53)*a53 + d(x,54)*a54 + d(x,55)*a55 + d(x,56)*a56 + d(x,57)*a57 + d(x,58)*a58 + d(x,59)*a59 + d(x,60)*a60 + d(x,61)*a61 + d(x,62)*a62 + d(x,63)*a63 + 1.000000*d(x,64)*(a1-a0) + 1.000000*d(x,65)*(a2-a1) + 1.000000*d(x,66)*(a3-a2) + 1.000000*d(x,67)*(a4-a3) + 1.000000*d(x,68)*(a5-a4) + 1.000000*d(x,69)*(a6-a5) + 1.000000*d(x,70)*(a7-a6) + 1.000000*d(x,71)*(a8-a7) + 1.000000*d(x,72)*(a9-a8) + 1.000000*d(x,73)*(a10-a9) + 1.000000*d(x,74)*(a11-a10) + 1.000000*d(x,75)*(a12-a11) + 1.000000*d(x,76)*(a13-a12) + 1.000000*d(x,77)*(a14-a13) + 1.000000*d(x,78)*(a15-a14) + 1.000000*d(x,79)*(a16-a15) + 1.000000*d(x,80)*(a17-a16) + 1.000000*d(x,81)*(a18-a17) + 1.000000*d(x,82)*(a19-a18) + 1.000000*d(x,83)*(a20-a19) + 1.000000*d(x,84)*(a21-a20) + 1.000000*d(x,85)*(a22-a21) + 1.000000*d(x,86)*(a23-a22) + 1.000000*d(x,87)*(a24-a23) + 1.000000*d(x,88)*(a25-a24) + 1.000000*d(x,89)*(a26-a25) + 1.000000*d(x,90)*(a27-a26) + 1.000000*d(x,91)*(a28-a27) + 1.000000*d(x,92)*(a29-a28) + 1.000000*d(x,93)*(a30-a29) + 1.000000*d(x,94)*(a31-a30) + 1.000000*d(x,95)*(a32-a31) + 1.000000*d(x,96)*(a33-a32) + 1.000000*d(x,97)*(a34-a33) + 1.000000*d(x,98)*(a35-a34) + 1.000000*d(x,99)*(a36-a35) + 1.000000*d(x,100)*(a37-a36) + 1.000000*d(x,101)*(a38-a37) + 1.000000*d(x,102)*(a39-a38) + 1.000000*d(x,103)*(a40-a39) + 1.000000*d(x,104)*(a41-a40) + 1.000000*d(x,105)*(a42-a41) + 1.000000*d(x,106)*(a43-a42) + 1.000000*d(x,107)*(a44-a43) + 1.000000*d(x,108)*(a45-a44) + 1.000000*d(x,109)*(a46-a45) + 1.000000*d(x,110)*(a47-a46) + 1.000000*d(x,111)*(a48-a47) + 1.000000*d(x,112)*(a49-a48) + 1.000000*d(x,113)*(a50-a49) + 1.000000*d(x,114)*(a51-a50) + 1.000000*d(x,115)*(a52-a51) + 1.000000*d(x,116)*(a53-a52) + 1.000000*d(x,117)*(a54-a53) + 1.000000*d(x,118)*(a55-a54) + 1.000000*d(x,119)*(a56-a55) + 1.000000*d(x,120)*(a57-a56) + 1.000000*d(x,121)*(a58-a57) + 1.000000*d(x,122)*(a59-a58) + 1.000000*d(x,123)*(a60-a59) + 1.000000*d(x,124)*(a61-a60) + 1.000000*d(x,125)*(a62-a61) + 1.000000*d(x,126)*(a63-a62) "data.txt" via a0 , a1 , a2 , a3 , a4 , a5 , a6 , a7 , a8 , a9 , a10 , a11 , a12 , a13 , a14 , a15 , a16 , a17 , a18 , a19 , a20 , a21 , a22 , a23 , a24 , a25 , a26 , a27 , a28 , a29 , a30 , a31 , a32 , a33 , a34 , a35 , a36 , a37 , a38 , a39 , a40 , a41 , a42 , a43 , a44 , a45 , a46 , a47 , a48 , a49 , a50 , a51 , a52 , a53 , a54 , a55 , a56 , a57 , a58 , a59 , a60 , a61 , a62 , a63 

 

このコマンドを実行すればGnuplotがフィッティングを実行する。

さて、結果を見てみよう。

まず、最初は理想データとノイズ入りデータである。

ガタガタしている。

λ=1でやってみよう。

少し滑らかになったのがわかるだろうか。

さらにλ=10ならばより滑らかになる。

目的関数の値もきちんと減少していることを確認した。

 

どうだろう、一応それなりに実践的なことも出来ることを示せたわけだが。面白くない?

役に立つかというと役には立たないわけだけど。

コード書かなくても最適化問題が解けるわけだが、単純な例なら今はWolfram alpha等便利なものがいくらでもあるし、複雑なものになると結局コードでコマンド生成しないと辛いものが・・・。

でも面白いからいいんだ。

ジャジャガッチ | 数学 | 22:42 | comments(0) | trackbacks(0) |

スポンサーサイト

スポンサードリンク | - | 22:42 | - | - |
Comment









Trackback
URL:

07
--
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
--
>>
<<
--
PR
RECOMMEND
RECENT COMMENT
MOBILE
qrcode
OTHERS
Since 2013/09/17
LATEST ENTRY
CATEGORY
ARCHIVE
LINKS
PROFILE
SEARCH