文献読みメモ: Causal inference in economics and marketing(その1)

以下のペーパーを読みます。

Hal R. Varian. Causal inference in economics and marketing. Proceedings of the National Academy of Sciences, 113 (27) 7310-7315, 2016. https://www.pnas.org/content/113/27/7310
※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらご指摘ください。
次回:まだ
f:id:cookie-box:20200101101603p:plain:w60

あらゆる因果分析のキモは「反実仮想」、つまり、「もしその処置がなされていなかったら何が起きていただろうかと推測すること」ですか。


冒頭の例は、以下のようなものですよね。
あるサーフィン映画のテレビ広告の効果を測ろうとして、一人当たり興行収入 y_c を一人当たり広告費 x_cy_c = b x_c + e_c と回帰しても、これは上手くいかない。なぜなら、ハワイ州でのデータを用いて広告の効果が b = 10 と推定されたとしても、ノースダコタ州では人々はハワイ州の人々ほどサーフィンに興味があるとは思えず、同じ効果を見込めるとは思えない。実際ノースダコタ州のデータで回帰すると b=0.1 になっている。これは「サーフィンへの興味」という説明変数が足りなかったのだ。
…うーん、なんだか、この例で「都市によって広告の効果 b の値が食い違った → 説明変数が足りなかったからだ!」といわれても、「都市によって広告の効果 b が違った → 本当に都市によって広告の効果は違うのだ」と思う方が自然ではないですか? 実際、「テレビ広告の効果」って在宅の主婦や高齢者が多く住んでいるかとか土地柄とかに左右されそうですし…まあどうでもいいですが…。
それで、じゃあどうして b が食い違ってしまうのか数式で考えると、b={\rm cov}(x, y)/{\rm cov}(x,x) なので…そうなんでしたっけ?

f:id:cookie-box:20200101101614p:plain:w60

そうだよ。厳密には、切片がないと考えた場合は以下の左下のスライドの通り  a = (X^\top X)^{-1} X^\top Y になるし(下のスライドでは ba になってるけど; X の定義は右下のスライドだけど説明変数が1次元なら X はデータの個数の長さの縦ベクトルだね)、切片 a_0 があると考えた場合は X にもう一列すべて1であるような2列目を補えばよくて、

 a = \left( \begin{array}{cc} X^\top X & \sum X_i \\ \sum X_i & N \end{array} \right)^{-1} \left( \begin{array}{c} X^\top Y \\ \sum Y_i \end{array} \right) = \displaystyle \frac{1}{N X^\top X - (\sum X_i )^2} \left( \begin{array}{cc} N & -\sum X_i \\ -\sum X_i & X^\top X \end{array} \right) \left( \begin{array}{c} X^\top Y \\ \sum Y_i \end{array} \right)
これの1つ目の成分が「傾き」だから、
 \displaystyle \frac{N X^\top Y - (\sum X_i )(\sum Y_i )}{N X^\top X - (\sum X_i )^2}
これは {\rm cov}(x, y)/{\rm cov}(x,x) に他ならないよね。それで、このペーパーでは切片を考えてないのに  b = (X^\top X)^{-1} X^\top Y じゃなくて b={\rm cov}(x, y)/{\rm cov}(x,x) になってるけど、これは xy も平均がゼロになるように中心化してるから後者が前者と同じなんだよね(だから切片を考えてないんだけど)。

f:id:cookie-box:20200101101603p:plain:w60

ああ、確かに b={\rm cov}(x, y)/{\rm cov}(x,x) ですね。それで、以下より、b にバイアスがないのは {\rm cov}(x, e) = 0 のときのみであると…。

 \displaystyle b = \frac{{\rm cov}(x, xb + e)}{{\rm cov}(x, x)} = \frac{{\rm cov}(x, xb)}{{\rm cov}(x, x)} + \frac{{\rm cov}(x, e)}{{\rm cov}(x, x)} = b +  \frac{{\rm cov}(x, e)}{{\rm cov}(x, x)}
{\rm cov}(x, e) = 0 はそもそも b を求めるときの仮定ですよね。ゼロであるとしたものをひねり出すこの式変形には違和感があるんですが…まあいいです。7310ページ左列の一番最後には因果推論のテキストでよくいわれるあれがかいてありますね。つまり、広告の打ち方に介入せず、単に広告費から売り上げを推測したいだけならさっきのシンプルな回帰でも構わない。しかし、いま知りたいのは「広告の打ち方を変化させたときに売り上げがどう変わるか」なのでそれではいけないということです。…そうですね、極端な話、配給会社は面白い映画にはその面白さの度合いだけ広告に熱を入れるかもしれません。そして、人々は面白い映画にはその面白さの度合いだけ映画館に足を運ぶかもしれません。このとき、広告費と興行収入は比例してみえますが、実は人々は配給会社が打った広告とは全く関係なく映画を見に行っていただけかもしれません。この場合、広告費を2倍にしたところで興行収入が2倍になるとは思えません。

f:id:cookie-box:20200101101614p:plain:w60

その「映画の面白さ」みたいに、広告費にも売り上げにも影響するのが交絡変数(confounding variable)だね。これが存在するから上辺の関係だけみていても効果を測ることはできない。

f:id:cookie-box:20200101101603p:plain:w60

そしてその問題は「人のふるまい」を分析するときにはいつもついてまわる…ですか。先の例だと、投じられる広告費はマーケティング担当の何らかの意思決定に基づくはずだが、アナリストはそれをわからず、意思決定時のファクターがエラーターム e_c になってしまう、だから x_ce_c が無関係でなくなってしまう、と。アナリストを雇ったのであればもう少しコミュニケーションしてくれませんかね…。しかし、現実にはマーケティング担当だって把握していない要素があるでしょうが…。ともあれ、通常は(マーケティング担当が incompetent でない限り)効果が出やすい広告の打ち方というのをしますから、本質的に交絡変数が存在するということですね。交絡変数が存在する例として、他にも以下が挙げられています。

コントロールできるものコントロールしたいもの
肥料の量農作物の収穫量
教育収入
健康管理収入
なんというか、2番目と3番目の「教育」や「健康管理」というのも曖昧ですが…これらは丸めていってしまえば、「『コントロールできるもの』をたくさん費やされるサンプルは、元々恵まれてるんじゃないの?(だから単純に回帰すると overstate になる)」という例ですよね。逆の例ってありませんか?

f:id:cookie-box:20200101101614p:plain:w60

あくまで例だけど、ある病院で、個々のお医者さんが手術した患者の生存率のデータがあるとする。でも、腕のいい人ほど重篤な患者さんを任されるかもしれない。だから、「お医者さんの腕のよさ」と「患者の生存率」は思ったほど正の相関をもたないかもしれないし、下手すると相関が負になってくるかもしれない、みたいなのはどうかな。

f:id:cookie-box:20200101101603p:plain:w60

なるほど。「運転しない人ほど事故を起こさない」「コードをかかない人ほどバグを出さない」ってやつですね?

f:id:cookie-box:20200101101614p:plain:w60

そこまでいってないよ!?

f:id:cookie-box:20200101101603p:plain:w60

先の例の後に、ラテン語のフレーズが含まれた文が出てきますね。

"We want an answer to a ceteris paribus question, but our data were generated mutatis mutandis."
他の事情が同じ条件での処置の効果を測りたいのに、実データにおいては処置にしかるべき調整が加えられている、という意味でしょうが、この文章は英語がわかる人には小気味よい文章なんでしょうか? わかりません。それで、この続きでは、まず因果効果を測る黄金の方法としての controlled experiments を紹介しますが、これは現実的ではないので、経済学で観測データから因果効果を推定するのに用いられる以下の4つの手法を紹介するということですね。
  1. natural experiments
  2. instrumental variables
  3. regression discontinuity
  4. difference in differences
…話が長くて面倒なので Difference in Differences の節までとぶと、ここでは下の表のような状況で、キャンペーンの効果 y_{TA} - y_{TF} を測りたいのですね。ここで y_{TF} は、「キャンペーンを受けた処理群がキャンペーンを受けなかった世界線での、キャンペーン後の時点での値」とでもいえばいいでしょうか。
キャンペーン前キャンペーン後
対照群 y_{CB} y_{CA}
処理群 y_{TB} y_{TF}(観測不可)
 y_{TA}
それで、ここでは y_{TF} - y_{TB} = y_{CA} - y_{CB} という仮定を置きます。処理群の人たちがもしキャンペーンを受けなかったとしたら、変化は対照群の人たちと同じであったはずだ、という仮定ですね。これが正しい状況かは慎重に判断しないといけませんが…。ただ、その仮定さえ認めれば、キャンペーンの効果は y_{TA} - y_{TF} = y_{TA} - y_{TB} - (y_{CA} - y_{CB}) と、観測データのみから計算することができます。何のことはない、処理群の成長度から対照群の成長度を差し引いただけですね。だから「差の差」というわけです。変化の大きさではなく変化率が同じであったという仮定 y_{TF}/y_{TB} = y_{CA}/y_{CB} を置いてもいいです。その場合、キャンペーンの効果は y_{TA}/y_{TF} = (y_{TA}/y_{TB})/(y_{CA}/y_{CB}) になりますね。処理群の成長率を対照群の成長率で割ることになります。…しかし、何か共変量 x_{td} があるなら、 y_{TF} の推定に機械学習的な手法を用いるとよいかもしれないというような提案がありますね。対照群でモデルを学習して処理群に適用すればよいと。

(その2があれば)つづく

雑記: t分布の話

キャラクターは架空のものです。お気付きの点がありましたらご指摘いただけますと幸いです。
参考文献

  1. 24-3. 2標本t検定とは | 統計学の時間 | 統計WEB ― 2標本t検定の話があります。
  2. Student's t-distribution - Wikipedia ― t分布を発表された方の顔写真があります。
  3. http://www012.upp.so-net.ne.jp/doi/biostat/CT39/ttest.pdf ― 「t分布はどのようにして生まれたか」というのと「分散をプールするとはどういうことか」というのがわかりやすかったです。この記事のt分布の導入の流れはこの参考文献の流れに近いです。
  4. 意外と難しいT分布についての証明 - Qiita ― t分布の導出について全面的に参考にしました。
  5. 不偏分散と自由度n-1のカイ二乗分布 | 高校数学の美しい物語 ― 上の記事にもリンクがありますが不偏分散の分布の導出について全面的に参考にしました。

まとめ

  • 対応のない2標本(分散は同一で未知)の平均に差があるか検定するには自由度 n_0 + n_1 - 2 のt分布をつかう。
    • まず N(\mu_0, \sigma_0^2) から独立に n_0 個のデータを生成したら  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがう。
    • 未知の \sigma_0^2 を不偏標本分散で置き換えた  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} は自由度 n_0 - 1 のt分布にしたがう。
      •  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうことと、それとは独立に (n_0 - 1)S_0^2/\sigma_0^2 が自由度 n_0 - 1 のカイ2乗分布にしたがうことから導かれる。
    • 対応のない2標本があるとき、 \bigl( \bar{X_1} - \bar{X_0} - (\mu_1 - \mu_0) \bigr) / \sqrt{(1/n_0 + 1/n_1)\sigma_0^2}N(0,1) にしたがう。他方、(n_0 - 1)S_0^2/\sigma_0^2 + (n_1 - 1)S_1^2/\sigma_0^2 は自由度  n_0 + n_1 - 2 のカイ2乗分布にしたがう(カイ2乗分布の再生性)。これらを組み合わせると、自由度 n_0 + n_1 - 2 のt分布にしたがう検定統計量が構成される。

f:id:cookie-box:20200101101603p:plain:w60

日常生活の中で、手元にある対応のない2標本の平均に差があるか気になって、「標本0は N(\mu_0, \sigma_0^2) から独立に生成されていて、標本1は N(\mu_1, \sigma_0^2) から独立に生成されているだろう(\sigma_0 は不明)」と考えて、以下の検定統計量 T を用いて2標本t検定をすることがありますよね。但し、標本0のサイズを n_0, 標本平均を \bar{X_0}, 不偏標本分散を S_0^2 とし、標本1のサイズを n_1, 標本平均を \bar{X_1}, 不偏標本分散を S_1^2 とします。

 \displaystyle  T = \frac{\bar{X_1} - \bar{X_0}}{\sqrt{S^2 \left( \frac{1}{n_0} + \frac{1}{n_1}\right)}}, \quad S^2 = \frac{(n_0 - 1) S_0^2 + (n_1 - 1) S_1^2}{n_0 + n_1 - 2}
この T の実現値 t が棄却域にあるかどうかを調べればいいわけです。帰無仮説に応じた棄却域は以下です。ここで t_{a\%}(n) は自由度 n のt分布の a パーセント点です。
帰無仮説対立仮説棄却域(有意水準5%の場合)
 \mu_0 = \mu_1 \mu_0 \neq \mu_1 t < t_{2.5\%}(n_0 + n_1 - 2) または t_{97.5\%}(n_0 + n_1 - 2) < t
 \mu_0 \leqq \mu_1 \mu_0 > \mu_1t_{95\%}(n_0 + n_1 - 2) < t
自由度 n のt分布とは以下のような確率密度関数をもつ分布ですね。
 \displaystyle f(t) = \frac{\Gamma \left( \frac{n+1}{2} \right)}{\sqrt{n \pi} \,\Gamma \left( \frac{n}{2} \right)} \left( 1 + \frac{t^2}{n}\right)^{-\frac{n + 1}{2}}
でもよく考えると、どうしてこれで検定できるんでしょう? そもそもt分布というのもよく知らない分布ですし。どこからきた分布なんです?

f:id:cookie-box:20200101101614p:plain:w60

どこからきたというか、そこからきたんだよ。まず標本0にだけ注目しよう。いま標本0の母平均 \mu_0 の信頼区間を知りたいとする。標本0が N(\mu_0, \sigma_0^2) から独立に生成されていることより、 (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうから、もし母分散 \sigma_0^2 がわかっていれば、適当な信頼係数で「 (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0} がとりうる区間」を決め打って、それを母平均について解けばいい。信頼係数 0.95 なら、信頼区間は標準正規分布の97.5%点 z_{97.5\%} を用いて以下のようになるね。

 \displaystyle \bar{X_0}- z_{97.5\%} \sqrt{\frac{\sigma_0^2}{n_0}} \leqq \mu_0 \leqq \bar{X_0}+ z_{97.5\%} \sqrt{\frac{\sigma_0^2}{n_0}}
でも母分散 \sigma_0^2 がわからなかったらどうする?

f:id:cookie-box:20200101101603p:plain:w60

母分散 \sigma_0^2 がわからなかったら、ですか…。先ほどと同様に  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうことを利用しようとしても、\sigma_0^2 がわからないので区間が定まりませんね。…であれば、\sigma_0^2 を不偏標本分散 S_0^2 などで置き換えて  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} としてはいけないでしょうか? いえ、しかし、不偏標本分散 S_0^2 は標本の出方によって真の母分散 \sigma_0^2 よりも大きく出たり小さく出たりしますから、やはりそのまま置き換えても N(0, 1) にはしたがわないでしょう。不偏標本分散 S_0^2 が大きく出たり小さく出たりという分布をも考慮して、 (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} の分布の形を考え直す必要があるでしょうね…。

f:id:cookie-box:20200101101614p:plain:w60

それでいいよ。まず S_0^2 が大きく出たり小さく出たりという分布を出してみよう。

f:id:cookie-box:20200101101603p:plain:w60

えっ、S_0^2 のしたがう分布ですか? うーん…とりあえず展開しましょう。標本0内の i 番目のデータを X_{0, i} とかくことにしましょう( i = 1, \cdots, n_0 )。

 \displaystyle  \begin{split} S_0^2 &= \frac{1}{n_0 -1}\sum_{i=1}^{n_0} (X_{0, i} - \bar{X_0})^2 \\ &= \frac{1}{n_0-1}\sum_{i=1}^{n_0} X_{0, i}^2 -\frac{2 \bar{X_0}}{n_0 -1}\sum_{i=1}^{n_0} X_{0, i} + \frac{n_0}{n_0 -1} \bar{X_0}^2 \\ &= \frac{1}{n_0 -1}\sum_{i=1}^{n_0} X_{0, i}^2 - \frac{n_0}{n_0 -1} \bar{X_0}^2 \end{split}
第1項(第1項の中に n_0 項分あるわけですが)の分布と第2項の分布は2乗する前は正規分布ですから、変数変換すればこれらの分布も求ま…あ、第1項と第2項は独立ではないですね、負けました。

f:id:cookie-box:20200101101614p:plain:w60

急に投了しないで! …まあ確かに第2項は厄介だ。以下にもトリッキーってかいてあるしね。

つまり、座標変換で第2項を消すことにする。まず X_{0, i} がしたがう分布が N(\mu_0, \sigma_0^2) ではなくて N(0,1) の場合を考える。このとき Y_0 = Q X_0 という変換をする( X_0Y_0 は縦に n_0 個並んだベクトルのイメージね)。あ、この Q は直交行列にしたいのね。変換前後で2乗和を保ちたいから。じゃあどんな Q を選べばいいかというと、変換後の第1成分が  Y_{0,1} = (1/\sqrt{n_0}) \sum_{i=1}^{n_0}  X_{0, i} となるようにしたい。つまり、Q の1行目の成分をすべて 1/\sqrt{n_0} にする。だってそうすれば、
 \displaystyle  S_0^2 = \frac{1}{n_0 -1}\sum_{i=1}^{n_0} X_{0, i}^2 - \frac{n_0}{n_0 -1} \bar{X_0}^2 = \frac{1}{n_0 -1}\sum_{i=1}^{n_0} Y_{0, i}^2 - \frac{1}{n_0 -1} \sum_{i=1}^{n_0} \bar{X_{0,i}}^2 = \frac{1}{n_0 -1}\sum_{i=2}^{n_0} Y_{0, i}^2
と、厄介だった第2項が消えるからね。実際にこんな直交行列 Q を手に入れるには、単位行列の1行目の成分をすべて 1/\sqrt{n_0} にしてからグラムシュミットの直交化法をすればいいだろう。

f:id:cookie-box:20200101101603p:plain:w60

なんと…第1項は正規分布にしたがう確率変数の2乗和だから、直交変換をしても構わない。適当な直交変換を選べば、厄介な第2項が消える…。そして変換後の Y_0 = Q X_0 のしたがう分布はというと、多変量標準正規分布にしたがう確率ベクトルを直交変換しても多変量標準正規分布にしたがいますね。であれば (n_0 -1)S_0^2 = \sum_{i=2}^{n_0} Y_{0, i}^2 のしたがう分布は正規分布の二乗和がしたがう分布であり、つまり、自由度 n_0 -1 のカイ2乗分布ですか(これも導出すべきですが尺の都合で割愛します)。

f:id:cookie-box:20200101101614p:plain:w60

いまは X_{0, i} がしたがう分布が N(0,1) の場合を考えたけど、一般に X_{0, i}N(\mu_0, \sigma_0^2) にしたがう場合は (n_0 -1)S_0^2 / \sigma_0^2 が自由度 n_0 -1 のカイ2乗分布にしたがうね。

f:id:cookie-box:20200101101603p:plain:w60

…いま  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} がしたがう分布を知ろうとして、S_0^2 のしたがう分布がわかった、という状態ですよね。では次はどうすればいいんでしょう?

f:id:cookie-box:20200101101614p:plain:w60

じゃあ分子の  \bar{X_0} のしたがう分布は?

f:id:cookie-box:20200101101603p:plain:w60

\bar{X_0} 自体は、 (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうのでしょう? まず母分散 \sigma_0^2 が既知の場合はといって副部長自身が言っていましたよ。しかし、\bar{X_0}S_0^2 がしたがう分布がわかったからといって、 (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} のしたがう分布はわかりませんよね。結局、\bar{X_0} がどんな値をとるときは S_0^2 がどんな値をとりうるという相関構造がわからなければ…\bar{X_0}S_0^2 が独立でもない限り…。

f:id:cookie-box:20200101101614p:plain:w60

\bar{X_0}S_0^2 は独立だよ?

f:id:cookie-box:20200101101603p:plain:w60

え? 直感的には標本平均と標本分散が独立には思えな…。

f:id:cookie-box:20200101101614p:plain:w60

S_0^2 のしたがう分布を導出したときに、Y_0 = Q X_0 という変換をしたけど、変換後の Y_0 も多変量標準正規分布にしたがうから各成分は独立だ。そして、変換後の第1成分は  Y_{0,1} = (1/\sqrt{n_0}) \sum_{i=1}^{n_0} X_{0, i} = \sqrt{n_0} \bar{X_0} で、変換後の第2~n成分の2乗和は \sum_{i=2}^{n_0} Y_{0, i}^2 = (n_0 -1)S_0^2 だった。

f:id:cookie-box:20200101101603p:plain:w60

あっ…\bar{X_0}S_0^2 は独立ということになりますね…。Y_0 = Q X_0 という変換が標本平均と標本分散を振り分けていたとは…。であれば、いま知りたい  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} のしたがう分布は計算すればよいです。

 \begin{split} \displaystyle P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{S_0^2/n_0}} < t \right) &= P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{\sigma_0^2/n_0}} \cdot \sqrt{\frac{\sigma_0^2}{(n_0-1) S_0^2}}< \frac{t}{\sqrt{n_0 - 1}} \right) \\ &= \iint_{\frac{a}{\sqrt{b}} < \frac{t}{\sqrt{n_0 - 1}} } z(a) \chi (b ; n_0 -1) da \, db \\ &= \iint_{\frac{a}{\sqrt{b}} < \frac{t}{\sqrt{n_0 - 1}} } \frac{e^{- a^2/2}}{\sqrt{2 \pi}} \frac{b^{(n_0 - 1)/2 - 1}e^{-b/2}}{2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} da \, db \end{split}
この変数を (a, b) から  (b, c=a/\sqrt{b}) に変換します。統計検定1級の教科書の1章でこれっぽい変数変換をちょいちょいやってるので何とかなるでしょう。この変換のヤコビアンは以下です。
 \begin{cases} a = h_0 (b, c) = \sqrt{b} c & \\ b = h_1 (b, c) = b & \end{cases}
 \left| \begin{array}{cc} \partial_b h_0(b, c) & \partial_c h_0(b, c) \\ \partial_b h_1(b, c) & \partial_c h_1(b, c) \end{array} \right| = \left| \begin{array}{cc} \frac{1}{2} b^{-1/2} c & \sqrt{b} \\ 1 & 0 \end{array} \right| = - \sqrt{b}
なので、このヤコビアンの絶対値をかければ変数変換できますね。
 \begin{split} \displaystyle P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{S_0^2/n_0}} < t \right) &= \int_{b=0}^{\infty} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{e^{- b c^2/2}}{\sqrt{2 \pi}} \frac{b^{(n_0 - 1)/2 - 1}e^{-b/2}}{2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \sqrt{b} \, db \, dc \\ &= \int_{b=0}^{\infty} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{b^{n_0/2 - 1}e^{- b c^2/2-b/2}}{ \sqrt{2 \pi} \cdot 2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \, db \, dc \end{split}
なんだこれ…どうしろと…この記事を参照するとガンマ関数をつくるために、エクスポネンシャルの肩に注目して (b, c) から  (s=b c^2/2 + b/2, c) と変数変換するのですか。ではまたヤコビアンの計算です。
 \begin{cases} b = h_0 (s, c) = s/(c^2/2 + 1/2) & \\ c = h_1 (s, c) = c & \end{cases}
 \left| \begin{array}{cc} \partial_s h_0(s, c) & \partial_c h_0(s, c) \\ \partial_s h_1(s, c) & \partial_c h_1(s, c) \end{array} \right| = \left| \begin{array}{cc}  1/(c^2/2 + 1/2) & \ast \\ 0 & 1 \end{array} \right| = \displaystyle \left( \frac{c^2}{2} + \frac{1}{2} \right)^{-1}
これを適用すると…なるほど s積分でガンマ関数がつくれますね。
 \begin{split} \displaystyle P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{S_0^2/n_0}} < t \right) &= \int_{s=0}^{\infty} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{s^{n_0/2 - 1}e^{- s}}{ \sqrt{2 \pi} \cdot 2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \frac{ \left( \frac{c^2}{2} + \frac{1}{2} \right)^{-1} }{\left( \frac{c^2}{2} + \frac{1}{2} \right)^{n_0/2 - 1}}\, ds \, dc \\ &= \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{ \Gamma \bigl( n_0/2 \bigr) }{ \sqrt{\pi} \cdot 2^{n_0/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \left( \frac{c^2}{2} + \frac{1}{2} \right)^{ - n_0/2} \, dc \\ &= \frac{\Gamma \bigl( n_0/2 \bigr)}{\sqrt{\pi} \, \Gamma \bigl( (n_0 - 1)/2 \bigr)} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \left( c^2 + 1 \right)^{ - n_0/2 } \, dc \\ &= \frac{\Gamma \bigl( n_0/2 \bigr)}{\sqrt{(n_0 - 1) \pi} \, \Gamma \bigl( (n_0 - 1)/2 \bigr)} \int_{c' < t} \left( \frac{{c'}^2}{n_0 - 1} + 1 \right)^{ - n_0/2} \, dc' \end{split}
最後 c' = \sqrt{n_0 - 1} \, c とおきました。そしてこれは自由度 n_0 - 1 のt分布ですね…。

f:id:cookie-box:20200101101614p:plain:w60

ここまでたどり着いたことになるね。

定理1.標本0( サイズ n_0, 標本平均 \bar{X_0}, 不偏標本分散 S_0^2 )があるとき、以下の T_0 は自由度 n_0 - 1 のt分布にしたがう。
 \displaystyle T_0 = \frac{\bar{X_0} - \mu_0}{\sqrt{S_0^2 / n_0}}

f:id:cookie-box:20200101101603p:plain:w60

t分布が何者かはわかりました。まさしくその T_0 がしたがう分布なのですね。しかし、本当に検定したいのは標本0と標本1の平均が同じかどうかです。

f:id:cookie-box:20200101101614p:plain:w60

だったら調べるべき \mu_1 - \mu_0 をつくろう。いま以下がわかっていた。

  •  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがう。
  •  (\bar{X_1}-\mu_1) / \sqrt{\sigma_0^2/n_1}N(0,1) にしたがう。
ここから \mu_1 - \mu_0 をつくるには以下のように変形して正規分布の再生性をつかえばいい。
  •  \bar{X_0}N(\mu_0,\sigma_0^2/n_0) にしたがう。
  •  \bar{X_1}N(\mu_1,\sigma_0^2/n_1) にしたがう。
  •  \bar{X_1} - \bar{X_0}N(\mu_1 - \mu_0, \sigma_0^2/n_0 + \sigma_0^2/n_1) にしたがう。
  •  \bigl( \bar{X_1} - \bar{X_0} - (\mu_1 - \mu_0) \bigr) / \sqrt{(1/n_0 + 1/n_1)\sigma_0^2}N(0, 1) にしたがう。
もし \sigma_0^2 が既知なら、この信頼区間はただちに求まる。
 \displaystyle \bar{X_1} - \bar{X_0} - z_{97.5\%} \sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \sigma_0^2} \leqq \mu_1 - \mu_0 \leqq \bar{X_1} - \bar{X_0} + z_{97.5\%} \sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \sigma_0^2}
でも実際には \sigma_0^2 はわからない。だから S_0^2, S_1^2 で代用するしかない。いま上の議論より (n_0 - 1)S_0^2/\sigma_0^2 は自由度 n_0 - 1 のカイ2乗分布にしたがい、同様に (n_1 - 1)S_1^2/\sigma_0^2 も自由度 n_1 - 1 のカイ2乗分布にしたがう。カイ2乗分布もまた再生性をもつから、 (n_0 - 1)S_0^2/\sigma_0^2 + (n_1 - 1)S_1^2/\sigma_0^2 は自由度  n_0 + n_1 - 2 のカイ2乗分布にしたがう。このシチュエーションでさっきの議論と同様にt分布にたどり着くようにしたいなら以下のようにすればいい。今回はt分布の自由度も  n_0 + n_1 - 2 になるね。
 \displaystyle P \left( \frac{\bar{X_1} - \bar{X_0}-(\mu_1 - \mu_0)}{\sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \sigma_0^2}} \cdot \sqrt{\frac{\sigma_0^2}{(n_0-1) S_0^2 + (n_1-1) S_1^2}}< \frac{t}{\sqrt{n_0 + n_1 - 2}} \right)
これを少し整理すればこうだ。帰無仮説 \mu_0 = \mu_1 のもとで2標本t検定の検定統計量に一致するよね。
 \displaystyle P \left( \frac{\bar{X_1} - \bar{X_0}-(\mu_1 - \mu_0)}{\sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \frac{(n_0-1) S_0^2 + (n_1-1) S_1^2}{n_0 + n_1 - 2} }}< t \right)

つづかない

雑記: numpy.fft の話

キャラクターは架空のものです。お気付きの点がありましたらご指摘いただけますと幸いです。全体的に離散フーリエ変換の備忘メモであり、高速フーリエ変換の原理の話は全くありません。
参考文献

  1. Discrete Fourier Transform (numpy.fft) — NumPy v1.17 Manual

まとめ

  • numpy.fft.fft は、長さ n の入力系列 \{a_m\}_{m=0}^{n-1} を渡すと、これを (4) 式のように分解したときの係数、つまり (4) 式の右辺の各列ベクトルの成分の和で表したときの係数 \{A_k\}_{k=0}^{n-1} を返してくれる。それが実は (3) 式で計算できるので実際は (3) 式をやってくれている(高速なアルゴリズムで)。
    • 入力系列 \{a_m\}_{m=0}^{n-1}numpy.fft.fft に渡して出た出力系列 \{A_k\}_{k=0}^{n-1}numpy.fft.ifft に渡すと元の入力系列 \{a_m\}_{m=0}^{n-1} に戻る。単に (4) 式をやってくれているだけである。元に戻ることは計算すればわかる。
    • 入力系列 \{a_m\}_{m=0}^{n-1}numpy.fft.fft に渡して出た出力系列 \{A_k\}_{k=0}^{n-1}(4) 式の右辺の各列ベクトルの成分の和で表したときの係数なので、周波数 0, 1/n, \cdots, (n-1)/n の波の振幅と解釈できる。ただし、整数でしか値を取らない波なので、周波数には整数を足し引きしても値は変わらない。なので、n/2 以上の周波数成分からは 1 を差し引いて周波数の取りうる範囲を -n/2, \cdots, -1,0,1, \cdots, n/2-1 (これは n が偶数のとき)と 0 を中心にしてもよい。こちらの流儀で numpy.fft.fft の出力系列がそれぞれどんな周波数の波の振幅と解釈できるかを取得できるのが numpy.fft.fftfreq である。入力系列のサイズ n を渡すだけである。
    • 入力系列 \{a_m\}_{m=0}^{n-1} が時間 \Delta t ごとにサンプリングされたものであるとき、これを反映した周波数にしたいならば、\exp(2 \pi i m k/n) = \exp(2 \pi i m \Delta t k/(n \Delta t )) より、m ごとの波を m \Delta t ごとの波にする代わりに周波数 k/nk/(n \Delta t) にすればよい。numpy.fft.fftfreq\Delta t も一緒に渡せばそれも踏まえて \{A_k\}_{k=0}^{n-1} はそれぞれどんな周波数の波の振幅なのかを返してくれる。もっとも、numpy.fft.fftnumpy.fft.ifft(3) 式と (4) 式をやるだけなので \Delta t を気にすることはない。\{a_m\}_{m=0}^{n-1}\{A_k\}_{k=0}^{n-1} をグラフにプロットしたときの横軸の目盛にのみ \Delta t が関わってくる。

f:id:cookie-box:20200101101603p:plain:w60

Discrete Fourier Transform (numpy.fft) — NumPy v1.17 Manual によると、numpy.fft における離散フーリエ変換と逆離散フーリエ変換の定義は以下であるとありますね。

\displaystyle A_k = \sum_{m=0}^{n-1} a_m \exp \left\{ - 2 \pi i \frac{mk}{n} \right\} \qquad k = 0, \cdots, n-1 \tag{1}
\displaystyle a_m = \frac{1}{n}\sum_{k=0}^{n-1} A_k \exp \left\{2 \pi i \frac{mk}{n} \right\} \qquad m = 0, \cdots, n-1 \tag{2}
行列でかくと以下でしょうか。
 \left( \begin{array}{c} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \end{array} \right) = \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ 1 & e^{-2 \pi i 1 \cdot 1 / n} & \cdots & e^{-2 \pi i (n-1) \cdot 1/ n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-2 \pi i 1 \cdot (n-1)/ n} & \cdots & e^{-2 \pi i (n-1) \cdot (n-1)/ n} \end{array} \right) \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) \tag{3}
 \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) = \displaystyle \frac{1}{n} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ 1 & e^{2 \pi i 1 \cdot 1 / n} & \cdots & e^{2 \pi i (n-1) \cdot 1/ n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{2 \pi i 1 \cdot (n-1)/ n} & \cdots & e^{2 \pi i (n-1) \cdot (n-1)/ n} \end{array} \right) \left( \begin{array}{c} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \end{array} \right) \tag{4}
これらは、元の信号をサンプリングした系列 \{a_m\}_{m=0}^{n-1} を周波数成分系列 \{A_k\}_{k=0}^{n-1} にする式と、その逆向きの式ということですよね…これって (1) で変換して (2) で逆変換したら元に戻ってくるんですか?

\exp(\pm 2 \pi i mk /n)kk-n としても値は変わらないので、k の動く範囲を k = 0, \cdots, n-1 とする代わりに、n/2 以上については n を差し引いて k = -n/2, \cdots, -1,0,1, \cdots, n/2-1 としてもよい(但しこれは n が偶数のとき;n が奇数の場合は k = -(n-1)/2, \cdots, (n-1)/2 とできる)。このとき (1), (2), (3), (4) は以下のようになる。A_{-n/2}ナイキスト周波数(元の信号をサンプリングした周波数の半分の周波数)の周波数成分に対応する振幅である。
\displaystyle A_k = \sum_{m=0}^{n-1} a_m \exp \left\{ - 2 \pi i \frac{mk}{n} \right\} \qquad k = -n/2, \cdots, -1,0,1, \cdots, n/2-1 \tag{1'}
\displaystyle a_m = \frac{1}{n}\sum_{k=-n/2}^{n/2-1} A_k \exp \left\{2 \pi i \frac{mk}{n} \right\} \qquad m = 0, \cdots, n-1 \tag{2'}
 \left( \begin{array}{c} A_{-n/2} \\ \vdots \\ A_{-1} \\ A_0 \\ A_1 \\ \vdots \\ A_{n/2-1} \end{array} \right) = \left( \begin{array}{cccc} 1 & e^{-2 \pi i 1 \cdot (-n/2) / n} & \cdots & e^{-2 \pi i (n-1) \cdot (-n/2) / n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-2 \pi i 1 \cdot (-1) / n} & \cdots & e^{-2 \pi i (n-1) \cdot (-1)/ n} \\ 1 & 1 & \cdots & 1 \\ 1 & e^{-2 \pi i 1 \cdot 1 / n} & \cdots & e^{-2 \pi i (n-1) \cdot 1/ n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-2 \pi i 1 \cdot (n/2-1)/ n} & \cdots & e^{-2 \pi i (n-1) \cdot (n/2-1)/ n} \end{array} \right) \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) \tag{3'}
 \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) = \displaystyle \frac{1}{n} \left( \begin{array}{ccccccc} 1 & \cdots & 1 & 1 & 1 & \cdots & 1 \\ e^{2 \pi i 1 \cdot (-n/2) / n} & \cdots & e^{2 \pi i 1 \cdot (-1) / n} & 1 & e^{2 \pi i 1 \cdot 1 / n} & \cdots & e^{2 \pi i 1 \cdot (n/2 - 1)/ n} \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ e^{2 \pi i (n-1) \cdot (-n/2) / n} & \cdots & e^{2 \pi i (n-1) \cdot (-1) / n} & 1 & e^{2 \pi i (n-1) \cdot 1 / n} & \cdots & e^{2 \pi i (n-1) \cdot (n/2-1)/ n} \end{array} \right) \left( \begin{array}{c} A_{-n/2} \\ \vdots \\ A_{-1} \\ A_0 \\ A_1 \\ \vdots \\ A_{n/2-1} \end{array} \right) \tag{4'}
f:id:cookie-box:20200101101614p:plain:w60

(1)(2) に代入したら戻ってくるね。

\displaystyle \begin{split} a_m &= \frac{1}{n}\sum_{k=0}^{n-1} \left( \sum_{m'=0}^{n-1} a_{m'} \exp \left\{ - 2 \pi i \frac{m' k}{n} \right\} \right) \exp \left\{2 \pi i \frac{mk}{n} \right\} \\ &= \frac{1}{n} \sum_{k=0}^{n-1} \left( a_0 + a_1 e^{-2\pi i k/n}  + a_2 e^{-2\pi i (2k)/n} + \cdots + a_{n-1} e^{-2\pi i (n-1)k/n} \right) e^{2 \pi i mk/n} \\ &= \frac{1}{n} \sum_{k=0}^{n-1} \left( a_0 e^{-2\pi i (-m)k/n} + a_1 e^{-2\pi i (1-m)k/n} + a_2 e^{-2\pi i (2-m)k/n} + \cdots + a_{n-1} e^{-2\pi i (n-1-m)k/n} \right) \\ &= \frac{a_0}{n} \frac{1 - e^{-2\pi i (-m)}}{1 - e^{-2\pi i (-m)/n}} + \frac{a_1}{n} \frac{1 - e^{-2\pi i (1-m)}}{1 - e^{-2\pi i (1-m)/n}} + \cdots + \frac{a_m}{n} n + \cdots + \frac{a_1}{n} \frac{1 - e^{-2\pi i (n-1-m)}}{1 - e^{-2\pi i (n-1-m)/n}} \\ &= a_m \end{split}
というか、(4) の左辺のベクトル(入力系列)を、(4) の右辺の行列の各列の基底ベクトルで表現したときの係数が  A_0, \, A_1, \, \cdots, \, A_{n-1} だからね。(1) をして (2) をしたら元に戻ったというより、(2) のような分解の係数を与えるのが (1) なんだよね。

f:id:cookie-box:20200101101603p:plain:w60

「元に戻ったというより元々そういう係数を出していたのだ」といわれても、(1) で係数が出てくること自体が非自明なんですがそれは…。まあ御託はいいです。とりあえず numpy.fft を使ってみましょう。とりあえず正弦波からサンプリングした系列を numpy.fft.fft して、再び numpy.fft.ifft してみればよいでしょう。

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 9, 3
rcParams['font.family']='Ume Hy Gothic O5'
rcParams['font.size'] = 15

T = 4
f = 1 / T
dt = 0.5

x = np.arange(-2, 2, dt)
y = np.sin(2 * np.pi * f * x)

X = np.fft.fftfreq(y.size, d=dt)
Y = np.fft.fft(y)

y_ = np.fft.ifft(Y)

print('周期   {}'.format(T))
print('周波数 {}'.format(f))
print('サンプル数             {}'.format(y.size))
print('サンプリング時間間隔   {}'.format(dt))
print('サンプリングレート     {}'.format(1/dt))

fig, ax = plt.subplots(5, 1, figsize=(7, 11))

ax[0].scatter(x, y)
ax[0].set_ylabel('元の信号')

ax[1].scatter(X, Y.real)
ax[1].set_ylabel('fftの実部')

ax[2].scatter(X, Y.imag)
ax[2].set_ylabel('fftの虚部')

ax[3].scatter(x, y_.real)
ax[3].set_ylabel('ifftの実部')

ax[4].scatter(x, y_.imag)
ax[4].set_ylabel('ifftの虚部')

plt.subplots_adjust(hspace=0.4)
plt.show()
周期   4
周波数 0.25
サンプル数             8
サンプリング時間間隔   0.5
サンプリングレート     2.0

f:id:cookie-box:20200128123638p:plain:w400
…確かに ifft の実部に元の正弦波が復元されたようにみえますね。みえますが、しかし、何が起きたのかよくわかりません…。numpy.fft さんは -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツの周波数成分の振幅を出してくださったようですが、そもそも私はそのようにお願いした覚えがなく…。

f:id:cookie-box:20200101101614p:plain:w60

みづらいけど、-0.25 ヘルツの振幅が -4i、+0.25 ヘルツの振幅が 4i になっているね。確かに周期4の正弦波が復元されるね。

 \displaystyle  \begin{split} & \frac{1}{8} \left( -4i e^{2 \pi i m (-0.25)} + 4i e^{2 \pi i m (0.25)}\right) \\ &= \frac{1}{8} \left( -4i \cos \left(-\frac{m \pi}{2} \right) + 4 \sin \left(-\frac{m \pi}{2} \right) + 4i \cos \left(\frac{m \pi}{2} \right) - 4 \sin \left(\frac{m \pi}{2} \right) \right) \\ &= - \sin \left(\frac{m \pi}{2} \right) \end{split}
部長が最初につくった正弦波と符号がずれているけど、部長は元の信号の左端を -2 にしているから、部長の座標の取り方に直せばこれを2だけマイナス側にずらして \displaystyle - \sin \left(\frac{(m+2) \pi}{2} \right) = \sin \left(\frac{m \pi}{2} \right) だね。

f:id:cookie-box:20200101101603p:plain:w60

ああ、numpy.fft さんは私が渡した系列の1つ目を時刻ゼロにおける信号の値と受けとっているのですか? …ん? それでは、2つ目の値は何であると受け取っているのでしょう? numpy.fft.fft にはサンプリング系列のみしか渡しておらず、どの時刻の値であるかなどは渡していません。

f:id:cookie-box:20200101101614p:plain:w60

…numpy.fft.fft は、入力系列のベクトルを (4) の右辺の行列の各列の成分ごとに分解するだけだよ。だから、1つ目の値はベクトルの第1成分、2つ目の値はベクトルの第2成分、…くらいにしか思ってないだろう。それで出てくる系列は (4) の右辺の行列の各列ベクトルの成分の係数になる。各列ベクトルは周波数が 0, 1/n, \cdots, (n-1)/n である波の m = 0, \cdots, n-1 での値になっているから、出力系列は周波数 0, 1/n, \cdots, (n-1)/n の周波数成分の振幅と解釈してもいい。

f:id:cookie-box:20200101101603p:plain:w60

では、上の例は n=8 ですから、周波数 0, 1/8=0.125, …, 7/8=0.875 ヘルツの成分の振幅になって…いませんね。出ているのは -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツの成分の振幅ですよ?

f:id:cookie-box:20200101101614p:plain:w60

2点あって、まず、m が整数しか取らない世界では \exp(2 \pi i m k/n)k/n に整数を足し引きしても値は変わらない。だから、0, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8 ヘルツの成分を、4/8, 5/8, 6/8, 7/8, 0, 1/8, 2/8, 3/8 ヘルツと並び替えて、さらに -4/8, -3/8, -2/8, -1/8, 0, 1/8, 2/8, 3/8 ヘルツだと考えていい。描画するならこう並べ替えた方が低周波数成分が真ん中に集まってみやすいだろう。加えてもう1点は、もし元の入力系列に対して「これは時間 \Delta t ごとにサンプリングしたものだ」という思いを込めたいなら、\exp(2 \pi i m \Delta t k/(n \Delta t )) と、周波数側にも \Delta t を押し付けることができる。部長は \Delta t = 0.5 にしたから、これを反映して周波数を1/0.5倍=2倍にする。つまり、-8/8, -6/8, -4/8, -2/8, 0, 2/8, 4/8, 6/8 ヘルツになる。-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツと一致するね。これまで単位時間に起きていたと思っていたことを「やっぱり時間 0.5 の間に起きていたことにして」としたんだから、周波数は2倍になるよね。でも、それらは人間がグラフの軸目盛にかく数字の問題でしかない。入力の点数が同じ8点なら、\Delta t がどうあったって基底ベクトルは同じだからね(記事の下部の2つのケース;横軸の目盛と何ヘルツと解釈するかが違うだけで、絵は全く同じ)。まとめると、numpy.fft.fft の返り値は素朴には 0, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8 ヘルツに対応する振幅と受け取れるんだけど、この2点を反映すると何ヘルツに対応する振幅と受け取るべきかを取得できるのが、部長が上のスクリプトでよんでいる np.fft.fftfreq だね。

つづかない

numpy.fft.fft に8点の系列を渡したときの基底(左列が実部、右列が虚部)

各行が (4') の各列ベクトルに対応。

x = np.linspace(0.0, 7.0, 8)
x_ = np.arange(0.0, 8.0, 0.01) # 点だけだとよくわからないのでcos, sinカーブを描画する用
vec_f = np.sort(np.fft.fftfreq(8, d=1)) # 8点の系列を渡したときの周波数成分

fig, ax = plt.subplots(8, 2, figsize=(11, 13))
for i, f in enumerate(vec_f):
    ax[i, 0].plot(x_, np.cos(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 1].plot(x_, np.sin(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 0].scatter(x, np.cos(2 * np.pi * f * x))
    ax[i, 1].scatter(x, np.sin(2 * np.pi * f * x))
    ax[i, 0].text(-0.5, 0.75, '{} Hz'.format(f), transform=ax[i, 0].transAxes)
plt.subplots_adjust(hspace=0.4)
plt.show()

f:id:cookie-box:20200127231000p:plain:w720

numpy.fft.fft に8点の系列を渡してΔt=0.5と指定したときの基底(左列が実部、右列が虚部)

各行が (4') の各列ベクトルに対応。

x = np.linspace(0.0, 3.5, 8)
x_ = np.arange(0.0, 4.0, 0.01) # 点だけだとよくわからないのでcos, sinカーブを描画する用
vec_f = np.sort(np.fft.fftfreq(8, d=0.5)) # 8点の系列を渡してかつΔt=0.5としたときの周波数成分

fig, ax = plt.subplots(8, 2, figsize=(11, 13))
for i, f in enumerate(vec_f):
    ax[i, 0].plot(x_, np.cos(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 1].plot(x_, np.sin(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 0].scatter(x, np.cos(2 * np.pi * f * x))
    ax[i, 1].scatter(x, np.sin(2 * np.pi * f * x))
    ax[i, 0].text(-0.5, 0.75, '{} Hz'.format(f), transform=ax[i, 0].transAxes)
plt.subplots_adjust(hspace=0.4)
plt.show()

f:id:cookie-box:20200128124655p:plain:w720

ベイズ統計の理論と方法: ノート5.2節その1(平均場近似)

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらコメント等でご指摘いただけますと幸いです。

ベイズ統計の理論と方法

ベイズ統計の理論と方法

書籍サポートページ
f:id:cookie-box:20200101101603p:plain:w60

…データから混合指数型分布 p(x|w) = \sum_{k=1}^K a_k v(x) \exp \bigr( f (b_k) \cdot g(x) \bigl) を推定するのに、

  • k の目が出る確率が a_k v(x) \exp \bigr( f (b_k) \cdot g(x) \bigl) であるようなサイコロ」
を仮想的に考えて、
  • n 個のデータ x^n を受け取ったとき、パラメータ a_{1:K}, b_{1:K} は何であったろうか」という問題を
  • n 個のデータ x^n を受け取ったとき、パラメータ a_{1:K}, b_{1:K} と各データのサイコロの目  y^n は何であったろうか」という問題に
読み換えたのですよね? そうすることで事後分布 P(y^n, a_{1:K}, b_{1:K} | x^n) を、独立な分布の積である試行分布 q(y^n) r(a_{1:K}, b_{1:K}) で平均場近似することが実現できるというのはわかります。わかりますが、そもそもそんなサイコロを考えなかったらどうなるんです?

f:id:cookie-box:20200101101614p:plain:w60

…まずサイコロを導入した場合はどうなったの?

f:id:cookie-box:20200101101603p:plain:w60

えっと、そもそも、p(w_1, w_2)= \exp \bigl( -\beta H(w_1, w_2) \bigr) / Z(\beta) という形で表される目標分布を、独立な分布の積である試行分布 q(w_1, w_2) = q_1(w_1)q_2(w_2) で(カルバック・ライブラー情報量  D(q || p) の意味で)最もよく近似したとき、q_1(w_1)q_2(w_2)p(w_1, w_2) の平均場近似とよぶのですね。


それではどのような q_1(w_1)q_2(w_2)p(w_1, w_2) の平均場近似なのか考えると、これらのカルバック・ライブラー情報量は具体的に、
 D(q || p) = \int q_1(w_1)q_2(w_2) \log \bigl( q_1(w_1)q_2(w_2)/p(w_1, w_2) \bigr) dw
 \quad = \int q_1(w_1)q_2(w_2) \log q_1(w_1)q_2(w_2) dw - \int q_1(w_1)q_2(w_2) \log p(w_1, w_2) dw
 \quad = \int q_1(w_1)q_2(w_2) \log q_1(w_1)q_2(w_2) dw + \int q_1(w_1)q_2(w_2) \bigl( \beta H(w_1, w_2) + \log Z(\beta) \bigr) dw
 \quad = \int q_1(w_1)q_2(w_2) \log q_1(w_1)q_2(w_2) dw + \beta \int q_1(w_1)q_2(w_2) H(w_1, w_2) dw - \beta F(\beta)
 \quad = \int q_1(w_1) \log q_1(w_1) dw_1 + \int q_2(w_2) \log q_2(w_2) dw_2 + \beta \int q_1(w_1)q_2(w_2) H(w_1, w_2) dw - \beta F(\beta)
 \quad \geqq 0
となりますから(ただし、 F(\beta) = - \beta^{-1} \log Z(\beta) は自由エネルギーです)、これを最小化する q_1(w_1)q_2(w_2)p(w_1, w_2) の平均場近似ということになります。等号成立は p(w_1, w_2) = q_1(w_1)q_2(w_2) のときですが、いま、試行確率分布は成分どうしが独立という制約がありますし、この等号が達成できる保証はありません。ところで、上式の下から2行目の、試行分布に依存しない - \beta F(\beta) を除いた部分を、平均場自由エネルギー \beta \hat{F}(\beta) とよぶとのことです。もし試行分布が目標分布を実現できれば真の自由エネルギーに等しくなる量ですが、一般的には真の自由エネルギー以上になる量ですね。最小化すべき目的関数の別表現といった感じがします。それで、もし \beta \hat{F}(\beta) が最小値になっているならば、この最小化問題のラグランジュ関数
 \int q_1(w_1) \log q_1(w_1) dw_1 + \int q_2(w_2) \log q_2(w_2) dw_2 + \beta \int q_1(w_1)q_2(w_2) H(w_1, w_2) dw_1 dw_2
\quad + \lambda_1 ( \int q_1(w_1)dw_1 - 1 ) + \lambda_2 ( \int q(w_2)dw_2 - 1 )
q_1, \, q_2 に関する1次の変化分が恒等的に 0 になっているはずです。このことから、
 \quad q_1(w_1) \propto \exp \bigl( - \beta \int q_2(w_2) H(w_1, w_2) dw_2 \bigr) \propto \exp \bigl(\int q_2(w_2) \log p(w_1, w_2) dw_2 \bigr)
 \quad q_2(w_2) \propto \exp \bigl( - \beta \int q_2(w_1) H(w_1, w_2) dw_1 \bigr) \propto \exp \bigl( \int q_2(w_1) \log p(w_1, w_2) dw_1 \bigr)
が要請され、これを自己無矛盾条件とよぶと。それで自己無矛盾条件を満たす試行分布をみつけよということですが…一般にラグランジュ関数の停留点は最小点である保証はないのであくまで最小点であることの必要条件になると思うんですが、この場合も極大点であることはないんでしょうか。そもそも必ず最小点があるんでしょうか。テキストの雰囲気的には「自己無矛盾条件を満たす試行分布の中には必ず平均場近似はある」といった感じですが。
あと、平均場自由エネルギー \beta \hat{F}(\beta) も何がエネルギーなのかわかりませんね。いうほど通常の自由エネルギーも何のことかわかりませんが。

f:id:cookie-box:20200101101614p:plain:w60

自由エネルギー \beta F(\beta) = -\log Z(\beta) って、確率分布が p(w_1, w_2)= \exp \bigl( -\beta H(w_1, w_2) \bigr) / Z(\beta) の形でかかれていることを前提としていて、このとき点 (w_1, w_2) の選択情報量(本当は  p(w_1, w_2) は確率じゃなくて確率密度だから選択情報量とはいわないけど)は  \log \bigl( 1/p(w_1, w_2) \bigr) だけど、「もし正規化を忘れて選択情報量を  \log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right) と考えてしまうと \beta F(\beta) だけ底上げされているよ! \beta F(\beta) だけ差し引く必要があるよ!」という量なのかな。w_1, w_2 によらずいつも  \log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right) から \beta F(\beta) が差し引かれているよね。


 \displaystyle \log \left( \frac{1}{p(w_1, w_2)} \right)  = \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) +  \log Z(\beta)= \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) - \beta F(\beta)
選択情報量を「その出来事を知らされたときの驚きの量」みたいに捉えるなら、「必ず起きることを知らされたときの驚きの量をゼロにするために差し引くべきオフセット」が自由エネルギーだといえるかもしれない。じゃあ確率分布  q_1(w_1)q_2(w_2) によるエンコードで同じ点 (w_1, w_2) の選択情報量を測ったものはというと、 \log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right) に比べて  \log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right) + \log q_1(w_1)q_2(w_2) だけ差し引かれたものになっているよね。トートロジカルだけど。
 \displaystyle \log \left( \frac{1}{q_1(w_1)q_2(w_2)} \right)  = \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) - \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) - \log q_1(w_1)q_2(w_2)
この差し引くオフセットは w_1, w_2 に依存しているけど、確率分布  q_1(w_1)q_2(w_2) で平均すると平均場自由エネルギー \beta \hat{F}(\beta) になるね。だからやっぱりこれも「驚きの量から差し引くべきオフセットの平均」になっている。このときの驚きの量のエンコーディングは最適ではないかもしれないけど。
\beta \hat{F}(\beta) = \int q_1(w_1)q_2(w_2) \log \Bigl( q_1(w_1)q_2(w_2)/ \exp \bigl( -\beta H(w_1, w_2) \bigr) \Bigr) dw_1 dw_2
 \quad \displaystyle =\int q_1(w_1)q_2(w_2) \left( \log \left( \frac{ 1}{ \exp \bigl( -\beta H(w_1, w_2) \bigr)} \right) + \log q_1(w_1)q_2(w_2) \right) dw_1 dw_2
ただ「平均場自由エネルギー」を「驚きの量から差し引くべき平均的なオフセット」と捉えても、「なぜ平均場自由エネルギーを最小化しなければならないのか」に示唆は得られないな…。いまは「カルバック・ライブラー情報量」つまり「真の驚きの量からの平均的なムダ」を最小化したいから、それと定数 \beta F(\beta) ずれているだけの「(仮の)驚きの量から差し引くべきオフセットの平均」を最小化したいはずだ。エンコードがずれているほど差し引くべきオフセットは平均的にかさむからね。

f:id:cookie-box:20200101101603p:plain:w60

だからその説明で誰かに伝わりますかね…? まあそれで冒頭の推定をするんですが、つまり、パラメータ a_{1:K}, \, b_{1:K} の分布 \varphi(a_{1:K}, b_{1:K}) を更新していくわけです。

  • パラメータ a_{1:K} の分布はディリクレ分布 {\rm Dir}(a_{1:K}|\phi) ですね。ハイパーパラメータ \phi \in \mathbb{R}_{+}^K を更新していくことになります。ディリクレ分布は K 面あるサイコロの各目が出る割合を生成する確率分布ですね。ある割合における確率密度は「各目が \phi_k - 1 回ずつ出たときにその割合であった確率密度」といえます。
  • パラメータ b_{1:K} の分布は  \prod_{k=1}^K \bigl(1 / z(\eta_k) \bigr) \exp \bigl( \eta_k \cdot f(b_k) \bigr) です。12ページでやった指数型分布の共役な事前分布ですね。ハイパーパラメータ \eta \in \mathbb{R}^{J \times K} を更新していくことになると思います。
12ページの指数型分布のときにはもう事後分布を直接計算しにいきましたが、ここでは謎の確率変数 Y を導入します。これは実質的には 1, \cdots, K のいずれかの値を取る確率変数で、元々の確率変数 XY の同時分布が混合分布の k 番目の山のみになるようなものです。言い換えると、YK 面サイコロではあるのですが、\mathbb{R}^N 空間内の点によってサイコロの各面の出る割合が変動します。その割合は、その点における混合分布の 1, \cdots, K 番目の山の密度の比になります。a_{1:K}, \, b_{1:K} から山たちが決まり、山たちから「各点における密度」「各点における各山の密度の比」が決まるので、x^n, y^n, a_{1:K}, b_{1:K} の同時分布はこうです。
\displaystyle P(x^n, y^n, a_{1:K}, b_{1:K}) = \varphi(a_{1:K}, b_{1:K}) \prod_{i=1}^n p(x_i, y_i |a_{1:K}, b_{1:K})
これをさらに観測データで条件付けた P(y^n, a_{1:K}, b_{1:K} | x^n)q(y^n) r(a_{1:K}, b_{1:K}) で近似したいのですが、 P(y^n, a_{1:K}, b_{1:K} | x^n) \propto P(x^n, y^n, a_{1:K}, b_{1:K}) なのでこのまま自己無矛盾条件に代入して差し支えありません。そうすると結局、q(y^n) は、各データの K 面サイコロの目があるパターンのときだけ 1 をとる確率分布になります(これは誤り)。r(a_{1:K}, b_{1:K}) の方はハイパーパラメータを更新した形になります。それで「各データのサイコロの目のパターン」と「ハイパーパラメータ」を交互に更新していけばよかろうという算段のようなんです。
…この手続きで、「謎の確率変数 Y 」と「交互の更新」がしれっと導入されましたが、そもそもこんなことをしてよいのでしょうか? 「謎の確率変数 Y 」は、いわば各データの実家が何番目の山かというようなものになっていますが、本当はどの山が実家とかないでしょう? そこには山の重なりのような分布があるだけで、各データがどの山からきたのかなんて人間の妄想に過ぎないはずなんです。それに、その「各データがどの山からきたのか」と「ハイパーパラメータ」を交互に更新するというのも、直感的には解に向かっていきそうですが、何の根拠があってそんなことをしているのか…。

2020-01-04 追記
f:id:cookie-box:20200101101614p:plain:w60

「こんなことをしてもよいのか」どうかは、「何をするためにそうしたのか」によるだろう。今は何をしたいんだっけ。

f:id:cookie-box:20200101101603p:plain:w60

何をしたいんだっけって、データから混合分布を推定したいんですよ。

f:id:cookie-box:20200101101614p:plain:w60

推定するって何?

f:id:cookie-box:20200101101603p:plain:w60

推定するって何って…この本では、4ページで定義した事後分布によって確率モデルを平均したものこそが推定分布でしたよね。だからこの定義に沿って事後分布を得ようとするのが筋のはずです。

\displaystyle \begin{split} p(a_{1:K}, b_{1:K}|X^n) &= \frac{1}{Z_n(1)} \varphi(a_{1:K}, b_{1:K}) \prod_{i=1}^n p(X_i |a_{1:K}, b_{1:K} ) \\  &= \frac{1}{Z_n(1)} \varphi(a_{1:K}, b_{1:K}) \prod_{i=1}^n \sum_{k=1}^K a_k v(X_i) \exp \bigr(f (b_k) \cdot g(X_i) \bigl) \end{split}
…山が1つのみのときは共役な事前分布がありましたが、山が複数あると山どうしのロスタームが出てきてしまいますね…少なくとも山の個数だけ共役な事前分布を用意すればよいということにはなりません…。

f:id:cookie-box:20200101101614p:plain:w60

まあ解析的に求まる例ではないんだよね。だからモンテカルロ法で解くとか、別の方法で解っぽい何かを探すことになる。5.2節の「平均場近似」はその「解っぽい何かを探す別の方法」だから、「4ページで定義した事後分布で確率モデルを得ること」はもはや放棄されている。「こんなことをしてもよいのか」というよりは「こうすることにした」んだね。…といっても、「だからどうしてそうすることにしたのか」の回答にはなっていないな。いま解析的には解けないけど「パラメータ  a_{1:K}, b_{1:K} の分布をよりよさそうなものに更新する」ことを達成したい。だから次善の策として平均場近似では別の分布でこう近似することにした。

  • 本当の事後分布とのカルバック・ライブラー距離の極小点を目指すことにする。
  • 極小点が満たすべき条件を更新式として利用できるようにするために、パラメータを2グループ(以上)に分けて、こちらのグループの分布を固定してあちらのグループの分布を更新する → 今度はさっき固定していたグループの分布を更新する → …という繰り返しができるようにする。この目的のため、パラメータの分布を独立な分布の積とする。
…実際に平均場近似が編み出された歴史的経緯はさておき、目的からブレークダウンするなら平均場近似というフレームワークは上のように解釈できるだろう。

f:id:cookie-box:20200101101603p:plain:w60

なるほど、だから「謎の確率変数」が導入された…説明になっていないんですが。「交互の更新」の方にはそれで説明が与えられましたが。

f:id:cookie-box:20200101101614p:plain:w60

うん、上のフレームワーク自体は「謎の確率変数の導入」を要請しないからね。ただ今回は、157ページの \log P は「謎の確率変数」を導入しなくては「あるパラメータたちの分布を固定して別のパラメータたちの分布のハイパーパラメータを更新する」という形に持ち込めないはず。「謎の確率変数」を導入しないと結局山どうしのロスタームが出ちゃうしね。

f:id:cookie-box:20200101101603p:plain:w60

「謎の確率変数」の導入はそういう背景が…しかし、「じゃあこういう確率変数を導入しよう」というのも突飛な気がしますが…。

f:id:cookie-box:20200101101614p:plain:w60

そんな突飛かな? 「山どうしのロスターム消えてほしい」→「各データの実家の山は一つまで」って思えば自然に導入できるんじゃない?

f:id:cookie-box:20200101101603p:plain:w60

いうほど自然に導入できますかね?? それに、クロスタームを消してしまうということはそれだけ正しいやり方を損ねているでしょう? 「各データの実家の山は一つまで」というのも結構大胆な仮定にみえるんですが、それが心配というか…。

もう少しつづく

雑記: サンプリング定理の話(その2)

キャラクターは架空のものです。お気付きの点がありましたらご指摘いただけますと幸いです。
参考文献

  1. Nyquist–Shannon sampling theorem - Wikipedia
  2. Fourier transform - Wikipedia

f:id:cookie-box:20200101101614p:plain:w60

前回は、関数 x(t)f_s/2 以上の周波数成分を含まないならば、サンプリング周波数 f_s でサンプリングした x(n/f_s), \, n \in \mathbb{Z}フーリエ変換[-f_s/2, f_s/2] の範囲を切り取って逆フーリエ変換することで元の x(t) を復元できる、という話をしたね。じゃあぴったり f_s/2 の周波数成分についてはどうなるのかを調べると、余弦波は復元できて正弦波は復元できなかった。これは、正弦波はそもそもサンプリングする点にいつも「波」がないのだから捉えられないのも当然といえる。

f:id:cookie-box:20200101101603p:plain:w60

それで前回の最後に副部長は「余弦波と正弦波をフーリエ変換してみてもわかる」といっていましたね。では余弦波からフーリエ変換してみますか。周波数 f_0余弦波は x(t)=\cos(2 \pi f_0 t) ですから、フーリエ変換の定義にしたがって、

\displaystyle \begin{split} X(f) &= \int_{-\infty}^\infty x(t) e^{-2 \pi i t f} dt = \int_{-\infty}^\infty \cos(2 \pi f_0 t) e^{-2 \pi i t f} dt = \int_{-\infty}^\infty \frac{e^{2 \pi i f_0 t} + e^{-2 \pi i f_0 t}}{2} e^{-2 \pi i t f} dt \\ &= \int_{-\infty}^\infty \frac{e^{-2 \pi i (f - f_0) t} + e^{-2 \pi i (f + f_0) t}}{2} dt = \frac{1}{2} \left[ \frac{e^{-2 \pi i (f - f_0) t}}{- 2 \pi i (f - f_0)}\right]_{t=-\infty}^{t=\infty} + \frac{1}{2} \left[ \frac{e^{-2 \pi i (f + f_0) t}}{- 2 \pi i (f + f_0)}\right]_{t=-\infty}^{t=\infty} \end{split}
…負けました。副部長、余弦波にはフーリエ変換はありません。

f:id:cookie-box:20200101101614p:plain:w60

投了しないで! 余弦波にもフーリエ変換はあるよ!!

f:id:cookie-box:20200101101603p:plain:w60

だって、上の積分は値をもちませんよ。フーリエ変換の定義は  X(f) = \int_{-\infty}^\infty x(t) e^{-2 \pi i t f} dt ですから、この積分ができない x(t)=\cos(2 \pi f_0 t)フーリエ変換をもたないと結論付けざるをえません。

f:id:cookie-box:20200101101614p:plain:w60

そのフーリエ変換の定義は x \in L^1(\mathbb{R}^n) である関数に対する定義だよ! L^1(\mathbb{R}^n) に属さない関数をフーリエ変換したいときは、定義を拡張するんだよ。つまり、その可積分関数に対するフーリエ変換の定義は、x(t), y(t) が2乗可積分でもあるならば以下のパーセバルの定理を満たす。「フーリエ変換L^2 -内積を保つ」といっても同じだ。

\displaystyle \int_{-\infty}^\infty x(t) \overline{y(t)} dt = \int_{-\infty}^\infty X(f) \overline{Y(f)} df
これをフーリエ変換の新しい定義にすることができる。x(t) が可積分でなくても、可積分y(t) との内積が定義できる関数なら(というか関数でなくても y(t) に作用してからそれを積分することができるようなものなら)、上式を満たす X(f)x(t)フーリエ変換とする(はず)。例えば  x(t) = e^{2 \pi i f_0 t} を左辺に代入してみよう。
\displaystyle \begin{split} \int_{-\infty}^\infty x(t) \overline{y(t)} dt &= \int_{-\infty}^\infty \overline{ y(t) } e^{2 \pi i f_0 t} dt = \int_{-\infty}^\infty \overline{ y(t)  e^{-2 \pi i f_0 t} } dt = \overline{ \int_{-\infty}^\infty y(t)  e^{-2 \pi i f_0 t}  dt } = \overline{Y(f_0)} \end{split}
よって \overline{Y(f_0)} = \int_{-\infty}^\infty X(f) \overline{Y(f)} df を満たす X(f)x(t)フーリエ変換だけど、この X(f) は「積分することで隣の関数のある点での値を抜いてくるもの」、つまりデルタ関数だ。よって X(f) = \delta(f - f_0) に他ならない(まあ以下の303番にもそうかいてあるし…)。

f:id:cookie-box:20200101101603p:plain:w60

 x(t) = e^{2 \pi i f_0 t}フーリエ変換X(f) = \delta(f - f_0) …となると、x(t) = \cos(2 \pi f_0 t)フーリエ変換

\displaystyle X_{\cos}(f) = \mathcal{F}\bigl[ \cos(2 \pi f_0 t) \bigr] = \frac{1}{2}  \mathcal{F}\bigl[e^{2 \pi i f_0 t}\bigr]+ \frac{1}{2}\bigl[e^{-2 \pi i f_0 t}\bigr] = \frac{1}{2}\delta(f - f_0) + \frac{1}{2}\delta(f + f_0)
ですね…。x(t) = \sin(2 \pi f_0 t)フーリエ変換も同様です。
\displaystyle X_{\sin}(f) = \mathcal{F}\bigl[ \sin(2 \pi f_0 t) \bigr] = \frac{1}{2i}  \mathcal{F}\bigl[e^{2 \pi i f_0 t}\bigr]- \frac{1}{2i}\bigl[e^{-2 \pi i f_0 t}\bigr] = \frac{1}{2i}\delta(f - f_0) - \frac{1}{2i}\delta(f + f_0)
まあ上のリンクの304番と305番にもこうかいてありますし…。しかし、余弦波のフーリエ変換も、正弦波のフーリエ変換も、2箇所に立つデルタ関数でできていて、どちらも似たような感じにみえますが。

f:id:cookie-box:20200101101614p:plain:w60

余弦波と正弦波の周波数 f_0 をサンプリング周波数の半分 f_s/2 にして、 X_{\cos}(f) X_{\sin}(f)f_s ずつずらして自身の「かげぶんしん」を重ねていくとどうなるかな。

f:id:cookie-box:20200101101603p:plain:w60

ああ、いまは余弦波/正弦波からのサンプルを離散フーリエ変換することを考えているのでしたね。余弦波はこうです。

\displaystyle \begin{split} \sum_{k = - \infty}^\infty X_{\cos}(f - k f_s) &= \frac{1}{2} \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s / 2) + \frac{1}{2} \sum_{k = - \infty}^\infty \delta(f - k f_s + f_s / 2)  \\ &= \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s/ 2) \end{split}
これ、2つの項がそれぞれ 1/2 の高さ(?)のデルタ関数を櫛のように立たせていますが、立たせている場所が同じなので、結局、f=f_s/2 を起点に f_s おきに高さ 1デルタ関数が立っていることになりますね。正弦波の方は、
\displaystyle \begin{split} \sum_{k = - \infty}^\infty X_{\sin}(f - k f_s) &= \frac{1}{2i} \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s / 2) - \frac{1}{2i} \sum_{k = - \infty}^\infty \delta(f - k f_s + f_s / 2)  \\ &= 0 \end{split}
2つの項が高さ 1/2iデルタ関数の櫛と高さ -1/2iデルタ関数の櫛をちょうど同じ場所に建設しているので、相殺して消えますね…なんという負の奇跡…。

f:id:cookie-box:20200101101614p:plain:w60

うん、だから余弦波の離散フーリエ変換は、区間 [-f_s/2, f_s/2] を切り取って、区間の端は 1/2 倍すれば余弦波の離散ではないフーリエ変換の結果に一致する。正弦波の離散フーリエ変換は、区間 [-f_s/2, f_s/2] を切り取る以前に 0 だ。これじゃどうしようもない。


じゃあ  x(t) = \exp(\pi i f_s t) という関数だとどうなるかな? これは高くなったり低くなったりを繰り返す波ではなくて、時間 2/f_s の間に複素数平面の単位円上を反時計回りに一周することを繰り返す、「回転する振り子」のようなものだね。ちなみにこれのフーリエ変換X(f) = \delta(f - f_s/2) だとさっきやった。

f:id:cookie-box:20200101101603p:plain:w60

それなら離散フーリエ変換はこうですね。

\displaystyle \sum_{k = - \infty}^\infty X_{\exp}(f - k f_s) = \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s/ 2)
ん? これ余弦波の離散フーリエ変換と一緒ですね…。複素数平面上を回転するものを時間 1/f_s ごとに観測した結果と、実軸上を往復するものを時間 1/f_s ごとに観測した結果とで、抽出した周波数成分が同じになってしまう…? いやでも、時刻 n/f_s における座標がどちらも \bigl( \cos(n \pi), \, 0 \bigr) で一致しているのでそれも道理な気はしますが。しかし、今回は正弦波のときのような相殺が起きたわけではありませんよね?? それなのに情報が欠落してしまうのでしょうか…?

f:id:cookie-box:20200101101614p:plain:w60

サンプリング定理の手続きには、暗黙裡に「ちょうど f_s/2 の周波数成分については、『反時計回りに回転する成分』と『時計回りに回転する成分』が等しい(区間 [-f_s/2, f_s/2] の端は 1/2 倍して切り取る)」という仮定が置かれている。だから、ちょうど 2/f_s の周期で複素数平面上を反時計回りに回転する振動子は、振幅の半分を「時計回りに回転する成分」にもぎとられてしまう。その結果、実軸上を往復する振動子に同一視されてしまう…って感じだと思うけどな。区間が重なる端っこの取り扱いはなんか決め打たないといけないし、対称だと考えるのが自然だろうしね。


結局、以下の関数をサンプリング周波数 f_s でサンプリングして元の信号を復元するとこうなるね。正弦波は「消える」。
  1.  x(t) = \cos(\pi f_s t) x(t) = \cos(\pi f_s t)
  2.  x(t) = \sin(\pi f_s t) x(t) = 0
  3.  x(t) = \exp(\pi i f_s t) = \cos(\pi f_s t) + i \sin(\pi f_s t) x(t) = \cos(\pi f_s t)
サンプリング周波数のちょうど半分 f_s/2 の周波数成分については 、「時刻 0複素数平面上のどこにいるか」の方向に射影した1次元の振動の成分しか生き残らない、とまとめることができるんじゃないかな。

f:id:cookie-box:20200101101603p:plain:w60

しかし、サンプリング周波数のちょうど半分の周波数の正弦波が消えてしまうとはややこしいですね。いったいどうするのがいいのか…。

f:id:cookie-box:20200101101614p:plain:w60

サンプリング周波数大きく取ればいいんじゃないかな。

f:id:cookie-box:20200101101603p:plain:w60

はい。

おわり