prml11.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 第11章
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{chapter}{10}
\chapter{「サンプリング法」のための物理学}

PRML下巻のp.263には11.5.1の理解に「物理学の背景知識は必要としない」と明記されている．しかし学生時代に物理学を専攻した私としては，むしろ物理学の背景知識を前提に説明されたほうが何倍も理解しやすい．このテキストは，私が入社間もない頃の社内勉強会でこの節を担当した際に，「ああ要するにこの文章は統計力学の○○のことか」などと理解したときのメモを再構成した文章を，さらに今回見苦しくない程度に再構成したものである．
\begin{flushright}
サイボウズ・ラボ  川合 秀実（特別寄稿）
\end{flushright}
\vspace{0pt}

\section{「11.5.1 力学系」のところを違う方法で説明する試み}
%ハミルトン力学によるハイブリッドモンテカルロアルゴリズムの解説

さて，PRMLでは式(11.54)が十分な説明もなくいきなり現れる．もちろんいきなり出てきても説明や証明には何の問題もないのだが，その背景を分かったほうが理解が深まると私は思う．私は個人的な趣味によりこの節の紹介を物理学的知識全開で行うが（笑），ぎりぎり高校物理の範囲で収まっていると思う．またサンプリング法に関する知識は式(11.2)しか必要としていない．
PRML下巻のp.239--263での議論を引用することはない．

\section{統計力学}

統計力学とは，分子運動の振る舞いの集合の結果として気体・液体・固体などを解釈し，それによりそれまで経験的に知られていた熱力学の法則を説明・再構築する，物理学の一分野である．統計物理学，統計熱力学とも呼ばれる．
統計力学の概論をここで扱うようなことはしないが，この節の理解に役立つことについては紹介する．

\section{ボルツマン因子}

統計力学で頻出の概念に「ボルツマン因子」というものがある．これは結局は次の法則を数式で表したものに過ぎない．

たとえば箱に収められた気体を考えよう．
この中には多数の気体分子が存在し，それぞれがさまざまな運動エネルギーを持って運動している．言い換えれば，速いものもあれば遅いものもある．このときあるエネルギー$E_1$を持つ分子の個数を数えたとする．それを$N_1$個としよう．
またエネルギー$E_2$を持つ分子の個数を数えてみて，それを$N_2$個とする\footnote{%
実際にはちょうど$E_1$に等しいエネルギーを持っている分子はいないかもしれない．
その場合は，$E_1$～$(E_1+dE)$の範囲のエネルギーを持っている分子の数，みたいなもので代用すればより正確な議論になる．
$E_2$についても同様にして，そしてこの$dE$を十分に小さくしていけば，やはり$N_2$と$N_1$の比は上式の値に収束する．}．
ここではとりあえず $0 < E_1 < E_2$ としておく．このとき，以下のような関係が成り立つ．
%
$$\frac{N_2}{N_1} = \exp\left( -\frac{E_2 - E_1}{k T} \right).$$
%
ここで$k$はボルツマン定数で，
SI単位系では $1.38064852(79) \times 10^{-23}$\xu{}J/K に
なる\footnote{2014年CODATA推奨値．}．
$T$は絶対温度．

この式の意味するところは，エネルギーが $k T$ だけ増えた状態の分子の個数は，約$0.36788$倍だろうと予想できるということである（この数値は $1/e$ である）．
そしてこの個数の比はエネルギー差にのみ依存する．エネルギーの絶対値は関係ない．ほかの要素も関係ない．

箱の中の多数の分子同士は常に衝突しあってエネルギーをやりとりしている．そうしているうちに運よく大きな運動エネルギーを得る分子もいるが，そんな偶然はめったに起きず，たいていは速くなった分子は周囲の遅い分子とぶつかってしまって，平均的な速度になってしまう．上記の関係式はそういう事情を反映している．

なお，統計力学においては，この関係式により系の温度を考える．つまりそれなりに多数の分子が存在してかつそのエネルギー分布がわからなければ温度は定義されない．
分布がわかれば，分子数の比が$0.36788$倍になるために必要なエネルギー差を観測し，それをボルツマン定数で割れば温度がわかる．

（余談：統計力学では分布により系の温度を定義するため，レーザー発振などのために逆転分布を作ると，それを「負の温度」として表すこともある．また，もしどのエネルギーの分子も等しい個数だけいるような系があったとしたら，それは温度無限大と解釈される\footnote{%
つまり$N_1=N_2$なので，$\exp$の中は$0$にならなければいけない．そのためには，$T$を無限大にするしかない．ということでこの状態の温度は無限大と定義される．
}．逆に最低値のエネルギーを持つ分子しかいない状態になったら，たとえその最低値がどんなに大きな値であろうと，温度は$0$である\footnote{%
$E_1$を最低値のエネルギー値，$E_2$をそれ以上の何らかのエネルギー値とおけば，$N_2$は常に$0$で，$N_1$は常に$0$ではない．となると，$\exp$の中はマイナス無限大でなければならず，そのためには$T$が$0$になるしかない．ということでこの状態の温度は$0$と定義される．}．
）

ボルツマン因子とは，上の式から自明にわかる次のことを指す．
%
$$p(E) = C \cdot \exp{ \left( - \frac{E}{k T} \right)}.$$

つまり系の中でエネルギー$E$を持つ分子の存在確率は，$\exp(-E/kT)$に比例する．これをボルツマン因子という．$C$は正規化定数で，$p(E)$の合計が$1$になることを保証するためのものでしかないと思ってよい（実際は，$C$の逆数は分配関数$Z$とも呼ばれ，熱力学と統計力学を関連付ける式の中で頻出する重要な因子である）．

なお，ここまでは話を「箱の中の気体分子」に限定してきたが，実は温度が定義できるような系ならば（つまりエネルギーが多数の要素に分配されうる系なら），いつでも成り立つ一般論である．

このボルツマン因子をより基本的な統計力学の定理から導くこともできるが，それは本筋からそれるので省略する．

\section{ポテンシャルエネルギー}

上記の話では，$E$は気体分子の運動エネルギーについてのみ議論したが，分子は運動エネルギーに加えて位置エネルギー（ポテンシャルエネルギー）も持っている（無重力とかなら話は別だが）．そしてポテンシャルエネルギーについてもボルツマン因子は有効である．たとえば地面を無限に広い平面で近似し，各分子がポテンシャルエネルギー $E = m \cdot g \cdot h$ を持つとする（$h$は地面からの高さ）．
そうすると高いところの分子は低いところに比べて高いポテンシャルエネルギーを持っているといえる．当然のことながら，ボルツマン因子によれば高いところまでのぼれる分子はそう多くない．
\pagebreak
高さに対して指数関数的に少なくなる．つまり上空ほど空気は薄くなる．これは私たちの常識と一致している．

ちなみに温度300\xu{}Kのときに窒素分子(分子量28)の存在確率$p$を
ボルツマン因子で計算すると，標高0\xu{}mと標高2000\xu{}mでの$p$の
比が0.80倍と出る．温度300\xu{}Kというのは$27^\circ\text{C}$に
相当し，窒素は空気の大半を占めることを考えれば，
これはそう悪くない大気の近似である．
そして標高2000\xu{}mでの一般的な気圧をGoogleで調べたら0.80気圧のようだ．
まあ上空に行くほど気温が下がるし，地球は平板ではないので，
この近似はいろいろと問題があるが，しかし参考にはなるだろう．

\section{サンプリングへの応用}

  さてやっと本題に入ろう．私たちは $p(z)$ を考えている．ここでの$p$はボルツマン因子の$p$ではなくて，サンプリングしたい対象の確率の$p$である．この$p(z)$の分布にそった出現頻度の$z$の集合がほしい．それならば，この$p$から対応するポテンシャル$E$を構成し，その中で物理学的なシミュレーションをすれば，その（仮想的な）粒子の位置$z$のログは，サンプリングにふさわしいものになりそうではないか\footnote{$z$が位置を表す……という表現がよくわからなければ，
$z$は仮想粒子の地面からの高さを表す，とでも思ってほしい．$z$がスカラーであるときは，このたとえは悪くないと思う．}．
これが基本的なアイデアである．
$E$をうまく作れれば，（そして精密に物理学を再現しさえすれば）期待通りの確率で$z$がサンプルできる．

ということで，ついに(11.54)に似た式が登場する（この式はもちろん統計力学のボルツマン因子に由来している）．
%
    $$p(z) = \frac{1}{Z} \exp{\left( - \frac{E(z)}{k T} \right)}.\label{Z}$$

教科書とは違い，私はまだこの$kT$を$1$に置きかえない．この式だと，$E$から$p$を導く式のように見えておかしいので，これを$E$について解いておこう．
%
    $$E(z) = - k T \cdot \ln{\left( p(z) \cdot Z \right)}.$$

与えられた確率分布$p$に対して，こういうポテンシャルエネルギー$E$を持つ系の中での
仮想的な分子の動きをシミュレーションすればいい．つまりはそういうことだ．

  PRMLではここから解析力学に突入するのだが，これは難易度が上がってしまうので
私は違う方法を選んだ．大丈夫，そんな小難しい理論を使わなくても，十分に説明でき
る．物理は私のようなものでも理解できるほどに簡単なのだ．ということで，私は高校
物理でおなじみのニュートン力学で進める．

  さて，この仮想分子にはきっと質量があるだろう．これをとりあえず$m$としよう．この
仮想分子は常にポテンシャルエネルギーから力を受けているのだが，その力は，
ポテンシャルエネルギーを$z$で偏微分し，$-1$を乗じれば求められる．
%
    $$F = - \frac{dE}{dz} = k T \cdot \frac{1}{p} \cdot \frac{dp}{dz}.$$

  この仮想分子の加速度を$a$とおけば，分子の質量$m$は時間変動しないので， $F = m \cdot a$となり，
%
    $$a = \frac{k T}{m} \cdot \frac{1}{p} \cdot \frac{dp}{dz}$$
%
となる．最初の因子 $k T / m$ は適当に$1$にしてしまってもいいかもしれない．
%分配関数は分子の分布を決めるうえで重要なのだが，温度$T$や質量$m$と一緒にしか出てこないので，結局それほど気にする必要はなさそうである．
また$p$も$a$に対してこの形でしか現れないので，つまり微分した関数との比だけが重要
なので，$p$を正規化し忘れていても$a$の値は変わらない．

  これでこの仮想粒子をこの加速度に沿って動かしていくシミュレーションをして$z$の
ログをとればいいだけなのだが，少し問題というか注意点がある．それについて，次の節に書こう．

\section{注意点}

  私たちは，統計力学という「多数の分子が衝突しながら乱雑に動く」物理学を使って
いる．ということは，このシミュレーションにおいて分子はどのくらいの個数を扱わな
ければいけないのだろう．100個くらいでいいのか？1万個か？100万個くらいだろうか．
……いやいや欲を言って厳密を期するのであれば1モルくらいはほしいかもしれない．
つまりは$10^{23}$個程度である．

  もちろんそんなことはやっていられない（メモリの消費量が尋常ではない）．という
ことで1個の分子の動きをシミュレーションしていくだけで同等の結果を得る方法を
考えよう．これはいわば本質ではなくただの技巧（テクニック）である．

  やることは簡単で，まずシミュレーション時間を十分に長くとることだ．そうすれば
分子は（たとえ一つであったとしても）さまざまなポテンシャルの場所を探検してくれる．
そうすれば一分子ながら $E(z)$ 全体を十分に反映した$z$のログができて，
サンプリングが
できる．……おっと言い忘れたが，もちろんシミュレーション内に複数の分子を配し
（つまり$z$や速度$v$の初期値が違う），それぞれシミュレーションしてもいい．それは
並列化が有効な手法である．

  またシミュレーションの最初のころの$z$は信用できないとしてログからは捨てることを
推奨する．しばらくシミュレーションしていると仮想分子は$E(z)$を反映した場所を
うろつくようになるが，それまでは初期値の選び方の影響を強く受けてしまい，サンプリングに際して有効とは思えない．
  なお，初期値に恵まれないと（ポテンシャルエネルギー的な意味合いで）どこかのく
ぼみにハマってなかなか出られない，なんてこともありうるだろう．こうなると$E(z)$全体を十分に反映しているとは言い難い．
そこで，たまに$z$を乱数で初期化したらより良いと私は思う．
つまり分子はワープするのだ．
ただワープ後しばらくはやはり$z$は$E(z)$を反映していないので，サンプリング用のログとしては使わずに捨てるのを推奨する．

  次は温度の問題を論じよう．私たちは先ほど気楽に$kT/m$という係数を1と置いたが，こ
れは$T$を定数として扱っていることになる．なぜなら，$k$も$m$も定数だからだ．つまり
このシミュレーションは温度が一定の仮想世界を考えているということである．

  しかし普通のニュートン力学でシミュレーションをしているだけだと，これは達成で
きない．なぜなら，普通のニュートン力学だけのシミュレーションでは，系の
エネルギーが（計算誤差を除いて）一定値をとるからだ．つまり，ポテンシャルエネル
ギーが低い地点では運動エネルギーが多くあって，高い地点では運動エネルギーが少な
いことになる．しかし自然界の温度一定系では，そうなってはいない．分子はほかの
分子との衝突やもしくは輻射のやりとりで運動エネルギーを得たり失ったりしており，
結果的にポテンシャルエネルギーの値とは独立に運動エネルギーを得ている．運動エネ
ルギーの平均値は（周囲の）温度にしか依存しない．つまり，温度一定とエネルギー
一定は一般には両立しないのだ．

  ということで，運動エネルギーをたまに調整してやらねばならない．これは分子同士
の衝突現象の代用である．つまりは速度を適当に決めなおすということだ．これをやら
ないと$p(z)$の分布は実現しないので重要である（実は当初私はこれをミスっておかしなサンプリング結果になり泣かされた）．
  速度ベクトルの各成分は乱数で決めればいいだろう．運動エネルギーもボルツマン
因子に従うはずなので，正規分布な乱数を使えばいいだろう．このとき，温度や質量の
関係が加速度$a$の算出に使ったものと矛盾しないようにすべきだろう．私がここで言わん
としていることは，速度を決めなおすときに乱数を適当に使うわけだが，その速度に
よる運動エネルギーの期待値が，$(kT/2) \times (\text{$z$の次元数})$ くらいになるように，
乱数の分布に気を配ってほしいということである．

  もし，$z$が4次元以上のベクトルであれば，それはもはや物理学からは離れてしまうが，
ニュートン力学は第4の座標変数があったとしても自然に拡張可能であり（誰でも容易に
類推できる），それでおそらく問題はないであろう．

\subsection{はみだしコラム「分配関数$Z$は$z$の関数なのか？」}

  分配関数$Z$は，ボルツマン因子を確率として正規化するための定数である．
%
$$ p(z) = \frac{1}{Z} \exp \left( - \frac{E(z)}{k T} \right)
\hseq \text{（\pageref{Z}ページ参照）．} $$

  この$Z$は，系が取りうるエネルギーのついてのボルツマン因子をすべて計算し，それを足し合わせることで求められる．これは確率の正規化の基本を思い起こしてもらえれば自明である．つまり，すべての確率を足したら$1$にならなければいけないので，正規化前のものをとにかく全部足して，それで割ってやればよいということだ．……正規化定数はその名の通り定数なのであるが，しかし別の視点から見ると関数でもある．それゆえに$Z$は分配「関数」などという別名があるのだ．その話を私はしたい．

  いろいろと理屈をこねてもいいのだが，とにかく一回$Z$を計算してみようではないか．それが一番わかりやすいだろう．今ここに，エネルギーがとびとびの値しか取れない実験装置がある．階段のような地形でしかも物体がなんらかの理由で地面から離れられないと想定すればいいだろう．もしくは，量子力学的に量子化された状態だと思ってもらってもいい．とにかく，ここではエネルギー$E$が$0$, $1$, $2$, $3$, $4$, \ldots と整数値しかとりえない．そういう状況を考えてほしい\footnote{%
何か適当な定数$c$を考えて，$E=0$, $c$, $2c$, $3c$, $4c$, \ldots とすれば
より一般的になるが，私は定数と言えども文字を増やして話をややこしくしたく
なかったのであえて単なる整数とした．エネルギーの単位を適当に調整したと
思ってもらってもいい．いずれにせよ，結論は全く変わらない．
}．

  このケースで$Z$を計算してみよう．実に簡単である．
$$Z = 1 + \exp\left(-\frac{1}{kT}\right) + \exp\left(-\frac{2}{kT}\right) + \exp\left(-\frac{3}{kT}\right) + \cdots.$$
これは無限等比数列の和なので，簡単に整理できる．
$$Z = \frac{1}{1 - \exp\left(-\frac{1}{kT}\right)}.$$

  さてこの$Z$を見てほしい．これは何の関数だろうか．なんの定数だろうか．……まず，$Z$は$E$を一切含んでいない．それは当然だ，なぜなら$E$に値を代入して数列を作り，それを全部足したのだから．代入したのだから，式中に$E$はもう残っていない．だから$Z$は$E$に対しては定数である．また$Z$の式の中には$T$という値が残っている．つまり$Z$は$T$の関数なのだ．

  私たちの考えている$E$は，$p(z)$から構成したものなので，当然のことながら$z$の関数であった．しかしこの$E$はもう$Z$には残っていないので，$Z$は$z$の関数ではない．……おっとこれは言い過ぎかもしれない．もし温度が場所によって違うような系を考えているのなら，$T$が$z$の関数になるので$Z$は$z$の関数であると言えるだろう．しかし私はそういう複雑な状況を今回は想定していない．温度は系全体で共通な定数だと想定している．

\subsection{さらに註……というかもはや追記}

  私はこの説明において，あまりよく考えずにPRMLの流れに合わせて$kT/m$を$1$として
しまったが，温度は本当はそんなに軽く扱っていいものではない．いやそれを言ったら
質量だって適当にしてはいけないかもしれない．物理屋の意地があるので少々語ること
にする．

  まず温度だが，もし温度があまりに低いと分子はほとんどエネルギーを持てないので，
最寄りのポテンシャルエネルギーの低い場所にすぐに収まってしまって，そこから二度と
出てこない．もちろんzの初期値に恵まれて，そこから下る過程で高い運動エネルギーを
一時的には持てるかもしれないが，それも温度を考慮した速度$v$の取り直しによって，
結局失ってしまう．この場合，結局仮想粒子はごく狭い範囲をちまちま動くことしか
できず，それはつまり$z$の初期値に強く影響されたサンプリング，言い換えれば$p(z)$全体
をほとんど反映していないサンプリングとなる．これではもちろんいけないだろう．

  では温度を高くしてやればいいのか．確かにそうすればポテンシャルの丘が多少あって
も難なく飛び越えていくだろう．つまり仮想粒子は$p(z)$全体を十分に探検できるように
なる．……しかし話はそう単純ではないのだ．もし温度が不適切なほど高ければ，
もはや仮想分子は$p(z)$に影響されなくなる．なぜなら，自分が温度由来で与えられている
速度vに対して，ポテンシャルから与えられた加速度$a$が小さく，もはやノイズ程度にしか
ならないからだ．こうなってしまうと，どの$z$に対しても$p(z)$は同じ，みたいな系の
サンプリングをしたような結果にしかならず，これも不本意だろう．

  今度は質量$m$について考えてみよう．質量はポテンシャルエネルギーの傾きがどのくらい
加速度に変換されるかの比例定数である．分かりやすくするために極端な例を考えよう．
もし質量が無限大だったらどうだろう．そうとも，加速度$a$は常にゼロとなり，仮想粒子は
$p(z)$に影響されることなく，温度をベースに与えられた$v$のまま等速直進運動をすること
になる．これは温度が高すぎた場合によく似ている．これはダメだ！……では質量が$0$
にかなり近い小さい値だったらどうか．これは温度が$0$だった場合のようにふるまうことに
なり，これもダメだ．

  ……と言いたいところなのだが，速度の初期値や取り直しの際に，私のアドバイスに
従い運動エネルギー $\frac{1}{2} m v^2$ の期待値が $\frac{1}{2} kT \times (zの次元数)$ となるようにとって
いるのであれば，$m$が大きいときには$v$が小さくなり，$m$が小さいときには$v$が大きくなる
ので，加速度$a$のスケールと自動的に同じになる．ということで，質量$m$をどのくらいの
値にすべきかについては，あまり深刻に悩まずともよさそうである．

  さらに温度についてもいいニュースがある．私たちはポテンシャルエネルギーを定義
するときに $E(z) = - kT \cdot \ln{( p(z) \cdot Z )}$ としている．つまり温度が高い場合は，
相応にポテンシャルの起伏も激しくなるのだ．温度が低い場合は，それに合わせて
ポテンシャルの起伏もなだらかになるのだ．だからたぶん温度についてもそこまで
神経をとがらす必要はない．

  あとは$v$の取り直しの頻度についても考えてみよう．
$v$の取り直しが頻繁に起こる場合，これは仮想世界中の粒子数が非常に多くて過密であることを意味する．だからしょっちゅうぶつかっているわけだ．また，取り直しの頻度が少ない場合，これはかなり希薄な気体をシミュレーションしているということになる．

  まず頻度が低すぎるというのはよろしくない．なぜなら温度が安定しなくなるからだ．
そもそもそんなに希薄な気体では総分子数は相当に少ないということになるだろうが，
そんな系では統計力学的な温度の定義が通用しなくなる．統計力学は十分に多数の分子が
エネルギーをやり取りしているような状況を前提に組まれているのだ．だからこそ
やっかい極まりない揺らぎの問題を解消できている．それが通用しなくなるほど希薄だと
するとさまざまな前提が崩れて，今までの説明通りにはいかなくなってしまう恐れが
出てくる．……別の言い方をするなら，この場合の仮想粒子は，運動エネルギーと
ポテンシャルエネルギーの和が一定になるような運動を過剰に長く続けてしまう．これは
先に書いたように，温度一定とは異なる挙動である．

  ……では頻度をうんと上げるのはどうだろう．今度は$p(z)$がほとんど反映されなく
なってしまう．ポテンシャルエネルギーから加速度を決めてようやくその向きに進み
始めたところで，$v$がリセットされてしまえば，結局$p(z)$は仮想粒子の運動にほとんど
影響できなかったことになる．ということは上記で温度が大きすぎた場合の考察のような，
残念な結果になるだろう．


\section{余談}

  PRML下巻のp.264では\textgt{ハミルトン}を紹介しているが，私ならその場所に，
\textgt{ルートヴィッヒ・ボルツマン}を置くだろう．
興味があれば Wikipedia で彼について調べてみてほしい．

\section{Verlet法}

  ここでは私の個人的な趣味により，蛙跳び法ではなくてVerlet法を紹介したい．

  もっとも単純に加速度$a$から速度$v$と位置$z$を更新すると，こうなるだろう．ここでは
時間刻みを$h$としている（PRMLではイプシロンを使っている）．
\begin{eqnarray}
\label{kawai1}z &=& z + v \cdot h + \frac{1}{2} \cdot a \cdot h^2,\\
\label{kawai2}v &=& v + a \cdot h.
\end{eqnarray}

  これでも$h$が十分に小さければ悪くない精度が出せる．しかし$h$が小さいとそれだけシ
ミュレーション内の時間経過が遅くなる．そこで，次のようなテクニックがある．

  まず$z(t+h)$を以下のように書き下すことができる．
%
    $$z(t+h) = z(t) + v(t)\cdot h + \frac{1}{2} \cdot a(t) \cdot h^2 + \frac{1}{6} \cdot e(t) \cdot h^3 + \cdots.$$

ここで$e(t)$は，$z$の3階微分である．加速度の1階微分（ここでは時間微分）であると
いってもよい\footnote{%
$z$の3階微分や4階微分などには，速度や加速度といったような物理学的な名前が
ついていないがために，これらはゼロであると決めつける誤解があるが，そんなことは
断じてない．
もし$z$の3階微分が常にゼロなら，加速度は時間変化しないと言っているの
と同値であるが，これはおかしい．かかる力が位置によって変化し，その位置は運動に
よって時間変化しているのであれば，当然加速度も時間変化するに決まっている．した
がって一般性を失わないためには，3階微分や4階微分などを仮定しておく必要がある．
}%
．上記の$h$を$-h$に置き換えれば，以下の式が得られる．
$$z(t-h) = z(t) - v(t) \cdot h + \frac{1}{2} \cdot a(t) \cdot h^{2} - \frac{1}{6} \cdot e(t) \cdot h^{3} + \cdots.$$

そしてこの2式を辺々加える．
\begin{equation}\label{kawai3}
    z(t+h) + z(t-h) = 2 \cdot z(t) + a(t) \cdot h^2.
\end{equation}
整理すれば，以下を得られる．
\begin{equation}\label{kawai4}
    z(t+h) = 2 \cdot z(t) - z(t-h) + a(t) \cdot h^2.
\end{equation}

  この式\eqref{kawai4}は$z(t)$, $z(t-h)$, $a(t)$が既知であれば$z(t+h)$が
計算可能であることを
示している．つまり現在の座標，1ステップ前の座標，現在の加速度が分かれば，
1ステップ後の座標がわかるということだ．しかもその$z(t+h)$には $h^4$ オーダーの
誤差しか含んでいない．式\eqref{kawai1}では$h^3$オーダーの誤差を含んでいたので，これは格段
に精度がよい．また特筆すべきは，$z$を次々と計算していくにあたって，$v$の計算を必要
としていない．つまり$z$だけでどんどん計算を進めていくことができる．私たちは$z$が
ほしいのであって，$v$のログには興味がないから，これは（計算量的に）好都合である．

問題は，最初に$z(t)$のみならず$z(t-h)$が必要になることである．この$z(t-h)$は
$v$と$a$を使って計算していくしかない．具体的な方法は複数考えられるが，例えば，
最初に$h$をかなり小さくとって式\eqref{kawai1}と式\eqref{kawai2}で計算を10回くらい進める．そうすれば
かなり精度のよい$z(t=0)$と$z(t=10h)$が手に入るだろう．そのあとは$h$を10倍して
Verlet法を使っていくことができる．なおここで10倍という数値を挙げたが，実際の
ケースではどのくらいが適当なのかは，各種の条件によるだろう．たぶん$1/10$という
小さな準備ステップが必要になることはまずない．たいていは$1/2$とか$1/4$程度で十分に
足りる．

  ちなみにVerlet法では以下の式も同じように導出できる．
\begin{equation}\label{kawai5}
    v(t) = \frac{z(t+h) - z(t-h)}{2h}.
\end{equation}

これを使えば，$z(t+h)$と$z(t-h)$から$v(t)$を計算することができる．ただし式\eqref{kawai5}の精度は
特に高いというわけではない．誤差は$h^2$オーダーである．つまりこれは式\eqref{kawai2}と
同程度ということになる．しかし私たちの目的においては正確な$v$は必要ではないので
これはあまり気にする必要はない……と私は思う．

  もしどうしても高い精度の$v$が必要ならば，こういう式も作れる．
\begin{equation}\label{kawai6}
    v(t+h) = v(t-h) + 2 \cdot a(t) \cdot h.
\end{equation}

この式には$h^3$オーダーの誤差しか残ってないが，この式でやるためには$v$もステップ
ごとに計算していかなければいけないのでVerlet法のメリットは少なくなってくる．

  ちなみに高い精度の$v$が必要になるかもしれないのは，このシミュレーションがどの
くらい正確なのかを見積もらなければいけない場合である．一般にこの手のシミュレー
ションで精度を確認するときには運動エネルギーと位置エネルギーの和が保存している
かどうかを見る．それがうまくいっていないときは，たいていどこかがバグっている．
エネルギーを計算するには運動エネルギーを求める必要があり，その際には高精度
な$v$がほしくなるだろう．