統計力学入門-化学の視点から-

本

学部時代に統計力学の講義を選択しなかったことが何となく心残りだったところに、たまたまamazonがこの本を押してきたのでノリで購入して読んでみた。

統計力学入門:化学の視点から

統計力学入門:化学の視点から

基本的には、エネルギーや粒子数を一定とするなどの仮定の下で可能な粒子の状態(微視的状態)のうち、組み合わせ数(微視的状態の重み)が最大となるものをLagrangeの未定乗数法で解いて、得られた式の解釈と熱力学における知識と照らし合わせるという流れであった。
あとは統計力学上での仮定が異なる場合や古典力学ではなく量子力学が必要な場合など、様々な条件のもとでの結果が端的にまとめられており、200ページくらいでかつ密度もそこまでないので気軽に読めた。
個人的には機械学習をいくばくかかじっていたことが功を奏して、基本的な考え方がすんなり理解できたのではないかと思う。

なかなか生物学へ応用のしがいのありそうで、大変参考になった。
気が向いたらもうちょっと専門的な統計力学の本も読んでみたい。

ドリフト項付きブラウン運動と線形回帰の関係

機械学習 確率過程

最近、いろいろな事象を確率過程の観点から捉えられないかと思考するうちに、不意に線形回帰との関係性に気付いたのでまとめてみる。

単純なブラウン運動は面白くないので、以下のようなドリフト項つきブラウン運動を考えていた。

{ \displaystyle N(x_t|x_{t-1}+t\mu, t\sigma^2) }

ここでμが固定だといまいちバリエーションに貧しいので、これも時間変化するような場合を考えた。
ただしなにも制限を加えないとよくわからないことになるので、ここでは以下のように何らかの値の線形結合で表現されるとした。

{ \displaystyle \mu = \sum w_i \theta_{i,t} }

ここでwは重みパラメータで時間に依存せず一定、Θは時刻依存的に変化するがその値は観測されるとする。

イメージとしては、xが何らかの農作物に関わる株価を表し、Θが天気や気温とかを表しているとする。
要はxの変動がランダムなブラウン運動的な挙動を示しながらも、Θに依存してある程度の方向性を示すような場合をモデル化し、重みパラメータwを推定することでどの要素がxの挙動にどのように影響を与えるかを知ることはできるか試したかったのだ。

それで、ドリフト項付きブラウン運動に戻ると、以下の式となる。

{ \displaystyle N(x_t|x_{t-1}+t\sum w_i \theta_{i,t}, t\sigma^2) }

さてここで株価のようにxを一定間隔の時間で観測できるような場合を考えると、

{ \displaystyle tw_i=w'_i\\ t\sigma^2=\sigma'^2 }

とでき、

{ \displaystyle N(x_t|x_{t-1}+\sum w'_i \theta_{i,t}, \sigma'^2) }

となる。


さて、上記式をちょっと移行して以下の通りに変形すると

{ \displaystyle N(x_t-x_{t-1}|\sum w'_i \theta_{i,t}, \sigma'^2) }

なんということでしょう、これはxの時刻差分がガウスノイズの下でΘの線形回帰と同じことではないか。

見方によっては面白い結果かもしれないが、確率過程で俺つえーしたかった身としては、面白くない結果となってしまった。

フィールドデータによる統計モデリングとAIC

本

統計数理研究所によるシリーズ本（ISMシリーズ:進化する統計数理）の第二巻である「フィールドデータによる統計モデリングとAIC」を読んだ。

フィールドデータによる統計モデリングとAIC (ISMシリーズ:進化する統計数理)

フィールドデータによる統計モデリングとAIC (ISMシリーズ:進化する統計数理)

「進化する統計数理」を掲げているように、本書では古典的な仮説検定型の統計学ではなく、70年代に提唱され今では非常に多くの分野で利用されているAICが、いかに有用な規準でありその実は何を意味するかを、計算機の発展によっていっそう強力になってきている現状を踏まえ、統計モデリングを実際のフィールドデータに適用させながら説明している。

1〜3章で簡単な統計モデルでAICの利用の仕方と、その有用性を示す。
4章ではAICの、パラメータ数を足す、という結果に至るまでの過程を数式のみではなく、どのような思考のもとでどのような変形がなされているかを非常に分かりやすく説明している。
5章からは実データに対し研究レベルで統計モデルを構築した中からそれぞれ変わった例をトピックとして書いたという感じである。（なので説明不足感はやや感じる。）

実データに対し、どのような思想からどのように統計モデルを構築したかを堅苦しくなく書いており、これからモデルに基づくデータ解析を行いたいという入門者にうってつけの本ではないか思う。
これまで、実験系からバイオインフォマティクスに入ってくる人にどのような本を薦めるべきか分からず有名どころで後述のPRMLを押していたが、今後はその一歩手前として本書を押すことにする。


なお、本書ではパラメータフィッティングに関してはおおよそ計算機上で最適化したという記述のみにとどまることが多かった。
また、何らかのデータに対しモデルを構築したい人にはモデルの種類やテクニック等の情報が不足しているだろう。
上記二点に関してさらなる追求をしたい人には、下記のPRMLをとりあえずお薦めする。

パターン認識と機械学習 上

パターン認識と機械学習 上

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

指数でパラメータを束縛する

機械学習

何らかのモデルに対するパラメータ最適化において、パラメータが正に限られるなどの場合に多々遭遇する。（例えばガウス分布のσ^2とか）

そのような正に制限されたパラメータを数値的にでも解析的にでもいいので最適化する場合、大抵は（私が出くわした）特に考慮せずとも勝手に正のパラメータに落ち着いていた。
ただし複数のパラメータが絡み合っていた場合、不意に負になってしまうこともあった。
そういうときはif文で負の場合のときの処理を書いて、最終的にそれっぽい解が見つかったからといってごまかしていたけど、ちゃんとした方法があったのでまとめてみる。

まとめるというほどまとめることのない方法であるが、例えば正のパラメータαを以下のように変換してβに関して最適化すれば、もとの値が負になることはない。

{ \displaystyle \alpha=\exp(\beta) }

なんてことはない簡単な話だけど、教えてもらうまで気付かなかったのは筆者に数学的センスが欠如しているからだろう。。。
ただし解析的に解く場合は式が複雑になる可能性があるので、ご注意を。


例として以下の式で表せるCox回帰モデル（比例ハザードモデル）を挙げてみる。

{ \displaystyle h(t, z_1, z_2, …, z_k)=h_0(t)\times\exp(\beta_1z_1+\beta_2z_2+…+\beta_kz_k) }

ここでβがパラメータでありこれを最適化したいのだが、expの肩に乗っているのは全体で見たとき必ず「正」であり都合がいいからである。

SRA Toolkit

バイオインフォ

自分が利用したオプションのメモ。
あくまで自分が利用した範囲にとどまるので、一般的な使い方は他をあたった方が無難だと思われる。

1. fastq-dumpを用いたsraからfastqへの変換

"--spot-group" or "-G"
バーコード付きの配列が1つのsraファイルに混在しているとき、自動的にバーコードを認識してバーコードごとに分割されたfastqを作成するオプション

"--split-files"
ペアエンドの配列が1つのsraファイルに混在しているとき、対応づけてhoge_1.fastq, hoge_2.fastqに変換するオプション

ある確率変数に対する複数のガウス分布の積

確率論

独立な確率変数に対するガウス分布の積は見かけるが、特定の確率変数に対しガウス分布が複数存在する場合の計算を発見できなかったのでまとめてみた。

動的計画法の計算過程で必要となったが、一般的な確率モデルで必要となることはなさそう？

ある確率変数に対するガウス分布の積が以下であるとき

{ \displaystyle \prod_{i=1}^NN(x|\mu_i,\sigma_i^2)=\prod_{i=1}^N\frac{1}{\sqrt{2\pi\sigma_i^2}}\times\exp\left( -\frac{1}{2}\sum_{i=1}^N \frac{1}{\sigma_i^2}(x-\mu_i)^2 \right) }

exponentialの中身を展開すると、
f:id:yrfwbioinfo:20140408183105p:plain

ここで、後半2項をまとめると、
f:id:yrfwbioinfo:20140408183129p:plain

以上から、正規化係数をひとまず無視すると、もとのガウス分布の積は以下のガウス分布の積となる。
f:id:yrfwbioinfo:20140408183145p:plain

正規化係数が不要な場合には、一つ目のガウス分布のみが重要で、右側のガウス分布の積はスケーリングに影響するのみである。
上記ガウス分布の積の正規化係数は、
f:id:yrfwbioinfo:20140408183201p:plain

よって、最後に正規化係数の帳尻を合わせるには以下の係数を掛ける必要がある。
f:id:yrfwbioinfo:20140408183220p:plain


ちょっと計算が雑多でどこかで計算間違いをしている危険性はあるが、プログラム上でN=5の場合で適当に値を設定していずれも数値計算が一致したのでおそらく正しいだろう。

memory-less propertyと指数分布のゆるくない関係

確率過程 確率論

時間を経ても次の起こりやすさが変化しないという性質は、無記憶性（memory-less property, loss of memory property）と呼ばれる。

ここである時刻までにある現象が起きる確率の累積分布関数を考えたとき、それが無記憶性を持つとは以下の式が成り立つ場合であり、これが成立するのは指数分布のみとされる。

{ \displaystyle Pr(T>t_1+t_2|T>t_1)=Pr(T>t_2) }


というのは、累積分布関数を{ \displaystyle Pr(T>t)=g(t) } とすると、上記定義から
{ \displaystyle g(T>t_1+t_2)=g(T>t_1+t_2,T>t_1)=g(T>t_1+t_2|T>t_1)g(T>t_1)=g(T>t_2)g(T>t_1) }
となる。

ついで{ \displaystyle \ln P(T>t)=h(t) } とすると、上記式から{ \displaystyle h(t_1+t_2)=h(t_1)+h(t_2) } が成立する。

ここで、微小時間{ \displaystyle dt=t_2 } を考えと、以下の式が導かれる。

{ \displaystyle h(t_1+dt)-h(t_1)=h(dt)=d(0)+dth'(0)+O(dt^2)\\ }

ここで確率の定義から{ \displaystyle Pr(T>0)=1 } すなわち{ \displaystyle g(0)=1 } であるので、{ \displaystyle h(0)=0 } となることを合わせて、

{ \displaystyle \frac{h(t_1+dt)-h(t_1)}{dt}=dth'(0)+O(dt^2)\\ }

ここで、{ \displaystyle h'(0)=-\lambda } とすると、上記式から{ \displaystyle h(t)=-\lambda t } が導かれる。

よって累積分布案数は
{ \displaystyle g(T>t)=e^{h(t)}=e^{-\lambda t}\\ }
となり

確率密度関数{ \displaystyle f(T=t) } は
{ \displaystyle f(T=t)=\frac{g(t)-g(t+dt)}{dt}=-g'(t)=\lambda e^{-\lambda t}\\ }
となる。

以上から、無記憶性を満たす分布は指数分布のみとなる。