論文読みメモ: Temporal regularized matrix factorization for high-dimensional time series prediction(その1)

以下の論文を読みます。

Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon. Temporal regularized matrix factorization for high-dimensional time series prediction. In Advances in neural information processing systems, pages 847–855, 2016. https://papers.nips.cc/paper/6160-temporal-regularized-matrix-factorization-for-high-dimensional-time-series-prediction
※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらご指摘ください。
f:id:cookie-box:20200101101614p:plain:w60

前回の論文の引用文献の話が長くなりそうだったから記事をスピンアウトするね。つまり、高次元時系列予測モデルの先行研究であって、ニューラルネットを利用するものであって、グローバルなパターンも学習するものの2つ目の方ね。論文の図のまんまだけど、時系列を下図のように  Y \approx FX と行列分解(MF)するんだね。F の各行は i 番目の商品の特徴で、X の各列はその時刻の特徴って感じだね。

f:id:cookie-box:20200216104950p:plain:w480
こんな行列分解を求めるには以下を解くんだけど、2番目の項と3番目の項は正則化項だね。
 \displaystyle \underset{F, X}{\rm min} \sum_{(i, t) \in \Omega} (Y_{it} - f_i^\top x_t)^2 + \lambda_f \mathcal{R}_f(F) + \lambda_x \mathcal{R}_x(X)
でも通常(とは)の行列分解みたいに \mathcal{R}_x(X) = \| X \| _F \| \cdot \| _F は後で出てくるけどフロベニウスノルムだね…最初  Y \approx FXF と関係ある何らかのノルムなのかと思った…)なる正則化は適切ではないといっているね。もし Y がアイテム×ユーザという行列だったら、各ユーザがどのアイテムをより好きかを表現できればいいけど、いまはアイテム×時刻という行列だから、「各時刻にどのアイテムがより強い(?)」だけでは不十分で、次の時刻にその強さがどう変わるかも必要だもんね。
多変量時系列分析に行列分解を用いた先行研究ではこういう時間方向の依存性をグラフベースで取り扱ってラプラス正則化を適用しているみたいだけど、これだと2時点間の負の相関を考慮できなかったり、そもそも時間依存性が明示的に手に入っていない場合は推測しないといけなかったり、欠損値には強いものの予測は不得手だったりするって。
でも、この論文で提案する temporal regularized matrix factorization framework (TRMF) は「時間的正則化項」を導入することでデータドリブンで時間依存性を学習できて予測性能もあると。MF の特徴であるスケーラビリティと欠損値への強さはそのままに。あと先行のグラフベースアプローチとの統一的な見方も与えると。そうすることで既存のソルバーが活用できるって。

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

MFを用いる先行手法にはかなり制約があるのですね…? 本論文の2節には先行手法とその限界についてまとめてありますね。

もう少しつづく

論文読みメモ: Think Globally, Act Locally: A Deep Neural Network Approach to High-Dimensional Time Series Forecasting(その1)

以下の論文を読みます。

Rajat Sen, Hsiang-Fu Yu, Inderjit S. Dhillon. Think Globally, Act Locally: A Deep Neural Network Approach to High-Dimensional Time Series Forecasting. In Advances in Neural Information Processing Systems 32, 2019.
https://papers.nips.cc/paper/8730-think-globally-act-locally-a-deep-neural-network-approach-to-high-dimensional-time-series-forecasting
※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらご指摘ください。
次回:まだ
f:id:cookie-box:20200101101603p:plain:w60

近年の時系列データには、一系列が一個人に対応し、何千もの系列が互いに相関するようなものがありうるとありますね。なのでそれらの系列からグローバルなパターンを搾り取った上で、個別データに対するキャリブレーションをするとよい予測ができるのではないかといっています。そこでこの論文で提案するモデル、DeepGLO は「グローバルな」Matrix Factorizartion モデルと、個別の系列の性質及びそれらの相関を捉えるネットワークを組み合わせるということです。


冒頭で時系列予測は様々な分野で重要な課題だといって挙げられている例は、1番上のは本論文と同じ Amazon の方々の論文で、2番目は割と古いですね。3番目も結構古い本です。最近のデータセットは何千もの相互に相関する時系列  Y \in \mathbb{R}^{n \times t} が含まれているということですが、ECサイトでいう「各カテゴリの需要」や「各商品の需要」というのを指しているのでしょうか。需要が高まりそうな商品については契約している業者に在庫の補充を促したりするのですかね?
それで、時系列予測は伝統的には個々の系列や少数の系列たちに対して AR、ARIMA、指数平滑法、Box-Jenkins法、状態空間モデルを適用する方法がとられてきたとありますね。このうち、状態空間モデルについては複数の系列を一緒に取り扱えますよね。しかし、膨大な数の系列にはスケールしないということですか。まあそんな気はしますが…。そして、すべての系列が共通にもっているパターンを上手く利用できないともあります。これは状態空間モデルでもそうなんでしょうか。よくわかりません。それらに対してニューラルネット手法として RNN、LSTM、WaveNet の話があって、これらは大規模な時系列も学習できるものの、2つの欠点があると指摘していますね。
  • 個々の時系列の規模が異なると上手くいかない。例えば個々の製品の需要予測の場合、人気商品の需要の規模はニッチな商品のそれよりも大きい。このような場合、上記のようなニューラルネット手法では正規化が必要になる。しかし、正規化のパラメータ調整は難しい。
  • モデルは予測時に個別のデータの過去データしか利用していない。しかし、データはすべての時系列にわたるグローバルなパターンをもっている。
…グローバルなパターンを活用しようとしたモデルについてこう簡単にかかれていますが、具体的にどのようにグローバルなパターンを捉えているのですかね?

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

1つ目の方は論文内の図がわかりやすいけど以下のような感じだね。まず2次元畳込みをして、その出力を普通にGRUする層と、スキップしてGRUする層もつくって(例えば1時間毎といった周期性があるデータなら1時間スキップして1時間毎のレベルの変化を学習するみたいな)、それらの出力を全結合してニューラルネットの予測結果にする。それを、通常のARモデルの結果と足し合わせるみたい。

f:id:cookie-box:20200208002515p:plain:w480

2つ目はスケーラブルで欠損値があっても大丈夫らしい。Wal-mart E-commerce datasets でよい予測ができたと。

もう少しつづく

文献読みメモ: 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