Rで微分方程式

微分方程式を解かず（積分せず）に＆楽に解（結果）を知りたい。大学で習うのは簡単に積分できる簡単な微分方程式ばかりですが、興味深い微分方程式の多くはそうした簡単に積分できるものばかりではありません。むしろ興味深い性質を持つのは簡単に積分できないからということも多く、そうした場合にどうするかがいつも頭の痛い問題でした。もちろんMathematicaを用いるとか、お金がないならOctaveを使うとか、微分方程式をコンピュータに解かせることはできるのでしょうが、それを学生にやらせるにはいろいろな問題があります。例えば、フリーでないので、学生の自習用に同じ環境を学生のパソコンに用意するのが難しいとか、数値計算でしか使わないソフトを習得させるためのコストを学生は払う価値があるのか、なのです。

そこで、フリーで自由にインストールして自習環境が容易に準備可能＆数値計算専用ソフトと遜色ない数値解析能力＆他の用途にも使用可能という3つの条件を満たすものとして、RとdeSolveというパッケージに着目し、3年前期の選択科目「コンピュータ数値計法」で扱うことにしました。2年後半の「統計データ解析」でRを使ってデータの統計解析（記述統計、回帰分析、母数の推定・検定など）を扱っているので、同じRの環境で数値計算を学んだほうが学生のコストを抑えられるという考えもあります。

deSolveは微分方程式を数値積分するRのパッケージで2010年に公開されたもの。もともとRにはodesolveという常微分方程式を数値積分するパッケージ(2001)があったのですが、それを徹底的に強化し、常微分方程式、拘束条件のある常微分方程式から、遅延のある常微分方程式、偏微分方程式を数値積分する環境としたものです。


このdeSolveのチュートリアルとして書かれたのが、K.Soetaert,J.Cash and F. MazziaによるSolving Differential Equations in Rというテキスト。K.SoetaertはdeSolveを開発者の一人です。



テキストは数値積分のアルゴリズムとdeSolveで用いるアルゴリズム、計算誤差の議論から始まり、さまざまな微分方程式をdeSolveで解く方法が解説されています。

例えば、左に示した二つの図は、第8章の偏微分方程式の差分解法のあとの第9章で扱っている熱伝導方程式、波動方程式、ラプラス方程式、１次元および2次元のブッルッセレータ方程式、サインゴルドン方程式、非線形シュレーディンガー方程式のうちの2次元ブルッセレーター方程式（上の図）と非線形シュレーディンガー方程式（下の図）を数値積分した結果を示したものです。

ソースコードは初期条件や格子の差分化の指定などを含めても10数行の短いもの。それでこうした結果を得られるのは素晴らしいの一言。時間が出来たら2次元ナウ゛ィエストークス方程式を数値積分して対流から乱流になる様子をローレンツ方程式と比較してみたいと思います。示したいのは波長の短い撹乱が長波長側には影響せず、バタフライ効果がないこと。つまり、ブラジルで蝶が羽ばたいてニューヨークでハリケーンが起きるなんてバカなことはないことを数値計算で検証したいだけなのですが、deSolveを使えば楽に出来そうで楽しみです。

下は2次元サインゴルドン方程式の解のスナップショット。

情報カスケードとドミノ効果

ヒトが順番に二択のクイズに回答する。自分より前に回答した人たちがどのように選択したかを教えたときどうなるのか？簡単な問題であれば、順番が前の人々がたまたま間違った選択肢を選んだとしても、後の人で正解を知っている人が正解を選び、いづれ正解を知っている人が多数派になります。そうなると正解を知らない人は多数派を選べば正解を選ぶことになり、正答率は上昇します。このような、たまたま最初に間違っても、いづれ多数派が正解を選ぶ状況を自己修正的であるといいます。この「自己」というのは、系の性質（問題が簡単なので、正解を知っている人の比率が高い）により系が自力で自動的に修正するためこう呼ばれています。また、正解を知らない人は多数派を選択ぶ傾向が強いため正答率が上昇しますが、このことを「集合知効果」と呼びます。

一方、難しい問題の場合、こうした修正機能や集合知効果がうまく働く保証はありません。正解を知っている人が少ないため、一旦多数派が不正解となってしまうと、それを修正するだけの正解を知った人が現れる前に、正解を知らない人が誤った多数派にどんどん流れ込んでくるため、修正されなくなるのです。こうした性質を「非自己修正的」と呼びます。このとき、正解を知らない人は多数派を選んでも正解とは限らないため、集合知効果はないか、あったとしても非常に弱くなります。

では、系が自己修正的か、非自己修正的かをどう判断すればいいのでしょう。自己修正的の場合、系の正答率は1/2より大きなある値に収束します。一方、非自己修正的な場合、正答率は1/2未満のある値か、1/2より大きなある値の二つの値のどちらかに収束します。つまり、ヒトが次々回答するときに正答率がどの値に収束するかで判断するわけです。実際、最初に集団実験で検証したときは、ノーヒントで自分の知識だけで回答したときの正答率が70%,80%の簡単な問題の場合、正答率は90%まで上昇して自己修正的。一方、50%,60%の難しい問題の場合、正答率は80%から90%に収束することもあれば、20%に収束することもあって集合知効果は弱く、非自己修正的になっていました。

左の図が、過去4回の集団実験の結果を示したものです。（実験の名前のカッコの中の数字は、クイズに回答した人数を表しています。）横軸はノーヒントでの正答率、縦軸は過去の回答を参考にしたときの正答率を表します。ノーヒントの正答率が70％以上なら、概ね過去を参照した回答の正答率は上昇し80％から100%になっています。一方、ノーヒントの正答率が50%、６０％の難しい問題の場合、過去を参照した回答の正答率はノーヒントの正答率から上がることもあれば下がることもあることが分かります。

これは直感的でわかりやすいのですが、一番の問題は「カッチョ悪い」こと。要するに、見たまんま。また、正答率の分散がゼロか有限か、と言っても、実験で高々60人程度が回答するだけなので、本当にゼロに収束することはありません。では、もっと統計物理学者として「カッチョいい」「正統的な」方法はないのか？それが、ドミノ効果の測定でした。

ドミノ効果とは、ドミノ倒しの様子から連想されるような、一つの現象が次々と連鎖的に別の現象を引き起こす共同的な効果のことです。二択のクイズの場合は、最初の人が間違うかどうかが後半に回答する人への影響のことです。この影響を、最初の人と、順番がｔ離れたｔ＋１番目に回答する人の選択の相関関数で計測しC(t)と書くことにします。するとC(t)/C(0)という規格化された相関関数は、最初の人が正しいか間違うかでのt+1番目の人の正答率の差を表すようになります。この規格化された相関関数の振る舞いから、系が非自己修正的かどうかが判定できます。もし、C(t)/C(0)が十分大きなｔに対して正の値をとるなら非自己修正的となります。

つまり、ドミノ効果が無限に続くなら、最初に間違うと延々と続くので非自己修正的と言えるわけです。一方、C(t)/C(0)がゼロになっても自己修正的とは言えません。最初の人の影響が無視できるほど小さくなっても、後半の人々が集団的に間違う可能性は排除できないからです。

この規格化された相関関数の振る舞い（十分大きなｔに対してC(t)/C(0)>0）が系の非自己修正性の十分条件を与えることをもとに、データ解析（ECCS'14のプロシーディングとして投稿中）を行ってみました。 上の正答率の散布図に示した実験のうち実験３(EXP3)のデータを用いて解析したものです。ノーヒントでの正答率が70％、80%の簡単な問題の場合、ｔが大きくなるとC(t)/C(0)はほぼゼロになります。一方、ノーヒントの正答率が50%,60%の難しい問題の場合、正の値にとどまっているように見えます。これを用いてｔが無限大のときのC(t)/C(0)の値を推定した結果
 が水平線で示されています。その値は大体30%前後なので、最初の人が間違えることにより、後ろの人が間違える確率が30％も下がることが分かります（難しい問題の場合）。

ドミノ効果の測定が統計力学的に「カッチョいい」「正統的」であると言っている理由は、C(t)/C(0)のｔが無限大での値が相転移での「秩序変数」になっているからです。このことを用いると、図のようにC(t)/C(0)のゆらぎが大きくても、ｔが無限大での値を推測可能になったりと、データ解析を行う上でも重要な点です（つまり、「統計物理学者の自己満足ではない」といいたいわけです）。

この解析方法を用いて他の情報カスケード実験のデータ解析も行っているのですが、ヒトの持っている情報の精度が低い場合、非自己修正的という結論になりそうです（現在進行形）。ヒトのハード（群れる）の相互作用がいかに強いかを示していると考えられますが、逆にどういう仕組みにすれば自己修正的になるのか気になります。選択肢にオッズをつけるとかやってみましたが、うまくいきませんでした。

現代芸術鑑賞

国立新美術館のオルセー美術館展にいって、イマイチ楽しめなかったので（最後の展示室の円熟機のマネの「ロシュフォールの逃亡」の海の青はよかったので、「すいつくような赤」とは違うマネの魅力を知ることができましたが）、家にもどってネットで現代芸術を検索して引っかかったのが、「現代芸術のハードコアは実は世界の宝である展」。金色の「ミニチュアのヴィーナス」（マーク・クイン）の画像を見て興味を持ち、調べてみると台湾の電子部品メーカーヤゲオコーポレーションのCEOピエール・チャンのコレクションの展覧会とのこと。8月24日までということなので、さっそく行ってきました。





国立近代美術館に行ってみると、門のところで 「ミニチュアのヴィーナス」のモデルであるケイト・モスがヨガのポーズでお出迎え。チケットの行列でうんざりしたオルセー展とは違い、人も少なく、アートを楽しめる雰囲気。


展示は、マン・レイのジュリエット、サン・ユウ（常玉）の裸体から始まり、いままで名前も知らなかった作家の作品に触れることが出来、楽しめるものでした。サン・ユウの「アヒルとボート」、ホワン・ミンチャン（黄銘昌）のスーパーリアルな森林や田園の風景画、マーク・タンジーの「サントビクトワール山」、杉本氏の海の写真、「最後の晩餐」、（史上最高額の写真家らしい）アンドレアス・グルスキーの「V＆R」、どれもこれもスゴイ。もちろん、現代美術なので見ただけでは意味不明、説明読んでもよく分からない作品もありましたが、印象派の絵を見て楽しむのとは異なる、現代の美とはどうマーケットで定義されているのかを感じることのできる展覧会でした。


この展覧会のおもしろいところは美術史的な位置づけと同時にマーケットでの位置づけも行うこと、つまり、一つ一つの作品の価格やサン・ユウやマン・レイ、村上隆の落札価格の比較などを説明したりなど、マーケットでの評価を前面に出していることです。最近、美術品の落札価格が高騰しているということをあるブログで知ったのですが、それを裏付けるものでした。いい勉強になりました。

もっとも、、一番クールだと思ったのは、ピエール・チャンの家の写真。白壁に真っ赤なサントビクトアール山、美しい！

統計力学

佐宗先生の量子統計力学に関する著書。「おわりに」では、「大胆な本を書いてしまった，と我ながら思う．」とあります。その意味は、統計力学の基礎づけを量子力学に置く点にあります。これは難しい問題で、統計力学が適用可能なのは、量子力学の方程式が計算できない（解けない）系の場合であることにあります。この場合、系の物理量の時間発展を量子力学で計算することはできません。つまり、「系が解けないからこそ統計力学は適用可能なのですが、解けない以上は、平衡状態とはどういう状態か分からず、その物理量＝平衡値も計算できず、統計力学の定式化には量子力学をまともに使うことができない」という構造になっているわけです。

このテキストでは、量子多体系の波動関数の時間変化から出発して、マクロな物理量を観測する時間Tとマクロな物理量の関係を議論して密度行列が出てくることを解説し、マクロな系の典型的な状態がマクロな系の状態のほぼすべてであることをもとにミクロカノニカル分布を導く点に特色があります。それを系が孤立している場合と、系以外の環境と接触している場合で行い、後者の例として、観測での波束の収縮も解説しています。観測時間Tは、十分短いことを前提としていて、従来の教科書によくあった「Tはほぼ無限大で、エルゴード性から長時間平均=状態での平均」といった、奇妙な定式化とは一線を画しています。

 田崎先生の統計力学が「基礎づけには異論がある」「教科書には不向き」と「参考書」のところで書かれています。個人的には、量子多体系での典型的な状態が位相空間で占める割合がシステムサイズが無限大の熱力学極限で１になること、平衡状態の性質とはその典型的な状態の性質に他ならないこと、の2点に統計力学の基礎（ミクロカノニカル分布）をおく田崎先生のテキストのほうが、観測時間Tを導入しなくてすむ（一個のミクロな状態でOKなので、Tは無限小でよい）ため簡潔で理解はしやすいと思います。実際、「観測時間Tでの平均＝エネルギーがEの状態での平均＝密度行列ではエネルギーがEの状態の成分のみ」を導いたあとは同じ論理展開なので、熱力学で統計力学を基礎づける場合は、観測時間Tの議論がないほうが楽です。ただ、統計力学は熱力学極限の場合だけ適用できればいいというものではないので、有限の系＆観測時間Tでの物理量の計算を定式化し、実験で検証することは重要な問題です。

その他、イジング模型の平均場理論での比熱や自由エネルギーの計算が丁寧に書かれている点も参考になります。アマゾンの書評では「誤植が多い＆売り物以前」などと酷評されていますが、著者の書きたかったことを理解するのには障害とはなりません。ただ、熱力学を田崎先生や佐々先生の本などでしっかり学んでから読まないと、量子論だけで統計力学はできている訳ではないので、理解不十分となる点は注意が必要です。

更新

前回の更新は8月末なので、7ヶ月ぶりのWEBページの更新です。昨年は、集団実験の実施やオッズの情報カスケード実験の論文などと平行して、ある論文をまとめるのに時間がかかりました。その論文は、０、１のデジタルの時系列において相関が非常に強く、べき則を示したり、相転移する場合の臨界的な振る舞いを相関距離とスケーリングという視点で理解するというものです。

 タイトルは、

Finite-size scaling analysis of binary stochastic processes and universality classes of information cascade phase transition

 詳しい解説はまた別の機会に行いたいと思いますが、この仕事が情報カスケードは相転移なのかどうか、Two-Peak相というself-correctingでない相が存在するのかどうかを考える際のベースとなるものです。高々数十名規模の集団実験データから、被験者数無限大の極限をどう推測するのか。しょせん推測であって、その極限である相が存在するかどうかを証明することは不可能なのですが、その証明不可能な問題に対してどこまで主張できるのか?

10/06追記：どうもレフェリーに次々と辞退されているようです。ジャーナルに投稿してからかなり経ちますが、エディターの作業経過を見てみると、レフェリーをお願いしたら、すぐ辞退か、お願いしてしばらくして催促したら辞退か、ということが繰り返されています。確かに、相関のあるRWの相転移をスケーリングで扱った論文が過去にない、などレフェリーが難しいのも事実。確率過程として解析するならあるのですが。１年以上かけた力作なので、また、実験データの解析を行う上で不可欠な論文なので、建設的でも批判的でも統計物理学者のコメントを頂きたいものですが、どうなるのか非常に心配です。

10/25追記：辞退(declined)のカスケードが止まりませんね。投稿して２ヶ月になるのに、誰も引き受けてくれない。相関のあるRWでの相関距離（時間）を定義して、過去に調べられている情報カスケードに関連するRWのスケーリング則と相転移の普遍クラスを理論的に求めただけなのに。情報カスケードの論文でdeclinedのカスケードが起きるとは。


2015/01/10追記：いまだに何にもなし。ジャーナルでは無理か？この論文で扱ったモデルを一般化した非線形ポリア壺の相関関数をスケーリングで解析する論文が完成したので、この論文はタイトルを短く「universality classes of information cascade phase transition」とし、イントロも簡潔にして別の雑誌に投稿し直すのがいいかもしれませんね。

2015/01/15追記：やっとレフェリーコメント到着。英語なので、日本人ではレフェリーを探せなかったのでしょうね。修正しないといけませんが、とにかくよかった。

2015/4/10追記：レフェリーとのやりとり2回でパス。レフェリーのdeclinedの連鎖が続いたとき（たしか7回か8回）はどうなることかと思いましたが。ジャーナルの編集者の方に感謝します。最初に投稿したPhys.Rev.Eではレフェリー１から、「お前たちが過去にやって実験の論文と全く同じじゃないか」という理由でrejectされたときは、「相転移の普遍クラスを分類する厳密な数理の論文がなぜ実験論文と同じと言われないといけないのか」とムカつくと同時に「この論文はどうなることか」と思ったものですが。次は非線形版の論文ですが、これはこれで苦戦しそうです。