49.9.1 大きな行列の固有値のグラフの表示
N = 20と実数gに対して行列
Q+(g) :=
1 √
1g 0 0 0 0 0 0
√1g 0 √
2g 0 0 0 0 0
0 √
2g 3 √
3g 0 0 0 0
0 0 √
3g 2 √
4g 0 0 0
0 0 0 √
2g 5 √
5g 0 0
0 0 0 0 √
5g 4 . ..
0 0 0 0 0 . .. . .. √
N g
0 0 0 0 0 √
N g N+ (−1)N
+g2I
Q−(g) :=
−1 g 0 0 0 0 0 0
g 2 √
2g 0 0 0 0 0
0 √
2g 1 √
3g 0 0 0 0
0 0 √
3g 4 √
4g 0 0 0
0 0 0 √
4g 3 √
5g 0 0
0 0 0 0 √
5g 6 . ..
0 0 0 0 0 . .. . .. √
N g
0 0 0 0 0 √
N g N−(−1)N
+g2I
を定義する。Iは単位行列。Q±(g)の固有値を小さい順からそれぞれλ±1(g), λ±2(g), λ±3(g),· · · , λ±N+1(g)と する。固有値のリスト
L1 = (λ+1(j/30))90j=−90, L2 = (λ+2(j/30))90j=−90, L3 = (λ+3(j/30))90j=−90, L4 = (λ+4(j/30))90j=−90, L5 = (λ+5(j/30))90j=−90, K1 = (λ−1(j/30))90j=−90, K2 = (λ−2(j/30))90j=−90, K3 = (λ−3(j/30))90j=−90, K4 = (λ−4(j/30))90j=−90, K5 = (λ−5(j/30))90j=−90,
を作りlist plotでそれぞれのリストをプロットせよ。グラフはオプションplotjoined=Trueにより隣り
合う点同士を線で繋ぐこと。すべてのグラフを足し合わせて一つのグラフにする(固有値ごとに色を変えてお くとよい)*19。
49.9.2 ランダム行列の固有値の分布
ランダムな数値を要素に持つ行列をランダム行列という。ランダム行列は,原子核の研究から現れたもので あるが,最近になってその数学的研究は飛躍的に発展した。また,ランダムなポテンシャルを持つシュレディ
*19注:このような行列は2準位を持つ原子とレーザー光との相互作用を表すモデルのハミルトニアンを記述しておりRabiモデルと 呼ばれる。対応する固有値はその量子系のエネルギー準位を表している。
49.9 練習問題 49 行列計算
ンガー作用素の固有値や固有ベクトルの挙動は絶縁体のメカニズムの数学的な説明を与える。
————– *** —————
問題 −1/2から1/2までのランダムな数を成分に持つ行列に持つ行列を作り,その固有値を複素平面上にプ ロットせよ。行列のサイズは800×800とし,縦横比は1(aspect_ratio=1)とする*20。
そのような行列の固有値はある半径の円の中に均一に分布する事が知られている(circular law)。わかりや すい現象だが,これを数学的に証明するのは難しい。
-5 5
-5 5
ちなみに,Sageではrandom()で(0,1)の中の数をランダムに返す:
1 for i in r a n g e(4 ) : 2 p r i n t r a n d o m () ,
実行結果の例 0.329694823321 0.831787904723 0.214856702912 0.615524329254
したがって,random()-1/2で(−1/2,1/2)の一様な乱数を得ることができる。
1 for i in r a n g e(4 ) :
2 p r i n t r a n d o m () - 1/2 ,
実行結果の例 -0.0357016806704 0.0999712876727 -0.169974222529 0.39728018763
行列のbase ringは何も指定しないか(その場合はRDF),もしくはCDFにするとよい。グラフの描画は
list plotを使うのがよい。
*20800×800のサイズの行列の計算ではエラーが起こることがあるが,そのような場合には,400×400のように行列のサイズを減ら して計算してください。
50 乱数
乱数とはさいころの目のようにでたらめ*21に表れる数字のことです。乱数は,乱雑*22な現象を記述する場 合だけでなく,決定的な量の近似値を計算するためにもしばしば用いられるます。ここではSageにおける乱 数の基本的な使い方と簡単な応用を紹介します。以下のプログラムはSage notebookで動作する事を確認し ていますが,plot等のSageの命令を除けば,純粋なPythonプログラムとしても動作します。)
Pythonで使用される乱数はメルセンヌ・ツイスター(Mersenne twister)と呼ばれるアルゴリズムによって
生成されます。これはある決まった手続きによって数を生成してゆくものですが,これは実用上十分なランダ ム性をもち非常に高速に生成する事ができるアルゴリズムです。ちなみにメルセンヌ・ツイスターは1996年 頃日本人(松本眞氏・西村拓士氏)によって開発されました。これは巨大なメルセンヌ素数219937−1が素数 であることを利用して作られるます。
コンピューターで作る乱数は,ある決定的な手続きで生成されます。したがって,それ自体にはパターンが あるので本当の意味での乱数ではありません。決定的な手続きによって作られる,乱数を近似するような数の 列を擬似乱数といいます。どのような数列なら真の乱数といえるのかというのは難しい問題*23です。ひとつ の擬似乱数は,ある種の問題を考えるときには,十分乱雑だけど,他の問題を考える場合には,乱雑さが不十 分かもしれません。ここではどのようにして良い乱数を生成するのかという問題には触れずに,乱数の使い方 と,その簡単な応用例を紹介したいと思います。
50.1 乱数の基本的な使い方
Pythonにはrandomという標準ライブラリがあり,様々なタイプの乱数を生成することができます。
まずは次のプログラムを実行してみましょう:
1 i m p o r t r a n d o m # ラ イ ブ ラ リr a n d o mを 呼 び 出 す 2 r a n d o m . s e e d (0) # 乱 数 種 を0と す る
3 p r i n t " [0 ,1)の 中 の 実 数 を ラ ン ダ ム に 生 成 す る"
4 for i in r a n g e(5 ) :
5 p r i n t r a n d o m . r a n d o m () # [0 ,1)の 実 数 を ラ ン ダ ム に プ リ ン ト す る
実行結果の例 [0,1)の中の実数をランダムに生成する
0.844421851525 0.75795440294 0.420571580831 0.258916750293 0.511274721369
random.seed(0)は乱数の種の設定ですが,これについては後で説明します。random.random()で区間[0,1) の実数をランダムに作ることができます。*24
さて,上のプログラムをもう一度実行して見ましょう。実行結果は先ほどと 全く同じ になります。また隣 の席の人とも同じになっているはずです。この数列は一見すると乱雑ですが,ある決定的な手続きで作られて います。このようにある種の数学的な手続きで作られる数列だけど,一見するとパターンがなく乱数のように 見えるものを疑似乱数といいます。
上の乱数には再現性がありますが,数値計算では乱数の再現性がない事にはあまり意味はなく,むしろ再現 性のある方が有益な場合が多くあります。たとえば,ある数値計算をしていて予想外の結果が出たときに,そ
*21出鱈目
*22乱雑の英語約はrandom(ランダム)だけど,二文字まで同じなのが不思議
*23数学的なというより思想上の問題
*24[0.0,1.0)の中にある浮動小数点数が無作為に抽出されます。
50.1 乱数の基本的な使い方 50 乱数
れが偶然なのか必然なのかの原因を知りたいと思っても,同じ状況を再現できなければ,原因を突き止める事 が難しくなるでしょう*25。
乱数の列は乱数種(random seed)ごとに決まっており,乱数種を変えることにより異なる乱数を使うこと ができます。例えば,上のプログラムでrandom.seed(1)のように乱数種0代わりに1を使えば,異なる乱 数が得られます。乱数種を変えてみましょう。
1 i m p o r t r a n d o m 2 r a n d o m . s e e d (1) 3 for i in r a n g e(5 ) :
4 p r i n t r a n d o m . r a n d o m ()
実行結果の例 0.134364244112
0.847433736937 0.763774618977 0.255069025739 0.495435087092
パスワードの生成などの時には,いつ実行しても同じ結果が出てしまうのは困るかもしれません。そういう 場合には,そのときの時刻や乱雑なキー入力,マウスの動きから乱数を作ればよいでしょう。放射性物質の崩 壊時間を利用すれば,真の乱数を作ることが出来ます*26。
プログラムを実行するたびに異なる乱数を使いたければ乱数種の指定をrandom.seed()とします。この場 合,プログラムの実行時刻から乱数種が決まります:
1 i m p o r t r a n d o m
2 r a n d o m . s e e d () # 乱 数 種 を 時 間 か ら 定 め る 3 p r i n t " [0 ,1)の 中 の 実 数 を ラ ン ダ ム に 生 成 す る"
4 for i in r a n g e(4 ) :
5 p r i n t r a n d o m . r a n d o m ()
実行結果の例 [0,1)の中の実数をランダムに生成する
0.0104079725999 0.605201486317 0.0631615144758 0.592872489434
このプログラムは毎回異なる乱数を生成します。このプログラムを何度か実行して結果が異なることを確認 してみましょう。このプログラムの実行する時刻が全く同じでない限り,同じ結果は得られません。
[0,1)の一様分布だけでなく,さまざまな乱数が用意されています。代表的なものを紹介します。
乱数の種類
• random() :[0,1)の中の実数をランダムに返す。
• randint(a,b):a≤N ≤bであるようなランダムな整数Nを返す
• choice(list):listの中からランダムに要素を取り出す。
• uniform(a,b):[a, b]の中の実数をランダムに返す。
• gauss(mu,sigma) :平均mu,標準偏差sigmaのガウス分布に従うランダムな実数を返します。
乱数生成の命令はすべてrandom.で始まり,乱数の種類に応じてその後の命令を指定します。
たとえばrandom.choice(....)はライブラリの中のchoiceという命令を呼び出すことを意味します。た
とえばリスト[’a’,’b’,’c’]からランダムに要素を取り出すにはrandom.choice([’a’,’b’,’c’])のよ うにします。上記以外にも,多くの種類の乱数が用意されていますが,必要に応じて調べて下さい。
*25幽霊の目撃情報のように
*26量子力学によれば,放射性物質の崩壊時間は真に確率的であり,それを正確に予測することは理論上不可能です
次のプログラムは1から8までの整数1,2,3,4,5,6,7,8をランダムに20個表示します:
1 i m p o r t r a n d o m 2 r a n d o m . s e e d (0)
3 for i in r a n g e( 1 0 ) : # 1か ら8ま で の 整 数 を ラ ン ダ ム に 抽 出 す る 4 p r i n t r a n d o m . r a n d i n t (1 ,8) ,
実行結果の例 7 7 4 3 5 4 7 3 4 5 8 5 3 7 5 3 8 8 7 8
次のプログラムではリストx=["red", "green", "blue", "black", "white", "yellow"]の中から要 素がランダムに抽出されます:
1 x = [" red ", " g r e e n ", " bl u e ", " b l a c k ", " w h i t e ", " y e l l o w "]
2 for i in r a n g e( 1 0 ) : # リ ス ト xの 要 素 を ラ ン ダ ム に1 0個 と り だ す
3 p r i n t r a n d o m . c h o i c e ( x ) ,
実行結果の例 green white yellow white blue red blue black yellow yellow
次はrandom.gauss(0,1)は平均0,分散1のガウス分布に従う乱数です:
1 # 平 均0, 分 散1の 正 規 分 布 に 従 う 乱 数 2 for i in r a n g e( 1 0 ) :
3 p r i n t r a n d o m . g a u s s (0 ,1) ,
実行結果の例 -1.98153307955 0.288243745581 -0.119123311085 1.80432993192
-0.160362179057 -0.0506597136487 -0.190873889601 -0.990606238348 0.673029984025 -1.32408246447
一般に平均µ,分散σ2のガウス分布は
P(x) := 1 σ√
2πe(x−µ)2/(2σ2)
で定義され,実数値をとる乱数がガウス分布に従うとは,その乱数が[a, b]に入る確率が
∫ b a
P(x)dx
である事と定義されます。この乱数が,本当にガウス分布であるかどうかは,ヒストグラムを書いてみるとよ く分かります。
1 i m p o r t r a n d o m 2 r a n d o m . s e e d (0)
3 lis = [ r a n d o m . g a u s s (0 ,1) for j in r a n g e( 2 0 0 0 0 ) ] # ガ ウ ス 分 布 に 従 う 乱 数 の リ ス ト 4 H = h i s t o g r a m ( lis , b i n s =40 , n o r m e d =True, c o l o r =’ l i n e n ’) # l i sの ヒ ス ト グ ラ ム 5 p ( x ) = 1/ s q r t (2* pi )* e ^( - x ^ 2 / 2 ) # g a u s s分 布
6 P = p l o t ( p ( x ) , ( x , -5 ,5) ) # p ( x )を プ ロ ッ ト 7 ( H + P ). s h o w ()