Diary October, 2007

 

01日Mon.

10月。

新二年生の歓迎会。 小腹がすいていたので、おやつ代わりに頂きました。

近年はお酒が出ないのがさびしい。 公には彼らは未成年な訳ですので、 当然といえば当然ですが。

しかし日本の慣習として「大学生は飲酒可能」な流れがあるのは確か。 皆、新入生コンパで酒飲んでる訳で。

近年飲酒運転はかなり厳しく取り締まっている(当然ではあるが。駄目、絶対。飲酒運転!)が、 こと年齢のことは無法地帯のような気がする。

世界に目を向けると、二十歳から飲酒可能というのは少数派で、 かなりの国は18歳以上からである。

厳しく取り締まるのであれば、そうすればよし、 黙認するのであればいっそ法律を変えてしまったほうがいいのでは、と思える。 なぜなら、「破ってもいい法律」なんてものは本当はあってはいけないからである。 そんなものがあったら「法」そのものの権威を失うことになる。

Gibbs samplerは計算ライブラリで計算することが出来る分布の場合、 採択率が1となるメトロポリス・ヘイスティング法の特別な場合だと理解できる。 例えば、ある変数がガウシアンに従う場合、 そこからのサンプルすることは容易である。

P(s|d,C_l) \propto \exp\left[ -\frac{1}{2} ^T(d-s)N^{-1}(d-s)               \right]     \exp\left[  -\frac{1}{2} ^Ts_lS_l^{-1}s     \right]

の様にデータ回帰、事前分布を分けて書いた場合、 これからsをサンプルするのは難しい。 メトロポリス・ヘイスティングを使えば出来なくも無い。

一方上式を変形して、

P(s|d,C_l) \propto \exp\left[ -\frac{1}{2} ^T(s-\hat{s})(S^{-1}+N^{-1})(s-\hat{s})           \right] ,\quad \hat{s} = (S^{-1} +N^{-1})^{-1}N^{-1} d

と変形することで、sをガウシアンからサンプルすることが可能である。Gibbs samplerである。 ここでは、今注目しているs以外の項を適当に追加しても、分布の形が変化しないことを用いている。

しかし、一方で問題がある。 それは試行後とに、逆行列(S^{-1}+N^{-1})^{-1}を計算し無ければならないことである。 この計算オーダーはO(N_{\rm pixel}^3)であるので、 現実的ではない。CMBが一様等方で揺らぎオンリー、 ノイズが平均零、分散ということから、最適化が可能である。 \hat{s} = (S^{-1} +N^{-1})^{-1}N^{-1} dはWeiner filterと呼ばれるもので、 信号処理等の分野でよく用いられる(統計でも出てくる)。

この枠組みで考えると、O(N_{\rm pixel}^2), O(N_{\rm pixel}^3/2)のオーダーになる。

この計算オーダーと、 ガウシアンに変形する前の分布からのメトロポリス・ヘイスティングでのサンプリングでの必要な計算オーダーの兼ね合いで、 どちらを使うかが決まるのだと思う。

ただし、gibbs samplerは「確実」な分布なので、メトロポリス・ヘイスティングにはぶが悪いのだろうか?

どうもGibbsでできるなら、Gibbsでやった方がよいといった風潮があるように見える。 特に宇宙論では。

上手く分布作る方法は、何もここで述べた二つだけではなく、 いくつか存在する。 MHの拡張のようなものである。

これをうまく使えば、 Gibbs samplerにこだわる必要も無い気がするが・・・・。

ude.jpg

中間発表・・・・。

02日Tue.

今、component separationの簡単な例題として単一ピクセルでの差し引きを行っている。 オフセット、振幅とそれのスペクトルインデックス。 パラメータ数は3、 データ周波数は4。

d_\nu=m\delta_{\nu,\nu_0}+A\leftfrac{\nu}{\nu_0}\right)^\beta

パラメータ自身が縮退しており、 そもそも自由度も小さい。

分布は左右非対称な歪な形をとる。

単一ピクセル、少ないパラメータならばグリッドを切ってすべての場合を計算できる。 この場合、分布は歪でも必ずもとまる。

一方これをMCMCで行おうとすると、 原理的には可能だが、なかなか難しい。 他の2つのパラメータについてマージナライズした分布は再現できるが、 一つだけをマージナライズしたcontourの再現性は悪い。

これはそもそも問題の設定が悪いのであって、 どうしようもない。 MCMCを使えばなんでも解ける、というのは間違いである。

現実の場合は、多ピクセルである。 各点各点で上と同じことをするのではなく、 ピクセルで縮約をとった場合、 大数の原理から、分布は正規分布に近づくことが予想される。

全体の確率を計算してchainを設計する際、 どのように提案を行うかが問題となる。

一回に何個変化させるのか、グループ化の問題(少ないと計算回数が増える。マルコフ確率場の問題もあるような)、 皆同時に変化させるのか(採択率悪そう)。 と、色々と問題がある。

ただし、空間統計をうまく使えば、 上記の不都合は上手くいく方向に進むことが予想される。

今後はテストを行いながら、 実際の観測に適応できるよう、 ツールの使い方(ツールといってもI/Oライブラリ集)を練習するつもりである。

kun.jpg

03日Wed.

規格化定数が合わないと思ったら、 P(d)を計算していなかった。

foregroundのモデルを記述するパラメータはほとんど線形である。 理想を言ってしまえば、すべて線形である(あるピクセル、ある周波数での、成分分だけパラメータがある)。

しかし、特徴的な周波数依存から、spectral indexで記述することがある(できる)。 これでパラメータの数が減る。

一方の弊害として、非線形であるためサンプリングが困難になる。 (温度というパラメータに対して、温度の周波数依存を決めてしまうパラメータは、一つ上の階層の様に思われる。 spectral indexが決まれば、それに従う領域で、温度が決まってしまう)

縮退を解く方法として、 事前分布の導入がある。 これはベイズからの自然な帰結である。

ジェフリーズ事前分布はその一つである。 そのほかに、滑らかな制約条件・・・・とかを考えているのですが、 どうでしょう。まだ試せていません。

あかりも新しいデータが届きました。でも修論には間に合いません。

tarako.jpg

04日Thu.

AKARIの解析に多くの時間をとられている現状を考えると、 AKARIの解析も修論の一部に含めるべきなんだと考えるようになった。 つまり中間発表でも話すと。

CMB, AKARIとも論文にできるようなレベルに達していない。 AKARIの方は外的要因だが、 CMBはそうではない。

でも今取り組んでる簡単な場合の再現ですら、 そのコーディングは手間が掛かる。

Eriksenの論文を全部再現するには、 最低でも一年ぐらいこつこつとコードを書いていかないといけない気がする。 一度かければ、これを改良して色々できる。 現にEriksenはここ数年で同じような論文(悪い意味ではない)をいくつも出している。

再現をさせながらコード開発、 並行して自分のアイディアの定式化、 実験をするのが吉。

取り組み方に失敗したと痛感。

AKARIは変な結果出るし、 Reply-Toかえるのめんどいし・・・。

今日は徹夜になりそうです。

pu.jpg

07日Sun.

とりあえず毎日大学へ。 休日関係なし。

お仕事をこなす。 が、進まない。

data2Clは現実的に考えて修論では出来ない。

  • data2Cl
    • 線形パラメータのGibbs samplerの実装
      • Conjugate Gradient (CG): 共役勾配法の実装
    • 非線形パラメータのsampler
      • MH法の実装
      • inversion sampler
    • 独自のforeground modelの実装
      • 場所に依存する線形パラメータの空間方向について、事前情報を加える。synchrotron amplitudeが候補。これが変化することで、spectral indexにも影響が?
      • 線形なので、解析解が求まる。
    • 実空間から波数空間への変換の実装、またその逆の実装
      f(\gamma)=\sum_{l=0}^{l_{\rm max}}\sum_{m=-l}^{l} a_{lm}Y_{lm}(\gamma) \quad\leftrightarrow \quad \hat{a}_{lm} = \frac{4\pi}{N_{\rm pix}} \sum_{p=0}^{N_{\rm pix}-1} Y_{lm}^*(\gamma_p)f(\gamma_p)
    • power spectramの計算
      \hat{C}_l = \frac{1}{2l+1}\sum_{m} |\hat{a}_{lm}|^2

独自の点以外は、 既存のライブラリを使えばいい。車輪の再発見は無意味である。 しかし組み合わせるのが困難である。

CGに関しては、よく研究されている線形・行列計算のライブラリを使えばよい。 FortranかCだろうか? mapに関するのとの大半はHealpixで解決できる。 これはfortran, IDL, C, C++で用意されている。 C++を触ってみたところ、上のことはほぼすべてできるようである。 しかし、template, classを総動員で書かれている(当然か)ので、 まだなれない。

これらのこと関しては、こつこつとやっていくしかない。 今現在出来ることは、独自の事前情報の導入についてである。 実のところこれが一番重要なことではあるが。 これを行うことで、改善されるという確かな根拠を構築しさえすれば、 結果はそのうちということでどうにかなる。

この空間方向の事前情報は、 多くの天文の人にとってなじみの無いはずであるから、 納得させるのは一苦労だと思う。 実際のところ、例を提示して効果を見せたいところだが。 (CMB解析の問題点を例示するにはどっかの論文を持ってくればいいわけだが・・・。)

少なくとも、この点に関する自分の結果は見せたい。

すぐにhealpixをマスターするのは現実的に無理なので、 これを使わず、平面でモデルを構築、テストをやってみる。

click.gif
click.gif
click.gif
click.gif

09日Tue.

Cat's eye nebulaの単純平均のイメージが何とかできた。 データのミスだと思っていたことは気のせいだと思う。

中間発表まで三週間を切った訳だが、 それそれトラペを各準備しましょうか。 図とかグラフは可能ならば自分で作ったものを入れたい。 その方が理解も完璧だし。 論文から引っ張ってきた図は、 理解した気でも質問されると「そういえばどういうことだ?」ということが間々ある。

自分の成果が出ていない現状はどうしようもないので、 それでもまともな発表が出来るよう、周到に準備できればと思う。

早起きしなければ!

ne.jpg

10日Wed.

ここ一週間ばかり懸案となっていたことが解決した・・・と思う。 ようやく次のステップに移れる。

desck.JPG

今日も早起きできなかった。

なので、実働12時間ちょい。 それでもなんだか、目がちかちかするのは、 パソコンばっかりいじっているからであろうか。

me.jpg

11日Thu.

昨日解決したことをまとめて、ひと段落。 やはり上手くいくと気分がいい。

早起きも出来た。

CMB, AKARIとも前振り(の前振り)が終わったので、 ようやく本番に入っていきましょう。

umi.jpg

13日Sat.

Healpixが使えるようになってきた。 色々といじっている。 a_{lm}C_lの計算なんかがすぐに(=命令一つで)出来る。 FITSのI/Oも実装されている。 pixel <-> (\theta,\phi)の変換も出来る。

以上の機能は、 実空間とY_{lm}空間を行ったり来たりする際に大変便利である。 実空間ではノイズは(ほぼ)独立(自己相関しかない)であり、 Y_{lm}空間だと、&mimetex({\bf C})は厳密に対角成分しか存在しない。

これらの機能を手で実装するのは、結構面倒である。 ソフトの開発者に感謝。

発表まで約二週間

ko.jpg

朝晩はめっきり寒くなって


15日Mon.

なぞの点源発見。 他のInfraredの観測では見えない。 (alpha,delta)で固定されているように見える。 でもtotal averageでは消える。 見えるのは各rotation毎の絵。

cross talkとなるbeamパターンをはじめて書いた。 無駄に自分のプログラムで実装してしまった。 最初はさらに無駄なアルゴリズムで・・・。

FITS化した絵の解析はIRAFを使った方がいいと悟った。 個別の複雑な問題(blind deconvolution)に対応するためには、自前で書く必要があるが、 差和etc等のよく使う操作をいちいち自前で書いていたらきりが無い。

CMBの方は手が進まない。 Healpixが動くようになったのはいいが、 そもそも作業が出来ない現状。 二つのことを同時進行させるのは大変ですね。

中間発表も真近に迫ってきているので、 無駄な作業をしている暇は無い。

それにしても、皆研究がまとまっていて偉いと思う。 つまりは、自分が偉くないということだが・・・。

01.jpg
02.jpg

16日Tue.

FITSを使った解析は様々なツールとの連携が大事である。 自分がプロットした画像と、 外部カタログとを重ねたり、 参照したりすることがしばしばある。 これをいちいち自前で行っていたら、 大変である。

車輪の再発見はしないのが賢いやり方である。

今までFITSを見るのにds9を使っていたが、 外部カタログとの連携、特にIRASとの連携が上手くいかないのが悩みの種であった(自分の無知ゆえかも知れないが)。 色々いじった結果、web browserとの連携が優れている、javaベースのoasisが気に入った。 コマンドラインベースのツールも便利だが、 画像を扱うので、目に見えて操作できるのはやはり便利である。

22日Mon.

テラヤバス

・・・・・・・・・・・・・・

パソコンのメモリが足りないと感じ始めた。 CMB,AKARI共に。

そのうち、大規模な計算をやるときはお金払ってどっかのPCクラスターを使うつもりではいる。 こことか。

Eriksenたちには敵わないが、 国内ではここが限界。 他のところも似たり寄ったり。 大体十個のオーダー。

因みにEriksenは2000CPUでやってます。

今までのやり方ではバイアスがかかった結果となる、 という思惑通りの結果がでた。

やる前から分かっていたが、 synchrotronが無いところでは、 indexは決まり用が無いので、分布が変になる。 ノイズが10あるのに、30,20,10,10の四つのデータから、 spectral indexを求めるのは、そもそも不可能である。 indexが決まらないというのがむしろ正しい結果である。

下限値を設定しておかないと、とんでもないことになる。 こればoffsetにも影響を与えるので、注意が必要である。

amplitudeはそれなりに決まる。 推定されたamplitudeが小さいところがbetaがよく決まらないところである。

次は、事前情報を加えた場合、 どのようにバイアスが変化するか焦点となる。

上述の通り、amplitude,indexに事前情報を組み込むことで、 データそれ自身からは制限がつかないことについて、 制限をつけましょうというものである。

近頃流行の・・・(注意、音がでます)で現実逃避。 耳に残ります。 気に入ったならパソコンの時計をこちらにしてはいかがでしょうか・・・・。

L先生、申し訳ないです。 アブストは明日の朝に書きます。

hon.jpg

25日Thu.

終わらんがな。

「逃げちゃだめだ、逃げちゃだめだ」とシンジ君は言ってましたが、 「逃げた方が辛いことを知っているから逃げない」とも言っていました。

hasa.jpg

31日Wed.

修論中間発表は無事終わりました。 取りでした。 時間オーバーしました。

直ぐ他の仕事です。

学振が届きません。