計算情報数学 2008 (第 11 回) 種々の数学ソフトウェア (その2)

計算情報数学 2008 (第 11 回)
種々の数学ソフトウェア (その 2)
野呂 正行
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 1
Q & A : 消去法の説明が不十分です
変数集合 X = Y ∪ Z, Y ∩ Z = ∅ とする
< が (Y, Z) に関する消去順序とは
t ∈ K[X] \ K[Y ], s ∈ K[Y ] ⇒ s < t
がなりたつことをいう
Y が Z に優先する辞書式順序は消去順序
しかし, 計算時間がかかる (かかりすぎる)
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 2
ブロック順序
<Y , <Z を K[Y ], K[Z] の項順序とする (たとえば全
次数逆辞書式順序).
t = tY tZ , s = sY sZ (tY , sY ∈ K[Y ], tZ , sZ ∈ K[Z]) に
対し
s < t ⇔ sY <Y tY または (sY = tY かつ sZ <Z tZ )
と定めれば, < は消去順序
Singular
ring r=0,(y1,y2,z1,z2,z3),(dp(2),dp(3));
Asir
Ord=[[0,2],[0,3]];
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 3
消去イデアル
定理
イデアル I ⊂ K[X], Z ⊂ X, Y = X − Z, < : (Y, Z)
に関する消去順序, <Z を < の K[Z] への制限とする.
G を I の < に関するグレブナー基底とすれば,
G ∩ K[Z] は イデアル I ∩ K[Z] の <Z に関するグレ
ブナー基底となる.
意味
消去したい変数集合を Y とする.
消去順序によるグレブナー基底の元のうち, Y を含ま
ない多項式のみを集めたものが, I から Y を消去し
て得られる全ての多項式を生成する.
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 4
Risa/Asir での典型的な計算例
P を入力イデアルとし, Y を消去したい変数とする.
load("primdec_mod")$
X=append(Y,Z)$
G0=nd_gr_trace(P,X,1,1,0)$
Ord=[[0,length(Y)],[0,length(Z)]]$
G1=nd_gr_trace(G0,X,1,1,Ord)$
G=elimination(G1,Z)$
Z が 1 変数 (z とする) なら
M=minipoly(G0,X,0,z,t)$
により, z が満たす方程式が t の多項式として得ら
れる
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 5
Q & A : いろいろ
Singular の option(redSB), option(prot)
って?
redSB : 簡約グレブナー基底 (お互いに割れるだ
け割る) を計算
prot : 計算の経過を表示
Maxima で張り付けができない
張り付けメカニズムが違うソフトがしばしばあり
ます.
一旦, ファイルにセーブして, ファイルを読み込め
ば確実です.
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 6
Q & A : いろいろ
DVD 起動で USB メモリが認識されない
やっぱり, 実物を見せてもらわないとわからない
です (すみません).
多項式集合の共通零点を, なぜ拡大体上で考える
のか?
グレブナー基底計算アルゴリズム自体は任意体で
動くので, 基礎体を一般の体にしてありますが, 零
点集合については基礎体に依存する性質があるの
で, 一般に拡大体上で考える設定にしてあります.
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 7
統計解析ソフト – R
統計解析ソフトいろいろ
S (S-PLUS)
商用ソフト – 30 万以上? (学生用は安いが)
R
S と同じ言語仕様を実装した, フリーソフト
Microsoft Excel
ご存知表計算ソフトだが, 統計機能もついている
⇒ 実は危険!! (だそうです)
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 8
R の特徴
一般的な数値計算言語として使える
オブジェクト指向言語だが, ベクトル, 行列, リス
トが使える C ライクな言語と思ってよい.
行列操作が便利
逆行列, 一次方程式系求解, 固有値, 固有ベクトル,
特異値計算などがコマンド一発でできる
大規模データ処理にも便利
データが並んでいるテキストファイルを用意すれ
ば, 簡単にデータを入力できる.
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 9
R の特徴 つづき
グラフィックス
plot でグラフ, hist でヒストグラムなど, 簡単なコ
マンドで簡単に描ける
統計解析のための豊富な機能
もちろん, 確率分布, 乱数発生, 推定, 検定, モデル
の当てはめなどが, 組み込み, あるいはパッケージ
として多数提供されている
ユーザコミュニティの充実
日本にもユーザが多数. 情報も豊富 (RjpWiki
参照)
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 10
R の使い方
起動
√
KNOPPIX/Math で, x → GNU R
orange のシェルから R
コマンドライン入力が基本
プログラムファイルを作り, 読み込み, 実行
⇒ Singular, Asir と同じ
マニュアル
書籍, フリーマニュアル多数. RjpWiki, あるいは
http://aoki2.si.gunma-u.ac.jp/R
/usr/local/Math-ja/manuals に一部入れてある
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 11
R におけるデータ
基本データ (スカラー)
数値または文字列 ( "abc" etc.)
複合的データ
ベクトル, 行列, データフレーム, リスト
変数 (オブジェクト)
なんでも入る入れ物
代入は x <- 1
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 12
ベクトル
ベクトル = 同種データの一次元配列
c(...) で生成できる
x <- c(1,2,3)
y <- c(x,x,x)
後者 : ベクトル x の成分を 3 つ分並べたベクトル
x <- 1:10
(1,2,...,10) が生成される
ベクトルの要素
x[1] . . . x[10] : インデックスは 1 から始まる
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 13
行列
行列 = 同種データの二次元配列
A <- matrix(0,2,3) ⇒ 2 × 3 の零行列
A <- matrix(1:6,2,3) ⇒ A =
1 3
5
2 4
6
!
第一成分はベクトルで, 列優先 で成分の初期化に
使われる.
成分, 部分行列の取り出し
A[1,2] : (1, 2) 成分
A[1,] : 第 1 行, A[,2] : 第 2 列
A[,c(1,3)] : 第 1, 3 列からなる部分行列
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 14
行列操作
対角行列 : diag(c(1,2,3))
単位行列 : diag(3) (サイズを指定)
行列積 : A %*% B
A*B は成分どうしの積
転置 : t(A)
逆行列 : solve(A)
Ax = b の求解 : solve(A,b) (b はベクトル)
固有値, 固有ベクトル : ev<-eigen(A)
固有値 : r<-ev$values
固有ベクトル : v<-ev$vectors
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 15
繰り返し, 条件分岐
for 文
for ( i in 1:1000) { s <- s+i }
if, if-else : Asir と基本的に同じ
break : 繰り返しを抜け出す
found <- 0
for ( i in 1:n )
if ( v[i] != 0 ) {
found <- 1
break
}
found が 1 なら, i に 0 でない最も左の成分の位
置が入る
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 16
その他
関数定義
f <- function(x,y,...) { ... }
ファイルから入力 : source("filename")
ヘルプ
help(関数名) または ?関数名
less が起動される.
less : テキストファイルを見るコマンド
j, k で上下, q で終了
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 17
分布の数表としての利用
例 : 正規分布 N (µ, σ 2 ) : 平均 µ, 分散 σ 2
確率密度関数 : f (x) =
√1 e
2πσ
(x−µ)2
− 2σ2
X が N (µ, σ 2 ) に従う ⇒ Z =
X−µ
σ
は N (0, 12 ) に従う
dnorm(x,m,s)
正規分布 N (m, s2 ) の密度関数 f に対し, 値 f (x)
を返す
pnorm(x,m,s)
Rx
P (X ≤ x) = −∞ f (t)dt (CDF) を返す
rnorm(n,m,s)
N (m, s2 ) に従う乱数を n 個生成し, ベクトルとし
て返す
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 18
種々の分布
関数名の作り方は共通
xxx : 分布名
⇒ dxxx : 密度関数, pxxx : CDF, rxxx : 乱数
分布
ベータ
2項
正規
カイ自乗
ポアソン
一様
R における名前
beta
binom
norm
chisq
pois
unif
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 19
使用例 : runif によるモンテカルロ法
半径 1 の 4 分の 1 円の面積 (= π4 ) を乱数により計算:
(x, y) ∈ [0, 1] × [0, 1] を N 個 ランダム生成して,
x2 + y 2 ≤ 1 となるものが n 個 ⇒ Nn はほぼ π4 だろう
N<-10000
set.seed(123)
x<-runif(N)
set.seed(456)
y<-runif(N)
n<-0
for ( i in 1:N )
if ( x[i]ˆ2+y[i]ˆ2 <= 1 ) n<-n+1
print(4*n/N)
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 20
(疑似) 乱数について
コンピュータによる乱数発生
通常, 何らかのアルゴリズムにより発生
seed (種) が同じなら, 同じ乱数を発生
周期をもち, 計算法によっては規則性を持つ
R の乱数発生器
さまざまなアルゴリズムが選べる
デフォルト : Mersenne Twister (日本製!!)
set.seed(s)
乱数生成器を seed s で初期化する
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 21
グラフィックス機能
点のプロット
v が配列 ⇒ plot(v) で (i,v[i]) のプロット
x, y が配列 ⇒ plot(x,y) で (x[i],y[i]) の
プロット
plot(f,min,max) で, 関数のプロット
ヒストグラム
v が配列 ⇒ hist(v) でヒストグラムを描画
多数のオプション
type="l" で折れ線 etc.
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 22
例 : ランダムウォーク
n<-10000
x<-runif(n)
x[1]<-0
for ( i in 2:n ) {
if ( x[i] > 0.5 ) x[i] <- x[i-1]+1
else x[i] <- x[i-1]-1
}
plot(x,type="l")
abline(0,0,col="red")
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 23
レポート問題
いろいろな乱数生成法について調べる
検定法を知っている人は, R を使って検定して
みる
ランダムウォークのシミュレーション
0 に戻る回数の分布
+ 側にいる割合の分布
2 次元版のランダムウォーク
その他, 教科書にあるいろいろな話題のシミュ
レーション
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 24
余談 : Excel について (伝聞)
http://aoki2.si.gunma-u.ac.jp/Hanasi/excel/
Excel のあら探しが詳しく載っている
Excel を統計解析に使うな
特に, 乱数については, Excel に限らず, MS は昔
からいいかげん
財務関数もだめ
ソースが古いので直っているかも
演算誤差
2 進-10 進の演算誤差がでないような仕組みをい
れているらしいが, どうもいいかげんらしい
計算情報数学 2008 (第 11 回)種々の数学ソフトウェア (その 2) – p. 25