古代の立方根計算


平方根の計算については、"ニュートン法"と等価なアルゴリズムがBabylonian methodという名前で呼ばれることもある(ニュートン法という名前は適切かどうか分からないけど、代替となる呼称がないので仕方ない)。このアルゴリズムが、本当に古代バビロニアで使われていたかは、疑問があるようだけど、YBC 7289という古い時代の粘土板に、2の平方根が書いてあるらしいし、二次方程式も解いてたという話なので、平方根を計算する何らかの方法は持ってたのかもしれない。以下では、そのことに留意した上で、Babylonian methodという呼称を使うことにする。

Babylonian methodを明確に記述した最古の文献の一つは、ヘロン(現在は、1世紀頃の人という説が有力っぽい)のMetricaという本っぽい。そういうわけで、Heron's methodと呼ばれることもあったそうだ。とはいえ、ヘロンのMetricaが発見されたのが19世紀後半らしいから、その呼称も20世紀に生まれたものだろう。別にヘロンは、自分で考えたとは書いてないし、ヘロンが考案者だと見なす理由はない。

もう一つの文献は、1983〜4年に中国で出土した竹簡『算数書』で、『九章算術』より古く、前漢代に成立した可能性もある。この『算数書』(検索しづらい!)に記述されている平方根の計算法は、Babylonian methodと等価と解釈できなくもない(後述)。そういう解釈が出来るなら、ヘロンのMetricaより早いかもしれない。


平方根に比べて、立方根を計算する機会は少ない。けど、歴史上の多くの算術書には、立方根の計算も記載されている。フィボナッチの"Liber Abaci"にも載ってるようだし、江戸時代初期の和算書「塵劫記」にも記述があるそうだ。多分、体積が与えられた立方体の一辺の長さや球の半径を知りたいことが時にはあったのだろう。で、立方根の近似計算は、どこまで遡れるのか。

立方根については、古代エジプトや古代バビロニアアッシリア、アケメネス朝で、何らかの計算が行われたのか、よく分からない。1930年に書かれた"The History of the Solution of the Cubic Equation"という論文には、"The Egyptians considered the solution impossible,"と書かれているけど、古い文献だし、根拠や出典もないので、本当かどうかは分からない。


大雑把には、紀元前600年前後(±200年)に、中国の諸子百家タレスを筆頭とする(ことになっている)古代ギリシャの哲学者、インドのヤージュニャヴァルキヤ(ウパニシャッド哲学のエライ人)や釈迦と六師外道(六師外道は悪口だけど、かっこいいので良しとする)など、ユーラシア各地に思想家が現れて、現代まで名前が残るようになっている。

散発的なものを別にすれば、系統的な日食記録が出現するのも、この時期と言える。中国では、『春秋』の中で、日付と共に日食記録を残していて、最古のものは、BC720年に起きた日食だと推定されている。期間が長いので、一人の人間が全てを実際に見たわけではないものの、記録数が多いので、そこそこ信用度は高いと思われる。この慣習は日本にも受け継がれたのか、『日本書紀』には、推古天皇の時代から日付とセットの日食記録がある。

プトレマイオス天文学書『アルマゲスト』には、古代の月食記録が記載されていて、新しいものは、プトレマイオスの時代に観測されたものだが、古いバビロニアの記録を引用してるらしく、その中で最古の記録は、BC721年に起こったものだそうだ。バビロニアに於ける系統的な天文記録は、バビロン第8王朝の王Nabonassar(在位BC747~BC734)の時代に始まるそうで、プトレマイオスも、この王の即位の年を暦計算の起点に使っている。

古代インドでは、日付の記載がない日食記録はrig vedaにあるらしいけど、紀元前の体系的な日食記録などはないっぽい。

現代まで継承されている"数学書"は、それから更に数百年以上後のものになる。古代エジプトの数学パピルスなどは遥かに古いけど、これらは19世紀以降に発掘されたものなので、それ自体が継承されていたとは言えないだろう。そのへんの時代から西暦600年くらいまでを対象として適当に文献を眺める。



ヘレニズム時代の地中海・オリエント地域

ヘレニズム時代は、一般的に、BC323年〜BC30年を指すそうで、この時代の初期に、アレクサンドリアが建設されて、学術都市として機能するようになったと言われている。古代ギリシャ語文献は、Loeb Classical Libraryとして、ハーバード大学出版局から英訳と共に出版されているらしい。残念ながら、数学関連文献は扱われてないものも多いので、このシリーズ以外から探さないといけない。

古代ギリシャ語圏の人々は、どういうわけか、立方根の計算を、作図することで解きたいと思ったらしい。この問題は、定規とコンパスのみでは原理的に解けないけど、古代ギリシャ人が、"定規とコンパスのみ"という制約を課したかは知らない。

彼らは、むしろ、立方根の計算を"作図"によって行うための数学的器具を開発するという謎の方向に進んだたらしい。19世紀〜20世紀初頭の欧米では、graphical method(図式解法)という数値計算法の一種が、それなりに研究されてたっぽいのだけど、発想としては、同じようなものかもしれない。graphical methodは、現在では、淘汰されてしまったので、数学史・科学史・技術史のいずれでも取り上げられることは少ない。

アスカロンのエウトキオス(Eutocius, CE480頃〜CE530頃の人とされる)が、アルキメデスの著書『球と円柱について』のcommentaryを残していて、そこで多くの方法を記載している。


この注釈書の(古典)ギリシャ語とラテン語の対訳は、1881年に出版された
Archimedis Opera omnia, 第3巻
https://archive.org/details/archimedisopera05eutogoog/page/58/mode/2up
https://books.google.co.jp/books/about/Archimedis_Opera_omnia.html?id=Va02AAAAMAAJ

に含まれている。英語訳は以下の本にある。

The Works of Archimedes: Volume 1, The Two Books On the Sphere and the Cylinder: Translation and Commentary (English Edition)
https://www.amazon.co.jp/dp/B001CJRQ28

この本には、誰が考えた方法かも、それぞれ書いてあるけど、古い話なので、エウトキオスの書いたもの以外、証拠がないものも多いと思われる。また、考案者が正しかったとしても、この本に記載されている各原理に基づいて、実際に何らかの器具が作成されたのかは分からない。全ての方法が実用的というわけでもないと思われる。ただ、この方法の中には、ルネサンス期に紹介されて、実際に、数学器具として実現したものもあるっぽい。


デューラーの 「幾何学世界」 について
https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/80780

によると、デューラーのUnderweysung der messung(1525,1534)に、立方根を"作図"する方法が3つ書かれているそうだ。

塑像などで小縮尺の模型を作成してから本体の構築をするが, その際, 鋳造材料の量の評価が必要であり, 倍積立方体の手法はこのときに不可欠な知識になる. しかし, 学者は秘儀として来たので, ここで, 始めて職人の言葉, つまり, ドイツ語で公開するとデューラーは強調している

とある。このような方法は、計算尺と同様、それほど高い精度は期待できない。

例えば、計測によって、2の平方根を決めようと思って、1辺が1メートルの正方形を書いて、対角線の長さを測ったとする。現代で最も普及している1ミリ刻みの定規を使って首尾よく行けば、1414〜1415ミリという結論が得られるはず。しかし、これは、古代バビロニア初期に知られていた2の平方根の近似値より精度が低い。これは理想的な場合の話で、例えば、直交すべき2辺の角度が90度から1度ずれるだけで、対角線の長さは容易に10ミリ以上変動する。そう考えると、相応の工夫をしない限り、1410〜1420ミリという結論が得られれば、上出来と思われる。

立方根の近似値でも同じような問題が起きるだろう。"作図"による近似解法は、当時の水準でも天文学で使うような計算には全く向いてなかっただろうけど、当時の"工業"用途で使う分には、手っ取り早い方法と思われてたのかもしれない。結局の所、加工のバラツキより細かい精度の近似値は意味がない。

また、例えば、金属の熱膨張率は温度1度につき、10万分の1のオーダーで、日常的な温度変化が数十度の幅であるということを考えると、寸法が有効数字3桁程度で決まってれば十分ということは多かったと予想される。円周率も3.14と3桁までの数値で習うし、3桁程度で十分なことは多いのかもしれない。それでも、計算した方が早いのではと思わなくもないけど。


3つの方法は、Sporus(多分、ニカイアのSporusと呼ばれるCE3世紀頃の人)、プラトン、(アレクサンドリアの)ヘロンによるとされている。これらは、全てエウトキオスの本に載っている。デューラーは、"プラトンの方法"を利用した器具を考案しており、実際に使われたのは、この方法だったのだろうと書いてるが、多分、該当箇所は以下の部分。

Albrecht Dürer's Unterweisung der Messung
https://archive.org/details/albrechtdrersun01peltgoog/page/n164/mode/2up

もう一つ比較的有名な器具として、1558年出版のGioseffo Zarlinoという音楽家の著書"Le istitutioni harmoniche"には、エラトステネスが考案したとされるmesolabioなる器具が記載されている。

Zarlino "Le istitutioni harmoniche",1558
https://archive.org/details/imslp-istitutioni-harmoniche-zarlino-gioseffo/page/n127/mode/2up

エラトステネスが、mesolabioについて書いた著書が残ってるわけでなく、出典は、エウトキオスの本。エラトステネスがプトレマイオス3世に宛てた手紙というのが収録されているが、これには偽造説もある。同じ"古代の人"といっても、二人が生きた時代は700年くらい離れている。


1885年の論文
The Ancient Methods for the Duplication of the Cube
https://doi.org/10.1017/S0013091500037615
には、古代に考案されたとされる10以上の"作図法"が載っている。全部が、エウトキオスの本に載ってるのかはチェックしてない。



"エラトステネスの手紙"には、自分の発見が、カタパルトや投石機(καταπαλτικά και λιθοβολα όργανα)の強化に役立つと書いてある。この記述は、明らかに軍事利用を想定しているので、古代ギリシャでの軍事技術の位置付けを確認する必要がある。

ヘレニズム以前のギリシャの学者は、(現代まで名前が伝わっている範囲では)思想家や哲学者の類が多く、道具を作ったり技術的なことを評価したようには見えない。プラトンもヘレニズム以前の人なので、立方根作図器具を作ったことがあるかは疑問がある。尤も、そういう印象は、もっと後の時代に形成されたものだという可能性もないわけではないけど。

ヘレニズム期になると、機械や器具を作るような技術的問題の解決で名前を残す学者が出てくる。一番有名なのはアルキメデスだろうけど、他にも名前が残ってる人はいる。この分野でアルキメデスの次くらいには有名と思われるビザンティウムフィロンも、アルキメデスと同時代の人で、フィロンは、Mechanike syntaxisという著作を書き、以下の9つの章からなるらしい(全部が残ってるわけではないそうだ)

(1)Isagoge: ギリシャ語では、エイサゴーゲー(εἰσαγωγή)で"introduction"を意味する
(2)Mochlica: "梃子"(Levers)のことらしい
(3)Limenopoeica: "poeica"は、makingのような意味(?)らしく、港の建設法を指すらしい。
(4)Belopoeica: 同じく"poeica"が付いていて、直訳では、"arrow-making"という感じっぽい。バリスタやカタパルトを指すのかもしれない
(5)Pneumatica: pneumaticsは、現在、空気圧工学と訳されるが、かつては、気体に限らず(水車やポンプのような)流体機械全般を扱う分野だったっぽい
(6)Automatopoeica: 直訳すれば"automaton-making"だと思う
(7)Parasceuastica: "準備"とかを意味するっぽい(?)。攻城の準備らしい
(8)Poliorcetica: 英語のpoliorceticsと同根の語と思われ"攻城法"に関する章っぽい
(9)Peri Epistolon: 直訳すると"about letters"とかっぽい。機密文書に関する話らしく、暗号化とかの話(?)

このリストを見るに、フィロンは軍事技術者だったと言って良さそうに思える。他に、アレクサンドリアのムセイオンの初代館長とされるクテシビオス(Ctesibius,Ktesibios)は、発明家/技術者の類だったと見なされている。クテシビオスは著作が伝わっておらず、業績は分からないことが多いとされている。

アルキメデス、クテシビオス、ビザンティウムフィロンの3人は、ウィトルウイウスのDe Architecturaで、機械に関する本を書いた人として名前が挙げられている。
De architectura/Liber VII [14]
https://la.wikisource.org/wiki/De_architectura/Liber_VII

全部で12人の名前が確認できるが、ヘレニズム時代より前の人らしき名前も見られる(素性の分からない人もいる)。比較的有名なのは、アルキタスで、プラトンの友人だったらしい。アルキタスは立方体倍積問題でも名前がよく出てくるが、機械に関する本の題名や内容は分からない。偽アリストテレスの『機械学』という本(この本には、梃子の原理が述べられている)が知られていて、この本の真の著者は、アルキタスだと主張している人もいる。
The Mechanical Problems in the Corpus of Aristotle
https://digitalcommons.unl.edu/classicsfacpub/68/

Polydius(おそらく、現代では出身地付きでPolyidus of Thessalyと表記される人)とDiades(おそらくDiades of Pella)という人も名前が出ているが、この二人は、マケドニアの軍事技術者だったようだ。Diadesは、Polyidusの弟子で、アレクサンダー大王の遠征にも従ったとか。彼らも、著作などは見つかっていない。ヘレニズム時代の少し前から、ヘレニズム時代にかけて、ギリシャ地域は兵器開発に熱心だったようで、色んな"新兵器"(Helepolisとかpolybolosとか)の名前が残っている。

1946年の古い論文
The Origin of Greek and Roman Artillery
https://www.jstor.org/stable/3291885
によれば、シケリアのディオドロス(BC1世紀の歴史家)は、カタパルトはシラクサの潜主ディオニソス1世(BC400年前後)が集めた技術者によって、シラクサで発明されたと書いてるそうだが、紀元前8世紀頃の中東に起源があるだろうと考察している。どっちが正しいにしろ、多分、アルキメデスよりずっと以前の時代から、シラクサは軍事技術の開発に力を入れていたんだろう。


エラトステネスは、ヘレニズム時代の人で、アルキメデスとは知人だったらしい。エラトステネスの著作には、詩や"伝記"、"年代記"の類と思われるものが結構含まれている(現存せず、タイトルが伝わるのみのものが多いっぽい)。アレクサンドリア図書館(アレクサンドリア図書館はムセイオンの付属機関?)3代目館長だったとされるけど、当時は、そうした人文系の業績を評価されたのかもしれない(初期アレクサンドリア図書館館長は、人文系の学者が多いように見える)。現代に伝わるエラトステネスの業績は、エラトステネスの篩や、地球の大きさの計測などで、アーミラリ天球儀(中国の渾天儀と同一?)の発明者ともされている。

(現代まで伝わってる範囲では)エラトステネスの書いたものに技術書的なものはないけど、こうしたことを踏まえると、エラトネテネスが、軍事技術への応用を念頭に、何らかの立方根作図器具を作っていたとしても違和感はない。



エウトキオスが書いてる方法を全部チェックするのは怠い。これらの成果の大部分は、何の役にも立たなかったし、実用化されたものも含めて、全ての結果は存在しなくても特に問題はなかったと思う。けど、残念ながら存在してしまったので、mesolabioだけ見る。

立方根の作図問題は、2つの長さa,bを持った線分が与えられた時に、a^{1/3}b^{2/3}及びa^{2/3}b^{1/3}の長さの線分を作るという形で定式化できる(両方、長さの次元を持つ)。当時は、x=a^{2/3}b^{1/3},y=a^{1/3}b^{2/3}は、 a : x = x : y = y : bを満たす量と捉えられていた。

これは、ユークリッドの『原論』8巻で述べられてる考え方でもある。そして、エウトキオスが列挙した方法のどれでも同じ方針を採っている。このように問題を定式化したのは、"エラトステネスの手紙"には、キオスのヒポクラテスという人物だと述べられている。キオスのヒポクラテスは、医者のヒポクラテスとは別人で、ユークリッドより前の時代の幾何学者。

mesolabioの作り方(の一例)は、以下の図に従って、

f:id:m-a-o:20220217212236p:plain

・上下の枠線は平行で距離はa
・3つの赤い線分は互いに平行で、左右に平行移動することだけできる
・青い線分OPは、Oを中心にして平面上を回転することだけできる(Oは左の赤い線分と上の枠線の交点)
という条件を満たすようになっている。

そして、赤い線分と青い線分を動かして、試行錯誤によって、以下の条件を満たす配置を見つける。但し、A,Bは、(真ん中と右側の)赤い線分と青い線分の交点。
・AD,BE,CFは上下の枠線と直交する
・線分CFの長さはb

こういう配置が見つかると、線分ADとBEの長さがxとyになる。”エラトステネスの手紙"は、こういう配置を発見すればいいということだけ書いてるのであって、mesolabioの詳細は、色々な作り方があると思う。

正直、めんどくさいとしか思わないけど、"エラトステネスの手紙"には、自分の方法は実用的だと書いている。古代ギリシャでは、幾何学は厳密な方法だと思われてたかもしれないけど、そこから得られた方法は、少なくとも、精密さは欠片もなかった。

古代ローマ時代のオリエント地域

別に時代区分に意味はないのだけど、ヘレニズム時代は、紀元前30年までということになっている。その後、アレクサンドリアは、古代ローマ帝国の一都市になったけど、数世紀の間は、学術都市としての機能を果たしていたっぽい。

エウトキオスは5〜6世紀の人だし、この時代になっても、"立方根作図器具"を探求する人はいたようだが、立方根の近似計算法を書き残している人に、アレクサンドリアのヘロンがいる。探した限り、このようなアプローチを採る人は例外的だったように見えるけど、平方根の近似計算は、ヘレニズム時代にも見られたし、立方根の近似計算は後世に残らなかっただけかもしれない。


アレクサンドリアのヘロンは、数学ではヘロンの公式に名前が残っていて、著作に記述されたアイオロスの球は、世界最初の蒸気機関と見なされることもある。ヘロンの生没年は不詳っぽいけど、現在は、CE1世紀頃の人と見なされているっぽい。

ヘロンの著作(であることが確実)とされているものは、以下の通り。
・Pneumatica: 流体機械を扱ったっぽい
・Automatopoietica: "自動機械"の作り方
・Mechanica: タイトルは、mechanicsと同根の語だけど、現在のmechanicsは"力学"を指すが、以前は、機械学に近い分野だった。滑車や梃子のような比較的単純な機械が対象となってるっぽい
・Catoptrica: "反射光学"と訳される
・Metrica: タイトルは、metricsと同根の語で、測量に関する本らしい
・Dioptra: dioptraは天体観測や測量で使用する計測機器で、BC3世紀から使われていたらしい
・Belopoeica: "catapult-making","arrow-making","artillery"などと訳されてる

機械に関する本が多いけど、ヘロン自身が、機械を作成したり発明したことがあったのかは分からない。ビザンティウムフィロンに比べると、軍事色は、それほどない。Belopoeicaは、オープンアクセスになっているのは、Wikisourceにあるものしか見つけられなかった。
Βελοποιϊκά
https://el.wikisource.org/wiki/%CE%92%CE%B5%CE%BB%CE%BF%CF%80%CE%BF%CE%B9%CF%8A%CE%BA%CE%AC
ギリシャ語しかないし、これで全文なのかも、分からない。一応、E.W. Marsdenという人が1971年に英訳を出版しているらしいのだけど、入手自体が困難。


ヘロンのMetricaは、19世紀のオスマン帝国で、12世紀の写本が保存されてるのが見つかったそうだ。1903年に出版されたヘロン全集(?)の3巻前半に収録されている(後半は、Dioptra)。

Heronis Alexandrini Opera quae supersunt omnia vol III
https://archive.org/details/heronisalexandri03hero/mode/2up
https://books.google.co.jp/books?id=09c2AAAAMAAJ&hl=ja&source=gbs_book_other_versions

この本は、ギリシャ語とドイツ語の対訳になっている。英語訳は、おそらく存在しない。私は、今の所、ギリシャ語は数字くらいしか読めないので、ドイツ語訳を信じるしかない。

立方根の計算は、MetricaのIII巻(?)XX節にある。
https://archive.org/details/heronisalexandri03hero/page/178/mode/2up

ギリシャ語側の5〜16行目と、ドイツ語側の12〜20行目が、立方根の計算を説明している部分。ギリシャ数字の部分を手掛かりにすれば、ギリシャ語とドイツ語で、対応してる箇所は大体分かる。

Λαβὲ τὸν ἔγγιστα κύβον τοῦ ρ τόν τε ὑπερβάλλοντα καὶ τὸν ἐλλείποντα· ἔστι δὲ ὁ ρκε καὶ ὁ ξδ. καὶ ὅσα μὲν ὑπερβάλλει, μονάδες κε, ὅσα δὲ ἐλλείπει, μονάδες λϚ. καὶ ποίησον τὰ ε ἐπὶ τὰ λϚ· γίγνεται ρπ· καὶ τὰ ρ· γίγνεται σπ. ‹καὶ παράβαλε τὰ ρπ παρὰ τὰ σπ·› γίγνεται θ/ιδ΄. πρόσβαλε τῇ [κατὰ] τοῦ ἐλάσσονος κύβου πλευρᾷ, τουτέστι τῷ δ· γίγνεται μονάδες δ καὶ θ/ιδ΄. τοσούτων ἔσται ἡ τῶν ρ μονάδων κυβικὴ πλευρὰ ὡς ἔγγιστα.

ヘロンは、100(ギリシャ数字ではρ)の立方根を例にとって説明していて、一文目は、100に近い"立方数"を(κύβον)大きい側と小さい側で取ると言っている。二文目の意味は、"それは125(ρκε)と64(ξδ)である"。それ以後は、ドイツ語訳に書いてある式を、ひたすら言葉で書いてるだけ。答えは、"δ καὶ θ/ιδ΄"(4 and 9/14)となっている。単語から推測するに、100の立方根は、"体積100の立方体の一辺の長さ"という形で表現されてるっぽい。

ヘロンの近似計算法を一般化すると、Sに対して、a^3 \le S \lt b^3となるa,bに対して、
S^{1/3} \approx a + \dfrac{b(S-a^3)}{b(S-a^3) + a(b^3-S)} (b-a) = \dfrac{b^2(S-a^3)+a^2(b^3-S)}{b(S-a^3)+a(b^3-S)} = \dfrac{(a+b)S + a^2 b^2}{S + a^2b + ab^2}
という近似値を得る方法と解釈できる。例えば、S=2に対して、a=1,b=2とすると、2の近似立方根として、5/4を得る。真の値との差は1/100より小さい("5/4+1/100"の3乗は、暗算で容易に計算でき、2より少し大きい)ので、それほど悪くない。

ヘロンの記述からは、反復計算を行うという含意は読み取れない。近似値を得るのに、上限と下限の2つの値が必要となるが、得られる近似値は一つなので、そのままでは反復できない。

反復計算への拡張法は色々考えうるが、最も安直に済ますには、近似値xに対して、もう一つの近似値をS/x^2に取るということが考えられる。a=x,b=S/x^2として、上記近似式を計算すると
S^{1/3} \approx x \dfrac{x^3 + 2S}{2x^3 + S}
が得られる。これは、"ニュートン法"とは少し異なる反復計算を与える。
 x \dfrac{x^3 + 2S}{2x^3 + S} = x + x \dfrac{S - x^3}{2x^3 + S}
と書き直せば、違いがはっきりする。S \approx x^3とすれば、\dfrac{x}{2x^3 + S} \approx \dfrac{1}{3x^2}なので、ニュートン法と同じ式が得られる。

ニュートン法が、一次までの近似であるのに対して、ヘロンの式は、二次までの近似になっていて、反復した場合、ニュートン法より収束が早い。この方法は、Halley's methodと呼ばれる方法と等価になっている。このHalleyは、ハレー彗星の人で、ニュートンの知人でもあったことは有名。Halley自身は、一般的なアルゴリズムを与えたわけではなく、立方根の計算アルゴリズムを考えただけらしい。

Haley's Methods for Solving Equations
https://doi.org/10.2307/2303467

Halleyは、有理反復式と非有理反復式の2つを示していて、実際に使ったのは後者らしいけど、231の立方根を有効数字18桁まで小一時間で計算できると書いてるそうだ。初期値を6にして、231の立方根を有効数字18桁まで計算するのに、Halleyの有理反復式だと反復回数は3回、ニュートン法だと4回なので、今となっては大差ない気もする。


Metricaには平方根の計算も書いてあって、立方根と違って、反復計算をすれば近似精度を改善していけることが明示的に書かれている。ヘロンが書いてるのはBabylonian methodと等価な方法で、MetricaのI巻(?)VIII節にある。
https://archive.org/details/heronisalexandri03hero/page/18/mode/2up
これは、三角形の三辺の長さから面積を求めるヘロンの公式を説明している箇所。例として、3辺の長さが7,8,9の三角形の面積が720の平方根であることを述べて、その近似計算法を説明している。現代日本の中学/高校数学では、\sqrt{720}を答えにして終わりだけど、昔は、これでは十分ではなかった。

ヘロンの書いてる方法を一般化すれば、任意の正数Sに対して、S \approx x^2であるxを持ってきて
\sqrt{S} \approx \dfrac{1}{2} \left(x + \dfrac{S}{x} \right) = \dfrac{S + x^2}{2x}
ということになる(720の平方根の場合、x=27を使っている)。新たに得られた近似値を使って、同じ計算をすれば、もっと良い近似値を得られると述べている。

Metricaには、ヘロンの公式に対する証明はあるけど、平方根の計算については、何故うまくいくのか説明はない。平方根に収束することの証明がないのは仕方ないとして、ヘロンが、反復計算によって、近似精度の改善する理由を説明できたかも分からない。

平方根の計算では反復による改善が指摘され、立方根で、それがない理由は分からないけど、単に、反復できる形にする方法を思いつかなかったのかもしれない。


ヘロンより後の時代に、アレクサンドリアのテオンという人が平方根の近似計算を説明しているが、その方法は、ヘロンのものとは異なっている。アレクサンドリアのテオン(335?~405?)は、アレクサンドリア図書館最後の館長とされる人で、アルマゲストの注釈を書いて、その中で、\sqrt{4500}平方根の計算法を説明しているそうだ(テオンの著書の原文が、どこで読めるのかは知らない)

プトレマイオス(83?~168?)は、天文学アルマゲストを書いて、その中で、文献上確認できる最古の"三角関数表"を載せた。Johan Ludvig Heiberg(1854〜1928)という人が出版したらしいギリシャ語板のアルマゲストが以下にあり、
https://www.wilbourhall.org/pdfs/HeibergAlmagestComplete.pdf

の48ページから、表が掲載されている。περιφερειών(periferion)と書いてある第一列が、角度\thetaを表し、ευθειών(eutheion)と書いてある2,3,4列目が、120 \sin(\theta/2)を60進表記で表したものっぽい。

アルマゲストの写本は、一杯ありそうだけど、例えば、以下には、16世紀前半のものとされるギリシャ語写本がある。不完全らしいけど、"三角関数表"は確認できる。
Ptolemy's Almagest (Queens' College MS 32)
https://cudl.lib.cam.ac.uk/view/MS-QUEENS-00032/34


この表の値を決定するのに、平方根を計算する必要がある。例えば、Heiberg版51ページの3行目は、36度に対する値で、37;04,55(=36+4/60+55/3600)と書いてある。この数値は、30(\sqrt{5}-1)=\sqrt{4500}-30の近似値になっている。実際
\sqrt{4500} \approx 67 + \dfrac{4}{60} + \dfrac{55}{3600} + \dfrac{2}{21600}
は、容易に計算できる。

これを計算するためのテオンの方法は、60進法ベースで開平方をやっているのと同じ手続きになっている。テオンは、プトレマイオスよりは大分後の人なので、プトレマイオス三角関数表が、この方法で作られたかは分からない。Wikipediaに書いてある推定生没年では、ヘロンが死んだ10数年後に、プトレマイオスが誕生している。


古代中国

中国の古い算術書に「算経十書」というのがあって、周髀算経,九章算術,海嶋算経,孫子算経,五曹算経,夏候陽算経,張邱建算経,五経算術,緝古算経,数術記遺を指すらしい。清の時代に、孔继涵(1739〜1783)という人が、この10冊に、戴震(1729〜1777)の『勾股割圓記』という本を附録につけて、『算経十書』を出版したらしいから、一応、全ての本は、現代でも参照可能ということになる。一番新しいのが、緝古算経で7世紀の成立とされ、三次方程式の解法が記載されているそうだ。

算経十書(新日本古典籍総合データベース)
https://kotenseki.nijl.ac.jp/biblio/100257118


日本では、CE700年頃の大宝律令によって大学寮が設立され、養老律令令集解算経条で、算道の教科書は、『孫子』、『五曹』、『九章』、『海島』、『六章』、『綴術』、『三開重差』、『周髀』、『九司』の9冊と定められているらしい。「算経十書」に含まれない4冊は、現存しないようで、詳細は分からない。『綴術』は祖沖之(429〜500)の著書とされ、唐の時代には、「算経十書」の一冊とされていたのが、宋代に失われたため、代わりに、数術記遺を「算経十書」に加えたらしい。残りの3冊は、中国以外(朝鮮や日本)で編纂されたと推測してる人がいる。
大学寮算科の教科書
https://www.lab.twcu.ac.jp/~osada/math_textbooks_at_NU.pdf

古代中国の算術書の著者が何者だったのか分からないことが多いけど、祖沖之は、技術者でもあったとされる。天文や測量だけでなく、機械設計に取り組む"数学者"というのも、洋の東西問わず、存在したのでないかと思われる。


科学の名著〈2〉 中国天文学・数学集
https://www.amazon.co.jp/dp/B01MFGUY31

には、九章算術,海嶋算経、周髀算経の訳が含まれている。電子化されてて偉いけど、海嶋算経は、九章算術の附録という扱いで、目次には載ってない。『周髀算経』は、算術書というより、天文学書という趣が強いけど、題名に「算」の字が入ってるし、北宋代の漢籍目録『崇文総目』でも、算術類(31書ある)に分類されている。『開元占經』なども数学的内容を含むけど、こっちは、天文占書類に分類されている。

《崇文總目》
https://ctext.org/wiki.pl?if=en&res=285530



少し前までは、『周髀算経』と『九章算術』が知られている限り中国最古の算術書で、7世紀以前の算術書で現存するのは、算経十書のみだったようである。

1983~84年にかけて発掘された竹簡の中に『算数書』(<<筭數書>>)と表題される数学書が含まれていたそうである。発掘箇所が前漢時代の墓と推定されていることや、内容的にも、『九章算術』より前の算術書と考えられている。また、2006年に、雲夢睡虎地漢簡『算術』が発掘され、2007年には、岳麓書院蔵秦簡『数』という竹簡が見つかり、2010年に、北京大学蔵秦簡『算書』が北京大学に寄贈されて、古い時代の算術書が、見つかってるそうだ。他に、2008年に清華大学に寄贈された『清華大学蔵戦国竹簡』は、戦国時代中期以降の楚の竹簡とされ、『算表』と仮に名付けられた乗算表が含まれてるそうだ。2013年出版の清華大学蔵戦国竹簡(肆)に収録されているらしい。

これらの中で、『算数書』は、発見が少し前なので、一番情報が多い。偽造でないことを確認する術はないけど、そこは信用することにする。中国では、何冊かの書籍が出版されているようだけど、日本語では、以下のような解説がある。

『算数書』日本語訳
http://www.osaka-kyoiku.ac.jp/~jochi/jochi2001.pdf

張家山漢簡『算数書』訳注稿
http://pal.las.osaka-sandai.ac.jp/~suanshu/SSS/j/publications.html

訳注稿は、順番がバラバラで分かりづらいけど、手っ取り早く、原文を全部見るには、以下で公開されている英訳本の末尾を見るのがいいと思う(PDF123ページ以降)
The Suan shu shu 筭數書 'Writings on Reckoning'
https://www.nri.org.uk/suanshushu.html


『算数書』には、立方根の計算は含まれてないけど、平方根の計算が書かれていて、後の時代のものとは異なっている。訳注稿では(5)の最後にある。該当箇所の文言は、以下のようになっている。原文には句読点はないけど、読みやすさのために補われている。

方田
田一畝方以幾何步?曰、方十五步卅一分步十五。
术曰、方十五步不足十五步、方十六步有徐十六步。曰、并贏、不足以為法、不足子乘贏母、贏子乘不足母、并以為實。復之、如啟廣之术。

古代中国語(漢文)の文法書というものは殆ど存在してないけど、日本語を読める人なら、漢字の雰囲気で、大まかな意味は察することができると思う。一畝は、面積の単位で、時代によってサイズが変わるが、ここでは240平方歩に相当する。"方田"の一辺を15歩とすると、15平方歩足りず、一辺を16歩とすると、16平方歩余るという記述と辻褄が合う。

『算数書』は、他の古代中国の算術書の多くと同じく、例題と解法が列挙された問題集のような形式になっているけど、単に、数式を表現する一般的な方法がなかったためだろう。"方田の術"を一般的に解釈すると、正数Sに対して、 a^2 \le S \lt b^2となるa,bを見つけて
\sqrt{S} \approx \dfrac{(S-a^2)b + (b^2-S)a}{b^2 - a^2}
と近似しろと言ってるように読める。

元のSを、この新しい近似値で割ると、もう一つの新しい近似値が得られる。"復之、如啟廣之术"が言ってるのは、そういうことらしい。

"啟廣之术"は、『算数書』の別の箇所に、以下があるのを指している。

啟廣
田從卅步為啟廣幾何而為田一畝。曰啟八步。朮曰以卅步為法以二百卌步為實。啟從亦如此

縦が30歩で面積が一畝の(長方形)田のもう一辺の長さを問う問題で、240を30で割れと言ってる。ついでに、この記述からも、一畝が240平方歩だと分かる。


"方田の術"に従うと、a,bより良い上限と下限が得られるけど、この記述が反復計算まで示唆していると言えるのかは分からない。研究者たちの解釈では、これは"盈不足術"(古代中国の算術で、一種の線形補間法を指す方法と言っていいと思う)の一種で、反復計算せずに、終わるものだと見なされているっぽい。ただ、より良い上限と下限の両方が得られると示唆しているのは、単なる線形補間とも言えないように見える。

そして、反復計算を行えば、この方法は、Babylonian methodと等価になり、収束の早さも同じである。実際、a=x,b=\dfrac{S}{x}と置いて、ちょっと計算すると
S \dfrac{b^2-a^2}{(S-a^2)b+(b^2-S)a} = \dfrac{S+x^2}{2x} = x + \dfrac{S-x^2}{2x}
となる。

『算数書』の近似平方根と、ヘロンのMetricaにある近似立方根の式は、似ている。n乗根への一般化として、
S^{1/n} \approx \dfrac{b^{n-1}(S-a^n) + a^{n-1}(b^n - S)}{b^{n-2}(S - a^n) + a^{n-2}(b^n - S)}
を考えると、それぞれ特殊な場合と見ることが出来る。ヘロンも『算数書』も、一般式を与えているわけではなく、唯一つの例題を示しているだけなので、こんな一般化は、勝手に私が考えたものに過ぎないけど。



後代の算術書に見られる平方根の計算では、『算数書』に見られる方向性は失われて、面白みがなくなる。例えば、『孫子算経』(4〜5世紀頃の成立と推定されている)の近似平方根は、以下の通り。

今有積,二十三萬四千五百六十七步。問:為方幾何?答曰:四百八十四步九百六十八分步之三百一十一。

術曰:置積二十三萬四千五百六十七步,為實,次借一算為下法,步之超一位至百而止。上商置四百于實之上(案:上商原本脫上字,今補。),副置四萬于實之下。下法之商,名為方法;命上商四百除實,除訖,倍方法,方法一退(案:原本脫方法二字,今補。),下法再退,復置上商八十以次前商,副置八百于方法之下。下法之上,名為廉法;方廉各命上商八十以除實(案:原本脫實字,今補。),除訖(案:原本脫除字,今補。),倍廉法,從方法,方法一退,下法再退,復置上商四以次前,副置四于方法之下。下法之上,名曰隅法;方廉隅各命上商四以除實,除訖,倍隅法,從方法(案:原本訛此六字,今據術補。),上商得四百八十四,下法得九百六十八,不盡三百一十一,是為方四百八十四步九百六十八分步之三百一十一。

日本語訳と解説は、以下にある。
孫子算経』訳注稿(2)
http://id.nii.ac.jp/1338/00002164/

孫子算経』(原文)
https://zh.wikisource.org/wiki/%E5%AD%AB%E5%AD%90%E7%AE%97%E7%B6%93

234567の平方根を求める問題で、"術曰"の後に、ごちゃごちゃと書いてあるが、自然数S平方根を計算する場合、x^2 \le S \lt (x+1)^2なる自然数xを見つけたら
\sqrt{S} \approx x + \dfrac{S - x^2}{2x}
だと言っている。これは、Babylonian methodと同じ式ではあるけど、反復適用するという含意は全く読み取れない。

説明が長いのは、どうやって自然数xを見つけるか書いてるからで、この部分は、いわゆる開平方の手続きによる。

234567の平方根を求めるのに、例えば、400と500から初めて、『算数書』の手順を反復してもいいのだけど、そういう方向には発展しなかったらしい。有理数のまま計算すると、分子、分母は急速に大きくはなるけど、400と500に対して、『算数書』の方法で得られる近似平方根の下限は434567/900=482.85222...で、近似の精度は、かなりいい。

手計算(そろばん等を使ってもいいけど)の時には、桁数が大きくなると、計算量は増える。『算数書』の手法はスマートだと思うけど、反復計算すると、桁数が急激に大きくなるので、桁数の大きな計算が出てこないという点では、『孫子算経』の開平方は、優れていたのかもしれない。あるいは、『算数書』の式は、幾何学的に理解しやすいとは言えないので、それが理由で嫌われたのかもしれないけど。

孫子算経』の方法は、小さい数の平方根の精度は、やや悪い。例えば、2の平方根は、3/2となる。現代人は小数を小学校で習うから、開平方を小数点以下まで続けることに違和感はない。小数を知らなくても、アレクサンドリアのテオンは1/60や1/3600のオーダーまで計算を行っているので、古代中国人も同様のことに気付いたかもしれない。あるいは、100倍したら、平方根が10倍になることには気付いていたと思うので、小数点を使う代わりに、"10の偶数べき乗"倍の平方根を計算してから、適当な10の冪乗で割るという方法を思いついててもよさそうではある(証拠はないけど)。200の平方根なら、『孫子算経』の方法では、14+1/7=99/7になり、10で割れば、99/70なので、3/2よりは大分ましになる。



孫子算経』には立方根の計算は見当たらないので、立方根の計算が書いてある算術書を探して、順番に見ていく。「算経十書」で最古のものは『周髀算経』だとされている。
《周髀算經》
https://ctext.org/zhou-bi-suan-jing/zh

これには、平方根や立方根の計算は見当たらない。次に挙げられるのが『九章算術』で、平方根と立方根の計算を扱っているけど、有理数の平方、立方の時(以下、有理数の平方、立方になる数を平方数、立方数と呼ぶ)のみしか扱っておらず、非有理数の時の扱いは、曖昧。

『九章算術』訳注稿
http://pal.las.osaka-sandai.ac.jp/~suanshu/j/publications2.html

『九章算術』は題名通り、9つの章からなる。平方根、立方根の計算は、少広の章にあり、訳注稿では、(10)と(11)に平方根と立方根の計算がある。

本題とは関係ないけど、この訳注稿では、方程と句股の章は扱われてないっぽい(『算数書』とは関係しない内容だから?)。方程の章では、負の数の概念が導入され、正数と負数同士の加減算規則が書かれている。
九章算術《方程》
https://ctext.org/nine-chapters/fang-cheng

正負術曰:同名相除,異名相益,正無入負之,負無入正之。其異名相除,同名相益,正無入正之,負無入負之。

第一文は、減算についての記述で、「同符号では相殺しあい、異符号では相補しあい、無入から正数を引くと負数になり、負数を無入から引くと正数になる」ということらしい。無入は0に相当する用語。第二文は、加算について、同様に述べられている。

劉徽の注釈では、正数は赤い算木、負数は黒い算木の本数で表すという"モデル"で説明されている。正数から負数を引いた時、あるいは、負数から正数を引く時、結果として得られる算木の本数は、引く側の算木や引かれる側の算木より多くなるので、"相益"と言ってるんだろう。

ちょっと不正確というか簡潔すぎるように思うけど、整数演算のことだと思って読めば、そう読めなくはない。この章の主題は、連立一次方程式の解法なので、負の数同士の乗算、除算は見当たらない。また、負の数が解になるケースは扱われてないので、計算途中で使う方便という程度の認識だったかもしれない。




『張邱建算経』では、非平方数の平方根、非立方数の立方根が扱われていて、古代中国の算術書で、非立方数の立方根が見られるのは、これが最初っぽい(?)。

『張邱建算経』は、割と誤植が多いように思われる。例えば、平方根を計算する問題は3題あり、原文では中巻にある。以下で見ることが出来る(18a末尾あたりから)

張邱建算經 卷中
https://www.kanripo.org/text/KR3f0039/002


与えられた円と同じ面積の正方形の一辺の長さを求めるとか、その逆といった問題が扱われている。最初の一題は、整数の平方になってる問題なのでいいとして、残りの二題は、175692の平方根を、419+131/829とし、13068の平方根を114+72/229としている。著者の意図は、自然数Sに対して、x^2 \le S \lt (x+1)^2となる自然数xを見つけて
\sqrt{S} \approx x + \dfrac{S-x^2}{2x+1}
と近似すると解釈するのが尤もらしく思える。とすると、前者は、419+131/839のミスということになる。

『算数書』の近似平方根の式で、a=x,b=x+1の場合を考えると
\sqrt{S} \approx \dfrac{S(S-a^2)b + (b^2-S)a}{b^2-a^2} = x + \dfrac{S - x^2}{2x + 1}
が得られる。

でも、多分、2x+1にしたのは、 もっと短絡的に、x + \dfrac{S - x^2}{2x}だと、常に、真の値より大きくなるのを嫌っただけじゃないかと思う。『張丘建算経』では、立方根の計算でも、同じような補正が行われている。

つまり、自然数Sに対して、x^3 \le S \lt (x+1)^3となる自然数xを(開立方で)見つけて
S^{1/3} \approx x + \dfrac{S-x^3}{3 x^2+1}
としているように読める。

立方根に関する箇所は、原文の下巻にある。
https://www.kanripo.org/text/KR3f0039/003

与えられた球と同じ体積を持つ立方体の一辺の長さを計算するとか、その逆が扱われている。最初の問題は、1572864=96*96*96*16/9の立方根を計算している。

今有立方九十六尺欲為立圓 問徑幾何

答曰一百一十六尺四萬三百六十九分尺之一萬一千九百六十八

"立圓"(立円)は、球のこと。立方(体)と同じ用法になってる。このあとに"術"(解法)の説明が長々とあるけど、開立方の説明があるだけなので省略。


立方根の問題でも、ミスがあるっぽい。

今有立圓徑一百三十二尺 問為立方幾何

答曰二百八尺三萬四千九百九十三分尺之三萬四千二十

術曰令徑再自乘九之十六而一開立方除之得立方

草曰置徑一百三十二尺再自乘得二百二十九萬九千九百六十八又以九因之得二千六十九萬九千七百一十二
又以十六除之得一百二十九萬三千七百三十二以開立方法除之得合前

この問題は、1293732の立方根(108.96…)が、208+34020/34993だと書いている。これは、108+34020/34993のミスと考えないとおかしい。そして、34993=3*108*108+1になっている。

古代インド

サンスクリットの文字記録が出現するのは、中国語やギリシャ語に比べると、やや遅い。サンスクリット話者がインドにやってきてから長い間、口承で知識を伝えていたらしい。従って、"文献"の成立時期と、文字で記録される時期が数百年以上ずれてるみたいなことが、ザラにあるよう。

サンスクリット固有の文字というのもなく、比較的新しいサンスクリット文献は、デーヴァナーガリー文字で書かれるが、デーヴァナーガリーの前身ナーガリー文字が出来たのが、西暦700〜900年頃で、それ以前にも、サンスクリットは、複数の異なる文字で記述されることがあった。

サンスクリットの碑文が出現するのが西暦100年頃で、Spitzer写本というのが、現存する最古のサンスクリット写本らしい。放射性炭素年代測定では、西暦130年前後のものとされている。これは、ブラーフミー文字などで書かれてるそうだ。法隆寺貝葉写本は、般若心経のサンスクリット版らしいけど、シッダマートリカー文字というので書かれてるらしい。

Spitzer写本は、他のインド近辺の古写本と同様、貝葉を書写材料としている。サンスクリットに限らなければ、貝葉の使用は、紀元前500年くらいまで遡れるそうだ。現時点で最古のギリシャパピルス(Derveniパピルス)や、中国の竹簡(断片しかないものも多いし、どれが最古かは分からないが)も、時代的に大きくは変わらない(成立時期の推定が100年程度の誤差はあるので)。スペイン侵入以前のメソアメリカでは、アマテという書写材料が使われてたようで、起源は分からないけど、Spitzer写本と同時期に成立したと見られる断片が発見されてるらしい。最古の貝葉写本の一つは、タミル語文法書だそうで、サンスクリットパーニニ文法成立から、それほど隔たってない時期に成立したと推定されている。

古代に、サンスクリット以外で書かれた数学書があったかもしれないけど、写本も残ってないのだと思われる。まぁ、存在しても読めないし、とりあえずサンスクリット文献に限ることにする(サンスクリットも読めないけど)。

古いサンスクリット文献は、色々な文字で記録されたので、デーヴァナーガリーで読む意味は特にないと思う。新しい文字を覚えたくない人のため(?)に、サンスクリットのAscii文字によるtransliteration方法が定義されている。transliterationの方法は複数あるけど、よく使われてるのは、IASTというものっぽい。多分、IASTが一番古いんだと思う。IASTでは、デーヴァナーガリーは、devanāgarīという風に書く。割と新しいISO15919は、他のインドの言語でも使用可能なように作られてるっぽいけど、ISO15919とIASTは、ごく僅かな違いしかない。



古代インドの"数学書"の内、アーリヤバティヤには、平方根、立方根の計算が記述されており、バクシャーリー写本には、平方根の計算が記述されている。

アーリヤバティヤは、CE499年に成立したと考えられていて、その根拠は、アーリヤバティヤの中に、"現在、(今のユガの)3/4と3600年が経過し、アーリヤバタの出生から23年になる"という記述があって、これからアーリヤバタの出生年と、アーリヤバティヤの書かれた年が決定できるらしい。

バクシャーリー写本は、現物が存在していて、放射性炭素年代測定の結果によると、3~4世紀ごろ成立したとか、複数の異なる時代に成立したものが混ざってるとか書いてある。成立年代はよく分からないけど、デーヴァナーガリー以前の写本で、シャーラダー文字というもので書かれてるらしい。

バクシャーリー写本は、1927年に出版された

The Bakhshālī manuscript : a study in mediæval mathematics (Archaeological survey of India : new imperial series, v. 43)
https://archive.org/details/in.ernet.dli.2015.23309/mode/2up

に、サンクスリットのtransliterationと、現物の写真が掲載されているっぽい。

バクシャーリー写本には、平方根の計算があるらしく、原文を確認していないけど、検索して出てくる多くの解説を見ると、S平方根に近い数aに対して、以下の式で改善された平方根が計算できると述べてるらしい。
\sqrt{S} = \sqrt{a^2 + r} \approx a + \dfrac{r}{2a} - \dfrac{(r/2a)^2}{2(a + r/2a)}

この式は、Babylonian methodを2回繰り返したものに等しい。つまり
f(x) = \dfrac{S + x^2}{2x}
として、
\sqrt{S} = \sqrt{a^2 + r} \approx f(f(a))
とするのと等価。古代インド人は、別の方法で、これを見つけたのかもしれないけど、結果的には、そうなっている。

別の見方をすれば、
\sqrt{1+x} \approx \dfrac{1+x+\frac{x^2}{8}}{1+ \frac{1}{2}x}
という有理関数近似を使ってるとも理解できる。いずれにしろ、本質的には、Babylonian methodと大差ない代物。

ただ、バークシャリー写本が反復計算を示唆しているかは知らない。負の数がなかった時代には、r=S-a^2は正の数でないとダメだったかもしれない。古代インド人が、そういう風に考えてたなら、この式で得られる近似平方根は真の値より大きいので、そのままでは反復計算できない。ヘロンの式なんかは、そういう問題が起きない形になっている。

S \lt a^2であっても、\dfrac{S}{a}として、小さい近似平方根を得られることに気付いてたかもしれない。例えば、2の近似平方根として、3/2は、出発点としては、そこそこ良い。しかし、2の平方根より大きいので、代わりに、4/3を使うことになる。a=4/3として、2の平方根に、バークシャリー写本の近似を適用すると
\sqrt{2} \approx \dfrac{4}{3} + \dfrac{1}{12} - \dfrac{1}{408}
が出てくる。これは、577/408に等しい。


ところで、インドの非常に古いヴェーダ文献にsulba sutraというのがある。そこに載っている2の平方根の近似値は、この形で表現されている。

sulba sutraは、流派ごとに微妙に異なる内容を伝えているらしく、区別する場合は、Baudhayana sulba sutraとかApastamba sulba sutraなどと呼ばれることもある(ほとんどの流派は、失伝して、現存するのは僅かとのこと)。多分、本文であるサンヒターは基本的には同一で、解釈・注釈が違うんだと思うけど、サンヒターの部分でも、多少の違いが存在するっぽい(?)。全ての流派が同時に誕生したわけではないので、それぞれのバージョンで成立年代も微妙に異なるらしいが、大雑把には、BC400年前後に成立したとされている。

Apastamba sulba sutraでは、三平方の定理が述べられた後、2の平方根の近似値が、以下のように書かれている

चतुरश्रस्याक्ष्णया रज्जुर्द्विस्तवतीं भूमिं करोति | समस्य द्विकरणी ॥ ५ ॥
प्रमाणं तृतीयेन वर्धयेत् तच् चतुर्थेनात्मचतुस्त्रिंशोनेन सविशेष ॥ ६ ॥

caturaśrasyākṣṇayā rajjurdvistavatīṁ bhūmiṁ karōti | samasya dvikaraṇī ॥ 5 ॥
pramāṇaṃ tṛtīyena vardhayet tac caturthenātmacatustriṃśonena saviśeṣa ॥ 6 ॥

私のサンスクリット歴は3時間位なので、文章が合ってるのかどうかすら分からんけど、デーヴァナーガリーとIASTを併記しておく。縦棒一本が、文の切れ目で、縦棒2本で、段落の切れ目のような感じだそうだ。

日本語では""正方形の対角線は、2倍の面積の正方形を作る。(これが)正方形のdvikarani。基準とする(正方形の一辺の)長さを1/3増やし、増分の1/4倍を加えるのは1/34だけ過剰。"みたいな感じっぽい。

辞書によると、pramāṇaはprootfという意味と、measure/scale/standardという意味があるようで、ここでは後者の意味。サンスクリットの数字は、1,2,3,4が、eka,dvi,tri,caturで、34はcatustriṁśaなので、なんかそれらしい数字が出てることくらいは分かる。dvikaraniは、two-fold producerのような意味らしい。

caturaśraは、本来は単に四角形のことだけど、ここでは正方形を意味するらしい。samaは、same/identical/equalなどの意味があって、samasyaは、多分、属格(?)という格変化だけど、ここでは"正方形の"という意味になる模様。等四辺形の略みたいな感じなんだろうけど、等四辺形は、菱形で正方形とは限らない。

それで、文章全体として、前半部分は、正方形の対角線の長さが一辺の√2倍となることを言ってるだけ。後半は、dvikaraniの近似値の説明をしている。後半は、自然言語で書くと分かりにくいけど、数式で直訳すると
\sqrt{2} \approx 1 + \dfrac{1}{3} + \dfrac{1}{3} \times \dfrac{1}{4} \times \left( 1 - \dfrac{1}{34} \right)
ということが言いたいのだと思う。これの右辺は、577/408に等しい。

こんな回りくどい書き方をしてるのは、何らかの形で、Babylonian methodを知ってた可能性が高いように思える。バークシャリー写本と同じ方法を使った可能性もあると思うけど、sulba sutraでは、この近似値を、どうやって出したか説明がないので、確実なことは何も言えない。

Katyayana sulba sutraにも、同様の内容があるようで、英語訳と解説が以下にある。

Katyayana Sulba Sutra
https://archive.org/details/in.ernet.dli.2015.142145/page/n33/mode/2up

Apastamba sulba sutraの日本語訳は以下にある。残念ながら原文は掲載されてない。

科学の名著〈1〉 インド天文学・数学集
https://www.amazon.co.jp/dp/B01M5KFZZO

バクシャーリー写本にもsulba sutraにも、立方根の計算は、ないっぽいので、この話はこれまで。



アーリヤバティヤは、天文学書で、gītīkā、gaṇita("数学")、 kālakriyā("時の計算") 、gōla("天球論")の4つの章からなる。gaṇitaでは、天文学そのものではなく、数学的内容のみが扱われてる。最初のgītīkāは、韻文の形式を表すらしいけど、ここには、天文定数が色々書かれている。アーリアバタが、この数値を、どこから持ってきたかは分からない。"正弦三角関数表"も、ここに含まれている。この三角関数表の計算法は、gaṇitaに書かれている。この計算の説明は、三平方の定理より前に出現するという謎配置になってる。

個々のトピックの説明は、簡潔な短文からなっていて、この"本"単体で理解するものではないらしい。元々は、まず暗誦できるようになって、それから先生意味を教えてもらうみたいな類のものだったのかもしれない。後世のインドの数学者も、注釈書を作ったりしていて、現代に伝わってるものは、これらの注釈書に基づくっぽい。

アーリヤバティヤは、英語訳も出ている。1976年出版の

Aryabhatiya of Aryabhata: 3 volumes
https://archive.org/details/Aryabhatiya1976

が、オープンアクセスになっている。上の"インド天文学・数学集"にも、アーリヤバティヤの日本語訳が収録されている。


平方根計算の記述はgaṇitaの4段目にあり、以下のようなもの。これで、開平方と、ほぼ同一の手続きを述べているらしい。

भागं हरेद् अवर्गान् नित्यं द्विगुणेन वर्गमूलेन |
वर्गाद् वर्गे शुद्धे लब्धं स्थानान्तरे मूलम् ॥ ४ ॥

bhāgaṃ hared avargān nityaṃ dviguṇena vargamūlena |
vargād varge śuddhe labdhaṃ sthānāntare mūlam ॥ 4 ॥

単語の意味と注釈付きの訳を眺めると、何の補足もなければ、直訳は"非平方位では常に平方根の2倍で割る。商は根の次の位を与え、平方位では商を平方して引く"みたいな感じじゃないかと思われる。ギリシャ語や中国語文献と違って、例の一つもない。簡潔すぎて、サンスクリットを読める人でも、これだけでは何を言ってるのか分からなそう。

別に10進法に依存してはいないけど、アーリヤバティヤでは10進法が使われているので、以下10進法で書く。まず、1の位、100の位、10000の位、...などが平方位(つまり、"100=10の平方"のべき乗の位)で、10の位、1000の位、...が非平方位。vargaは、辞書によると「カテゴリー、クラス、グループ」みたいな意味と、「平方」の意味がある。接頭辞a-が付くと、否定語になるっぽい。語尾の変化はよくわからんけど、avargānは、接頭辞のa+varga+語幹anなんだと思う。

そして、アーリアバタの意図は以下のようなものらしい。自然数X=100A+Bで、Bは、0〜99の範囲にあるとする。そして、x^2 \le A \lt (x+1)^2となる自然数xが分かっているとする。その時、
100A+B \approx (10x+y)^2 = 100 x^2 + 20 xy + y^2
を満たす自然数yを探す。正確には、開平方と同様、
 (10x+y)^2 \le X \lt (10x+y+1)^2
を満たすyを探すと考えるべきなんだろう。これは、0〜9の範囲にある。もし、10以上なら
(10x+y)^2 \ge (10(x+1)))^2 = 100(x+1)^2 \gt 100A \ge X
となってしまうので、条件を満たさない。

アーリアバタの指示に文字通り従うなら、
y = \left[ \dfrac{100(A-x^2)+ B}{20 x} \right]
yを決めることになる。\left[ z \right]は、zを超えない最大の整数。これは正しくないこともあるが、一旦は気にしないことにする。

20で割ってるので、実際には、Bの下一桁は、yには影響しない。つまり
y = \left[ \dfrac{10(A-x^2) + \left[B/10 \right]}{2 x} \right]
と計算できる。これで、xは、Aの近似平方根で、その2倍で割ってるから前半部分の意味が理解できる。そして、商yは、Xの近似平方根の1の位になる。

Aの近似平方根xを見つける方法が説明されてないけど、Aが、100より小さい場合(従って、Xは4桁以下の自然数)は、xは、0〜9の範囲なので、容易に見つけられる。Aが2桁以上の時は、Xの近似平方根を探す手続きを、A再帰的に適用すればいい。

ところで、yの計算で、A-x^2を使ってるので、X-(10x+y)^2まで計算しておくと都合がいい。アーリヤバティヤの2行目の文は、この計算を述べていて
X-(10x+y)^2 = 100(A-x^2) + B - 20 x y - y^2
で、100(A-x^2) + B - 20 x yは、100(A-x^2) + Bを、20xで割った余りだと理解できる。そこから更にy^2を引くと、X-(10x+y)^2が計算できる。

Xが非平方数の場合に、どうするのかは何も指示してない。アーリヤバティヤには例が一つもないけど、アーリヤバテイヤと成立時期が近いと思われる天文学書Sûrya-Siddhântaには、平方根の計算法の説明は見当たらない代わりに、平方根の計算例はあるっぽい。(英訳を見ると)Sûrya-Siddhântaとアーリヤバティヤには、同じ(正弦)三角関数表があり、同じ漸化式による計算法が記載されているっぽい。

Sûrya-Siddhântaでは、8907173の平方根(2984.4887…)が2984だと書いてある。また、1238125の平方根(1112.7106…)は、1113と書いてある。一方で、782;48=782+48/60の平方根(27+58/60+43/3600+…)は、27;59=27+59/60と書いてある。小数点以下は、60進法になってるけど、有効数字4桁くらいの精度になるまで計算するという方針だったのかもしれない。

また、アーリヤバティヤの方法だと、問題が起きる場合もある。例えば、29の2乗=841の場合、A=8,B=41で、x=2と計算される。すると
y = \left[ \dfrac{441}{40} \right] = 11
となって、0〜9の範囲で収まらない。例外があることは知りつつ要点だけ述べたのかもしれない。



同じ調子で開立方も述べられてる。立方根計算は、5段目にある。

अघनाद् भजेद् द्वितीयात् त्रिगुणेन घनस्य मूलवर्गेण |
वर्गस् त्रिपूर्वगुणितः शोध्यः प्रथमाद् घनश् च घनात् ॥ ५ ॥

aghanād bhajed dvitīyāt triguṇena ghanasya mūlavargeṇa |
vargas tripūrvaguṇitaḥ śodhyaḥ prathamād ghanaś ca ghanāt ॥ 5 ॥

直訳は"2番目の非立方位を立方根の平方の3倍で割り、商の二乗と立方根の積の3倍を1番目の非立方位から引き、立方位からは商の立方を引く"とかいう感じだと思う。ghanaは立方を意味するらしい。
(10 x + y)^3 = 1000 x^3 + 300 x^2 y + 30 x y^2 + y^3
なのでX=1000A + Bに対して、Bは0〜999の範囲。x^3 \le A \lt (x+1)^3となる自然数xが既知の時、
 y = \left[ \dfrac{1000(A - x^3) + B}{300 x^2} \right]
のように計算するというのが、アーリアバタの指示。

先ほどと同じく
 y = \left[ \dfrac{10(A - x^3) + \left[B/100 \right]}{3 x^2} \right]
となるので、A-x^3の10倍とBの100の位の数を足したものを、 3 x^2で割ることになる。

立方根でも平方根と同様の問題は存在する。例えば、
2000 = (10+y)^3 = 1000 + 300 y + 30y^2 + y^3
とした場合、アーリアバタの手続きでは、y=3で、1000+300y=1900は2000より小さいけど、30y^2を更に足すと、2170になって、2000を超える。また、非立方数の時にどうするかについても、何も言ってない。


感想。中国語圏、ギリシャ語圏、サンスクリット圏の平方根、立方根の計算を見てみると、古い時代の方法は意外と面白い。時代が下ると、みんな開平方、開立方になってしまって、とても詰まらない。今の所、確認できた範囲では、開平方が現れるのは、紀元後っぽい。古代のギリシャ語圏で開立方に言及した文献は発見できてないけど、開平方と殆ど同じなので、必要とあれば、計算できたんじゃないかと思う。


Appendix:9世紀の"ケプラー方程式"

"古代"を逸脱してしまうのだけど、古い時代に使われていた超越方程式の数値解法の事例が面白い。

A Note on 'Kepler's Equation'
https://doi.org/10.1007/BF00376451

という論文を読むと、9世紀に、形式的にKepler方程式と等価な方程式が考えられていて、数値的に解く方法が知られていたらしい。Habash al-Hasib al-Marwaziという人のZij(アラビア語で天文表を指す言葉)に書かれてるらしいけど、写真とかが公開されてたりはしないよう。されてても、アラビア語だと思うので、読めないけど

当然、9世紀にケプラーの法則が知られてたわけではなく、Kepler方程式と同形なのは、たまたまらしい。詳細は知らないけど
t = \theta(t) - \epsilon \sin \theta(t)
を満たす解\theta(t)を見つけて\sin theta(t)の表を作ると、なんか便利だと思われたらしい。\epsilonは定数。

この方程式を解くのに、Habashが用いている方法は
\theta_{0}(t) = t
\theta_{n+1}(t) = t + \epsilon \sin \theta_{n}(t)
という逐次近似だそうだ。

|\theta - \theta_{n+1}| = |\epsilon||\sin \theta - \sin \theta_{n}| \lt |\epsilon||\theta - \theta_{n}|
なので、|\epsilon| \lt 1なら収束する。

この方法を考えたのは、Habashではなく、インドあたりで知られていた方法を持ってきただけではという推測がされているようだけど、証拠もないし、それ以上のことは分からない。

近代ヨーロッパで、Kepler方程式を同様の方法で解くことを最初に考えた一人は、Eulerだと書いてある(1740年)。Lagrangeは、1760年代に、Bessel関数を使って級数展開することを考えた。

John Couch Adamsの1882年の論文

On Newton's solution of Kepler's problem
https://adsabs.harvard.edu/pdf/1882MNRAS..43...43A

には、Kepler方程式を"ニュートン法"で解くことを最初に明示的に書いたのは、1740年出版のThomas Simpsonの本

Essays on Several Curious and Useful Subjects in Speculative and Mixed Mathematics
https://archive.org/details/essaysonseveralc00simp/page/n53/mode/2up
だとしている。

Simpsonの本では、”A new Method for the Solution of Equations in Numbers"と題された節が、全く別のページ(81ページ)にあって、Simpsonが"ニュートン法"の発見者だとされる時は、こっちが引かれる。2つのテーマで、ページが離れてる理由は謎。

Adamsは、"ニュートン法"によるKepler方程式の解法が、Principia第二版以降で、微分を使わずに説明されていると書いている。

常微分方程式のLie symmetryと宇宙の作者の気持ち

たまたま、以下の論文を読んだところ、2階の常微分方程式y''=0のLie symmetryというものが、sl(3,R)をなすと書いてあった。

Symmetries, integrals and solutions of ordinary differential equations of maximal symmetry
https://doi.org/10.1007/s12044-010-0001-8
[PDF] https://www.ias.ac.in/article/fulltext/pmsc/120/01/0113-0130

この結果自体は、19世紀に、Lieによって知られていたらしい。上の論文では、独立変数、従属変数がx,yだけど、ここでは、t,xを使うことにする。


常微分方程式x''(t)=0は誰でも解けるので、対称性とか面倒なことを考えてどうするんだという気もするけど、やってみる。
\hat{x}=\hat{x}(x,t) , \hat{t}=\hat{t}(x,t)という変換で、\dfrac{d^2 \hat{x}}{d \hat{t}^2} = 0を満たすものを探そうというのが基本的発想らしい。対称性という言葉は、何かを不変にする(可逆な)変換を指すけど、今の場合は、素朴に微分方程式を不変にする変換を考えている。

このような変換は一杯ある。例えば、定数Cに対して
\hat{x} = x + C t
\hat{t} = t
なども条件を満たす。

それで、このような変換を生成する無限小変換を全部決定しようというのが、Lieの発想だったらしい。具体的には、新しいパラメータaにも依存する\hat{x}=\hat{x}(x,t,a) , \hat{t}=\hat{t}(x,t,a)
\dfrac{d \hat{x}}{d a} = X(\hat{x},\hat{t})
\dfrac{d \hat{t}}{d a} = T(\hat{x},\hat{t})
\hat{x}(x,t,0) = x
\hat{t}(x,t,0) = t
を満たし、
\dfrac{d^2 \hat{x}}{d \hat{t}^2} = 0
となるようなX(x,t),T(x,t)を探すという問題になる。X(x,t)=t , T(x,t)=0なら、
\hat{x}(x,t,a) = x + a t
\hat{t}(x,t,a) = t
で、上に書いた変換が出てくる。

\hat{t}(x,t,\epsilon) = t + \epsilon T(x,t) + O(\epsilon^2)
\hat{x}(x,t,\epsilon) = x + \epsilon X(x,t) + O(\epsilon^2)
として
\dfrac{d \hat{x}}{d \hat{t}} = \dfrac{x' + \epsilon(X_{t} + x' X_{x}) + O(\epsilon^2)}{1 + \epsilon (T_{t} + x' T_{x}) + O(\epsilon^2)} = x' + \epsilon(X_{t} + x'(X_{x} - T_{t}) - (x')^{2})T_{x} + O(\epsilon^2)
だから、
\dfrac{d \hat{x}}{d \hat{t}} = x' + \epsilon U + O(\epsilon^2)
の形になる。

同様に、やや面倒な計算をすると
\dfrac{d^2\hat{x}}{d\hat{t}^2} = \dfrac{d}{d \hat{t}} (\dfrac{d \hat{x}}{d \hat{t}}) = \dfrac{d}{d\hat{t}} ( x' + \epsilon U + O(\epsilon^2)) = x'' + \epsilon A + O(\epsilon^2)
の形で書け、
A = X_{tt} + x'(2 X_{xt} - T_{tt}) + (x')^2(X_{xx} - 2 T_{xt}) - (x')^{3} T_{xx} + x''(X_{x} - 2T_{t}) - (3 x'x'')T_{x}
となる。条件からx''=0で、また、A=0になって欲しいので
X_{tt} + x'(2 X_{xt} - T_{tt}) + (x')^2(X_{xx} - 2 T_{xt}) - (x')^{3} T_{xx} = 0
が満たされればいい。X,Tは、x'とは独立なので、X_{tt} = 2 X_{xt} - T_{tt} = _{xx} - 2T_{xt} = T_{xx} = 0が満たされればいい。

一般解は、8個のパラメータa_{1},\cdots,a_{8}を使って
X(x,t) = a_{1} + a_{2} x + a_{3} t + a_{4} xt + a_{5} x^2
T(x,t) = a_{6} + a_{7} t + a_{8} x + a_{4} t^2 + a_{5} xt
になる。

無限小変換をベクトル場X(x,t) \dfrac{\partial}{\partial x} + T(x,t) \dfrac{\partial}{\partial t}の形で書くと、論文の式(2.4)にある通りの結果となる
 G_{1} = x \partial_{x}
 G_{2} = \partial_{x}
 G_{3} = t \partial_{x}
 G_{4} = \partial_{t}
 G_{5} = t \partial_{t} + \dfrac{1}{2} x \partial_{x}
 G_{6} = t^2 \partial_{t} + xt\partial_{x}
 G_{7} = x \partial_{t}
 G_{8} = x t \partial_{t} + x^2 \partial_{x}

x,tを、力学の空間変数と時間変数と思えば、x''(t)=0は、ニュートン力学自由粒子運動方程式で、G_{2},G_{4}は空間並進と時間並進に対応し、G_{3}は、Galilean boostに対応する。この3つが、空間1次元のガリレイ対称性を作る。

G_{1}は、空間変数のスケール変換で、G_{5}は、時間変数と空間変数をスケール変換する。G_{6}が生成する大域変換を計算すると
\hat{x}=\dfrac{x}{at + 1} . \hat{t} = \dfrac{t}{at + 1}
になる。

G_{4},G_{5},G_{6}は閉じたLie環をなし、sl(2)と同型になる。この3つの無限小変換で生成される大域変換は
\hat{t} = \dfrac{a t + b}{c t + d} , \hat{x} = \dfrac{x}{c t + d}
の形で、a d - b c=1を満たす。特に、時間変数に対しては、一次分数変換の形。

G_{7}は、Galilean boostG_{3}と、x,tの役割が逆になったもの。G_{8}は、G_{6}と、x,tの役割が逆になったもの。


G_{2},G_{3},G_{4}は、一次元ガリレイ代数の生成元だけど、G_{2},G_{3},G_{4},G_{5},G_{6}は、一次元シュレディンガー代数(シュレディンガー群という方が一般的だけど、Lie環を考えるので、"代数"としておく)で質量を0にしたものと同型になる。シュレディンガー代数は、自由粒子シュレディンガー方程式の対称性で、多分、1970年代初頭に、そう命名された。

The maximal kinematical invariance group of the free Schrödinger equation
https://doi.org/10.5169/seals-114417

シュレディンガー代数はガリレイ代数を部分群に含んでいて、量子力学で、ガリレイ対称性を考えるには、中心拡大が必要で、空間並進演算子とGalilean boost演算子は非可換になる。古典的には、[G_{2},G_{3}]=0で、Galilean boostと空間並進は可換なので、中心拡大項を0にしないと、上のLie symmetryと同型にならない。

The Maximal Kinematical Invariance Group of the Harmonic Oscillator
https://www.e-periodica.ch/cntmng?pid=hpa-001:1973:46::960

調和振動子シュレディンガー方程式も、同じ対称性を持つ。調和振動子自由粒子は、互いに変換可能なので、当たり前かもしれない(物理的には、散乱状態しかない自由粒子と、束縛状態しかない調和振動子が"同じ"というのは、不思議なことかもしれない)。熱方程式も同じ対称性を持つので、この対称性は、19世紀には知られていたという記述も見るけど、確認はしてない。


そんなわけで、8つのLie symmetry、3つのガリレイ対称性がある。G_{2},G_{3},G_{4},G_{5},G_{6}は、シュレディンガー方程式の対称性から来ているのだけど、ネーターの定理によって対応する保存量が存在する。なので、ネーター対称性とでも呼べばいいかもしれない。

Lie symmetryから発見的に、ネーター保存量を見つけるために、まず、ベクトル場G_{i}のprolongationを考える。
\hat{t}(x,t,\epsilon) = t + \epsilon T(x,t) + O(\epsilon^2)
\hat{x}(x,t,\epsilon) = x + \epsilon X(x,t) + O(\epsilon^2)
の時
\dfrac{d \hat{x}}{d \hat{t}} = x' + \epsilon U(x,t,x') + O(\epsilon^2)
だったので、u=x'を新しい変数と見て
G = X(x,t) \partial_{x} + T(x,t) \partial_{t}
に対して、
G^{[1]} = X(x,t) \partial_{x} + T(x,t) \partial_{t} + U(x,t,u) \partial_{u}
という(x,u,t)空間のベクトル場を考える。これは、ベクトル場を、Jet spaceに拡張してることに相当する。

Uは既に計算した通り
U(x,t,u) = X_{t} + u(X_{x} - T_{t}) - u^2 T_{x}
となる。

例えば
G_{1}^{[1]} = x \partial_{x} + u \partial_{u}
G_{2}^{[1]} = \partial_{x} + 0 \cdot \partial_{u}
G_{3}^{[1]} = t \partial_{x} + \partial_{u}
G_{4}^{[1]} = \partial_{t}
G_{5}^{[1]} = t \partial_{t} + \dfrac{1}{2} x \partial_{x} - \dfrac{1}{2} u \partial_{u}
G_{6}^{[1]} = t^2 \partial_{t} + xt \partial_{x} + (x - ut) \partial_{u}
など。

(x,u,t)空間上の二次微分形式として
\omega = du \wedge dx - u d u \wedge d t
を定義する。運動量p = m uとエネルギーH = \dfrac{ m u^2}{2}を導入すると
\omega = \dfrac{1}{m} (dp \wedge dx - d H \wedge dt)
であることに注意。extended phase spaceで、ハミルトニアン=エネルギーという拘束条件を課したものと思ってもいい。

(x,u,t)空間は、3次元空間なので、\omegaは、symplectic形式ではないけど、symplectic多様体のネーターの定理と同様に、(x,u,t)空間のベクトル場Zに対して
\omega(Z) = d I_{Z}
を満たすI_{Z}が存在すれば、ネーター保存量となる。保存量なので、I_{G}の時間微分は0になる。

今考えるのは、(存在すれば)I_{G_{n}^{[1]}}なので、これをI_nと略記することにする。

\omega(G_{2}^{[1]}) = du = dI_{2}
なので、I_{2} = u + \mathrm{const.}になる(積分定数は重要でないので、このあとは省略する)。

u=x'なので、\dfrac{d I_{2}}{dt} = u' = x'' = 0になる。G_{2}は空間並進で、運動量保存則が対応する。同様に計算していくと
I_{3} = u t - x
I_{4} = -\dfrac{u^2}{2}
で、重心位置の保存と、エネルギーの保存則が出る。この3つは、ガリレイ代数の対称性に対応する保存則なので、よく知られてる。シュレディンガー代数の生成元に対しては、
I_{5} = \dfrac{1}{2} u (x-ut)
I_{6} = -\dfrac{1}{2} (x-ut)^2
となる。これらの保存量に名前はなく、既知の保存量を使って書けてる。特に新しいものは出ないが、これも自由粒子の保存量ではある。

そもそも、自由粒子の解は、x(t)=A+ Btで、保存量はu = Bx - ut = Aの関数で書けるものしか出ない。


G_{1},G_{7},G_{8}については、対応する保存量が存在しないことを確認する。
\omega(G_{1}^{[1]}) = x du - u dx + u^2 dt
で、(x,u,t)空間は、幾何学的には、単なるEuclid空間なので、これが完全形式になることと、閉形式であることは同値。しかし
d \omega(G_{1}^{[1]}) = -2 (du \wedge dx - u du \wedge dx) = -2 \omega
なので、閉形式ではない。

同様に
d \omega(G_{7}^{[1]}) = 3 u \omega
d \omega(G_{8}^{[1]}) = (3ux-3u^2t-u^2)du \wedge dt + (2 - ut -u)dx \wedge du + u^2 dx \wedge dt
で、これらのどのような一次結合も、閉形式にはならない。


以上から、Lie symmetryの内、ネーター保存量を持つものが、シュレディンガー代数をなすということができる。Lie symmetryは、常微分方程式だけでなく、偏微分方程式でも考えることができて、自由粒子シュレディンガー方程式には、G_{7},G_{8}に対応するLie symmetryは出てこない。

自由粒子のシュデレィンガー方程式のLie symmetryは、シュレディンガー代数、スケール変換、あと謎の無限個の族からなる。スケール変換は、色んなところで出てくる"対称性"で、時間の一次分数変換にもスケール変換は含まれるけど、一般的には、ネーターの定理が成立するとは限らない。ネーターの定理を拡張して、スケール変換にも使える形にしようと試みる人もいる。

A Note on the scale symmetry and Noether current
https://arxiv.org/abs/hep-th/9807086

A generalized Noether theorem for scaling symmetry
https://arxiv.org/abs/1903.05070

まぁ、ネーターの定理自体は、純粋に数学上の定理なので、何かしら一般化は出来るかもしれないけど、それに意味があるかは別問題。ネーターの定理が成立しないLie symmetryは、スケール変換以外にも出てくることがある。

この宇宙を作ったクソ野郎が、どういうつもりだったのか分からないけど、ネーターの定理が成立するような対称性を重視してるように見える。



冒頭の論文では、3階の常微分方程式に対して、Lie symmetryを見ていて、これはこれで面白いと思う。上で見たような点変換以外に、点変換のprolongationでない接触変換が3つあって、それらを合わせると、Lie代数sp(4,R)をなすらしい。点変換だけ見ても、十分な対称性を得られないというのは、多分、一般的によくあることなんだろうと思う。

3階の微分方程式接触変換は、最近でも論文が書かれている。例えば、以下は、適当に検索して出てきた論文。

Third order ordinary differential equations and Legendre connections
https://doi.org/10.2969/jmsj/05040993

Geometry of Third-Order Ordinary Differential Equations and Its Applications in General Relativity
https://arxiv.org/abs/0810.2234

Geometry of third-order ODEs
https://arxiv.org/abs/0902.4129

偏微分方程式では、KdV方程式の対称性が、点変換だけ見ていたのでは十分でない別の例だと思う。KdV方程式を不変に保つ点変換を考えても、1次元ガリレイ群とスケール変換しか出てこない。KdV方程式の対称性を扱う枠組みとしては、Lie-Bäcklund symmetryなどがある。KdV方程式のLie-Bäcklund symmetryと保存則には、対応関係がある。KdV方程式が、Euler-Lagrange方程式になるようなLagrangian密度は知られてないが、KdV方程式と殆ど等価な方程式を与えるLagrangian密度は知られている(誰が考えたのか知らないけど、arXiv:9210226v1とかに書いてある)。


点変換でない対称性を調べるという問題は面白いと思うけど、置いておいて、ここでは、二次元空間の自由粒子についても、一次元空間と同様の議論が出来ることを概観しておく。空間変数をx,yに、時間変数をtにする。

微分方程式x''(t)=0,y''(t)=0を考える。やることは、さっきと同じだけど、計算は、やや長いので省略。結果として、Lie symmetryの生成子は15次元のLie環を作る。出てくるベクトル場は
\partial_{t} , \partial_{x} , \partial_{y}
t \partial_{x} , t\partial_{y} , t \partial_{t} , x\partial_{x} , x \partial_{y} , x\partial_{t} , y \partial_{x} , y\partial_{y} , y \partial_{t}
t(x\partial_{x} + y\partial_{y} + t \partial_{t}) , x(x\partial_{x} + y\partial_{y} + t \partial_{t}),y(x\partial_{x} + y\partial_{y} + t \partial_{t})
になる。一般に、空間の次元がnの場合は、(n+2)^2-1個のベクトル場が出ると思う。

自由粒子の場合、xとyの任意の線型結合も、2階微分が0になるので、単なる空間回転より多くの対称性がある。二次元シュレディンガー代数は9次元で、内一つは、中心拡大項。なので、空間変数の線形変換は殆どがネーター対称性ではない(と予想される)。

ベクトル場G = X \partial_{x} + Y \partial_{y} + T \partial_{t}のprolongationを考える。u=x',v=y'として
G^{[1]} = X \partial_{x} + Y \partial_{y} + T \partial_{t} + U \partial_{u} + V \partial_{v}
U = X_{t} + u(X_{x} - T_{t}) + v X_{y} - u^2 T_{x} - uv T_{y}
V = Y_{t} + + u Y_{x} + v(Y_{y} - T_{t}) - v^2 T_{y} - uv T_{x}
となる。

さっきと同じく、(x,y,u,v,t)空間上の二次微分形式を
\omega = du \wedge dx + dv \wedge dy - u du \wedge dt - v dv \wedge dt
で定義する。

例えば、G = y \partial_{x}に対して、prolongationは
G^{[1]} = y \partial_{x} + v \partial_{u}

\omega(G^{[1]}) = y du - v dx + uv dt
となる。これは閉形式ではない。




Lie symmetryが物理の問題で考えられるようになったのは、比較的最近のことらしい。多くの(簡単な)微分方程式に対して、Lie symmetryの計算が、1970年代頃から、行われるようになったっぽい。

(古典)一次元調和振動子も、自由粒子と同型なLie symmetryを持つことは、以下に書いてある。
The Lie group of Newton's and Lagrange's equations for the harmonic oscillator
https://doi.org/10.1088/0305-4470/9/4/007

ネーター保存量を持つものに限定すると、シュレディンガー代数と"同じ"Lie環が得られるという議論は、以下でされている。
Symmetry groups and conserved quantities for the harmonic oscillator
https://iopscience.iop.org/article/10.1088/0305-4470/11/2/005

調和振動子自由粒子は、互いに変換可能なので、同じ結果が成立するのは、当然ではある。


(古典)Kepler系のLie symmetryは、以下で計算されている。

On the Lie symmetries of the classical Kepler problem
https://doi.org/10.1088/0305-4470/14/3/009

この論文では、スケール変換とRunge-Lenzベクトルが対応するようなことを述べていて、何を言ってるのかよく分からないけど、多分正しくない。そもそも、対称性の(生成子の)数と保存量の個数が合ってない。

Noether’s theorem and Lie symmetries for time-dependent Hamilton-Lagrange systems
https://doi.org/10.1103/PhysRevE.66.066605
[PDF] https://web-docs.gsi.de/~struck/hp/noether/noether.pdf

は、流し読みする限り、(二次元Kepler系における)Runge-Lenzベクトルに対応するLie symmetryとして
(q_{1}\dot{q_{1}} + q_{2}\dot{q_{2}}) \dfrac{\partial}{\partial q_{n}}
を考えればいいと述べている。これも正しくないように思える。これらの生成子が、どういうLie環を作るのか、prolongationを計算しないといけないので、やってないけど、正準形式でのPoisson括弧と同じ交換関係を満たさないのでないかと思われる。

正しい答えは、正準形式で書いて、Runge-Lenzベクトルに対応する無限小生成子の計算をすれば分かる(Poisson括弧の計算をするだけ)はず。多分、以下が答えだと思う。
K_{j} = \sum_{i=1}^{n} (2 \dot{q_i} q_{j} - q_{i} \dot{q_{j}} - (\sum_{k=1}^{n} q_{k}\dot{q_{k}}) \delta_{ij})\dfrac{\partial}{\partial q_{i}}
これは、点変換だけ見ていても、ネーター対称性を得るには十分でないことがあるという一つの例と思われる(解析力学の教科書とかで、一般化座標の微小変換を考える時、どのようなクラスの変換を考えてるか明記されてないことが多いようなので、注意する必要がある)。

前者の論文のAbstractで、Runge-Lenzベクトルが、"“non-Noether invariant”(この表現は、後者の論文にある)だとされているけど、正準形式で見れば、"Noether invariant”と見る方が自然。




spectrum generating algebraとの関係について。

調和振動子の対称性とspectrum generating algebra
https://m-a-o.hatenablog.com/entry/20161117/p1

では、調和振動子の"対称性"として、symplectic代数の作用を書いた。spectrum generating algebraが何を不変にしてるのかは、よく分からない(古典的には、symplectic形式を不変にしてはいるけど...)

symplectic群Sp(2n,\mathbf{R})の部分群として、SL(2,\mathbf{R}) \times O(n)があって、自由粒子の場合、SL(2,\mathbf{R})が、時間の一次分数変換、O(n)は空間回転を与える。自由粒子の場合、空間並進とGalilean boostの無限小生成元、ガリレイ代数の中心拡大項(質量に相当する)は閉じたLie環を作り、Heisenberg代数と同型になる(古典的な場合は、中心拡大項はいらないので、群の次元は一つ小さくなる)。sl(2,\mathbf{R}) \oplus o(n)とHeisenberg代数の半直積が、シュレディンガー代数と同型になる。

こうして、Heisenberg代数が、シュレディンガー代数の部分代数になっている。普通、量子力学におけるHeisenberg代数は、運動量演算子と位置演算子CCR関係から出てくるものだけど、ここのHeisenberg代数は、ガリレイ代数の部分代数に過ぎず、幾分異なる出自を持っている。運動量演算子は、そのままだけど、対になるのは、Galilean boostで、その非可換性は、ガリレイ代数の中心拡大の帰結。

Galilean boost演算子は、(自由粒子の)ハミルトニアンとは可換でないけど、i \hbar \dfrac{\partial}{\partial t} + \dfrac{\hbar^2}{2m} \Deltaとは可換になる。

確率微分方程式のカオス展開による近似解法

適当な参考文献。

The asymptotic error of chaos expansion approximations for stochastic differential equations
https://arxiv.org/abs/1906.01209

Dynamical Polynomial Chaos Expansions and Long Time Evolution of Differential Equations with Random Forcing
https://arxiv.org/abs/1505.00047

Wiener Chaos Expansion and Numerical Solutions of Stochastic Partial Differential Equations
https://thesis.library.caltech.edu/1861/




多項式カオス展開と呼ばれる方法によって、確率常微分方程式を等価な無限個の常微分方程式に置き換えることが出来る。同様に、確率偏微分方程式も無限個の偏微分方程式に置き換えることが出来る。これらの置き換えは厳密に正しいが、無限個の微分方程式から、適当に有限個の方程式を選んで、近似計算に使うことができる。

確率微分方程式数値計算への応用が現れたのは、1980年代か1990年代初頭あたりっぽい(ちゃんと調べてない)けど、多項式カオス自体は、Wienerの1938年の論文で導入された。
The Homogeneous Chaos
http://users.dimi.uniud.it/~giacomo.dellariccia/Table%20of%20contents/Wiener1938.pdf

論文タイトルは、"homogeneous chaos"だけど、本文中で、polynomial chaosと呼んでいる。現在、"カオス"という名前は、カオス理論で使われる方が有名になった。カオス理論と多項式カオス展開は、関係ない。ただ、Wienerは、上論文中で、

Physical theories of chaos, such as that of turbulence, or of the statistical theory of a gas or a liquid, may or may not be theories of equilibrium.

と書いているので、動機自体は、カオス理論と近い所にあったと思われる。個人的に、カオス理論は別に好きじゃないけど、カオス展開は面白いので、知名度が上がってほしい。



Wienerの論文は長い上に、少なくとも現代の人にとっては、あんまり読みやすくない。CameronとMartinの1947年の論文は短いし、現代の人でも読みやすい。

The Orthogonal Development of Non-Linear Functionals in Series of Fourier-Hermite Functionals
https://doi.org/10.2307/1969178

CameronとMartinの1947年の論文では、いわゆるCameron-Martin空間の正規直交基底がimplicitに導入されている。


標準Wiener過程の適当な説明。Wiener過程は、ホワイトノイズの積分として「定義」されるガウス過程であるという性質と、推移確率密度が熱核で書けるマルコフ過程であるという性質を持つ。2つの性質は、見かけ上異なるけど、どっちを出発点としてもいい。WienerがWiener過程を考えた時、まだマルコフ過程ガウス過程も定義されてなかったけど、Wienerの出発点は、どっちかというと後者に基づいていたっぽい(論文を見てないので、二次情報からの推測)。今となっては、前者の方が、Wiener過程が何故基本的なのか理解しやすいと思う。

確率解析でWiener過程が重要な理由も前者に由来すると思う。物理の教科書では、ホワイトノイズZ_tを、E[Z_t=0]かつE[Z_s Z_t]=\delta(s-t)を満たす"定常ガウス過程"もどきと見なして扱うのが一般的に行われている。数学では、自己相関がデルタ関数であるために、ホワイトノイズを直接扱いづらいので、ホワイトノイズの積分であるWiener過程をベースに理論を作るのが一般的な処方となった。

W_tを標準Wiener過程とすると、形式的に
W_{t} = \displaystyle \int_{0}^{t} Z_{s} ds
と考えられる。Wiener過程の標本経路が微分不可能であることを思うと、ホワイトノイズのサンプルは、少なくとも、普通の関数であると考えることはできない。この積分を、正当化する方法は分からないけど、無視して進む。

E[\cdot]で確率変数の期待値を表す時、期待値を取る操作が、和を取る操作と可換であることを、積分にも適用できると考えると
E[W_{t}] = E \left[ \displaystyle \int_{0}^{t} Z_{s} ds \right] = \displaystyle \int_{0}^{t} E[Z_s] ds = 0

E[W_sW_t] = E \left[\displaystyle \int_{0}^{s} \int_{0}^{t} Z_u Z_v du dv \right] = \displaystyle \int_{0}^{s} \int_{0}^{t} E[Z_uZ_v] du dv = \displaystyle \int_{0}^{s} \int_{0}^{t} \delta(u-v) du dv = \mathrm{min}(s,t)
が成立してほしい。この議論を正当化することはできない(か、出来ても却って準備が面倒になる)ので、この性質で特徴付けられる(非定常)ガウス過程というのが、標準Wiener過程の一つの定義になる。

E[W_s W_t] = \mathrm{min}(s,t)が成立するなら、E[(W_t-W_s)^2]= t + s -2 \mathrm{min}(s,t)=|t-s|などから、増分が分散|t-s|正規分布に従うと期待される。また、t_1 \lt t_2 \lt t_3に対して、W_{t_3}-W_{t_2}W_{t_2}-W_{t_1}は、
E[(W_{t_3} - W_{t_2})(W_{t_2}-W_{t_1})] = 0
で互いに直交し、正規分布の性質から、重なりのない区間の増分は互いに独立であることも分かる。推移確率密度が熱核で書けるマルコフ過程になっているという出発点からも、同じ要請が出る。


h \in L^2[0,1]に対して、\displaystyle \int_{0}^{1} h(t) dW_tが、平均0で分散がhのノルムの二乗で与えられる正規分布に従うことは、積分の定義に戻れば、すぐに理解できる(hは確定的な関数なので、特に難しいことはない)

L^2[0,1]の正規直交基底e_n(n=1,2,\cdots)を適当に固定して、\chi_I(t)は、 t \in Iの時に1、それ以外では0となる関数とする。0 \le s \le 1に対して、
c_n = \displaystyle \int_{0}^{1} \chi_{[0,s]}(t) e_n(t) = \displaystyle \int_{0}^{s} e_n(t)
を定義して、
W_s = \displaystyle \int_{0}^{s} dW_{t} = \displaystyle \int_{0}^{1} \chi_{[0,s]}(t)dW_{t} = \displaystyle \int_{0}^{1} \sum_{n=1}^{\infty} c_n e_n(t) dW_{t}
で、形式的には
W_s = \displaystyle \sum_{n=1}^{\infty} \left(\displaystyle \int_{0}^{1} e_n(t) dW_{t} \right) c_n = \displaystyle \sum_{n=1}^{\infty} \left(\displaystyle \int_{0}^{1} e_n(t) d W_{t} \right) \displaystyle \int_{0}^{s} e_n(t)
という式変形ができる。実用上は、どうせ有限和しか取らないけど、収束するかどうかについては、もう少し議論がいる。

証明は、例えば、1966年のSheppという人の論文に書かれていて(Theorem4)、”This result appears to be new except for two special cases of Wiener and Levy"と述べられている。
Radon-Nikodym derivatives of gaussian measures
https://doi.org/10.1214/aoms/1177699516

どうやって証明するのか知らないけど、Kolmogorovのthree-series theoremというものがあるらしく、それを使っている。Kolmmogorovが1925年に証明して、確率論の人には、よく知られた定理らしい。原理的には、勝手な正規直交基底を取っていいとしても、数値計算する時には、有限項で打ち切ることになるので、誤差評価ができることが望ましい。


一般のGauss過程に対しては、類似の定理として、Karhunen–Loèveの定理と呼ばれるものがある(Wiener過程に対しては、基底の形が制限されるので特殊ケースしか示せないけど)。この定理は、1947年に、定常ガウス過程に対して、KacとSiegertによって述べられている。

An explicit representation of stationary gaussian process
https://doi.org/10.1214/aoms/1177730391

KacとSiegertが定常性を仮定した理由は、分からないけど、(多分ガウス過程を最初に導入したと思われる)1934年のKhintchineの論文でも、定常性が仮定されているので、それに従っただけかもしれない。

Korrelationstheorie der stationären stochastischen Prozesse
https://doi.org/10.1007/BF01449156

KarhunenとLoèveも、1947年前後、独立に、なにか論文を書いたらしいけど、電子化されたものを発見できなかったので内容の詳細は不明。1951年のCramérの論説に2人への言及が見られる。
A Contribution to the Theory of Stochastic Processes
https://projecteuclid.org/proceedings/berkeley-symposium-on-mathematical-statistics-and-probability/Proceedings-of-the-Second-Berkeley-Symposium-on-Mathematical-Statistics-and/Chapter/A-Contribution-to-the-Theory-of-Stochastic-Processes/bsmsp/1200500237




何を典拠にするかはともかく、平均0、分散1の正規分布に従う無限個の独立な確率変数\xi_i(i=1,2,\cdots)とL^2[0,1]の正規直交基底e_nを適当に取って、標準Wiener過程を
W_t= \displaystyle \sum_{n} \xi_n \displaystyle \int_{0}^{t} e_{n}(u)du
と表示することができる。上の話は全部忘れて、これを定義のように思っても、差し支えない。

\xi_1,\cdots,\xi_nは互いに独立で、標準正規分布に従うので、
E[f(\xi_1 , \cdots , \xi_n)] = (2 \pi)^{-n/2} \displaystyle \int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty} f(u_1,\cdots ,u_n) e^{-(u_1^2 + \cdots + u_n^2)/2} du_1 \cdots du_n
が成り立つのは、直感的に明らかで、\xi_1,\cdots,\xi_nの関数を、Hermite多項式の直積を基底にして展開したくなる。これが、polynomial chaosの"polynomial"の起源。

多重指数\alpha=(\alpha_1,\alpha_2,\cdots)(\alpha_j自然数\sum_{n} \alpha_n \lt \infty)に対して
T_{\alpha}(\mathbf{\xi}) = \displaystyle \prod_{n=1}^{\infty} H_{\alpha_n}(\xi_n)
を定義する。H_k(x)は正規化したHermite多項式

一般の標準Wiener過程の汎関数
f[W_t] = \displaystyle \sum_{\alpha} c_{\alpha} T_{\alpha}(\mathbf{\xi})
という形で展開できる。この展開は、L^2[0,1]の正規直交基底e_nの取り方に依存していることは忘れるべきでない。

ここで、
c = \displaystyle \sum_{\alpha} T_{\alpha}(\mathbf{\xi})
に対して
E[c T_{\beta}(\mathbf{\xi})] = c_{\beta}
なので、c=0という確率論的な条件は、T_{\alpha}の係数が全部0という条件と同じ。

確率常微分方程式の場合は、係数c_{\alpha}が時間に依存し、確率偏微分方程式なら位置(や時間)に依存すると考える。係数自体は、普通の(確率的でない)微分方程式に従い、全ての係数を決定すれば、確率微分方程式が解けたことになる。


簡単な例として、Ornstein-Uhlenbeck過程の従う確率常微分方程式
d x_{t} = (a - b x_t) dt + c \cdot dW_{t}
x_{t=0} = x_{0}
を考える。これは、Langevin方程式を少し一般化したものと見なせる。x_tは、Langevin方程式の文脈では、時間tにおける粒子の速度と解釈される。物理の教科書では、ホワイトノイズZ_{t}が定義できるとして、dW_{t} = Z_{t} dtであるかのように書かれてるのが普通。

Langevinは1908年に、(フランス語タイトルを英訳すると)"On the Theory of the Brownian Motion"という論文を書き、UhlenbeckとOrnsteinは1930年に、"On the Theory of the Brownian Motion"という論文を書いている。
Sur la théorie du mouvement brownien
https://fr.m.wikisource.org/wiki/Sur_la_th%C3%A9orie_du_mouvement_brownien

On the Theory of the Brownian Motion
https://doi.org/10.1103/PhysRev.36.823

論文を読むと、方程式をLangevin方程式と呼び、その解をOrnstein-Uhlenbeck過程と呼ぶのは、一応、妥当に思える(他に、貢献者がいないかどうかは分からないが)。Doobは1942年の論文で、Ornstein-Uhlenbeck過程を、(非自明な)定常Gauss-Markov過程として特徴付けた。
The Brownian Movement and Stochastic Equations
https://doi.org/10.2307/1968873

重要というわけではないけど、以下の論文は、読み物として面白い。
Laplace and the origin of the Ornstein-Uhlenbeck process
https://doi.org/10.2307/3318524


で、変数x_t
x_t = \displaystyle \sum_{\alpha} x_{\alpha}(t) T_{\alpha}(\mathbf{\xi})
と展開する。

形式的に計算して
d x_{t} = \displaystyle \sum_{\alpha} x_{\alpha}'(t) T_{\alpha}(\mathbf{\xi}) dt
 (a - b x_{t})dt = (a - b \displaystyle \sum_{\alpha} x_{\alpha}(t) T_{\alpha}(\mathbf{\xi})) dt
 c \cdot dW_{t} = \displaystyle \sum_{\alpha} \xi_{n} e_{n}(t) dt
と思うと、無限個の微分方程式
x_{\alpha}'(t) = (a I_{\{\alpha=0\}} - b x_{\alpha}(t)) + c \displaystyle \sum_{n} e_{n}(t) I_{\{\alpha=\epsilon_n\}}
x_{\alpha}(0) = I_{\{\alpha=0\}} x_0
に書き直すことが出来る。

\epsilon_nは、k番目だけ1で、他は0という形の多重指数((\epsilon_n)_{k}=\delta_{nk})。従って、T_{\epsilon_n}(\mathbf{\xi}) = H_1(\xi_n) = \xi_nとなる。I_{\{\alpha=\beta\}}は、\alpha=\betaの時、1で、それ以外の時は0となる定数。

一般には、確率変数の積を普通の積と同じように扱って、F \cdot dW_{t} = (F \cdot Z_{t})) dtみたいな形式的計算をやると、間違った結果が得られる。Ornstein-Uhlenbeck過程の場合には、確率変数の積が出てこないので問題なかった。しかし、例えば、幾何ブラウン運動の満たす確率微分方程式だと、失敗する。それでも、確率変数の積をWick積で解釈してやると、うまくいくっぽい(f(t, x_t)dtみたいな式を計算する時は、確率変数の積も、普通の積として計算していい)。何でWick積にするといいのか、直感的な説明が欲しいけど分からない。

cf)幾何ブラウン運動で、うまくいかない理由を検討した時に、
chaos expansion methods in malliavin calculus. a survey of recent results
https://sites.dmi.uns.ac.rs/nsjom/Papers/45_1/NSJOM_45_1_045_103.pdf
のRemark5.9とかを見た。他に、
ホワイトノイズ解析の新展開
https://www.math.is.tohoku.ac.jp/~obata/research/file/2008-IIAS-report.pdf
の5.3節とか。しかし、何で、Wick積にするのかは、(私にとっては)magicのままである。


得られた常微分方程式は、3パターンに場合分けできて、\alpha=\mathbf{0}の時は
x_{\mathbf{0}}'(t) = a - b x_{\mathbf{0}}(t)
x_{\mathbf{0}}(0) = x_0
となり、解は
x_{\mathbf{0}}(t) = x_0 e^{-bt} + \dfrac{a}{b}(1 - e^{-bt})

第2のパターンは、\alpha=\epsilon_nの時で
x_{\epsilon_n}'(t) = -b x_{\epsilon_n}(t) + c \cdot e_n(t)
x_{\epsilon_n}(0)=0
となり、解は
x_{\epsilon_n}(t) = c \cdot e^{-bt} \displaystyle \int_{0}^{t} e_n(u) e^{bu} du

第3のパターンは、|\alpha| = \sum_{n=1}^{\infty} \alpha_n \ge 2の時で
x_{\alpha}'(t) = -b x_{\alpha}(t)
x_{\alpha}(0) = 0
で、解は
x_{\alpha}(t)=0

一方、確率積分を使って解を書けば
x_t = x_0 e^{-bt} + \dfrac{a}{b}(1-e^{-bt}) + c \cdot e^{-bt} \displaystyle \int_{0}^{t} e^{bs} \cdot dW_s
で、両者は同じ結果を与える。


Ornstein-Uhlenbeck過程に対しては、全ての係数が満たす微分方程式を解けるので、カオス展開は、厳密解を与える。一般的には、カオス展開で与えられる微分方程式が解けるとは限らない。けど、係数の式の形が単純な場合は、微分方程式そのものは書くことができる。なので、(確率的でない)微分方程式を数値的に解くことで近似解が得られる。そこでは、伝統的な微分方程式数値計算アルゴリズムが、自由に使える。

カオス展開で得られた微分方程式の解が分かれば、適当な個数の標準正規分布に従う乱数を振って、エルミート多項式に代入して、
 \displaystyle \sum x_{\alpha}(t) T_{\alpha} (\mathbf{\xi})
を計算するだけで、一つのpath(の近似)をサンプリングできる。


確率常微分方程式の数値解法として、Euler-丸山法というのが、よく出てくる。これは、常微分方程式のEuler法を、確率常備分方程式に、そのまま拡張したもので、一番簡単なアルゴリズムと思う。
Euler–Maruyama method
https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method

Ornstein-Uhlenbeck過程を例にして、カオス展開による数値解と、Euler-丸山法による解と比較することを考える。カオス展開での具体的な計算には、L^2[0,1]の正規直交基底を決める必要があるので、
e_{1}(t) = 1
e_{n}(t) = \sqrt{2} \cos( (n-1)\pi t) (n=2,3,\cdots)
と取る。これは、Wienerが使ったらしい基底。


Kを適当に選んだ自然数として、正規乱数\xi_1,\cdots,\xi_{K}に対して、近似Wiener過程を、以下の有限和で生成する。
W(t) = \displaystyle \sum_{n=1}^{K} \xi_n \displaystyle \int_{0}^{t} e_n(u)du

実験した限り、t \gt sに対して、W(t)-W(s)は、t-sがある程度大きい時には、分散t-s正規分布に従ってそうな感じがする。t-sを、小さくすると、期待される分布からの乖離は大きくなって、分散の大きさが期待される値と10倍位違うということも起こる。

t-sをどこまで小さく出来るかは、Kの大きさに依存しているように思われる。標準正規乱数Xを使って、W(t)=W(S) + \sqrt{|t-s|}Xによって逐次的にWiener過程を生成する場合でも、|t-s|より小さい時間幅では線形補間されるだけので、これが悪いわけではない。結局、近似の代償は、同じ部分に出るということらしい。

近似Wiener過程のサンプルを繰り返し作って、それをW^{(n)}(t)として、
E[W_s W_t] \approx \dfrac{1}{N} \displaystyle \sum_{n=1}^{N} W^{(n)}(s) W^{(n)}(t)
を評価すると、min(s,t)に近い値が得られる。こうした実験から、カオス展開によるWiener過程の近似は、悪くなさそうだと分かる。


近似Wiener過程を使って、Euler-丸山法では、Ornstein-Uhlenbeck過程は
Y(0) = x_0
Y(t_{n+1}) = Y(t_{n}) + (a - b \cdot Y(t_n))(t_{n+1} - t_n) + c(W(t_{n+1}) - W(t_{n})
で計算される。近似Wiener過程のサンプルごとに、近似Ornstein-Uhlenbeck過程のサンプルが一つ作れる。

通常は、正規乱数Xを使って、逐次的にW(t_{n+1}) = W(t_{n}) + \sqrt{t_{n+1} - t_{n}} XとWiener過程のサンプリングを行っていく。カオス展開で作った近似Wiener過程を使う場合は、t_{n+1}-t_{n}が小さすぎると、W(t_{n+1}) = W(t_{n})の分布が期待したものから大きく外れてるので、時間幅を小さく取りすぎないように注意する必要があると思われる。

一方、カオス展開では、Ornstein-Uhelenbeck過程は、
X(t) = x_0 e^{-bt} + \dfrac{a}{b}(1-e^{-bt}) + c \displaystyle \sum_{n=1}^{K} \xi_{n} \displaystyle \int_{0}^{t} e_n(u) e^{b(u-t)} du
で近似される。近似Wiener過程のサンプルを作るのに使用したのと同じ乱数を使って、近似Ornstein-Uhlenbeck過程のサンプルを作ると、Euler-丸山法による近似解の比較が出来る。

積分の計算が怠いけど、数値計算には必要なのでサボらずにやれば以下の通り。
X(t) = x_0 e^{-bt} + \dfrac{a}{b}(1-e^{-bt}) + c \xi_{1} \dfrac{1 - e^{-bt}}{b} + \sqrt{2} c \displaystyle \sum_{n=1}^{K-1} \xi_{n+1} \dfrac{(b \cos(\pi n t) + \pi n \sin(\pi n t)) - b e^{-bt}}{b^2+\pi^2n^2}

カオス展開、dt=1/50のEuler-丸山法、dt=1/500のEuler-丸山法の3つで比較したのが以下のグラフ。タイムステップを細かくしたEuler-丸山法と、Kを大きくしたカオス展開の結果は、ほぼ一致している。まぁ、そもそも、Kを有限に取ったら、確率的でない微分方程式を考えるのと、何ら変わらないので、当たり前という気もするけど、とりあえず、積分の計算が間違ってなさそうというくらいは分かる。

f:id:m-a-o:20211114233000p:plain

import math
import random
import numpy as np
import matplotlib.pyplot as plt

class WienerPath:
    def __init__(self, _rnds):
        self.zs = _rnds
    def __call__(self, t):
        assert(t>=0 and t<=1)
        K = len(self.zs)
        fvals = [t]
        fvals.extend([math.sqrt(2)*math.sin(k*math.pi*t)/(k*math.pi) for k in range(1,K)])
        ret = 0
        for (z,c) in zip(self.zs , fvals):
           ret += z*c
        return ret


#-- Ornstein-Uhlenbeck process by chaos expansion
class XPath:
    def __init__(self , _rnds , a=1.05 , b=0.7 , c=0.5 , y0=0):
        self.zs = _rnds
        self.a = a
        self.b = b
        self.c = c
        self.y0 = y0
    def __call__(self , t):
        assert(t>=0 and t<=1)
        K = len(self.zs)
        a,b,c=self.a,self.b,self.c
        fvals =  [c*(1-math.exp(-b*t))/b]
        ret = (self.y0)*math.exp(-b*t) + a*(1-math.exp(-b*t))/b
        fvals.extend([math.sqrt(2)*c*(b*math.cos(math.pi*n*t) + math.pi*n*math.sin(math.pi*n*t) -b *math.exp(-b*t))/(b*b+(math.pi*n)**2) for n in range(1,K)])
        for (z,c) in zip(self.zs , fvals):
           ret += z*c
        return ret


#-- Ornstein-Uhlenbeck process by Euler-Maruyama method
class YPath:
    def __init__(self , _W , a=1.05 ,b = 0.7 , c=0.5 , y0=0):
        self.W = _W
        self.a = a
        self.b = b
        self.c = c
        self.y0 = y0
    def trace(self , ts):
        t_n = 0
        y = self.y0
        ys = []
        a,b,c=self.a,self.b,self.c
        for t in ts:
           if t <= t_n:
              ys.append( y )
           elif t>1:
              pass
           else:
              y_next = y + (a-b*y)*(t - t_n) + c*(W(t) - W(t_n))
              t_n = t
              ys.append(y_next)
              y = y_next
        return np.array(ys)



if __name__=="__main__":
   random.seed(1934)   #-- for experiment
   W = WienerPath([random.gauss(0,1) for _ in range(5000)])
   X = XPath(W.zs)
   X50 = XPath(W.zs[:50])
   X500 = XPath(W.zs[:500])
   Y = YPath(W)
   ts = np.linspace(0,1,num=51)
   ts2 = np.linspace(0,1,num=501)
   ts3 = np.linspace(0,1,num=2001)
   ys = Y.trace(ts)
   plt.plot(ts3,np.vectorize(X)(ts3), label="chaos expansion (K=5000)" , linewidth=3 , alpha=0.4, color="blue")
   plt.plot(ts3,np.vectorize(X500)(ts3), label="chaos expansion (K=500)" , linewidth=3, color="pink" , alpha=0.8, linestyle="dashed")
   plt.plot(ts3,np.vectorize(X50)(ts3), label="chaos expansion (K=50)" , linewidth=3, color="black" , linestyle="dashed")
   plt.plot(ts,ys , label="Euler-Maruyama(dt=1/50)" , linestyle="dashed" , linewidth=3 , color="red")
   plt.plot(ts2,Y.trace(ts2), label="Euler-Maruyama(dt=1/500)" , linestyle="dotted" , linewidth=4  , color="grey")
   plt.legend()
   plt.show()

続・組合わせ範疇文法の漸進的構文解析

組合せ範疇文法の漸進的構文解析
https://m-a-o.hatenablog.com/entry/2021/05/05/085810
の補足。


type raisingとGeach ruleを拡張して
A|B => (X/(X\A))|B
A|B => (X\(X/A))|B
(A/B)|C => ( (A/T)/(B/T) )|C
(A\B)|C => ( (A\T)\(B\T) )|C
のような規則を考えることにする。ここで、記号"|"は、"/"か"\"のどちらか(どちらでもいい)。一般には
( (...(A|B1)|..)|B_n) => ( (...( (X/(X\A)) |B1)|...)|B_n)
のような形も許すべきだと思う。これらを、generalized type raising/generalized Geach ruleと(勝手に)名付けておく。これらの規則は、free pregroupにmapしても成立することは、容易に確認できる。semanticsも問題ないだろうから省略


以下は、人工的に作った例だけど、incrementalに構文解析することを考えなければ、

D/(C/(C\A))    A/B      B 
              --------------
                    A
              ---------------(type raising)
               (C/(C\A))
-----------------------------
            D

は、何も問題がない。しかし、これをincrementalにparseするには、上記のような規則を許すしかないように思える。generalized type raisingを使うと、以下のようになる。

D/(C/(C\A))           A/B                                 B
                  --------------(g-type raising)
                   (C/(C\A))/B   
--------------------------------
            D/B
---------------------------------------------------------------
                            D

generalized Geach ruleも必要である例を示す。前回、argument swapを導入したが、単純なargument swapは、以下のようにtype raisingによって実現できる。

  B                             (A\B)/C
---------(type raising)
X/(X\B)
---------------------------------------- (X:=A)
               A/C

しかし、以下のような形にすると、generalized Geach ruleが必要になる。

   B\D                                  (A\B)/C
----------- (type raising)         -----------------(g-Geach)
X/(X\(B\D))                          ((A\T)\(B\T))/C
-----------------------------------------------------(X:=A\D , T:=D)
                      (A\D)/C

この導出は、free pregrouprでも成立するものである。実際、冒頭に、Dを補って、以下が成立することに注意。

D     B\D          (A\B)/C
----------
    B
----------
 X/(X\B)
--------------------------(X:=A)
           A/C

主語と目的語を取る動詞の統語範疇は、(S\NP)/NPだったか、(S/NP)\NPだったかと一瞬でも悩む必要はない。どっちでも正しく構文解析できる。



argument swapで、
(A/B)/C => (A/C)/B
という形のものは、free pregroupでは成立しない規則であるが、同様にfree pregrouprで成立しない(が、CCGでは通常仮定される)crossing compositionがあれば実現できる。

(A/B)/C                B                             C
                 -------------(type raising)
                    A\(A/B)
------------------------------(xB<)
          A/C
---------------------------------------------------------
                        A


以下の例は、type raisingにgeneralized type raisingを重ねるなど、結構難しい。

       B                                        A\B                                            C\A                (C\B)\(C\B) 
------------------(type raising)  -----------------------------------(type raising)      -----------------(Geach)
    C/(C\B)                               (C\B)/((C\B)\(A\B))                                (C\B)\(A\B)
------------------(Geach)          ----------------------------------(g-type raising)
  (C/T)/((C\B)/T)                     (X/(X\(C\B)))/((C\B)\(A\B))
---------------------------------------------------------------------(T:=(C\B)\(C\B) , X:=C\B)
                  (C/((C\B)\(C\B)))/((C\B)\(A\B))
-------------------------------------------------------------------------------------------------------------
                              C/((C\B)\(C\B))
-------------------------------------------------------------------------------------------------------------------------------
                                                                     C

前回、incrementalにparseするのが難しい例として挙げられていた英文は、generalized type raisingとgeneralized Geach ruleを使って、以下のようにincrementalに導出木を作れる。

the   woman            that                          every                                    man                 saw
------------     ---------------           ---------------------------                     --------            ----------
    NP            (NP\NP)/(S/NP)                      NP/N                                    N                 (S\NP)/NP
------------
NP/(NP\NP)      
----------------------------------         ----------------------------(g-type raising)
       NP/(S/NP)                                 (S/(S\NP))/N
                                           ----------------------------(g-Geach)
                                               ((S/X)/((S\NP)/X))/N
----------------------------------(Geach)  ----------------------------(X:=NP)
      (NP/T)/((S/NP)/T)                        ((S/NP)/((S\NP)/NP))/N
-----------------------------------------------------------------------(T:=((S\NP)/NP))
                 (NP/((S\NP)/NP))/N
-----------------------------------------------------------------------------------------------------          --------------
                                  NP/((S\NP)/NP)                                                                  (S\NP)/NP
------------------------------------------------------------------------------------------------------------------------------
                                                            NP


Geach ruleには、次のようなcrossing型のものも考えられる。
A/B => (A\T)/(B\T)
A\B => (A/T)\(B/T)
これらのcrossing Geach rule(と例によって勝手に命名)を許すと、次のようにcrossing compositionが成立する。

  A/B                         B\C
--------------(crossing Geach)
  (A\T)/(B\T)
------------------------------------(T:=C)
              A\C
     B/C                   A\B
                     ----------------(crossing Geach)
                        (A/T)\(B/T)
--------------------------------------(T:=C)
                  A/C

crossing compositionは、通常のCCGで仮定される規則だけど、crossing Geach ruleの有無によっても制御できる。


application/composition/type raisingのみが許されるCCGの任意の導出木を、generalized type raisingとgeneralized Geach ruleを使って、(semanticsを変えることなく)incrementalな導出木に変形できることが示せると思う。基本的には、何か3つの統語範疇の並びがあって、後ろの2つが先に評価される部分を、前の2つが先に評価されるように変形していけばいい。上に挙げた例で複雑なものは、このプロセスを手で追って作った(これは、完全に機械的な作業で、この変形を行うアルゴリズムになってるはず)

更に、crosing compositionw持つ場合も、全く同様の方法で、対応するgeneralized crossing Geach ruleの追加によって、任意の導出木をincrementalな導出木に変形できる。但し、generalized crossing Geach ruleは、pregroup grammarからは逸脱した規則なので、generalized crossing Geach ruleの追加によって、文法の表現力が拡大しないことは何か別の手段がないと保証できない。


CCGには、他にも、色々な規則が仮定される場合があり、incrementalにparseするのを邪魔する可能性がある。例えば、substitutionと呼ばれる規則
f : (A/B)/C g : B/C => λc. f(c , g(c)) : A/C
は、変数cが非線形に出現するので、free pregroupでは成立しない。このような規則を追加したCCGをincrementalにparseできるのか、あるいは、incrementalにparseするのに、どのような規則を追加すればいいのか、私には予想が付かない。substitutionが、自然言語構文解析で、必要になる例があるのかは分からないので、とりあえず、気にしないことにする。


実装は次回。

Daniel Bernoulliの発想を継承した二項分布の正規近似の導出

Daniel Bernoulliの正規曲線の方程式の導出法
https://doi.org/10.34336/jhsj.25.158_89

に、Daniel Bernoulliが1770年頃に発表したらしい、二項分布の正規近似の導出法が書かれてて、ちょっと面白い議論だった。



二項分布の正規近似は、de Moivreの1733年の論文に起源があるとされる。de Moivreは、本質的には、現代の殆どの教科書で使われてるのと同様、Stirlingの公式を使った議論を行ったのだろう。このStirlingの公式も、発見したのは、de Moivreとされている。

Stirlingの公式に関するde Moivreの議論は、以下の論説に詳しい。

The Early History of the Factorial Function
https://doi.org/10.1007/BF00389433

de Moivreの元々の動機からして、二項係数に端を発していて、その途中で、de Moivreは、
\log (m-1)! = (m-\dfrac{1}{2})\log m - m + C + \dfrac{B_2}{1 \cdot 2} \dfrac{1}{m} + \dfrac{B_4}{3 \cdot 4} \dfrac{1}{m^3} + \dfrac{B_6}{5 \cdot 6} \dfrac{1}{m^5} + \cdots
C = 1 - \dfrac{B_2}{1\cdot2} - \dfrac{B_4}{3\cdot 4} - \dfrac{B_6}{5 \cdot 6} - \cdots
という形の式を証明したらしい。B_nはベルヌーイ数。この定数C\log(2 \pi)/2であることを指摘したのが、Stirlingだったそうだ。この公式が、Stirlingの公式と呼ばれるようになった経緯は知らない。

この論説の中にあるDaniel Bernoulliから、Goldbachへの1729年の手紙というのは

File:DanielBernoulliLetterToGoldbach-1729-10-06.jpg
https://commons.wikimedia.org/wiki/File%3ADanielBernoulliLetterToGoldbach-1729-10-06.jpg

で読める。簡潔で証明もないけど、ガンマ関数の無限積表示が与えられているっぽい。これはGoldbachが出した問題に対する解答っぽいけど、同じ問題に対するEulerの解答が、
De progressionibus transcendentibus seu quarum termini generales algebraice dari nequeunt
https://scholarlycommons.pacific.edu/euler-works/19/
にある(Euler archiveのE19)。この報告では、ガンマ関数の積分表示(9ページ目。式の形は、現在の標準的なものとは違い、ちょっとした変数変換がいる。あと、英訳版は、lnが誤植っぽい?)が見られる。何にせよ、彼らは、確率論とは異なる文脈で、ガンマ関数に到達したらしい。



正規近似に関するD.Bernoulliの論文は、

Mensura sortis ad fortuitam successionem rerum naturaliter contingentium applicata (1770)

Continuatio argumenti de mensura sortis ad fortuitam successionem rerum naturaliter contingentium applicata (1771)

と2本あったようだが、電子化された原文を見つけられなかった。


D.Bernoulliによる正規近似の導出は、冒頭の論文に書かれてるけど、少し記号を変えて再構成すると以下のような議論。

正数N自然数k \le N、0〜1の実数pに対して、
q(k) = \dfrac{N!}{k!(N-k)!} p^k(1-p)^{N-k}
とすると、
 q( k+1 ) - q( k ) = \displaystyle \frac{N-k}{k+1} \displaystyle \frac{1-p}{p} q(k) - q(k)
が成り立つ。M=pNk=M+x,f(x)=q(k)=q(M+x)と記号を導入すると、
f(x+1) - f(x) = -\dfrac{\frac{1}{1-p} x + 1}{M+x+1} f(x)
と変形できる。f(x)を実数上定義された微分可能な関数であれば
f(x+1)-f(x) = \displaystyle \int_{x}^{x+1} f'(t)dt
なので、f'(t)が、(x,x+1)で、ほぼ一定であれば、微分方程式
f'(x) = -\dfrac{\frac{1}{1-p} x + 1}{M+x+1} f(x)
を得る。これは初期条件f(0)=Qの元で解くことができる。


関数方程式を微分方程式で近似するタイプの議論は、他にもいくつか存在する。例えば、
F(x+1)-F(x) = x^k
の解は、ベルヌーイ多項式だけど、
F'(x) \sim x^k
と置き換えると
F(x) \sim \dfrac{x^{k+1}}{k+1} + C
と解けて、ベルヌーイ多項式の主要項が得られる。

もうひとつの例として
F(x+1)= x F(x)
なら、対数を取って
\log F(x+1) - \log F(x) = \log(x)
なので、同様の近似によって
\log F(x) \sim x \log(x) - x
という形で、de Moivre-Stirlingの公式の主要項が得られる。

そうは言っても、この議論は、ちょっと荒っぽい。後者の議論に関連して、大学一年レベルの解析学の教科書には
\displaystyle \int_{0}^{N} \log(x) dx \le \log(N!) \le \displaystyle \int_{1}^{N+1} dx \log(x)
という不等式が書いてある。正規近似の導出でも、同じようなことができそうか考える。


関数方程式に立ち戻ると、解きたいのは
f(x+1) - f(x) = -\dfrac{\frac{1}{1-p} x + 1}{M+x+1} f(x)
だけど
g(x) = 1 - \dfrac{\frac{1}{1-p} x + 1}{M+x+1} = \dfrac{M- \frac{p}{1-p}x}{M+x+1}
と置くと
f(x+1)=g(x) f(x)
という形になる。

x \lt \dfrac{1-p}{p} M = (1-p)Nの時、g(x) \gt 0なので対数を取って、
\log f(x+1) - \log f(x) = \log g(x)
になる。x \gt -M-1=-(pN+1)で、g(x)は、単調減少関数なので
\log g(k+1) \le \displaystyle \int_{k}^{k+1} \log g(x) dx \le \log g(k)
が成り立つ。従って
\displaystyle \sum_{k=1}^n \log g(k) \le \displaystyle \int_{0}^{n} \log g(x) dx \le \displaystyle \sum_{k=0}^{n-1} \log g(k)
なので
\log g(n) + \displaystyle \int_{0}^{n} \log g(x) dx \le \displaystyle \sum_{k=0}^{n} \log g(k) \le \log g(0) + \displaystyle \int_{0}^{n} \log g(x) dx
となる。

更に
\displaystyle \int_{n}^{n+1} \log g(x) dx \lt \log g(n)
\log g(0) \lt 0
を使うと、
\displaystyle \int_{0}^{n+1} \log g(x) dx \le \displaystyle \sum_{k=0}^{n} \log g(k) \le \displaystyle \int_{0}^{n} \log g(x) dx
が示せる。

上側の評価は、似たような議論で、もう少し改善できて、
\sum_{k=0}^{n} \log g(k) \le \displaystyle \int_{-1}^{n} \log g(x) dx \lt \displaystyle \int_{0}^{n} \log g(x) dx
を示すことが出来る。

x \lt 0の時も、同じような議論ができるが省略。以下、x \ge 0で考える。xは、実数値を取るとしてるが、例えば、f(x)が、x \gt 0で単調減少すると仮定すると

\log f(x) - \log f(0) = \log f(x) - \log f(x-[x]) + \log f(x-[x]) - \log f(0)

と分解([x]は、xを超えない最大の整数)して、|\log f(x-[x]) - \log f(0)| \lt | \log f(1) - \log f(0)|は、Nが大きくなると、0に近付く。ちゃんと計算すると、(f(x)が、x \gt 0で単調減少するという仮定のもとで)xが自然数の時と同じように、xが任意の実数の時でも同様の不等式の成立を示すことができる。


数値例を一つ見ておく。ベータ関数B(z,w)を使うと、関数方程式の一つの厳密解を書くことができ、これに対して、
E(x) = \log f(x) - \log f(0) = x \log(\dfrac{1-p}{p}) - \log B(M+x+1,\dfrac{1-p}{p}M - x + 1) + \log B(M+1,\dfrac{1-p}{p}M+1)
という関数を定義する。ここで考えた近似解に対して同様の関数を計算すると、
M(x) = \displaystyle \int_{0}^{x} \log g(x) dx = \dfrac{(w x -M) \log(\frac{M-wx}{M+x+1}) - (Mw+M+w)\log(\frac{M+x+1}{M+1})}{w} + \dfrac{M}{w} \log(\dfrac{M}{M+1})
w=\dfrac{p}{1-p}
を得る。また、ベルヌーイによる微分方程式の解からは
C(x) = -(1+w)x + (Mw+M+w)\log( \dfrac{M+x+1}{M+1} )
が得られる。

これらの関数をN=106,p=0.23の時、プロットすると、以下のようになる。

f:id:m-a-o:20211114225129p:plain

厳密解は、M(x)M(x-1)-M(-1)のほぼ中間を通り、\dfrac{M(x)+M(x-1)-M(-1)}{2}は、Nが比較的小さい時でも、厳密解とよく一致する。微分方程式解から得られるC(x)は、xが0の近くにある時には、M(x)と近い値を取るが、xが大きくなると、乖離が目立つ。

正規近似を得るには、M(x)をTaylor展開して、二次の項まで取ればいい。大体の計算は、素直に展開してけばいい。
x \log \left(\dfrac{M-wx}{M+x+1} \right) = x \log (\dfrac{M-wx}{M+1}) - x \log(\dfrac{M+x+1}{M+1})
という分解を使うといい。このへんの議論は、冒頭の論説と本質的には変わらないので省略。


グラフ作成コード

import math
from scipy.special import betaln
import numpy as np
import matplotlib.pyplot as plt

N,p=106,0.23
M=p*N

def exact(x):
   return x * math.log(p/(1-p)) - betaln(M+x+1 , (1-p)*M/p - x + 1) + 0*(M*math.log(p) + ((1-p)*M/p) * math.log(1-p) - math.log(N+1)) + betaln(M+1 , (1-p)*M/p+1)


def approx(x):
    w = p/(1-p)
    C = (M*math.log(M/(M+1)) + (M*w+M+w)*math.log(M+1))/w
    return ((w*x-M)*math.log((M-w*x)/(M+x+1)) - (M*w+M+w)*math.log(M+x+1))/w + C


def C(x):
  w = p/(1-p)
  return -(1+w)*x + (M*w+M+w)*math.log((M+x+1)/(M+1))


if __name__=="__main__":
   x = (np.arange(4000)/100 - 20)
   y1 = np.vectorize(exact)(x)
   y2 = np.vectorize(approx)(x)
   y3 = np.vectorize(approx)(x-1)-approx(-1)
   y4 = np.vectorize(C)(x)
   y5 = (y2+y3)/2
   plt.plot(x,y1 ,label="E(x)" , linewidth=3 ,alpha=0.6)
   plt.plot(x,y2 , label="M(x)",linestyle="dashed")
   plt.plot(x,y3 , label="M(x-1)-M(-1)",linestyle="dashed")
   plt.plot(x,y4 , alpha=0.8 , label="C(x)" , linewidth=5 , linestyle="dotted")
   plt.plot(x,y5 , alpha=1 , label="(M(x)+M(x-1)-M(-1))/2" , linewidth=5 , linestyle="dotted")
   plt.legend()
   plt.show()


パンがなければ合成すればいいじゃない

食糧合成は、食糧を農業によらず化学的に合成する技術を指す言葉で、1960年前後には、割と真剣に検討されていたっぽい。1950年代、1960年代に、日本の科学雑誌で、食糧合成を取り上げた記事を、いくつか見つけることができるし、日本国外にも、この手の検討を行った人はいる。以下は、その例。

(1955年)食料合成の可能性 https://doi.org/10.1295/kobunshi.4.292

(1956年)必須アミノ酸の合成 https://doi.org/10.5059/yukigoseikyokaishi.14.367

(1960年)Chemistry, Food, and Civilization https://www.jstor.org/stable/24534794

(1966年)合成食糧と2人の先達 https://doi.org/10.1271/kagakutoseibutsu1962.4.546

(1966年)Chemical and Biochemical Production of Food for Man and Animal https://doi.org/10.2527/jas1966.252575x

(1967年)食糧合成の方向 https://doi.org/10.1271/kagakutoseibutsu1962.5.534

(1967年)これからの合成食料(ベランダ) https://doi.org/10.14894/faruawpsj.3.9_639

農業で扶養できる人口の見積もり

人口増加に伴う食糧問題は、19世紀にも、色んな人が言及したテーマだった。

21世紀初頭の現在は、 (1)地球上の農業による食料生産で原理的に扶養可能な人口の上限が見積もれるようになった (2)食糧を化学的に合成するための知識と要素技術が、ある程度揃っている という点で昔とは違う。

(1)の人口の見積もりは、1925年のAlbrecht Penckによるもの以降、何度も行われてきた。現在は、Penckの時代より、光合成の理解も進んだので、もう少し、原理的な上限を考えることができる。

農業にしろ、狩猟採集にしろ、そのエネルギー源は、元を正せば、ほぼ太陽光に由来すると言っていい。狩猟採集では、生えている植物の大部分は、食用には適していない。農業の場合は、雑草などが生えてくるものの、農地に育った植物のかなりの部分を食べることができる。従って、農業は、狩猟採集より遥かに面積あたりの生産効率が高い。現在の地球人口で、大多数の人が、昔ながらの狩猟採集生活に戻ろうと判断した場合、食糧供給は破綻する。しかし、当然、農業でも扶養できる人口は無限ではない。

光合成のエネルギー変換効率から、単位面積当たり年間に収穫可能なエネルギー総量を計算でき、それに基づいて、単位面積あたりで扶養可能な人口の上限を試算できる。人間の食事には、様々な必須栄養素の確保といった目的もあり、エネルギー以外の制約条件もあるが、扶養可能人口は、これより小さい可能性はある。

Wikipediaには、(光合成によって)年間に地球上で貯蔵されるエネルギーは、1.0e18(kJ)と見積もられていると書いてある。出典は、ヴォートの『生化学』となっている。80億人の人間が、一日2500kcalを摂取する場合、年間3.05e16(kJ)が必要になる。単純計算で、現在、人類は、光合成で貯蔵されるエネルギーの1/30くらいを食糧として消費していることになる。これは、食糧廃棄や食糧以外での植物利用(木材や薪炭など)は考慮していない数値。光合成で貯蔵されるエネルギーは、他の生物も利用するものなので、その1/10くらいを人間が利用できるとすれば、扶養可能人口は240億人。

1.0e18(kJ)を陸地面積と一年間の秒数で割ると、大体、0.2(W/m2)程度で、太陽定数1370(W/m2)の1/1000にも満たない。勿論、砂漠や氷雪地帯では、殆ど光合成が行われないので、光合成の効率は、もう少しは大きいだろう。1.0e18(kJ)は、多分オーダーなので、数倍程度の食い違いは生じうる。もう少し、微視的観点からの試算として、

光合成と地球環境
http://hdl.handle.net/2261/53196

によると、日本の場合、単位面積当たり年間に収穫可能なエネルギー総量は、最大で、1.45(W/m2)≒1.0e8(kcal/ha/year)と見積もられている(表4)。これは、光合成のエネルギー変換効率から想定される数値で、割と本当の上限に近いものと思われる。これに、一期作の場合の米作りの期間が、一年の約半分(5ヶ月としている)であること、可食部は収穫物の一部のみである(約50%としている)ことを考慮すると、より小さな値となる。

表6(1992年の穀物生産高)で、米の数値は、おそらく、籾高(英語圏で、paddy riceと書いてある場合は、籾高を指してるっぽい)の数値であって、日本の統計では玄米収量、世界の統計では精米収量(milled rice)が用いられることが多いが、いずれにせよ、可食部重量は籾高より少なくなる(収量は、籾高>玄米>精米であり、精米収量は、籾高の6割程度になるっぽい)。ただ、オーダー計算としては致命的な誤りではない。

現在、世界の農地面積は、約14億haで、世界の陸地面積の1/10を占めていると見積もられている。世界全体の農地増加ペースは、数十年前から、かなり緩やかである。仮に、農地面積がこれ以上増えず、農業技術が物凄く向上して、世界平均で、農業のエネルギー利用効率が、0.5(W/m2)を達成して、かつ、100%が人間の食糧になっても(つまり、家畜飼料や廃棄は一切ないという理想的状況で)、一人一日消費カロリー2500kcal換算で養える人口は、500億人である。

Penckによる見積もりだと、上限は、約160億人とされていたそうだ(Penckが農地面積として、どの程度を想定していたかは知らない)。例えば、農業のエネルギー利用効率が年間平均0.25(W/m2)で、家畜飼料や廃棄によるロスが50%程度あると仮定すれば、農地拡大なしに扶養可能な人口は、125億人程度にまで減る。

現時点で、農作物には、備蓄に回される分、家畜の飼料になる分や廃棄など、非食糧用途の利用が、結構あると思われる。どれくらいあるか目安をつけるために、現在の世界人口を約80億人として、主要穀物の生産量(2020年)から、作物ごとに賄える一人一日あたりのカロリーを計算すると、以下の表のようになる。

生産量(Mt) カロリー換算係数(kcal/g) 一人一日カロリー(kcal/人・日)
とうもろこし 1125 0.92 354
小麦 775.8 3.5 930
精米 505 3.5 605
大麦 159.74 3.5 191
じゃがいも(2017年) 388.2 0.76 101

Mtは、メトリックトンじゃなく、メガトン。合計すると、2000kcalを超えるので、非食糧用途の消費が結構多いと思われる。一人一日の総摂取カロリーに占める穀物の割合は、国によって違うけど、800~1200caklくらいと見積もられているらしく、生産された穀物の半分以上は、直接、食糧にはなっていないということになる。穀物に関する限り、直接、食糧になるのが、生産量の50%という見積もりは、現在のところ悪いものではないと思われる。

20世紀には、色々な方法(肥料、農薬、品種改良など)で、面積あたりの農業の生産性を向上させてきたが、上文献の見積もりによれば、いくつかの地域と作物では、おそらく、上限に近い水準に達している。このことを確認する。

農林水産省が公開している作物統計によれば、2019年の水稲作付面積は1469000(ha)で、収穫量(1.70㎜のふるい目幅で選別された玄米の重量)は、7762000(t)で、反収は5.28(t/ha)となる。玄米の重量あたりカロリーを、3.5(kcal/g)とすれば、年間平均で、5280e33.5e34.185/(365864001.0e4)=0.245(W/m2)だったということになる。あるいは、米作りの期間が5ヶ月であるとすれば、この期間のエネルギー利用効率は、0.6(W/m2)に近い。可食部が50%とすれば、非可食部も含めたエネルギー利用効率は、1.2W/m2ということになるので、これは、上文献の限界値に近い。

一応、米の生産量について、他のデータと整合性をチェックしておく。7762000(t)を、(全部国内で消費し、かつ外国からの輸入を考慮しない場合)、人口と年間日数365で割ると、一人一日当たり177g程度の割り当てになる。精米にすると、もう少し減るけど、炊くと2倍ちょっとの重量になるそうなので、一日で茶碗2〜3杯程度の米を食べてる感じだろう。

また、玄米の相対取引価格(令和二年、米に関するマンスリーレポート)は、60kgあたり15000円くらいなので、収穫量7762000(t)は、(全部販売したとすれば)2兆円弱の売上になる。検索すると、米農家の戸数は、約93万9千戸(2015年)とあり、ざっくり100万戸として、一戸あたりでは、200万円程度の売上という計算になり、ここから諸経費等を引いたものが収入だろう。現代日本世帯年収200万円を下回るのは厳しいが、米農家は兼業農家が多いそうで、こんなものかもしれない。

他の作物についても見ておく。出典は、FAOSTATのCrops and live stock produtcsで、2019年の日本のデータ。

作物 英語名 反収(kg/ha) カロリー(kcal/g) エネルギー効率(W/m2)
じゃがいも Potatos 29869 0.76 0.30
小麦 Wheat 4900 3.5 0.23
大麦 Barley 3621 3.5 0.17
人参・カブ Carrots and turnips 33671 0.37 0.17
りんご Apples 19489 0.54 0.14
かぼちゃ Pumpkins, squash and gourds 11123 0.91 0.13
なす Eggplants 34879 0.22 0.10
キャベツ等 Cabbages and ... 41696 0.23 0.13
ぶどう Grapes 10404 0.59 0.081

主食になってる穀物類は、効率が高い傾向にあるが、概ね、近い水準にある。

農地拡大なしでも、他の条件次第で、扶養可能人口は結構幅はある。農地面積は増やせるかもしれないが、収穫逓減に対抗する技術が必要になるだろう。20世紀にも、かつては、農業に不向きだった土地を農地化したケースはある。例えば、鹿児島は、昔は、稲作に向いておらず、江戸時代には、サツマイモが栽培されるようになったが、20世紀半ば頃、灌漑施設が整備され、現在では鹿児島県の水稲の生産性は日本の他地域と大差ない。現在、砂漠化している陸地面積は、約36億haと言われてるので、緑化できれば、大幅な農地の増大はできるかもしれない。

今の所、2100年の人口予測は100億人前後となっているし、最近は、出生率の低下を心配してる国が増えてきた。たとえ農地がそれほど拡大しなくても、食糧不足について、差し当たって大きく心配することはないようにも思える。

このまま農業に依存した食料生産を続けても、当面(現在生きてる人が全員死ぬまでの期間程度)問題なさそうだけど、数十年後の人口予測なんか当てにならない。例えば、ありそうなシナリオとしては、老化の治療による死亡率の大幅な低下が考えられる。

それに、地球が寒冷化して農業が大打撃を受ける可能性が今後もないとは言えないし、農業を自動化するより化学プロセスを自動化する方が多分簡単。動植物、微生物関係なく、不殺生を実現できる可能性があるので、食糧合成技術は、真のジャイナ教徒にも優しいかもしれない(他の戒律に触れないかは知らない)。

食糧合成で扶養できる人口の見積もり

農業で利用可能なエネルギー源は太陽光だけど、食糧合成には、そういう制約はない。といっても、石油とかの正確な埋蔵量はわからないし、ここでは太陽光発電を利用した場合の食糧合成を考える。同じ太陽光依存なので、農業との比較も公正な気がする。

Googleが教えてくれる"設置面積あたりの太陽光発電の設置可能容量"の表には、設置面積30坪=99.17(m2)で、年間発電量が11400kWhとある。これは、13.1(W/m2)に相当する。上に挙げた文献「光合成と地球環境」では、日本が存在する中緯度地方で、"地表付近の実効太陽定数(※)"を145(W/m2)としているので、太陽光発電の効率が10%程度なら、妥当な数値だろう。

※)太陽定数は、大気上端で、太陽光線に垂直な面に入射するエネルギー量で、そのうち地表に届く光は、7割程度とされているが、更に、曇りや雨の日があったり、夜もある。そういった諸々の変動を平均化して、地表付近で、利用可能な太陽光エネルギーを考えている

実際の太陽光発電所を見ても、平均出力は、10(W/m2)前後になるようである。太陽光発電所は、定格出力を公開している場合もあるが、これを敷地面積で割っても、得られる数値は欲しいものではない。定格出力ではなく、実際の年間発電量を敷地面積で割るのがより妥当だろう。実際は、敷地をパネルで埋め尽くしているわけでもないだろうが、そこまで考慮しても、実用上の意味は薄いだろう。

株式会社関電エネルギーソリューションの発電所の実績・事例にある数値を抜粋すると、以下のようになっている。

年間発電量 敷地面積 面積出力 運転開始日
けいはんな太陽発電所 250万kWh 4ha 7.1(W/m2) 2013年12月
有田太陽光発電 3100万kWh 45ha 7.9(W/m2) 2015年10月
山崎太陽光発電 280万kWh 4ha 8.0(W/m2) 2016年11月
赤穂西浜太陽光発電 260万kWh 2.5ha 11.9(W/m2) 2018年6月
けいはんな第二太陽光発電 140万kWh 1.0ha 16.0(W/m2) 2018年9月

下に行くほど、運転開始日が新しく、面積あたりの発電量も大きいけど、技術的改善があったのかは定かでない。いずれにしても、単位面積から回収できるエネルギー量という点で見れば、太陽光発電は、農業の10倍以上効率的と言える。

人間は、今の所、直接電気を摂取できないので、人間が利用可能な形に変換する必要があって、その過程で、どれくらい損失が発生するかは何とも言えない。変換効率以外に、元素収支など、他の制約条件も考慮する必要がある。一方で、農地面積の拡大に比べれば、太陽光発電所の拡大には原理的な制約が少ない。

太陽光は、地球に到達しない分は、宇宙の彼方に捨てられてて、これを回収できれば、相当に余裕が生まれる。月面とかで発電して、無線送電するということも原理的にはできる。遠くに発電施設を作ると、メンテナンスが難しいという問題はあるけど、地球外で発電すれば、気象条件や時間帯によって発電量が大きく変動する問題は解消される。

100億人程度の人口を扶養するなら農業でも問題ないと思われるが、1000億人とかになってくると、農業では厳しい。全ての制約を詳細に検討してはいないが、食糧合成なら、それに対応できる可能性がある。

食糧合成小史

19世紀〜20世紀前半

いきなりパンの合成を考えるのは飛躍しすぎなので、最初は、当然、グルコースアミノ酸脂肪酸などの合成を考える。これらの分子が、栄養の基盤であると分かったのも20世紀のことだけど、有機化学の世界では、古く1850年には、ドイツの化学者Adolph Strecker(1822〜1871)によって、アラニンのストレッカー反応が報告されたとされている。また、19世紀末には、ドイツの化学者ヘルマン・エミール・フィッシャーらによって、グルコースやマンノースの化学合成が報告された。

1908年に、池田菊苗は、特許第14805号「グルタミン酸を主要成分とせる調味料製造法」を出願、登録し、1909年には「味の素」の販売が開始されたそうである(グルタミン酸自体は、1866年にドイツのKarl Ritthausenが単離)。この時代は、食品由来のタンパク質を加水分解して、製造していた。

当時は、まだタンパク質を構成するアミノ酸が全て特定されてなかった。ヘルマン・フィッシャーが、ペプチド合成法を確立して、タンパク質は、アミノ酸のポリマーであるという考えが出て、そんなに年月も経ていない。人間にとって、必要な栄養が何かということも、分子レベルで十分理解されていなかったので、食糧を化学的に合成する可能性を想像したとしても、実際に何を合成すべきか明らかではなかった。

1916年に、以下のような論文が出ていて、1910年代には、(主に動物実験によって)いくつかのアミノ酸が、成長や生存に必須であるという考え方は存在していたらしい。

cf) THE AMINO-ACID MINIMUM FOR MAINTENANCE AND GROWTH, AS EXEMPLIFIED BY FURTHER EXPERIMENTS WITH LYSINE AND TRYPTOPHANE87509-3)

cf) Experiments That Changed Nutritional Thinking

タンパク質を構成するアミノ酸20種が全部単離されたのは、1936年にスレオニン(アメリカのWilliam Roseによる)が発見された時で、スレオニン必須アミノ酸の一つである。1930年代には、必須脂肪酸があるということも、徐々に認められつつあった。以下の文献を見る限り、少なくとも、1947年には、タンパク質を構成するアミノ酸が20種類であること、必須アミノ酸が9種類あることなどは、確立してたようだ。

The Amino Acid Requirements of Man
https://doi.org/10.1016/S0065-3233(08)60081-9

また、1930年代には、ビタミンCの化学合成(1933年、Tadeus Reichstein)、ビタミンB1の化学合成(1935年、Robert R. Williams)が報告され、Hermann O.L. FischerとBaerは、1936年に、ジヒドロキシアセトンとグリセルアルデヒドの2つのC3糖から、フルクトース(とソルボース)の合成を報告した。

cf)Synthese von d‐Fructose und d‐Sorbose aus d‐Glycerinaldehyd, bzw. aus d‐Glycerinaldehyd und Dioxyaceton; über Aceton‐glycerinaldehyd III

このグルコースの合成は、生体内の糖新生と解糖系の反応、フルクトース-1,6ビスリン酸<->グリセルアルデヒド-3-リン酸+ジヒドロキシアセトンリン酸に似ていて、解糖系のエムデン・マイヤーホフ経路に名前が残っているGustav Emdenは、1932年に、解糖系の仮説的な過程を提案したそうだから、Fischer&Baerの研究も、それに触発されたのかもしれない。

どうでもいいけど、有機化学には、色んなフィッシャーさんがいる(英語で、フィッシャーというと、Fisherだと思うけど、ドイツでは、みんなFischerとなるっぽい?)。1902年にノーベル化学賞を受賞したヘルマン・エミール・フィッシャー(1852〜1919)は、一番有名なフィッシャー。次に、フィッシャー・トロプシュ法で知られるフィッシャーは、フランツ・フィッシャー(1877〜1947)。そして、グルコース合成を行った(恐らく一番無名の)このFischerは、また別の人っぽい。

1930〜40年代は、栄養の化学的基礎が概ね明らかになった時代と言える。

20世紀後半

食糧合成の検討が本格的に始まったのは1950年代になってからのようである。最初に挙げた「合成食糧と2人の先達」という文献では、

私たちは赤堀先生を中心として昭和27年,合成食糧の総合研究班を結成し,文部省より機関研究費を得て,まず従来の文献をたずねて1,500件の関係文献を集め,食品となり得るあらゆるアミノ酸,有機酸,アルコール類の化学的合成法と,微生物による前記諸化合物の生合成法を調べ,合成食糧研究資料を編纂し,さらにいっそう有効な新合成法の研究に努力したのである.それらの結果,赤堀先生,村上増雄教授らによって新しい各種アミノ酸の合成法が発見され,それらの新旧合成法のいくつかは直ちに工業生産に移され,今日のアミノ酸合成工業に発達したのである.

と、1952年(昭和27年)に食糧合成研究を開始したという記述がある。

アミノ酸発酵技術の系統化調査(PDF)によると、味の素では、1950年代から、アミノ酸の化学合成法の研究を開始したそうで、その後、以下のようにある。

アクリロニトリル一酸化炭素とシアン化水素を結合させるストレッカー反応を応用した方法を完成し、同社の東海工場で工業生産が行われた。この技術は、当時は無尽蔵と言われた石油製品から、食品であるグルタミン酸を作ると言う画期的な技術として称賛され、日本化学会技術賞(1964)、大河内記念生産賞(1965)を受賞した。 しかしながら、短時間の内に消費者の意識は大きく変わり、「食品を化学合成法で作る」というコンセプトは受け入れられない時代になった。ほぼ同時にスタートした発酵法の進展もあって、化学合成法によるMSGの合成工場は1973年に閉鎖された。

アミノ酸ビジネスと技術開発(味の素グループの100年史、第7章2節)には、次のように書いてある。

補足するなら、合成法のあらゆるメリットをもってしても、どうにもならない現実として残ったのは、6番目の環境要因にほかならない。1960年代前半には高い評価を得ていた技術が、わずか10年余りのうちに一転して受容されにくい技術になった。その背景には、化学に対する消費者の劇的な意識変化があった。化学合成技術を最先端技術としてもてはやす時代から、公害の元凶として排撃する時代に変わったのであれば、食品を化学合成で作るというコンセプトは受け入れられない

日本国内で四大公害裁判が始まったのが1967年あたりで、そのへんの出来事を指してるのだろう。

ちなみに、現在は、遺伝子組み換え細菌の利用によって生産効率の向上が図られているっぽい。2009年に開催された「第308回食品安全委員会会合」の議事録

第308食品安全委員会議事概要
http://www.fsc.go.jp/iinkai/i-dai308/dai308kai-gijigaiyou.html

に"遺伝子組換え食品等「GLU−No.2株を利用して生産されたL−グルタミン酸ナトリウム」に係る食品健康影響評価について"という記述が見られ、これは、「味の素」が開発・申請したものらしい。ベースになってるのは、Corynebacterium glutamicumという細菌っぽい。生物種は、安全性審査の手続を経た遺伝子組換え食品及び添加物一覧 (平成25年10月4日)の表を見ると分かる。

元々、この細菌は、協和醱酵工業株式会社の研究者によって、1956年に発見されたものっぽい(論文は1957年)。これは、上の引用文で、「発酵法」と書かれている方法。

Studies on the amino acid fermentation. Production of L-glutamic acid by various microorganism
https://doi.org/10.2323/jgam.3.193

GLU-No.2株というのは、遺伝子組み換え技術によって、生産効率を増大させたものなんだろう。

他に、

アミノ酸の製造について(1962年)
https://doi.org/10.5059/yukigoseikyokaishi.20.676

には、

また構造の簡単なアミノ酸であるグリシン ・DL-アラニンは合成によって製造され、タンパク質加水分解物から分離されるロイシンを主とした中性アミノ酸混合物と共に,合成酒用の用途がひらかれた。戦時中医薬用途の開発されたメチオニンはDL-態で合成利用され,さらに必須アミノ酸を主としたアミノ酸の混合水溶液に輸血代用の効果が認められて,1955年頃より市販が開始された。

のように書かれていて、1950年代には、いくつかのアミノ酸の化学合成が工業的に行われていたようである。メチオニンは、現在でも、家畜飼料用途としての需要が高いらしく、2019年の住友化学のプレスリリースでも、生産体制の強化を行う旨が書かれている。

飼料添加物メチオニン事業の競争力強化について
https://www.sumitomo-chem.co.jp/news/detail/20191001.html

コスト的な観点からすると、化学合成の場合は、炭素源として原油を利用することになるので、原油価格が高くなれば、化学合成は不利になる。化学合成は、原料価格以外に、精製コストの高さで、発酵法より不利になりがちなんだろうけど、実際の精製コストが、どんなものなのか詳細は知らない。

発酵法の場合、炭素源となるのは、遡れば、農作物に由来するだろう。酵母エキスとか使ってたとしても、酵母の培養には、食品から抽出した成分を含む培養液(サトウキビとか使うらしい)を使うものと思われる。なので、現時点で発酵法が安価だとしても、普遍的な話ではなく、食品価格が高騰することがあれば、化学合成法の方が低コストになるということは考えうる。発酵法の場合、人間が食べられない植物を利用する可能性はあるが、大規模に行われてる事例が実際にあるのかは不明。

別に、発酵法と化学合成法を併用していけない理由もなく、例えば、アミノ酸ではなく核酸だけど、

新規イノシン-グアノシンキナーゼを用いたイノシン及びグアノシンのリン酸化による核酸系調味料の製造法
https://ci.nii.ac.jp/naid/500000230297/

には、『(核酸系調味料の製造について)数量面で 、現在(2001年)主流となっている方法は 、イノシン酸の製造においても、グアニル酸の製造においても、発酵生産したヌクレオシドをPOCl3を用いてリン酸化する方法である』と書いてある。

環境調和型オリゴヌクレオチド合成を指向するリン酸とアルコールの触媒的脱水縮合反応の開発研究
https://ueharazaidan.yoshida-p.net/houkokushu/Vol.22/category/pharmacy_03.html

では、2008年の文章だけど、以下のように書かれている。

リン酸モノエステルの合成法は炭素や窒素の化学と比べると合成法の種類が限られており,優れた合成プロセスがないのが現状である.現在,リン酸モノエステルを合成する場合,リン酸化剤として塩化ホスホリルを用いてアルコールと縮合させる方法が汎用されている.しかし,塩化ホスホリルは反応性が高いため,等モル量のアルコールと反応させるとリン酸ジエステルやリン酸トリエステルが副生してしまうという問題がある

生物で、リン酸化反応が果たす重要性を考えると、合成法が限られてるのは意外な感じもする。生物では、特定の生成物だけが欲しいということはない点で、工業化学とは違ってて、生命誕生時には、リンが沢山あって、片っ端から反応したのだろうけど。

グルコースの非生物的・非酵素的合成は、学術研究としては当然あるが、産業利用で検討されたことは、あまりないっぽい。

エネルギー源と炭素源

エネルギー源

1950〜60年代の食糧合成では、エネルギー源も炭素源も、原油が想定されていたが、現在は、 (1)原油がいつ枯渇するのか不明なので、今後も、エネルギー源として原油を利用し続けられるかは不透明である (2)二酸化炭素の排出量を問題視する風潮(本当に問題なのかは別として)があって、原油から食糧合成した場合でも、最終的に体内で二酸化炭素になって放出されることになるので良くないかもしれない という2つの問題がある。

エネルギー源に関しては、太陽ですら、いつかは燃え尽きる(と考えられている)ので、枯渇しない資源はないと思うけど、太陽光が一番豊富な資源なのは間違いないだろう。五年とか十年程度の短期的なスパンでのことは知らないけど、地球人口が現在の何倍〜何十倍にもなるという状況を想定するなら、結局は太陽光発電で何とかするしかない状況が発生すると思われる。

そして、太陽光発電が、農業と比較して効率的なのは、もう確認したので、いいだろう。

炭素固定

炭素源として(原油であれ原油以外であれ)化石燃料を用いて食糧を合成した場合、地球温暖化を心配する人は、燃料にする場合と同様の懸念を持つことになる。地球温暖化を心配しない人も、枯渇時期を検討しておく必要はある。化石燃料の生成速度について、現時点で確実なことが分からないので、問題が顕在化するのが数十年後か数百年後か数万年後か全然分からないけど、不確実な事項が複数ある中で、ブラックボックスにしたまま消費し続けるより、制御可能な技術を開発した方が話は単純になる。

炭素源の候補として、農業残渣を利用するという方法も考えられるけど、食糧に含まれる炭素の数倍程度なので、そこまで多いわけではない。将来的には、現在は未知の技術によって安価に、元素を生成できることも起きるかもしれないけど、現時点では、何も言えない。

そんなわけで、長期的に利用可能な炭素源として、空気中の二酸化炭素を炭素固定する効率的な方法を確立しておくのは望ましいと考えられる。炭素源が空気中の二酸化炭素であれば、それを食糧として消費した後、二酸化炭素を放出しても、単に、二酸化炭素が循環するだけである。

炭素固定には、色んな反応が利用できるだろうけど、メタン生成菌は、水素ガスと二酸化炭素から、メタンを生成し、これによって得られるエネルギーを利用して、ATP合成を行う(資化できる物質は、他にも色々あるっぽいが、最も広く見られる反応は、これらしい)。収支だけ見れば、サバティエ反応と同じ反応っぽいけど、複数の酵素が関与しており、当然、素反応は異なる。

メタン生成菌の生理と利用
https://doi.org/10.1271/kagakutoseibutsu1962.30.537

メタン生成菌では、メタンは副産物であって、メタンから複雑な有機化合物を合成しているわけではないっぽいが、工業的には、メタンは、有機合成の出発点として単純であり有用である。とはいえ、メタンを分解する微生物もいないとメタンが蓄積していくばかりになる。メタンを主要な炭素源として利用しているメタン酸化古細菌もいるようであるが、代謝経路の詳細は、未だよく分かってないっぽい。

メタン生成と嫌気メタン酸化の酵素化学
https://doi.org/10.1271/kagakutoseibutsu.52.307

同じく、水素ガスと二酸化炭素から酢酸を作れる酢酸生成菌も、割と色んなとこにいるらしい(エタノールと酸素から酢酸を作る発酵とは全く別)。こっちも炭素固定兼エネルギー生成経路になってるようだ。酢酸の非酵素的/非生物的合成法としては、メタノール一酸化炭素から合成するものが工業的に使われることがあるそうだ。

人間にとってメタンは食べれないけど酢酸は食べれるという点で異なるものの、水素ガスと二酸化炭素から、メタンを生合成するのも、酢酸を生合成するのも、途中過程が割と似てるということで、どっちもWood–Ljungdah経路と呼ばれてたりする(狭義には、酢酸生成経路のみを指すっぽい)。メタンを生成するのは、古細菌に多く、酢酸を生成するのは、真正細菌に多いらしい。これらの代謝経路は、LUCA(現生生物の最も新しい共通祖先)が出現した時には既に存在していたという推測もある。

Beating the acetyl coenzyme A-pathway to the origin of life
https://royalsocietypublishing.org/doi/10.1098/rstb.2012.0258

窒素固定考

メタン生成菌や酢酸生成菌が、炭素固定を行うだけでなく、(エネルギーを消費するどころか)ATPを作り出せる理由は、水素ガスのおかげ。深海熱水中では、水素ガスが豊富に存在するという報告もあり、水素は、光合成より前の時代におけるエネルギー源だったのかもしれない。人間が、この反応を利用して炭素固定を行う場合、水素ガスを、どこから持ってくるかが当然問題となる。現代日本人なら、誰でも思いつく方法は、水の電気分解。けど、当然ながら、そんな方法で炭素固定しても、日本のように火力発電がメインの国では、天然ガス由来のメタンにコストで勝てない。

現在、産業の窒素源となっているハーバー・ボッシュ法でも、似たようなことはある。

アンモニア製造プラントの原料原単位
http://comtecquest.com/RD/rd011.html

では、水の電気分解で、水素ガスを生成した場合に、ハーバー・ボッシュ法で消費されるエネルギー量が書いてある。投入エネルギーの大部分(ここの計算では、94〜95%)は、水素ガスを作るのに使われてるっぽい。アンモニアの冷却・分離コストを無視すると、全体では、22(GJ/ton-NH3)=374(kJ/mol)相当という計算になってる。kJ/molと書いてるのは、アンモニア1モルを製造するのに、374kJのエネルギー消費ということ。

考慮されているのは、水の電気分解による水素生成エネルギー、深冷分離によって空気から窒素を分離するエネルギー、窒素と水素を高圧に圧縮するエネルギー。なんか微妙に途中計算の単位が奇妙(間違っているわけではないものの)だし、あと、水素と窒素の温度をあげる分のエネルギーが入ってない気もする。

後者について、現在のプラントは、摂氏500度くらいで運転してるそうなので、アンモニア1モルの製造について、窒素分子0.5モルと水素分子1.5モルを500度ほど加熱するエネルギーが必要になる。これは、気体定数をRとして、500R=4.157(kJ/mol)の7倍程度のオーダー(二原子分子気体の定圧モル熱容量は、並進、回転、振動自由度を考慮すれば、"高温"で約3.5R)で、ざっくり30(kJ/mol)くらいだろう。

更に、反応が一瞬で終わるのでない限り、反応が終わるまで、温度と圧力を維持する必要があるが、このために必要なエネルギーは、反応速度以外の条件にも依存するので一概には決められない(プラント内が冷える速度は、体積:表面積比で決まるだろう)。あと、これに、アンモニアの冷却・分離コストが加わると書いてある。

生物学的窒素固定では、アンモニア1分子を生成するのに、ATP8分子を使うようである。ATPの加水分解エネルギーを50(kJ/mol)とする単純計算だと、生物学的窒素固定は、400(kJ/mol)程度で機能してることになる。ATP加水分解の自由エネルギー変化50(kJ/mol)は、ヒトの細胞内環境の数値だと思う(標準Gibbs自由エネルギーは31kJ/mol)ので、窒素固定細菌の細胞内環境では、多少値が変わるかもしれないけど、生物の窒素固定は、水の電気分解によるハーバー・ボッシュ法の限界効率と同水準を達成しているのかもしれない。

日本では、メタンなどの炭化水素と水蒸気を反応させて、水素を作る方法(水蒸気改質。フィッシャー・トロプシュ法の逆の反応)が、よく使われてるらしい。この場合、副産物として、二酸化炭素が生じて、この二酸化炭素アンモニアと反応させて尿素が製造できるので、アンモニア製造と尿素製造を一体化できるっぽい。現在、アンモニアのかなりの部分は、尿素合成に使われるようなので、一緒にやってしまうと都合がいいということだろう。

エネルギー的にも、水蒸気改質(多分、高温の水蒸気を作るのにエネルギーの大部分を消費すると思われる)は、水の電気分解より、大分有利なようである。火力発電が主流な現在の日本では、水の電気分解による水素製造は、水蒸気改質による水素製造より経済的にも不利だけど、他の発電方式が主流な場合は、無条件に、水蒸気改質が有利とは言えない。

以下のスライドTABLE12には、現在のハーバー・ボッシュ法で、水素の生成法に応じて、80.5(g/kWh)〜127(g/kWh)程度の効率を実現しているとある(冷却・分離精製コストを含んでるのかは不明)。アンモニアの分子量17(g/mol)と、1kWh=3600kJから、480~760(kJ/mol)程度の範囲にあるという計算。

地球環境にやさしい社会のために、化石燃料から再生可能エネルギー
http://web.tuat.ac.jp/~doce/yumenabi.pdf

差分の280(kJ/mol)は、水の電気分解で水素1.5モルを製造するのに必要なエネルギーと、水蒸気改質で水素1.5モルを製造するのに必要なエネルギーの差に由来すると思われる。上の計算と違って、水の電気分解によって水素を生成する場合でも、水素の生成エネルギーは、半分未満ということになる。何にエネルギーを消費してるのか内訳までは書かれてない。

"経済産業省生産動態統計年報 化学工業統計編"を眺めると、近年のアンモニア国内販売価格は、6万円/t=1.02(円/mol)前後を、うろうろしてる。スライドの生成コストを信用すると、水蒸気改質による場合、エネルギー単価が、7.65(円/kWh)未満でないと利益が出ない計算。例えば、2010年代の日本のLNGとLPGの価格が、おおまかに50(円/kg)で、単純に発熱量が50(kJ/g)とすれば(メタン、エタン、プロパンの燃焼熱がこんなもの)、3.6(円/kWh)で辻褄は合ってそう(火力発電で、エネルギー変換効率が60%なら、燃料費は、6円/kWh)

水蒸気改質の場合、これとは別に、材料としても天然ガスが必要で、材料費は別計上。これも計算しておく。仮に、天然ガスが100%メタンである(LNGの90%以上はメタンらしい)とすれば、理想的なメタンの水蒸気改質では、メタン1モルあたり水素3モルが生成するから、アンモニア1モル当たりメタン0.5モルが材料として要求される。エネルギー源としては天然ガスを使わなくてもいいが、同じ天然ガスを使うとした場合、メタンの燃焼熱は、890(kJ/mol)なので、アンモニア1モルの製造にメタン0.54モルの燃焼が必要で、加えて材料として0.5モルのメタンが必要なので、合計1.04モルのメタンが要求される。従って、100%メタンの天然ガスを使う場合、50(円/kg)=0.8(円/CH4-mol)で、アンモニア1モル当たり0.83(円/mol)が製造にかかる費用となる。

プロパン100%と仮定して計算するなら、50(円/kg)=2.2(円/C3H8-mol)ということになるけど、燃焼熱2200(kJ/mol)で、プロパン1モルから水蒸気改質で7モルの水素が作れる。従って、メタンの計算で使った数値を流用すると、プロパン100%の天然ガスを使う場合、アンモニア1モルの製造に、プロパン0.43モルが必要な計算になり、アンモニア1モルの製造価格は、0.94(円/mol)となる。メタン100%、プロパン100%のいずれの場合でも、アンモニア国内販売価格を下回ってはいる。これに、人件費や装置の保守費、減価償却費などが上乗せされるはず。

水の電気分解によって、水素1.5モルを製造する場合、「アンモニア製造プラントの原料原単位」の数値を使うと、356kJのエネルギーが必要で、これが、メタンの水蒸気改質による水素1.5モルの材料費用0.4円(製造費用は、水の加熱が必要なので、もう少し高いが、簡単のため材料費を使う)を下回るためには、電力価格が4.04(円/kWh)未満になる必要がある。国内の業務用電力単価が10〜20(円/kWh)なので、現時点で、日本では水蒸気改質の方が圧倒的に有利だということになる。

電気料金の国際比較-2019年までのアップデート-
https://criepi.denken.or.jp/jp/serc/discussion/20010.html

の図2を見ると、産業用電気料金が(為替レートにも依存するが)4円/kWhを達成している国はないようである。アメリカの一部地域は電力価格も安いが、天然ガス価格も安いので、水蒸気改質が圧倒的に有利なのは変わらないだろう。

電力会社から電気を買わず、自社で原子力発電施設を構築できるなら、4円/kWhを下回ることはできるかもしれないが、なんか色々(法律的と技術的な)問題がありそうである。水蒸気改質の方は、天然ガスの採掘、冷却と輸送に掛かる費用が必須で、こっちも安くなる可能性はあるが、考え出すとキリがないのでpass。しかし、水素ガス製造について、水の電気分解のほうがコスト的に不利なのは、物理的な制約によるものではないので、将来的に、コストが逆転する可能性がないわけではない。

いつかは、水の電気分解による水素ガスを用いたハーバー・ボッシュ法が経済的に有利という時代が来るのかもしれない。それくらいになれば、炭素固定によるメタン生成も、視野に入ってくると予想される。窒素に比べると、大気中の二酸化炭素濃度は低いけど、それが、どの程度の問題になるのかは、よく分からない。

メタンや酢酸以外の有機物(メタノールやギ酸)をターゲットにした炭素固定も色々検討されてはいる。二酸化炭素アンモニア尿素を生成する反応も炭素固定と言えなくもないけど、そういう言い方をしないのは、尿素有機合成の出発点としては便利じゃないせいだろう。

窒素固定の時も、ハーバー・ボッシュ法以前に、Franck-Caro法という石灰窒素製造法が主流だった時代がある(日本でも、1908年に日本窒素肥料が特許実施権を購入)し、窒素と酸素を反応させて、硝酸を得るBirkeland–Eyde法というものも、1903年に特許が取られ、ノルウェーの会社で、暫くは使われていたっぽい。ちょっと面白い話として、「稲妻」は雷が豊作に結び付いてることから来る呼称だという話があるけど、自然界では、雷によって、この反応に基づく窒素固定が起こってて、豊作に結び付いてるのかもしれない。

炭素固定でも、暫くは、複数の方法が競合する時代が暫く続くのかもしれない。

高分子の合成

水と空気と太陽光から、グルコースアミノ酸を合成したとして、それだけでは、まだパンは作れない。パン生地やうどん生地の素になるらしいグルテンは、名前は聞いたことあるけど、実体は、グルテニンとグリアジンという名前のタンパク質らしい。

Glutenin
https://www.uniprot.org/uniprot/P08488

alpha-gliadin
https://www.uniprot.org/uniprot/Q9ZP09

配列を見ると、(タンパク質としては標準的だけど)結構でかい。これを、化学的手法でアミノ酸を繋げて合成するのは、コストが高く付きすぎる。将来的には、DNAの鋳型から転写、翻訳をin vitroで安価に行う系が作れる可能性もあるかもしれないけど、現時点では、生物に合成させるのが手っ取り早そう。

グルテンフリーのパンでは、片栗粉や米粉を使うらしい。片栗粉の実体は、デンプンで、これは、グルコースのポリマーである。"でんぷんうどん"というものも存在するらしいので、デンプンで、グルテンの代わりができるのかもしれない(ベトナム料理のフォーは、米粉から作るそうだし)。やったことないので本当に片栗粉オンリーでパンが作れるのか定かではないけど、検索すると、卵白、砂糖、片栗粉のみでパンを作ったという人もいる。砂糖は原理的には合成できるし、味付けのためで無くてもいいだろう(料理としては、マズイだろうが)。卵白は、多分、空気を含ませてフワフワした食感を作るためのものじゃないかと思うので、何かで代替はできるだろう。

卵白を含まない起泡性組成物、メレンゲ様起泡物及び食品
https://patents.google.com/patent/JP2011147357A/ja

とかを見ると、セルロースのある種の誘導体が代替として利用できるのかもしれない。用途に応じて、増粘剤、安定剤、ゲル化剤などと表記されるらしい。セルロースなので、多分消化できないと思うけど、食品としては、食物繊維と同系統だろう。

cf)加熱するとゲル,冷却するとゾル溶液となるメトローズⓇ
https://doi.org/10.5458/bag.6.1_64

グルコースポリマーには、セルロース、グリコーゲン、デンプン(アミロースとアミロペクチンの混合物)などあって、異なる性質を持つ。種類が色々あるせいで、樹脂のように化学合成するのは、一般に容易ではなさそうだが、酵素合成はできるかもしれない。

セルロースは、(フッ化β-D-セロビオシルという分子を重合することで)in vitroで酵素合成ができてるっぽい

Novel method for polysaccharide synthesis using an enzyme: the first in vitro synthesis of cellulose via a nonbiosynthetic path utilizing cellulase as catalyst
https://doi.org/10.1021/ja00008a042

酵素触媒重合によるセルロース・セロオリゴ糖の合成 糖鎖工学の新しい側面:糖鎖工学の新しい側面
https://doi.org/10.1271/kagakutoseibutsu1962.31.385

また、江崎グリコ株式会社は、アミロースとグリコーゲンの酵素合成技術を持っているっぽい。

グリコーゲンの製造方法(特願2006-537787)
https://www.j-platpat.inpit.go.jp/c1800/PU/JP-2006-035848/6ED54827558E11EA130EBA8777161A32BC4C3148B82AAAD8B4FB2FDCBFA00E98/19/ja

ジメチルスルホキシド中でのアミロース合成(特願2008-121691)
https://www.j-platpat.inpit.go.jp/c1800/PU/JP-2009-270012/61BCFD7D7AF23E382A6EB81CC3996BE94340E2C769A24B152F821763828605C6/11/ja

デンプンの酵素合成と利用
https://doi.org/10.3136/nskkk.56.551

調べて初めて知ったけど、社名のグリコは、グリコーゲンに由来するそうで、「グ・リ・コ・ゲ・ー・ン」「パ・イ・ナ・ッ・プ・ル」「チ・ョ・コ・レ・ー・ト」なら、全部6文字で平等だったということになる。創業当初は、カキから抽出してたそう。

アミロース酵素合成は1986年には論文があったっぽいけど、いずれにしろ、意外と新しい技術ということになる。

The enzymic utilization of sucrose in the synthesis of amylose and derivatives of amylose, using phosphorylases
https://doi.org/10.1016/0008-6215(86)85078-9

アミロペクチン酵素合成は、調べたけど見つからないので、ないっぽい(?)。有用性の有無に関わらず、試そうと考えなかったとは思えないので、なにか難しいことがあるのかもしれない。

駿河屋のHPの和菓子技術者になるための必須知識を見ると、デンプンを、どの植物から抽出するにしろ、アミロペクチンの含量が多い(最も少ない小麦で70%がアミロスペクチンとなっている)。

澱粉に水を加えて熱すると、全体が一様にコロイド状態となる」ことを糊化(α化)と呼ぶそうだが、アミロースは直鎖状で、アミロースだけでは糊化しないだろう(やったことないので推測だけど)。そうすると、デンプンでパン生地が作れるとしても、アミロペクチン含量の高いデンプンが必要そうで、現時点では、グルコースアミノ酸から、パンを合成するのは難しいかもしれない。

味の作り方

アミノ酸の味 with 食レポ

アミノ酸単体の味

現在の人類には、パンを化学合成する力はないという結論になってしまった。高分子の合成は大変だけど、そもそも、高分子自体には、一般的に味がないので、食感を作り出すためにしか寄与していない(分解産物は、味にも栄養にも寄与するけど)。デンプンだって、分解しなければ、別に甘味とかない。そんな高分子を作るのは無駄な気もする。

肉や魚を熟成させるのも、タンパク質を分解してアミノ酸にするためという説明は、よく見る。しかし、糖と違って、殆どのアミノ酸は、単体で摂取しても、あんまり美味しくない。味が殆ど感じられないアミノ酸も多い。逆に、アルギニンは、とても苦い上に、特有の臭いがある。アルギニン粉末を直に熱すると、臭いが拡散して、酷いことになる。グリシンは甘味があるが、砂糖ほど美味しい感じでもない。

参考までに、私の主観で、アミノ酸粉末を舐めて感じた味を、以下に記録しておく。主観なので、信用はできない。タンパク質を構成するアミノ酸20種類全部をコンプできなかった。D-アミノ酸に至っては、全く確認できてない。

甘味 旨味 苦味 酸味 塩味
グリシン +
DL-アラニ + +
L-グルタミン酸 +
L-アルギニン ++
L-プロリン + +
L-スレオニン +?
L-ヒスチジン +?
L-リジン +? +?
L-グルタミン +?
β-アラニ +?
タウリン
  • L-α-アラニンは入手できなかったので、DL-アラニンしか評価できなかった
  • β-アラニンは、α-アラニンの構造異性体
  • "+?"とあるのは、ほぼ無味に近く、微かに感じた味を記した
  • タウリンは、ほんとに完全に無味な気がする

以下の論文によると、マウスとヒトの旨味受容体は、アミノ酸選択性が全然違うらしい。

機能解析技術が明らかにした味覚受容体と食物成分のかかわり味の感じ方は生き物それぞれ
https://katosei.jsbba.or.jp/view_html.php?aid=1108

ヒトでは、グルタミン酸が一番強く結合し、アスパラギン酸も他のアミノ酸より強く結合する傾向があるが、それに比べると、残りのアミノ酸は、ほぼ横並びっぽい。一応、ヒトでも、どのアミノ酸も、旨味応答を示すっぽいので、全く旨味がないわけではないのかもしれない。私が旨味を感じなかったものは、旨味が非常に弱いか、他の味が強いかのどちらかかもしれない。

犬猫にも食べさせてみたいけど、飼ってないので試せていない。

ヒトの甘味受容体は、T1R2とT1R3という2つのタンパク質のヘテロ二量体らしい。2002年の論文

An amino-acid taste receptor
https://doi.org/10.1038/nature726

には、マウスのT1R2+T1R3受容体の、いくつかのアミノ酸への応答の強さが計測されている(Figure1(b))。ヒトT1R2+T1R3で、網羅的に、同様の計測をしたデータは発見できなかった。

Human receptors for sweet and umami taste
https://doi.org/10.1073/pnas.072090199

には、複数の物質に対するヒト甘味受容体の応答を計測したデータが載っており、グリシンは、弱い応答を示すようである(Figure2B)。また、ヒトの場合、D-トリプトファンは、とても強い甘味応答を示すようである。データを信じるなら、グルコース、フルクトース、スクロースよりも甘味が強いらしい。D-トリプトファンは、食品添加物としては販売されてなくて、簡単に入手できそうにないので、味を確認できていない。

タンパク質は、L-アミノ酸のみから構成されていることから予想される通り、一般的に、食品のD-アミノ酸含有量は多くないっぽい。マウスの場合、いくつかのD体アミノ酸が甘味応答を示し、L-体のアミノ酸の甘味応答は、極めて弱いらしい。ところで、ヒトである私がL-プロリンを舐めると、甘味と苦味があると感じられる。L-プロリンがヒトT1R2+T1R3受容体に結合するのか、他にも、甘味受容体が存在するのかは不明。

ズワイガニ風溶液の調整

味を決めるアミノ酸(生物工学基礎講座-バイオよもやま話-)
https://dl.ndl.go.jp/info:ndljp/pid/10517786

には、『グリシンやアラニンは,甘味が主体であるが弱いうま味も呈し ,ヌクレオチドとうま味の相乗効果も示す.アルギニン自体は弱い苦味を呈するアミノ酸であるが,ズワイガニでは “こく”と表現されるような複雑な感覚に関わっている.』とあり、他の物質と混ぜることが重要らしい。

このPDFの表2に、ズワイガニ脚肉エキスの成分含量というのがある。この成分はオミッションテストで特定されたと書いてあるが、この表を参考に、以下の成分を混ぜて、味を見た。甲殻類アレルギーの人でも大丈夫なはず(甲殻類アレルギーの原因は「トロポミオシン」タンパクのことが多いらしい)。本物のカニと違って、ボトルに入れて、数日間、室温放置しても腐ったり劣化した形跡はなかった。多分、微生物が繁殖するには栄養が不足しているのだろう。

成分 分量
ミネラルウォータ(SUNTORY天然水) 2000g
グリシン 12.43g
DL-アラニ 3.74g
グルタミン酸ナトリウム 0.44g
L-アルギニン 11.58g
核酸調味料 0.18g
食塩  5.18g
塩化カリウム 7.51g
重曹 (食品添加物) 6.51g
85%リン酸 3.1ml(5.2g)

重曹は、ナトリウムイオン濃度の調整とリン酸の中和という2つの役割を兼ねている。ミネラルウォータは製品によっては、高い濃度で、ナトリウムイオンやカリウムイオンを含んでいることがあり、成分バランスを大きく変化させる場合があるので注意が必要。

PDFの表によると、核酸で重要な成分とされているのは、AMPとGMPであるが、残念ながら、AMPは入手できなかった。GMPも単体では入手できなかったので核酸調味料で代用するしかなかった。核酸調味料は、通常、5'-IMPと5'-GMPのナトリウム塩らしいが、比率などは不明で、100g肉中に含まれるIMPとGMPは合計9mgなので、合計重量で調整した

リン酸は、完全に電離するわけでもなく、1価や2価のイオンとして存在するものも多いはずで、3価のリン酸イオンが重要なのだとしたら、もっと量を増やすべきなのかもしれない。しかし、計算が面倒なので、3価リン酸イオンに含まれるリン原子量を再現すればいいという方針で、85%リン酸の分量を決めている。

少し別の視点での試算も行っておく。ヒトのリン元素重量比は、体重の1%とされ、その85%は骨に存在するそうなので、骨以外の組織に存在するリン原子は、体重の0.1%程度の重量しか占めない。カニでも同水準と仮定して、肉100gに対して、リン原子が100mgくらいあると推測する。316mgのH3PO4に含まれるリン原子の質量は約100mgである。核酸調味料を無視すれば、水2000gに対して、7.4gの85%リン酸を加えることで、肉100gに対してリン原子が100mgという存在比と同等の配分が実現する計算。上記リン酸量は、それより少なめだが、許容範囲内ではあるだろう。

肝心の味に関して、カニなのかどうかは、よく分からなかった。カニだと言われて飲めば、そんな気もしたが、AMPが含まれてないので、カニとは微妙に違うのかもしれない。少なくとも、何も教えられずに摂取して、何の味か問われたとしても、カニと答えるのは難しいように思えた。上記スープを加熱して、ポン酢を加えて飲むと、カニの味だという感じがずっと強まった。どうやら、私の脳内では、カニの味は、ポン酢とセットで記憶されているらしい。

(ポン酢なしで飲んで)美味しいかどうかで言うと、私の味覚は全く当てにはならないが、不味くはないと言っていいと思う。ほぼ無臭なので、物足りなさがあるとすれば、それが理由かもしれない。ごま油と米のみで雑炊にしても美味しい。このスープが何だろうと、多量のアルギニンが含まれているにも関わらず、不味くないのは意外だった。

粉末を舐めた時は気付かなかったけど、アルギニン10gを水100ccに溶かして飲むと、苦いけど、飲み込んだ後、エビやカニっぽい後味が残る。そういや、カニ味噌も苦みがあるけど、あれはアルギニンの味なんだろうか。水100ccにアルギニン1gを混ぜた場合、苦味も薄まるけど、エビ・カニっぽさも減るので、上記分量は、カニらしさを感じるには、少ないのかもしれない。いずれにせよ、グルタミン酸核酸などの旨味物質をドバドバ入れればいいというものではないということらしい。

味とは直接関係ないけど、エビやカニが茹でると赤くなるのは、アスタキサンチンが遊離するためだそうだ。アスタキサンチンは、β-カロテンと非常に似た構造をしていて、吸収されると、βカロテンと同じく、一部はビタミンAに変化するらしい。β-カロテンは、赤い人参の主要色素だと言われている。植物は、複数の色素を持つことが一般的らしく、人参の赤とエビやカニの色が微妙に違うとすれば、そのへんの差に由来するのだろう。赤くない緑黄色野菜にもβカロテンは含まれてて、紫蘇なんかは人参より含有量が多いようなのだけど、アントシアニン色素(ぶどう、紫人参や紫キャベツなどに含まれる。pHによって色が変わるそうなので、食べられるpH指示薬として利用できるかもしれない)によって紫色に見えてるらしい。これらの色素に味があるのかは知らない。

日本食品標準成分表にあるズワイガニアミノ酸成分は、上記脚肉エキスと異なっている。多分、脚肉エキスの成分は、熱水抽出液で、日本食品標準成分表は、タンパク質の加水分解を行っているようなので、脚肉エキス成分は、主に遊離アミノ酸の量を表しているのだろう。実際の料理の味に近いのは、脚肉エキスの成分で、一方、日本食品標準成分表は、味より栄養を評価したものと考えられる。燻製などにすれば、(タンパクの分解が進んで)全体の成分量は、日本食品標準成分表に近付くのかもしれない

特に意図したわけではなかったけど、グリシン、アラニン、アルギニン、グルタミン酸(ナトリウム塩)は、アミノ酸の中では、個人レベルでも、比較的安価に入手できる。現在のところ、1kg当たり1000〜3000円くらいで買える(核酸調味料も1kg3000円程度で入手できる)。タンパク質のアトウォーター係数4kcal/gで概算すると、0.25〜0.75円/kcal相当。一日の摂取カロリーが2000kcalで、食費1500円とした場合、0.75円/kcalなので、一般的な食費と同水準と言っていいだろう。そういうわけで、上記カニ風エキスを作るのに一番金がかかる材料は、水になる。

かつお出汁風溶液の調整

マグロ粉末も作りたいと思ったが、上記のように、味の重要成分が特定されているケースは少ないようで、探しても見つからなかった。1973年の論説

魚貝類の味
https://doi.org/10.3136/nskkk1962.20.432

には、

魚類については,呈味に関係ありそうな個々の成分の分析はかなりよく行なわれているが,前述の無脊椎動物のように,多数の成分を詳しく分析して呈味との関連を解析した例がほとんどないようである。 と書いてあり、50年近く経った現在でも、この状況に大きな変化はないのかもしれない。

食品会社が、データを公開せずに持っているという可能性はないわけではないが、マグロについては、あまり情報がない。かつおは、マグロに近い種だが、かつお出汁について調べた報告が多数ある。上記報告でも、また、1989年に書かれた

かつお節のエキス成分
https://doi.org/10.3136/nskkk1962.36.67

でも、かつお節の遊離アミノ酸は殆どヒスチジンのみで、核酸ではIMPが主だと書いてあって、更に、この報告は、冒頭で

かつお節だし中に多量含まれる遊離ヒスチジンの呈味効果については,ほとんど効果が認められないという報告があるが,同じく遊離ヒスチジンを著量含む赤身魚の肉の味には深い関係をもつと推測する報告もあり,いまだ明解な結論は得られていない.

と述べている。omission testをやる予定と書いてあるが、続報は存在してない。

仕方ないので、重要成分を推理してみる。

まず、ズワイガニ脚肉エキスと比較すると、IMP量が際立って多く、AMPもズワイガニと同水準で存在している。また、ナトリウム、カリウム、塩素イオン濃度も高く、塩味もズワイガニ脚肉エキスに比べて、強そうである。グルタミン酸量は、ズワイガニ脚肉エキスと同水準なので、旨味には寄与して可能性が高い。IMPの量を勘案すれば、旨味もズワイガニより強いのかもしれない。グリシンズワイガニ脚肉エキスより顕著に少なく、数%しかないが、アラニンはズワイガニ脚肉エキスの1/3ほど存在しているので、弱いが甘味に関与してると推測される。

ヒスチジンは、わずかに苦味がある気がするので、大量にあれば、苦味に一定の貢献をしている可能性はある。マグネシムイオンも、そこそこあり、塩化マグネシウムは、にがりの主成分だから、これも苦味に寄与してるかもしれない。また、かつお出汁は、乳酸量が多い。かつお節を煮詰めすぎると、酸味が出るというのは、よく知られたことらしい。逆に、カニにポン酢を合わせるのは、足りない酸味を補っているのかもしれない。

通常、かつお出汁のpHは、5〜6程度らしい(自分で測定したわけではないが)。発酵調味料である味噌や醤油は、乳酸菌が産生した乳酸を含んでおり、pHは5.0を若干下回るそうだ。かつお出汁の酸味は、強いというほどでもないらしい。

ズワイガニ脚肉エキスでは、リン酸イオンが重要成分として挙げられているけど、その影響は不明。リン酸自体は、強い酸味を持つが、それは水素イオンの効果で、リン酸イオンの方に味があるのかは分からない。他に陽イオン成分として存在すると考えられ、栄養上も重要である鉄や亜鉛などは測定されてないが、これらの含有量は、ほぼすべての食品で、100g中10mg未満なので、ナトリウムやカリウムイオンなどに比べれば、相当に少ない量しか存在しないと予想される。

総合的に見ると、かつお出汁は、ズワイガニ脚肉エキスに比べて、旨味と塩味が強く、甘味と苦味は弱めで、酸味もあるというバランスだと推測される。

上の報告中のかつお節エキス成分には、100g中乳酸量が3.415gと書かれているが、水100gに乳酸3.415gを混ぜると、とても酸っぱい。計算上のpHは、2.1強になる。しかし、仮に、乳酸の味への寄与が、酸味のみだとすれば、乳酸の量は、酸味が殆ど感じられなくなる程度まで減らしてもいいだろうと予想される。

以上の推測のもと、以下の溶液を調整した。市販の乳酸には、二量体の無水乳酸が結構含まれてるそうだけど、もし、無水乳酸が全く含まれてなければ、以下の分量で、pHが6くらいになるはず(計算が間違ってなければ)。 | 成分 | 分量 | |:--------------------------:|:------------------:| | ミネラルウォータ | 550g | | DL-アラニン | 0.28g | | グルタミン酸ナトリウム | 0.13g | | L-ヒスチジン | 11.0g | | 核酸調味料 | 2.6g | | 食塩 | 6.1g | | 塩化カリウム | 7.2g | | にがり | 2.7g | | 90%乳酸 | 3.3g |

雑炊を作ったところ、美味しかった(雑炊ベンチマーク)。海産物感はあったけど、かつお出汁なのかは不明。そもそも、昆布出汁とセットでお出しされるのが一般的な気がするので、かつお出汁単体で摂取した経験が自分にあるのか分からない。世の中には、「にんべん 液体かつお節だし」という商品があるようで、これと比較すればよさそうだけど、1L単位で販売されてて、使いきれる気がしないので、未購入。

この溶液は、塩分と核酸調味料が比較的多く、塩水とは違うのは明らかだけど、単純に、核酸調味料の味なんじゃないかという疑惑はある。それを調べるために、上記分量比率を固定して、以下の成分のみを含む溶液も作った。

(1)水+核酸調味料のみ

(2)水+核酸調味料+グルタミン酸ナトリウム

(3)水+核酸調味料+グルタミン酸ナトリウム+L-ヒスチジン

(1)(2)(3)のいずれも、味が薄くて、出汁と言うより味の付いた水という感じだった。かつお出汁風溶液の味は、核酸調味料のみで再現できるものではなく、塩味は重要らしい。L-ヒスチジンは、現在のところ、比較的高価な方のアミノ酸なので、不要なら使わずに済ませたいところだが、(2)と(3)には、味の違いがあるように感じられた。

L-ヒスチジンが高価なので、現時点では、かつお節から出汁を取るほうが、上記溶液より安価だと思う。

自作合成醤油

Wikipediaによると、第二次大戦中には、既に代用醤油というのが製造されていて、人毛醤油も検討されたとか書いてある。

実験室でのアミノ酸醤油の作り方
http://www.eiyotoryoris.jp/archive/ER10_06_/ER10_06_010.jpg

は1944年の資料だそうだが、18%塩酸でタンパク質を加水分解しと後、重曹でpHを調整する的なことが書いてある。アミノ酸の比率は、あんまり関係ないのかもしれない。異世界転生して醤油が欲しくなった時には、有用な可能性もある。

醤油の組成も、Wikipediaに詳しい。アミノ酸では、グルタミン酸が多い。また、無機イオンでは、ナトリウムが突出している。これは、塩分の濃さを考えれば当然。醤油の色は、メイラード反応によるものと書いてあるが、本当かどうかは分からない。Wikipediaには、核酸の組成が記述されていない。

市販加工醤油の分析
https://agriknowledge.affrc.go.jp/RN/2010690655

の表3を見ると、昆布醤油という名称の商品では、イノシン酸グアニル酸は殆ど含まれていない。昆布醤油は、砂糖やぶどう糖も比較的少ない。

昆布醤油が成分的に一番意外性があったので、これを模倣することにして、以下の溶液を調整した。

成分 分量
ミネラルウォータ 100ml
食塩 14g
砂糖 3g
90%乳酸 1.1g
グルタミン酸ナトリウム 0.7g
グリシン 0.3g
DL-アラニ 0.3g
L-アルギニン 0.5g
L-プロリン 0.5g

混ぜながら、これただの塩水だろと思ったけど、完成品の味は、意外と、醤油感があった。検索すると、「もろみ」の食塩濃度は、16〜19%くらいに調整されると書いてあるので、食塩濃度は、こんなものでいいっぽい。

無色透明の醤油というのは、食欲が削がれるので、色を付けたいと思って、上記溶液で、メイラード反応が起きんかと、100度程度で10分ほど加熱したけど、色は付かなかった。

料理を科学する!~至高の食を目指して~
http://www.hikonehg-h.shiga-ec.ed.jp/blog/wp-content/uploads/2020/08/0765daf40c7ddb7dbe631260fa600100-2.pdf

を見ると、90〜95度では、スクロースグリシンorリジンと混合して、加熱しても着色しなかったとある。メイラード反応で着色するなら、グルコースでもリジンでも、フルクトースが一番いいようだ。pHは、塩基性条件下で反応が速かったとあり、酸性条件下では全然ダメっぽいので、乳酸は、着色後に追加する方がよさそう。いずれ、気が向いたら試してみようと思う。

メイラードの原論文は電子化されてないのか見つけられなかったけど、以下の論文に、メイラードがやったことが多少書かれている。

アミノーカルボニル反応のLouis Camille Maillard (1878~1936)
https://doi.org/10.1271/nogeikagaku1924.56.1199

スクロースグリシンは短時間加熱しても着色しないが、3時間加熱すると(恐らくスクロース加水分解されて)着色が始まるそうだ。また、キシロースは、ヘキソースよりも反応しやすいと書いてある。

戦時中のアミノ酸醤油は、カラメルで色を付けていたそうで、スクロースでやるなら、その方が簡単だろう。

醤油の香味成分HEMF
https://doi.org/10.6013/jbrewsocjapan1988.101.151

によると、香気成分は無数にあるが、酵母が生合成するHEMFという物質とエタノールが混ざると、醤油様の香りがすると書いてある。HEMF単体では、甘い香りがすると書いてある。

HEMFは、構造的には、フラネオールと呼ばれる物質に似ていて、フラオネールには、ストロベリーフラノンという呼称もある通り、苺の匂いを呈するのだと思われる。HEMFをフラオネールで代用できるかは分からない。

HEMFは、化学的な名称は、4-ヒドロキシ-2-エチル-5-メチル-3(2H)-フラノンで、ストロベリーフラノンは、4-ヒドロキシ-2,5-ジメチル-3(2H)-フラノンで、一箇所、エチル基かメチル基かという違い。

食肉

海産物は、出汁(熱水抽出液)の分析が多いけど、牛、豚、鶏なんかでは、ミンチにして、肉の遊離アミノ酸分析がされている。

食肉の遊離アミノ酸
https://agriknowledge.affrc.go.jp/RN/2010721305

では、2%スルホサリチル酸を加え、ホモゲナイズしたと書いてるので、ミンチにしてるっぽい。

ズワイガニ脚肉エキスに比べると、全体的に遊離アミノ酸量が少ない。その中では、牛、豚、鶏共に、グルタミンとタウリンが多い。牛ではグルタミンの方が多く、鶏ではタウリンの方が多い。豚は、その中間で、重量比で、ほぼ同量。ただ、グルタミンもタウリンも、単体では、殆ど味がないので、かつお節のヒスチジンと同じく、あんまり意味はないかもしれない。

アンセリンカルノシンのようなジペプチドも多いが、味に影響してるのかは分からない。牛や豚には、結構な量の尿素が含まれるが、鶏にはない。これも理由不明。

異なる品種間の鶏肉における遊離アミノ酸、ジペプチド、イノシン酸
http://www.naro.affrc.go.jp/org/tarc/to-noken/DB/DATA/063/063-073.pdf

でも、鶏肉をミンチにして分析したと書いてある。ズワイガニ脚肉エキスと比べると、アミノ酸量が少ないのは変わらないが、IMP量も測定されていて、これは、ズワイガニに比べて、とても多いが、かつお出汁よりは少ない。グルタミン酸量は、かつお節出汁やズワイガニエキスより少し多い(しかし、これは前者の報告では、もっと少なくなっている)。鶏肉の味は、ほとんど旨味なのかもしれない。なんか海産物に比べて、面白みがない。そして、牛肉、豚肉、チキンの味の差が何に起因してるのかも、よく分からない。そもそも、これらの肉の(匂いや食感を除いた)"純粋な味"に本当に差はあるのだろうか。

もっと多くの食材で、味の重要成分を特定してもらいたいものである。人肉の味とか化学的に作れれば面白いと思う。アミノ酸や無機イオンの定量分析は、高価な機械なしでも、比較的容易に家で実施できると思うので、自分でやるべきかもしれない。

全体的に、一発勝負で混合しただけの割には、そこそこの物ができたので、味を特徴付ける成分というのは、意外と少なく、味を作るのは、それほど難しくないのかもしれない。マッポーの世では、カチグミ・サラリマンが、"天然物は香りが違う"とドヤるのかもしれない。

アミノ酸以外の味

核酸単体の味については、1968年の古い報告

ヌクレオチド系呈味成分と調理
https://doi.org/10.11402/cookeryscience1968.1.1_12

によると、美味しさでは、5'-IMPと5'-GMP、5'-XMPが強いらしい。ATPの重要性を考えると、5'-AMPは、もっと美味しくてもいい気もするけど、単体では、それほどでもないらしい。ただ、前述のカニのように、食品の味を特徴付ける重要な成分としては、カウントされていることもある。

味の評価に於いて脂質の影響は考慮されてないことが多い。1970年に書かれた

油脂の味
https://doi.org/10.5650/jos1956.19.612

を見ても、油脂は、基本的には無味無臭だとか書いてある。

トリアシルグリセロールを分解するリパーゼは、ヒトの唾液にも含まれているらしく、そうすると脂肪酸が遊離することになる。教科書的には、リパーゼは、胃液や膵液での記述が中心なので、昔は、知られてなかったのかもしれない。炭素数の少ない脂肪酸は、水に溶けて酸味成分となり、また揮発しやすく、香りの元でもある(小さい例は酢酸で、場合によっては、臭い)。

脂肪酸の炭素数が増えるにつれて、揮発しなくなり、水にも溶けにくくなる。長鎖脂肪酸は、水に不溶とされている。21世紀に入って、(長鎖)脂肪酸が基本味だという報告もされている。

脂肪酸の美味しさ不味さの生体メカニズムの解明へ向けて
https://doi.org/10.5650/oleoscience.21.261

Free fatty acid receptor 1 (Mus musculus)
https://www.uniprot.org/uniprot/Q76JU9

Free fatty acid receptor 4 (Mus musculus)
https://www.uniprot.org/uniprot/Q7TMA4

Free fatty acid receptor 1 (Homo sapiens)
https://www.uniprot.org/uniprot/O14842

Free fatty acid receptor 4 (Homo sapiens)
https://www.uniprot.org/uniprot/Q5NUL3

マウスでは、FFAR1とFFAR4が、直接の受容体として機能しているそうだが、それぞれのアナログがヒトにも存在する。ヒトFFAR4は脂肪酸受容体として機能してるらしいが、ヒトFFAR1も脂肪酸受容体として機能してるかは分からない(そもそも、舌で発言してるかも)。Wikipediaによると、唾液中のリパーゼの至適pHは5前後であり、弱酸性食品だと、口内での脂肪酸分解が促進され、脂肪味も感じやすくなるかもしれない(つまり、ココアはpH5にして飲むべき?)

高度不飽和脂肪酸と鶏肉とのおいしさの関連性の解明(1)
https://agriknowledge.affrc.go.jp/RN/2010831952

は、鶏もも肉の遊離アミノ酸含有量と脂肪酸含有量を調べている。遊離アミノ酸の量やIMPの量に関しては、既に見た他の報告と同じ傾向にある。そして、意外でもないけど、脂肪酸組成は、中華調味料の鶏油とよく似ている。

匂いは、微量の揮発成分に由来し、アミノ酸脂肪酸の分解、糖とアミノ酸の反応などによって生じるものと思われる。未調理状態での匂いは、良いものばかりでもなく、例えば、肉や魚は、それぞれ特有の生臭さがある。こうした匂い成分も、それなりに特定されているようだ。

魚の生臭さとその抑臭
https://doi.org/10.5650/jos1956.29.469

揮発性成分を指標としたベニズワイの品質評価(第3報)
https://tiit.or.jp/userfiles/report2014/yanohara.pdf

では、63種類の揮発成分が列挙されている。中には、食品添加物として利用されてる成分もある。カニで揮発成分の一つとして挙げられているヘキサナールは、豆臭の主成分とされ、大豆に於いては、リノール酸が、酵素的に分解されて生じるという報告もある。ヘキサナールは、リノール酸が自然に酸化する過程でも、ある程度は生じると思われ、酵素的に作られるのが、どれくらい一般的かは知らない。他の揮発成分も、自然に生成したものも含まれてるだろう。

食べるためだけに高分子を合成する意味は、あまりないが、高分子を含まない食事を、何世代も続けていると、"進化"によって消化能力を失うかもしれない。尤も、その頃には、遺伝子改変技術も進歩してるだろうから、そうなったとしても、問題ない可能性もある。

信頼区間のオーバーラップで比率の差の検定をするのは、どれくらい悪いのか

ダメな統計学
https://id.fnshr.info/2014/12/28/stats-done-wrong-ja-pdf/
という本の第6章に、二群間の平均値や比率の差の検定をする代わりに信頼区間が重なってるかどうかで、差があるかどうかを判断している論文が存在すると書いてある。

平均値や比率の差の検定は、ほぼ全ての統計学の教科書に書いてあるし、"有意差"という言葉があるくらいに統計学で一番使われるパターンでもあるだろうから、普通に検定すればいいという気はする。けど、信頼区間の重なりを見る方法でも、十分多くのサンプル数を確保して、信頼度を1に近づけていけば、正しい判定ができる確率は高くなるはずだから、本当に無条件で悪いと言えるかは疑問がある。

信頼区間の重なりによる"差の検定"は、有意水準が自明には分からないという問題があるけど、これは計算可能なので、致命的な問題とは言えない。教科書的な仮説検定に優位性があるとすれば、より少ないサンプル数で、より良い判定が下せるということであるべきだろう。


というわけで、比率の差の検定について、信頼区間の重なりを見る方法が、どれくらい悪いか計算したい。仮説検定の失敗としては、本当は差がないのに差があると判定してしまう場合と、本当は差があるのに差がないと判定してしまう場合がある。それぞれの確率を評価すればいいだろう。

差の検定では、帰無仮説として「母比率が等しい」、対立仮説として「母比率が異なる」という仮説を採用することが一番多いと思うので、このケースを考えることにすると、本当は差がないのに差があると判定するのは、帰無仮説が正しいのに棄却してしまう第一種のエラーとかいうやつで、本当は差があるのに差がないと判定してしまうのは、対立仮説が正しいのに帰無仮説を採用してしまう第二種のエラーとかいうやつ。

この計算は純粋に数学上の問題に過ぎないけど、無味乾燥なので、具体的な問題として、ゲームで、クリティカル率がアップするというスキルがあって、スキル使用時と未使用時のクリティカル率に本当に差があるか(プレイヤーが)検証したいという状況があることにする。そんなスキルのあるゲームを思い出せないので、ゲーム例が念頭にあるわけではない。製品の不良率とかと違って、仮にバグでスキルが反映されてなかったりしたら、クリティカル率は完全に等しいだろう。

ドラクエでは、武闘家は、他に比べて、会心の一撃が出やすいらしいので、職業間比較でもいい。ただ、検索したところ、レベルがあがると会心の一撃の確率もあがると書いてたりする。データ収集中にレベルがあがると困る(セーブデータを何度もロードし直せばいいけど)。ドラクエ3では、レベル99武闘家は、確率99/256=38.6%で会心の一撃が出て、武闘家以外は、1/64=1.6%なんだそうだ。これくらい差があると、統計学に頼らなくても体感で分かるだろう。



スキル未使用時のクリティカル率がp_1で、使用時のクリティカル率がp_2だとする。大抵のゲームでは、クリティカル率が実際に何%かという情報は、プレイヤーには公開されないので、これらは未知の量である(公表されてても、開発者がバグを仕込んでる可能性もある)。スキル未使用時に、n_1回の攻撃をしてx_1回クリティカルが出て、同様に、スキル使用時にn_2回の攻撃をして、x_2回クリティカルが出たとする。

x_1x_2は、二項分布\mathrm{Bin}(n_1,p_1),\mathrm{Bin}(n_2,p_2)に従う。二項分布の正規近似によれば、信頼度1-\etaの信頼区間は、
(\dfrac{x_i}{n_i} - z_{\eta/2} \sqrt{\dfrac{p_i(1-p_i)}{n_i}} , \dfrac{x_i}{n_i} + z_{\eta/2} \sqrt{\dfrac{p_i(1-p_i)}{n_i}} )
になる(i=1,2)。

但し、z_{\eta/2}は、\mathrm{erf}(z_{\eta/2}/\sqrt{2}) = 1 - \etaを満たす定数。
\mathrm{erf}(x) = \dfrac{2}{\sqrt{\pi}} \displaystyle \int_{0}^{x} e^{-y^2} dy

実際に信頼区間を計算する時には、未知量p_iが入ってるので、この形では使えないけど、今は、信頼区間を計算すること自体が目的ではないので、このままにしておく。

本当は差がないのに差があると判断してしまうということは、p_1=p_2=pなのに、上記の信頼区間が重ならない場合である。その条件は
|\dfrac{x_1}{n_1} - \dfrac{x_2}{n_2}| \ge z_{\eta/2}( \sqrt{\dfrac{p(1-p)}{n_1}} +  \sqrt{\dfrac{p(1-p)}{n_2}} )
と書ける。簡単のためにn_1=n_2=nの場合を考えると
|\dfrac{x_1}{n} - \dfrac{x_2}{n}| \ge 2 z_{\eta/2} \sqrt{\dfrac{p(1-p)}{n}}
という条件になる。

二項分布の正規近似に基づくと、\dfrac{x_i}{n}は、平均p、分散\sigma^2 = \dfrac{p(1-p)}{n}正規分布に従うと近似でき、その差であるX=(x_1-x_2)/nは、平均0、分散s^2=2\sigma^2正規分布に従う。上の不等式は
|X| \ge \sqrt{2} z_{\eta/2} s
と書ける。x_0 = \sqrt{2} z_{\eta/2} sと置くと、本当は差がないのに差があると判断してしまう確率は
\displaystyle \int_{x_0}^{\infty} \dfrac{1}{\sqrt{2\pi s^2}} e^{-x^2/2 s^2} dx + \displaystyle \int_{-\infty}^{-x_0} \dfrac{1}{\sqrt{2\pi s^2}} e^{-x^2/2 s^2} dx = 2\displaystyle \int_{x_0}^{\infty} \dfrac{1}{\sqrt{2\pi s^2}} e^{-x^2/2 s^2} dx
で計算できる。

ちょっと計算すると、この確率は1 - \mathrm{erf}(z_{\eta/2})となる。一般的な信頼度95%の時(従って\eta=0.05)は、これは、大体0.56%くらいになるから、第一種のエラーは起きにくくなっている。ということは、第二種のエラーが起きやすくなっているということでもあるが、オーバーラップ法で検定する場合、個々の信頼区間の信頼度は、もっと低くていいとも言える。オーバーラップ法で第一種のエラーが起きる可能性を5%に調整するには
z_{\eta/2} = z_{0.025} / \sqrt{2}
を満たせばよくて、\eta = 1 - \mathrm{erf}(z_{\eta/2} / \sqrt{2}) = 1-\mathrm{erf}(z_{0.025}/2)は、0.166くらいになる。つまり、(n_1=n_2=nという条件で計算していることは忘れるべきでないが)信頼度83〜84%の信頼区間を使えば、オーバーラップ法で第一種のエラーが起きる可能性は、5%程度になる。

あるいは、信頼区間の重なりで判断する場合に、第一種のエラーを起こす確率を\alphaにしたければ、次の信頼区間を使えと言ってもいいかもしれない。
(\dfrac{x_i}{n_i} - \dfrac{z_{\alpha/2}}{\sqrt{2}} \sqrt{\dfrac{p_i(1-p_i)}{n_i}} , \dfrac{x_i}{n_i} + \dfrac{z_{\alpha/2}}{\sqrt{2}} \sqrt{\dfrac{p_i(1-p_i)}{n_i}} )



次に、本当は差があるのに差がないと判定してしまう確率は、どれくらいあるか考える。ちょっと想像すれば分かることだけど、(他の条件が同一なら)検出しようとする差が小さいほど、差がないと判断する確率は高くなるだろう。クリティカル率が、2%から2.01%になったところで、体感で差はないだろう(スキルの説明文には、具体的な数値は書いてないのが普通なので、仕様通りだと主張はできるかもしれないけど)し、多くのデータを集めなければ、差があると確信できないだろう。というわけで、どれくらいの差を検出したいか、最初に決めておく必要がある。

今、p_1=p_2+\Deltaであるとして、信頼区間(信頼度は先ほどと同様\gamma=1-\etaとしておく)にオーバーラップが生じる確率を計算する。クリティカル率の問題については、\Delta<0であるべきだが、先入観に囚われず、スキル使用時のクリティカル率の方が低い可能性も検討することとする(そういうバグもないとは言えないし)。\Deltaは、これより小さい差は、検出しそこなっても許容できる大きさに取る(あるいは、有意差がないことを示したい場合は、どれくらい小さい差なら無視できる差なのかを指定することにする)。勿論、これは一種の目安であって、これより大きい差は、より高い確率で検出できるし、これより小さい差も検出できる確率は低くなるが、0になるわけではない。

さっきと同様、二項分布の正規近似に基づくとx_1/n_1は平均p_1、分散\sigma_1^2 = \sqrt{\dfrac{p_1(1-p_1)}{n_1}}正規分布に近似的に従い、x_2/n_2は平均p_2、分散\sigma_2^2 = \sqrt{\dfrac{p_2(1-p_2)}{n_2}}正規分布に近似的に従い、X=\dfrac{x_1}{n_1} - \dfrac{x_2}{n_2}は、平均\Delta、分散\sigma^2=\sigma_1^2+\sigma_2^2正規分布に従うと近似できる。

信頼区間にオーバーラップがある条件は
|X| \lt z_{\eta/2}(\sigma_1 + \sigma_2)
と書ける。これが起きる確率は
\beta = \displaystyle \int_{-z_{\eta/2}(\sigma_1+\sigma_2)}^{z_{\eta/2}(\sigma_1+\sigma_2)} \dfrac{1}{\sqrt{2\pi\sigma^2}} e^{-(x-\Delta)^2/2\sigma^2} dx
となる。

計算すると、これは
\beta = (-\mathrm{erf}(\Delta - \dfrac{z_{\eta/2}(\sigma_1+\sigma_2)}{\sqrt{2}\sigma}) + \mathrm{erf}(\dfrac{z_{\eta/2}(\sigma_1+\sigma_2) + \Delta}{\sqrt{2}\sigma}))/2
になる。この式を変形して、\eta,\beta,p_1,p_2から、nを計算する式を得ることにする。

二通りに場合分けして、第一に
|\Delta| \lt z_{\eta/2}(\sigma_1+\sigma_2)
の時は、ほぼ100%、差がないと判断してしまう。

それ以外の場合の
|\Delta| \ge z_{\eta/2}(\sigma_1+\sigma_2)
という条件成立時の評価は、面倒くさいが、統計学の力が必要なのは、p_1p_2の差が小さい場合なので、(\sigma_1+ \sigma_2)^2 \approx 2(\sigma_1^2+\sigma_2^2)=2\sigma^2と近似できるとしよう。すると、第二種エラーを起こす確率は
\beta = (\mathrm{erf}(z_{\eta/2} - \dfrac{\Delta}{\sqrt{2}\sigma})  - \mathrm{erf}(-z_{\eta/2} - \dfrac{\Delta}{\sqrt{2}\sigma}))/2
と近似できる。厳密さが気になるなら、\sigma_1+ \sigma_2 \le \sqrt{2} \sigmaなので、第二種エラーを起こす確率は、この近似値を超えることはないと言ってもいい。

更に、\dfrac{\Delta}{\sqrt{2}\sigma} = \delta z_{\eta/2}によって、定数\deltaを導入すると
\beta = (\mathrm{erf}( (1-\delta)z_{\eta/2}) - \mathrm{erf}((-1-\delta)z_{\eta/2}) )/2
という形になる。上に出てきた
|\Delta| \ge z_{\eta/2}(\sigma_1+\sigma_2)
という条件は、大体、|\delta|>1に相当する。

信頼度が十分大きいとして、\delta \gt 1の場合、
-\mathrm{erf}((-1-\delta)z_{\eta/2}) \approx 1
で、\delta \lt -1の場合、
 -\mathrm{erf}((1-\delta)z_{\eta/2}) \approx 1
なので、以下のように近似できる。
1 - 2 \beta = \mathrm{erf}( (|\delta|-1)z_{\eta/2} )

この式は
z_{\beta} = \sqrt{2} (|\delta|-1)z_{\eta/2}
と同値なので、n_1=n_2=nとして、これを変形すると、
n = \dfrac{(\sqrt{2}z_{\eta/2} + z_{\beta})^2}{\Delta^2}  \left( p_1(1-p_1) + p_2(1-p_2) \right)
を得る。この式によって、\eta,\beta,p_1,p_2から試行回数を決定できる。

オーバーラップ法で、第一種のエラーを起こす確率を\alphaにするためには、信頼区間の信頼度1-\etaは、z_{\eta/2} = z_{\alpha/2}/\sqrt{2}と選べばよかった。すると、上の関係式は
n = \dfrac{(z_{\alpha/2} + z_{\beta})^2}{\Delta^2}  \left( p_1(1-p_1) + p_2(1-p_2) \right)
と書ける。

実際には、p_1,p_2は未知量なので、事前に試行回数を決めるためには、p_1(1-p_1) + p_2(1-p_2)の値は、何らかの当て推量を行う必要がある。何の予測もなく、データを集めるのに(時間や予算の)制限がないなら、とりあえず、想像される範囲で、最大値を取っておけばいいだろう。例えば、クリティカル率なら、p_1,p_2は共に10%は超えないと思って、2*0.1*0.9=0.18にするとか、そのような仮定をおかず、最大値の2*0.5*0.5=0.5を仮定しておくなど

こういう見積もりは、必要なサンプルサイズを過大に評価する可能性があり、検出したい差\Deltaを、1-\betaより高い確率で検出できるようになる(これは朗報だろうけど)。一方で、試行回数を増やせば、どんな小さな差も検出できるので、あんまり意味のない小さな差でも、有意差があると結論付ける場合がある。

クリティカル率を上昇するスキルが有用かどうかは、実際の上昇量や効果時間、使用コストなどの兼ね合いで決まるはずなので、有意差があるというだけでは不十分だろう。

時には、サンプルサイズの設計を行わずに、事前に集められたデータのみから結論を出したいこともあるかもしれない。その時には、試行回数は変更することができないので、\alpha,\beta,nを指定した上で、どのくらいの差なら検出できるのかという計算をすることになるのだと思う。



以上は、信頼区間の重なりを見るという方法で"検定"した時の性能評価なので、他の方法でどうなるかを調べて比較する必要がある。比率の差を検定する方法の一つに、z-unpooled testと呼ばれるものがある。これは、1985年に提案されたもの。
Exact Unconditional Sample Sizes for the 2 × 2 Binomial Trial
https://doi.org/10.2307/2981892

z-unpooled testで使用する検定統計量は、Wald信頼区間の計算に使われるものと同じ。似た方法にz-pooled testもあり、1986年の以下の論文に書かれている
An exact unconditional test for the 2 × 2 comparative trial
https://doi.org/10.1037/0033-2909.99.1.129

z-pooled testとz-unpooled testは、統合z検定、非統合z検定とでも訳すのが適当なんだろうけど、一般的な訳語を知らないので、英語のまま使う。実際の所、二群の比率の差の検定に関して、z-pooled testはカイ二乗検定と同値である(z-pooled testの検定統計量は、正規分布に従うが、その二乗は、自由度1のカイ二乗分布に従う)。

カイ二乗検定は、統計学のどの教科書にも載ってるけど、初出は1900年のKarl Pearsonの論文とされている
On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling
https://doi.org/10.1080/14786440009463897
Pearsonの論文を読むと、カイ二乗検定は、標本空間が離散集合である時の適合度検定として提案されているが、標本空間が離散集合である時には、有限回の独立試行は、多項分布で記述されるので、多項分布に対するパラメトリック検定と思っても同じことになる。


z-pooled testとz-unpooled testでは、比率の差の分散で、微妙に異なる計算をする。"pooled"という名前は、分散の計算時に、二群間の比率を個別に計算する代わりに、両群の統合比率(pooled proportion)を使うことに由来する。二群の母比率が等しい場合には、z-pooled testの方が良い推定値を与える。z-unpooled testでは、比率の差の分散を小さめに評価するが、サンプル数が多ければ、両者の性能に、あんまり違いは出ない(サンプル数が増えれば、分散は小さくなり、真の分散に収束するはずなので)。サンプル数が少ないと、二群の母比率が等しい時に、「二群の母比率に差がない」という帰無仮説を棄却する確率が、設定した有意水準よりも大きくなる。

z-pooled testもz-unpooled testも、二項分布の正規近似に基づいてるので、サンプル数が物凄く少ない時には、設定した有意水準と、帰無仮説を棄却する確率が食い違うということは起こるが、z-unpooled testの方が、食い違いが大きく、サンプル数を増やした時の改善が少し遅い。そんなわけで、「二群の母比率に差がない」という帰無仮説を検定する場合には、z-pooled testを使った方がいいと言われている。


信頼区間の重なりを見る方法は、二群間で個別に比率を出してるので、z-unpooled testに近いと予想できるが、実際に、検定の有意水準\alpha、検出力1-\beta、検出したい差\Delta、試行回数nの関係は、オーバーラップ法のそれと似た形で書ける。オーバーラップ法の第二種エラーの計算の時と同様に、真の比率p_1,p_2にはp_1=p_2+\Deltaという関係があるとしておく。

まず、\dfrac{x_1}{n_1} - \dfrac{x_2}{n_2}は、平均\Delta、分散\sigma^2=\sigma_1^2+\sigma_2^2正規分布に従うが、z-unpooled testでは検定統計量として、
z = \dfrac{p_1 - p_2}{\sigma}
を使う。これは、平均が\Delta/\sigmaで、分散が1の正規分布に従う。検出力と有意水準は、
\beta = \displaystyle \int _{-z_{\alpha/2}}^{z_{\alpha/2}} \dfrac{1}{\sqrt{2\pi}} e^{-(x-\dfrac{\Delta}{\sigma})^2/2} dx
という関係を満たす。y=\sigma xと変数変換をすると、この関係は
\beta = \displaystyle \int_{-z_{\alpha/2} \sigma}^{z_{\alpha/2} \sigma} \dfrac{1}{\sqrt{2\pi\sigma^2}} e^{-\dfrac{(y-\Delta)^2}{2\sigma^2}} dy
オーバーラップ法の時と積分区間が微妙に異なっている以外は、同じ形をしている。

従って、同様に計算していくと、
n = \dfrac{(z_{\alpha/2} + z_{\beta})^2}{\Delta^2}  \left( p_1(1-p_1) + p_2(1-p_2) \right)
という近似式が得られる。

つまり、適切に信頼区間の信頼度を設定すれば、信頼区間の重なりを見るのと、z-unpooled testの性能は、ほぼ同等になる。オーバーラップ法の時は、(\sigma_1+ \sigma_2)^2 \approx 2(\sigma_1^2+\sigma_2^2)という近似を使ったが、この近似がない場合、第二種エラーの確率を与える積分区間が微妙に違っている。帰無仮説が正しい時、つまり、2つの比率に差がない時には、\sigma_1=\sigma_2なので、両者の積分区間は厳密に等しいが、2つの比率が異なっている時には、z-unpooled testの方が、積分区間が大きいので、第二種エラーを起こす確率は高く、検出力は少しだけ小さくなるはず。\sigma_1\sigma_2の差が小さいほど、検出力の違いも小さくなる。

というわけで、信頼区間の重なりで差の有無を判定するのは、信頼区間の信頼度を適切に設定すれば、z-unpooled testより、僅かに良いということになる。但し、多分、実際に、信頼区間の重なりを見てる人たちは、信頼区間の信頼度を適切に設定してはいないと思われる。それに、信頼区間自体の幅自体が比率に依存してて、実際には、これも正しく推定できないという問題もあるので、そこまで考えると、性能がどうなるかは分からない。



z-unpooled testで必要なサンプル数nの具体的な数字を見ておく(n_1=n_2=nとしているので、二群のそれぞれについて、以下の試行回数を要することに注意)

p_1,p_2が、それぞれ、20%と40%で、有意水準5%、検出力80%とした場合、nは、78.5

p_1,p_2が、それぞれ、20%と40%で、有意水準0.5%、検出力95%とした場合、nは、198.2

p_1,p_2が、それぞれ、10%と20%で、有意水準5%、検出力80%とした場合、nは、196.2

p_1,p_2が、それぞれ、2%と10%で、有意水準5%、検出力80%とした場合、nは、134.4

p_1,p_2が、それぞれ、2%と10%で、有意水準0.5%、検出力95%とした場合、nは、339.4

p_1,p_2が、それぞれ、5%と10%で、有意水準5%、検出力80%とした場合、nは、431.7

p_1,p_2が、それぞれ、35%と40%で、有意水準5%、検出力80%とした場合、nは、1467.7

p_1,p_2が、それぞれ、2%と4%で、有意水準5%、検出力80%とした場合、nは、1138.1

p_1,p_2が、それぞれ、2%と4%で、有意水準0.5%、検出力95%とした場合、nは、2873.8

p_1,p_2が、それぞれ、12%と14%で、有意水準5%、検出力80%とした場合、nは、4434.6

p_1,p_2が、それぞれ、12%と14%で、有意水準0.5%、検出力95%とした場合、nは、11197.9

p_1,p_2が、それぞれ、0.1%と0.2%で、有意水準5%、検出力80%とした場合、nは、23507.4

計算式から明らかかもしれないけど、有意水準と検出力を、かなり厳しくしても、必要な試行回数は意外と増えないが、小さな差を検出しようとすると、必要な試行回数は急激に増える。また、検出しようとする差が小さい時には、比率の差だけでなく、絶対的な大きさも強い影響がある。1%と2%の差を検出するより、50%と51%の差を検出するのは、ずっと多くの試行を必要とする。

疫学、心理学、教育学などの(人間の個体を対象とする)諸分野では、サンプル数nが、100〜1000の規模であることが多いように思う。比率の差10%くらいを検出しようとしてるなら、n=1000くらいあれば概ね問題なさそう(有意水準、検出力、母比率の差を固定した時、nが最大になるのはp_1,p_2が45%と55%の場合で、有意水準0.5%、検出力95%でも、nは981.1)。数%の差を検出しようとしてるなら、n=1000でも不十分な場合があるという感じだと思う。

サンプル数100~1000くらいの実験が多い(ように思われる)のは、サンプルサイズを慎重に検討した結果というよりは、データ収集コストの兼ね合いで決まってることも多そうだけど、検出できる差には(確率的なブレがあるとはいえ)、サンプルサイズ(と有意水準、検出力)から決まる限界がある。簡単な計算で、z-unpooled testでは、A=(z_{\alpha/2}+z_{\beta})^2とする時、
|\Delta| \ge \dfrac{A}{n+A}
を示すことができる。Aの値は、有意水準5%かつ検出力80%なら、7.85で、有意水準5%かつ検出力90%なら、10.5になり、有意水準0.5%かつ検出力5%なら、19.82になる。

この式が言ってるのは、例えば、比率0%と1%の違いを調べたいなら、それぞれの群で、n=1000回くらい試す必要があるというだけのこと。nの値は、有意水準と検出力次第で変わるけど、標準的な設定なら、このオーダーになる。少なくとも、n=100では、比率1%の群で、観測したい事象が起きる回数の期待値は1回に過ぎず、一度も観測できない確率も30%以上ある(一度も観測できない確率が1%を下回るのは、nが459以上の時)のだから、結論を出すのに少なすぎることは直感的にも納得できる。勿論、偶然の偏りによって、小さいnでも有意差があると判定される可能性がないわけではない。けれど、基本的には、物理実験などと同じく、小さな差を検出したいなら、相応のコストを払う必要がある。



ゲームのクリティカル率の例題に即した話では、n=100で(有意水準5%とか1%とかで)有意差なしと結論が出たなら、クリティカル率が2%と20%みたいな大きな差があるということは、余りなさそう。けど、2%と10%くらいの差ということは、まだありそうなので、バグだと結論付けて、メーカーに苦情をいうには、試行回数が足りなすぎるだろう。

n=1000で有意差なしという結論が出ても、2%と4%の差を検出しそこなってる可能性は、そこそこある。2%と4%の差は、大きくはないかもしれないけど、2倍違うということもできるので、(ゲームデザインを勘案して)使えるスキルかどうか判断するには十分だとしても、バグを主張するには、ちょっと心もとない。

n=10000で有意差なしと結論が出たなら、2%と3%くらいの差は検出できてる可能性が高いので、控えめな姿勢で、メーカーにバグの可能性を報告をしてもいいかもしれない。しかし、一秒に一回攻撃できるとしても、2n=20000のデータを集めるのに、6時間くらいかかる。そんな検証より楽しいことは、多分沢山あるだろう。



信頼区間のオーバーラップの話に戻ると、以下の文献には、信頼区間の重なりで、差の有無を判定している文献が挙げられてる。
On Judging the Significance of Differences by Examining the Overlap Between Confidence Intervals
https://doi.org/10.1198/000313001317097960
ここで例として挙げられているのは、医療系の論文っぽいけど、標準的な検定を行った場合と、結果がどう変わるのかは書かれていない。信頼区間の信頼度を変更したら、どうなるかということも議論されてないっぽい

A brief note on overlapping confidence intervals
https://doi.org/10.1067/mva.2002.125015
も、内容的には同じ感じ。