ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2013-04-07 尤度関数を解く

ryamada222013-04-07

[][][][][][]Solving the Likelihood Euations

  • 昨日の記事の続き
  • 論文はこちら
  • 計算機代数統計的に尤度関数を解くのはこんな感じ
    • 尤度関数は必ずしも凸ではなく、そうすると複数の極大値があり得る
    • 探索すると最大値にたどり着き損ねるので、すべての微分0の点を列挙して
    • それらを比較して最大値を得ることとする
    • それができるか、というと
    • 複素数解を許せば、得られるべき解の数も決まるし、その数のことをMakimum likelihood degree (of the model) と呼ぶ
    • そのような解の列挙はできるので、複素数空間でひとまず解く
    • そのうえで、有効な解空間内にあるものを拾いだす
  • 確率の単体¥Delta_n = ¥{(p_0,...,p_n) ¥in ¥mathbf{R}^{n+1}: p_i  ¥ge 0 ¥text{ and } ¥sum_{i=0}^n p_i = 1¥}を考える
  • (p_i)が(対数)尤度関数を構成する
  • 微分が0の条件から、連立多項式ができる(斉次な多項式でないと射影空間で代数多様体として扱えないので斉次が必要)
  • ここで連立多項式の解を複素数の空間で漏れなく列挙するために、連立多項式を複素射影空間¥mathbf{P}^nの代数多様体として取り扱うこととする(PはProjectiveのP)。こうすることで、「無限遠点」が解である場合などもそれが1個の解として扱われるので、解の個数は条件によらず常に等しくなる
  • また、このような解の個数が構成要素のdegreeの積となるという条件も満たすことから、単純に計算できることも強み
  • そんな、こんなで、尤度関数を代数計算するアルゴリズム(Algorithm 6)と尤度関数の極大を与える点を代数計算するアルゴリズム(Algorithm 8)が作れます(とここまで至るのに、証明とアルゴリズムの説明が続きますが…追えません。けど、気にせず、先へ)
  • 尤度関数が書けて、それについての最尤推定をしてもよいし
  • さらに制約を入れれば、そのうえでの複素射影空間での解の数の列挙を経て最尤推定してもよい
    • たとえば、種間の配列比較における系統樹の最尤推定では、系統樹の形ごとに制約を入れることが可能で、樹ごとにパラメタセットの最尤推定ができる

2013-01-08 SNPを代数統計する

[][][]Hardy-Weinberg Equilibriumを代数統計する

  • こちらから
  • 2アレル多型の2つのアレルの頻度をそれぞれp1,p2とする
  • HWEの下での3ディプロタイプの頻度をg1,g2,g3とする
  • g1 = p1^2,g2=2p1¥times p2,g3=p2^2の関係にある
  • g1,g2,g3には制約関数があるが、それを代数アプリケーションSingularで求めてみる
ring p = 0,(p1,p2),lp;
ring g = 0,(g1,g2,g3),lp;
setring p;
map f = g,p1^2,2*p1*p2,p2^2;
ideal i0 = 0;
setring g;
preimage(p,f,i0);
  • "ring p = 0,(p1,p2),lp;"
    • 環pを定義する。これはアレル頻度の環
  • "ring g = 0,(g1,g2,g3),lp;"
    • 環gを定義する。これはディプロタイプ頻度の環
  • "setring p;"
    • アレル頻度の世界で考えることにする
  • "map f = g,p1^2,2*p1*p2,p2^2;"
    • アレル頻度p1,p2を使ってディプロタイプgの3変数g1,g2,g3を表す関数fを定義する
  • "ideal i0=0;"
    • アレル頻度の世界(setring p;がまだ生きている)でのイデアルを定める
  • "setring g;"
    • ディプロタイプ頻度の世界に移る(変数としてg1,g2,g3を使うことにする)
  • "preimage(p,f,i0);
    • p =(p1,p2)とg=(g1,g2,g3)の関係はfで与えられているとき、ディプロタイプの変数g=(g1,g2,g3)の多項式で、その値がi0=0になるような式を作れ。
  • 出来た式-4g1g3-g2^2g1=p1^2,g2=2p1¥times p2,g3=p2^2を代入すると、-4(p1^2)(p2^2)+(2p1¥times p2)^2=0となっている

[][][][]連鎖平衡を代数統計する

  • 2アレル多型が2つある。そのアレル頻度をそれぞれp1,p2,q1,q2とする
  • 連鎖平衡の下での4ハプロタイプの頻度をh1,h2,h3,h4とする
  • Singularを使ってやってみる
ring p=0,(p1,p2,q1,q2),lp;
ring h=0,(h1,h2,h3,h4),lp;
setring p;
map f = h,p1*q1,p1*q2,p2*q1,p2*q2;
ideal i0=0;
setring h;
preimage(p,f,i0);
    • -h1¥times h4 + h2¥times h3が出ます

[][][]分割表を代数統計する

  • 2x2表をSingular でやってみよう
  • x+y=p1,z+w=p2,x+z=q1,y+w=q2とする(p1,p2,q1,q2が周辺度数)
ring p = 0,(p1,p2,q1,q2),lp;
ring x = 0,(x,y,z,w),lp;
setring p;
map f = x,p1*q1,p1*q2,p2*q1,p2*q2;
ideal i0 = 0;
setring x;
preimage(p,f,i0);
  • これは前記事の連鎖平衡とまったく同じ
  • 2x3表をやってみよう
ring p =0,(p1,p2,q1,q2,q3),lp;
ring x =0,(x1,x2,x3,y1,y2,y3),lp;
setring p;
map f = x,p1*q1,p1*q2,p1*q3,p2*q1,p2*q2,p2*q3;
ideal i0=0;
setring x;
preimage(p,f,i0);
  • -x2y3+x3y2,-x1y3+x3y1,-x1y2+x2y1が返る

2013-01-01 Singular を使ってみる 代数統計に役立つソフトウェア

[][][]Singular 代数統計ソフトウェア

  • 代数統計はデータの構造(グラフや群環体)に注目してデータを把握しようとします
  • 便利なソフトウェアもあります
  • 多項式、可換環論、計算機代数幾何。代数統計にとても便利だ、というSingularを使ってみる
  • ダウンロードはこちらのcygwinバージョンで言われるがままに入れるだけ(Windowsの場合)
    • 以下のsetup.exeをダウンロードしてダブルクリックする
All versions of Singular and all setup routines are powered by Cygwin and their setup program. Since version 3.0.0, Singular has been in the official Cygwin repository. We wish to thank the developers of setup.exe and the maintainers of the Cygwin packages.
  • cygwinのデフォルトのダウンロードだと、singularに関するものが取れてこないので、上記のダウンロード解説サイト(こちら)のうち、Mathに関するところに含まれるsingular関係のだけを選ぶ点にだけ注意する

http://www.singular.uni-kl.de/WINDOWS/cygwin_install_screen_7.JPG

  • デスクトップアイコン"Singular"ができるので、ダブルクリックして立ち上げる

f:id:ryamada22:20121231112128p:image

  • 使ってみよう!
  • ひとまず使う
    • 入力プロンプトは">"で、コマンドの終りを";"で示す
    • 整数の足し算
> 37+5;
42
>
    • 変数に整数を付値して、それを表示させてみる
> int k = 2;
> k;
2
>
    • 変数の値の真偽を調べ、それを表示させる。kは2ですか?「はい=1」、kは2ではない、ですか?「いいえ=0」
> k == 2;
1
> k != 2;
0
>
    • すでに定義した変数を使って、新たな変数を定義する。j=k+1とすればk=2,j=3となる
      • kだけ値を3に変えると…kは変わってjは変わらず
      • j=k+1をやり直せば、k=3,j=4に
> int j = k+1;
> j;
3
> k = 3;
> k;
3
> j;
3
> j = k+1;
> j;
4
>
  • 整数行列はintmat型
    • 3行4列の行列mを作る
> intmat m[3][4] = 1,2,3,4,5,6,7,8,9,10,11,12;
> m;
1,2,3,4,
5,6,7,8,
9,10,11,12
>
    • Rで言えば
m <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),3,4,byrow=TRUE)
  • ヘルプ記事を立ち上げる
> help intmat;
  • とすると"intmat"に関するヘルプが立ち上がる

[][][]Singularで環 代数統計ソフトウェア

  • ソフトウェアのインストールが済んで、簡単なことができるようになったので、環を扱ってみよう
> ring myring = 32003,(x,y,z),(dp,C); 
  • のように"myring"という名前の環を定義する
  • 以下が定義の3要素 coefficients, names_of_ring_variables,ordering
ring name = ( coefficients ), ( names_of_ring_variables ), ( ordering ); 
  • 3変数 x,yからなって標数0で、変数の並べ方の順序をlpで取り扱うような環r1を作ってみる
    • グレブナー基底の計算には斉次逆辞書式をよく使う、とか、連立方程式を解くときには「辞書式順序の斉次逆辞書式順序にはない特殊性を使う」とか、使いまわしのルールもあるようだ(参考こちら)
ring r1 = 0,(x,y),lp;
    • 標数は、ぐるっと回って「0扱いする」もの
      • 例えば、標数3とすると…
> ring rr = 3,(x,y,z),dp;
> setring rr;
> poly f = x+2y+3z+4xy+5yz+6zx+7x2y+8y2z+9z2x;
> f;
x2y-y2z+xy-yz+x-y
      • となる。3は0のことなので、3z,6zx,9z2xの項は消え(係数が0)ている。また、3で割って余りが2の項は-1をその係数として取っている
  • 変数の並べ方の順序は色々なルールが入れられるが、lpは辞書式順序で、rpというのは逆辞書式順序なので、それをr2として作ってみる
ring r2 = 0,(x,y),rp;
  • この二つの体を「環境」に設定した上で、多項式を登録してみる
>setring r1;
>poly f = (x-y)^3;
>f;
x3-3x2y+3xy2-y3
  • (x-y)^3x^3-3x^2y+3xy^2-y^3に書き換えて登録された
  • 環境で使う体をr2に変えてやってみる
>setring r2;
poly f = (x-y)^3;
f;
-y3+3y^2x-3yx2+x3
  • x,yの優先順位が入れ替わって-y^3+3y^2x-3yx^2+x^3として登録されていることがわかる
  • さて、標数0の3変数の環をlpで作って環境の環に指定しよう
>ring r = 0,(x,y,z),lp;
>setring r;
  • x^2+y^2+z^2-1=0の球を考えて、さらに(x^2+y^2+z^2-1)(z-2)=0というような球とz=2という平面の合わさった面とを考える
>poly f = x2+y2+z2-1;
>poly g = f*(z-2);
>ideal I = f,g;
>I;
I[1]=x2+y2+z2-1
I[2]=x2z-2x2+y2z-2y2+z3-2z2-z+2
>
  • このイデアルを操作して点(0,1,0)がfとgの両方の交わり?上にあることを確認するのには次のようにするそうだ
    • イデアルIについて第一変数xの値が0の部分を取り出してイデアルJにする
>ideal J = subst(I,var(1),0)
>J;
J[1]=y2+z2-1
J[2]=y2z-2y2+z3-2z2-z+2
  • f,gについてx=0を代入したものになっている
  • 引き続いて第二変数yの値が1である部分を取り出す
>J = subst(J,var(2),1)
>J;
J[1]=z2
J[2]=z3-2z2
>
  • x=0,y=1の部分が取り出されている
  • さらにz=0を指定する
>J = subst(J,var(3),0);
>J;
J[1]=0
J[2]=0
>
  • どちらも0担っているので、(0,1,0)はf,gで定まる多様体の上、ということ(らしい)
  • コマンドをいつも手打ちするのは面倒
    • テキストファイルに書いたものを読み込ませるには
> < "test.txt"
  • 逆にオブジェクトを描き出すには(output.txtファイルにオブジェクトpを書き出す)
> write("output.txt",p)
exit;