2013-04-09

Singularを使って多項式を解く

代数代数統計 Singular 多項式

Singularのインストールなどはこちら
- 多項式のままの扱いや変数の変換やらはこちら
この記事は、「値(近似値)」を求める話
Singularにテキストファイルを読ませながら実行するのがよいので、そのやり方を。
- 普通に言われたとおりにインストールすると(Windows7の場合)、"C:\cygwin\home\myusername"がワーキングディレクトリになっている(ようだ)ので、そこにテキストファイルを保存して...
- 別のPCでは"C:\cygwin\bin"直下がワーキングディレクトリらしい

> <"test.txt"

- のように読み込ませる。冒頭の">"はプロンプト。"<"test.txt" "が打ち込むコマンド
まずは、Singularを使うまでもない１変数で１個の多項式の数値解をSingularを用いて求めてみる
- $f(x) = (x-1)*(x-2)*(x-3)=0$
- "test3.txt"の中身

LIB "solve.lib";
ring r = 0,(x),dp;
poly f = (x-1)*(x-2)*(x-3);
ideal If = f;
ideal gf = groebner(If);
solve(gf,6);

- 各行の解説
  - solve()関数を使うために"solve.lib"ライブラリを読み込み
  - ringを定義。変数はx、１個だけ
  - 多項式を定義
  - 多項式fのイデアルを作る
  - イデアルIfのグレブナー基底を作る
  - グレブナー基底を６桁精度で数値的に解く
- その実行コマンド

> <"test3.txt"

- 結果

> <"test3.txt"
// メッセージがたくさん…
// ** loaded /usr/share/Singular/LIB/solve.lib (13733,2010-12-06)
[1]:
   1
[2]:
   2
[3]:
   3

// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
// is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; SOL;
//   characteristic : 0 (complex:6 digits, additional 6 digits)
//   1 parameter    : i
//   minpoly        : (i^2+1)
//   number of vars : 1
//        block   1 : ordering lp
//                  : names    x
//        block   2 : ordering C

- - 解として1,2,3が得られている
多項式を微分して、それについて解きたくなることもあるだろう
- 多項式 fをjacob()関数に渡す(そこには微分された多項式が登録されている。オブジェクトとしてはイデアル)を作り、それについてグレブナー基底をとり、それの数値解を得る

ring r = 0,(x),dp;
poly f = (x-1)*(x-2)*(x-3);
ideal If = f;
LIB "solve.lib";
ideal gf = groebner(If);
solve(gf,6);
ideal j = jacob(f);
ideal gj = groebner(j);
solve(gj,6);

- 得られた値は、1.42265,2.57735。これは $f'(x) = 3x^2-12x+11$ に代入すると $f'(x)=0$ となることは確かめられる
じゃあ、２SNPの９ジェノタイプ観測度数について、HWEを仮定して尤度関数がわかればその微分がゼロになるようなハプロタイプ頻度は次のような処理で求まる？
- 求まりません…

LIB "solve.lib";

ring r=0,(h(1..4),G(1..9)),dp;
intmat g[3][3] = 1,2,3,4,5,6,7,8,9;

poly f1= h(1)^2-G(1);
poly f2= 2*h(1)*h(2)-G(2);
poly f3= h(2)^2-G(3);
poly f4= 2*h(1)*h(3)-G(4);
poly f5= 2*(h(1)*h(4)+h(2)*h(3))-G(5);
poly f6= 2*h(2)*h(4)-G(6);
poly f7= h(3)^2-G(7);
poly f8= 2*h(3)*h(4) - G(8);
poly f9= h(4)^2-G(9);
poly H = h(1)+h(2)+h(3)+h(4)-1;
poly L= G(1)^g[1,1]*G(2)^g[1,2]*G(3)^g[1,3]*G(4)^g[2,1]*G(5)^g[2,2]*G(6)^g[2,3]*G(7)^g[3,1]*G(8)^g[3,2]*G(9)^g[3,3];
ideal iL = jacob(L);
ideal iME = iL,f1,f2,f3,f4,f5,f6,f7,f8,f9,H;
ideal Gg = groebner(iME);
solve(Gg,6)