Pythonで多くの回帰を計算する最速の方法は？

Question

私はこれを達成する方法についてはかなり合理的な考えを持っていると思うが、私はすべてのステップで100％確実ではない。この質問は、私が最も効率的なやり方でこれをやっていること、そして私の統計知識が実際には健全であることを保証するための健全性チェックとして主に意図されています。

とにかく、私がやろうとしているかについて、いくつかの説明：

私は上のいくつかの線形回帰を実行したい時系列データがたくさんあります。特に、私は約500の異なる変数について2000回の観測を行っています。各変数について、2つの説明変数（約2000件の観測の2つの追加ベクトル）を使用して回帰を実行する必要があります。ですから、500種類のYのそれぞれについて、aとbを次の回帰​​で見つける必要があります。

これまで、私はregsを実行するためにstatsmodelsパッケージのOLS関数を使用してきました。しかし、私が知る限り、問題を解決するためにstatsmodelsパッケージを使用したいのであれば、それを何百回も呼び出さなければならないでしょう。これは一般的には効率が悪いようです。

代わりに、私は長い間本当に触れていない統計をいくつか再訪することにしました。私の知識がまだ間違っていれば、私はすべての観測値を約2000 x 500の1つの大きなYマトリックスに入れることができます。次に、説明変数を約2000 x 2のXマトリックスに張り、すべての結果(X'Y)/(X'X)を計算して私の回帰の500。私が基本的なnumpyのもの（matrix.Iを使った*と逆を使った行列乗算）を使ってこれを行うなら、それは何百ものstatsmodel OLS呼び出しよりもずっと速くなると思います。ここで

は、私が持っている質問です：

• は私がstatsmodels何度も呼び出して、以前の方法よりも高速をしていますnumpyのものですか？もしそうなら、私が欲しいものを達成するための最速/最も効率的な方法ですか？私はそれがあると仮定していますが、あなたがより良い方法を知っていれば、私はそれを聞いてうれしく思います。 （確かに、私はこのように多くの回帰を計算する必要がある最初の人ではありません）
• 行列に欠けているデータをどうやって扱いますか？私の時系列データは素敵で完全なものではなく、たまに価値が失われます。私がnumpyで規則的な行列乗算をしようとすると、NAの値が伝播し、私は最終的には最終的な結果としてほとんどがNAの行列になります。私が個々の回帰を単独で行う場合、私は回帰を実行する前にNAsを含む行を削除することができますが、これを大規模な2000 x 500行列で実行すると、実際の非NAデータを他の変数から削除します私は明らかにそれが起こることを望んでいません。
• 最初に行列に入れる前に、私の時系列データが実際に正しく並べられることを保証する最も効率的な方法は何ですか？私の観測の開始日と終了日は必ずしも同じではなく、一部のシリーズでは他の人がそうでない日があるかもしれません。もし私がこれを行うための方法を選ぶなら、私はすべての観測を日付で索引付けされたパンダのデータフレームに入れます。そうすれば、パンダは私のためにすべてを整理して作業を終わらせ、それが終わった後に基礎となるndarrayを抽出することができます。これは最善の方法ですか、またはパンダには別の方法でマトリックス構築を行うことで回避できるオーバーヘッドがありますか？

Answer 1

繰り返しstatsmodelsを呼び出すと、いくつかの簡単な答え

1）は、最速の方法ではありません。パラメータ、予測、残差が必要で、同じ説明変数が必要な場合は、通常はparams = pinv(x).dot(y)を使用します。ここでyは2次元であり、そこから残りを計算します。問題は、推論、信頼区間などが作業を必要とするため、速度が重要でなく、限定された結果セットしか必要とされない限り、統計モデルOLSはさらに便利です。

これは、すべてのyとxが同じ観測値インデックス、欠損値、およびギャップを持たない場合にのみ機能します。

脇に：セットアップは多変量線形モデルであり、今後は統計モデルによってサポートされる予定です。

2）および3）欠落した細胞があるか観察の完全な重複がない場合（指標）、ケース1）の高速単純線形代数は機能しません。アナログからパネルへのデータでは、最初のケースでは「バランスの取れた」パネルが必要であり、他のケースでは「アンバランス」データが含まれています。標準的な方法は、データをブロック対角形式で説明変数に積み重ねることです。これにより、メモリが大量に増加するため、疎行列と疎線形代数を使用するほうが優れています。それは個々のOLS回帰をループするよりもスパース問題の構築と解決が速いかどうかという特定のケースに依存します。

専門コード：（ジャスト思想）：ケース2では

）完全に重複しているか、またはcellwiseない欠損値を持つが、我々はまだ、すべてのyのためにすべてのX'X、およびX'Y行列を計算する必要がありますすなわちそれらのうち500個。 500 x 2 x 2の2つの回帰子しか持たないことを考えると、大きなメモリを必要としません。したがって、不完全マスクを交差積計算の重みとして使用することによって、パラメータ、予測および残差を計算することが可能である可能性があります。 numpyは私が知っている限り、linalg.invをベクトル化しています。だから、私はいくつかのベクトル化された計算でこれを行うことができると思います。