新年一発目に真面目な話。 numpy.powerのソースコードどこにあんねん。見つからんわ。ma.power.pyはあったけど、c言語の実装がほしいよね。
数百万から数億の波形を取得して2グループに分割し、二者の平均値に差があるか検定を行うことがある。分野ではTVLAというのだけど、詳しい話は割愛。サンプルポイントごとに検定を行うので、波形データをcsv形式とかで保存しておいて、numpyを使って一枚ずつ処理させると実装が楽。もちろん波形を取得しながらOn the flyでの計算も可能。
d-th orderの検定には、2*d次までのモーメントを計算しなければならない。numpy,narrayの各要素をk乗する方法はいくつかあるが、どれを使うのが良いだろうか?というのが今回の趣旨。
結論を先にいうと、TVLAのように配列の要素数が多いケースではnumpy,powerを使わずベタ書きするほうがいい。 つまり、np,power(x, 3)とするよりもx*x*xとしたほうがはるかに速くなる。 似たような話をstack overflowかredditで見つけたが、もう一度探すのは怠いので引用しない。自分でコードを書いて検証する。
third-order testを行いたいので6乗まで計算させる。ループ回数N=1000はやや小さい気もするが、これ以上大きくすると実行時間が長くなるのでやらない。
おっっっっっっっっっっっっっそ。お話になりませんわね。 M=300だと.powerのほうが10倍速かった。arrayのサイズで処理の切り替えとかやってんのか? サンプルポイント数はたいてい1万以上だし、powerは使わないほうがいいね。
けつ
今回は指数がたかだか10程度に留まっているが、50あたりからはnp.powerのほうが速くなると思う。やってないからわからん。25次の検定なんかやらないしな。ケースバイケースで使い分ければええんちゃう。
np.powerのほうも、y6 = np.power(y3, 2)みたいにしないとフェアじゃなくない?wowwowという話もある。うるせぇよ。
注意として、**dやnp,power(, d)の計算結果がmypowerと一致しない事がある。と言っても誤差レベル。IEEE754ﾑｽﾞｶｼｰ。
