スポンサーサイト

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

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

フィッティングで方程式を解いてみる(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) |

フィッティングで方程式を解いてみる(2)連立方程式

前回、f(x)=cをフィッティングで解いてみたわけだが、今回は連立方程式でやってみる。

基本は前回と同じである。

解きたい方程式を

¥[
¥sum_j A_{ij}x_j = b_i
¥]

とする。

次のようなデータを用意する。

0 b0

1 b1

2 b2

・・・

 

フィッティング関数は

¥[
f(x;¥{a¥}) = ¥sum_j A_{ij}a_j¥delta_{xi}
¥]

とする。aが推定するパラメータ、δはクロネッカーのデルタである。

フィッティングすれば

¥[
 ¥sum_i(¥sum_j A_{ij}a_j-b_i)^2
¥]

を最小とする解が求まる。これにより元の方程式の解が求まる。

具体的にやってみよう。

x+y=3

2x-y=0

を解いてみる。明らかに解はx=1, y=2である。

”data.txt”として

0 3
1 0

を用意する。

クロネッカーのデルタを定義する。

d(i,j) = (i==j) ? 1 : 0 

後は上で述べたようにフィッティング関数を定義してフィッティングすればよい。

 fit d(x,0)*a+d(x,0)*b+2*d(x,1)*a-d(x,1)*b "data.txt" via a,b

はい、求まった。

ジャジャガッチ | 数学 | 23:17 | comments(0) | trackbacks(0) |

フィッティングで方程式を解いてみる

実に半年ぶりのブログ更新。

Fallout4やったり、ウィッチャー3やったりしてた。

あんまり家でガッツリ趣味やる時間がなくて放置していた。

会社ではいろいろ面白いことやってるんだけどそれをここに書くわけにはいかないし。

趣味やるくらいなら英語やらねばというのもある。ここしばらく英語は勉強してなくて、3月に入ってから勉強を再開した。

 

本題。

今日は以前ふと思いついたことを実際にやってみることにした。

それは・・・、フィッティングを使って方程式を解くこと!

さて、どういうことをするかは後で説明することにして思い出話から。

高校の数学では、基本的に解析的に解けない方程式は出てこない。

しかし、大学に入ると解析的に解けない方程式のオンパレード。

ここで面食らうわけだが、大学1年生のときに先生にいいことを教わった。

グラフソフトでグラフ書けば数値解簡単に求められるよ、と。

例えば

x^2 -2x = -1

を考えよう。

この数値解はグラフソフトでy=x^2-2xとy=-1をプロットして交点の座標を読み取ればいい。

この例は解析的に解けるが、解析的に解けない方程式も同じだ。

 

さて、ここからが本題。

同じグラフソフトを使って解を求めるならフィッティング機能使えないかな?

グラフソフトにはフィッティング機能がついていることがある。

例えば僕が愛用しているGnuplotは

fit a*x "data.txt" via a 

と入力すれば直線フィットが出来る。

与えたデータに最も近くなるようにパラメータを推定してくれるわけだ。

これを利用して方程式を解けないだろうか。

 

フィッティングは次のようなことをしている。

対象のデータを{xi,yi}としよう。これをパラメータ{a}の関数fでフィッティングするとは

¥[
¥sum_i (y_i-f(x_i;¥{a¥}))^2
¥]

を最小化する{a}を求めることに他ならない。

 

それではこれを利用して方程式を解いてみよう。

ここでは

g(a) = c

をaについて解くことにする。cは定数である。

Gnuplotではデータは2点以上必要である。

次のようにデータを用意する。

0 0

1 c

フィッティング関数は

f(x;a) = g(a)*x

とする。

するとフィッティングすることにより求まるパラメータaは

¥[
(g(a)¥times 0-0)^2 + (g(a)-c)^2
¥]

を最小化するものとなる。これはg(a)=cである。

つまり数値解が求まる。

 

さて、x^2-2*x = -1をこの方法で解いてみよう。

"data.txt"として

0 0
1 -1

を用意する。

fit (a**2-2*a)*x "data.txt" via a

を実行すれば

はい、ばっちしa=1が得られました!

いや、だから何?なんだけど(笑)

役に立たなくても面白いものは面白いじゃない?

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

数学定数の分数表現探してみた

355/113的なやつ探してみた。

https://jajagacchi.github.io/pi2/pi.html

ジャジャガッチ | 数学 | 23:39 | comments(0) | trackbacks(0) |

三角関数の値を素朴に計算してみた

三角関数表ってどうやって作ってんだろうっていう疑問から素朴に三角関数の値計算してみた。

 

https://jajagacchi.github.io/tri/tri.html

ジャジャガッチ | 数学 | 23:05 | comments(0) | trackbacks(0) |

モンテカルロ法で円周率を計算してみる

誤差解析してみた。

https://jajagacchi.github.io/pi/pi.html

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

円周率は3.14か

ゆとり教育批判真っ只中の頃は円周率が3なんて!という批判がよく聞かれた。結局これはデマだったようだが。

円周率3が批判される理由のひとつとして、円周率が3だとその円に内接する正六角形と周長が等しくなってしまう、というものがある。

でもこれは3.14でも同じ話だ。円周率が3.14の場合は正57角形が3.1400*で少し円より周長が長くなる。

批判するならこういうことを踏まえた上で批判しましょう。

 

ゆとり教育の円周率3は、場合によっては3で代替するというくらいの話だったようだ。

これは賛成だ。概算という考え方は非常に重要だ。

買い物するときにこれは108円でこれは216円でみたいな足し算中々しないでしょう。

大体これが100円でこれが200円で大体いくらかな、という考え方をするでしょう。

例えば、次のような問題を考える。

 

問. A君は半径2の円の面積が30よりも大きいと言っています。真偽を答えよ。

 

愚直に3.14×2×2を筆算しちゃう子も出てくるんじゃないかと思うんだけど、これは円周率を3ないし安全策で4ととれば瞬時に判断できる。

こういう判断が出来るのは実生活でも科学でも重要じゃないでしょうか。

 

あと関係にないけど曲がった空間では円周率は3.1415…じゃないよ。みんなボールに円を描いて円周率考えてみよう!この場合、円周率は描く円の半径に依存する。面白いね。

 

ジャジャガッチ | 数学 | 21:49 | comments(0) | trackbacks(0) |

チューリングパターン

最近会社で英語(会話)を使わなくてはならなくなって刺激的な毎日。
毎夜英語勉強してます・・・。
暇な時間は、あれは英語で何て言うんだろう、って調べてみたりしてます。
最近はgoogleでワード検索したら英訳も勝手にしてくれるのね。
今日は"よっこいしょ"を調べてみた。
 

へー、英語ではBy Koisho!って言えばいいんだ!
・・・っって、んなわけあるかーい!
よっ、がbyに変換されることが驚き。
流石にこの場合に騙される人は多くないだろうけど、もっと普通の内容を英語で何て言うか調べた時に紛らわしい。
機械翻訳なのであんまり信用できないのに、これだけでかでかと表示されると正しい表現だと騙されてしまう人も結構いるんでないか。よろしくない。音声付きだし余計騙されそう。機械音声だけど。
計算とかグラフとかはいいんだけどこれはやめてほしいなあ。
僕は参考にもせず完全に無視することにしてます。
ちなみによっこいしょを英語で何というかは謎のまま。

/////////////////////////////
チューリングパターンというものを知っていますか。
この前所さんの目がテンでやっていたらしいのだが、ある種の動物の体の模様はチューリングが考えた偏微分方程式に従うそうだ。
チューリングというのはチューリングマシンで有名な、あのアラン・チューリングだ。
番組をぜひ見たかったのだが、妻が録画を消した後に教えてくれたため見れなかった。非常に残念。

で、この方程式がすごいのは模様なので空間変数が入っているのは当然なのだが、時間変数も入っていること。
つまり、模様が経時的にどのように変化するか予想することが出来る。
なんでも最近あるお魚の模様を経時的に追うとチューリングパターンで予測されるパターンに近かったという結果が出たらしい。
しかもこの方程式、いわゆる反応拡散方程式というやつでその名の通り、物質の反応、拡散を記述した方程式だ。
つまり、ある種の動物の模様が変化するのは物質の反応と拡散のようなことが起きた結果である可能性がある。
細胞間で何かやりとりしてるのかな。
上述のお魚の模様に関する論文も入手したのだがまだフォローできていない。
偏微分方程式を解くプログラムも書いてみたのだが、どうもあまりうまくいっていない。
パラメータが悪いか、バグがあるか、解法が悪いか・・・。
うまくいったら公開します。
ジャジャガッチ | 数学 | 22:27 | comments(0) | trackbacks(0) |

CADでルーローの三角形

年始に予告したとおり、夜は時間が出来たら英語の勉強をするようにしている。
残念だが今年は技術的な取り組みの時間が減るだろう。
と言いつつcadでルーローの三角形描いたり無駄なこともしてたりする。

そういえば、おととい大失敗した。
妻のスマホがソフト的に調子が悪いので夜中にいじっていたのだが、再起動した際にPINコードなるものの入力を求められた。
適当に入れれば当たるかな、と思って入力してたら3回失敗して今度はPUKコードを入れないと使えなくなってしまった。
ぐぐると申込書控えに書いてあるとか。
色々書類探したけど見つからない・・・。
結局妻は昨日ロックを解除してもらいに店舗まで行きましたとさ・・・。

ごめんよ・・・。しばらく頭が上がりません・・・。

/////////////////////////////////////////////////////////
電子工作をするようになってからたまーにCADも使うようになった。
使っているのはjw_cad。

今回はjw_cadでルーローの三角形を描いてみよう。
まず正三角形を描きます。
次に各頂点から正三角形の一辺を半径とする円を描きます。
次に円のいらない部分を削除します。
最後に正三角形を削除すればルーローの三角形の出来上がり。
じゃあついでに外接する正方形も描こう。面倒くさいし簡単だから解説なしでいいよね。



 
ジャジャガッチ | 数学 | 21:27 | comments(0) | trackbacks(0) |

スマホのゲームでガチャは何回まわす必要があるのか

最近は家族3人で近くの川に鴨を見に行くのが楽しみだ。
それでこの前妻が鴨の柄のキーケースをプレゼントしてくれた。これは嬉しい。癒される。
僕のお気に入りだ。
ちなみに僕の財布にははしっこにペンギンがいる。これも妻からのプレゼントだ。
30歳だけど・・・、いいよね!

///////////////////////////////////////////
巷ではずいぶん前からスマホのゲームが流行っている。
これらのゲームにはよくガチャというシステムが導入されているらしい。
これはある確率でレアアイテムが入手できるくじだ。
さて、一回あたりの当選確率pがわかっているとき、何回ガチャを回せば目的のアイテムが手に入るだろう。

時々勘違いされているのは例えば当選確率5パーセントなら20回引けば100パーセント当たる、という考え方。
これは大間違い。
計算してみると、当選確率が低い場合は当選確率によらず約63パーセントとなる。
上の例における20回(=1/p)を想定回数と呼ぶことにするとpが小さいとき、90パーセントの確率で当選するためには想定回数の2.3倍ガチャを回す必要がある。
この2.3は
-ln(1-0.9)~2.3からきている。
もっと安全に95パーセントで当てたければ
-ln(1-0.95)~3
より想定回数の3倍が必要だ。同様に99パーセントなら4.6倍が必要となる。
ガチャを回すときは想定している回数の2倍程度は回す覚悟が必要になるということだ。
詳しい計算はこちら→gacha.pdf
計算間違い、考え違いがあれば教えてください。
 
ジャジャガッチ | 数学 | 21:09 | comments(0) | trackbacks(0) |
1/8PAGES | >> |

06
--
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
--
>>
<<
--
PR
RECOMMEND
RECENT COMMENT
MOBILE
qrcode
OTHERS
Since 2013/09/17
LATEST ENTRY
CATEGORY
ARCHIVE
LINKS
PROFILE
SEARCH