今回は、要素間の足し算引き算を便利に行なってくれる関数2つを紹介します。

np.diff

まずは要素間の差分をとるnp.diff関数についてです。この関数は、ある時間ごとに観測したデータ値の差分を分析するときや微分などをとりたい場合に有効です。APIドキュメントを見てみましょう。

**np.diff(a, n=1, axis=-1, prepend=, append=)**

params:

パラメータ名 概要
a array_like
配列に相当するもの
差分をとりたい配列を指定します。
n int (省略可能)初期値1
差をとりたい微分の階数を指定します。
axis int (省略可能)初期値-1
どの軸方向に差分をとっていくのかを指定します。初期値では最も次元が低い方向になっています。
prepend,append array_like (省略可能)初期値no value。
差を計算する前に、指定した軸に沿ってaの前に追加するか後に追加する値。

returns:

差分のとられた配列が返されます。

np.diffの引数は5つ存在します。第一引数で差分をとりたい配列を指定し、第三引数でどの軸方向で差分をとるかを指定します。 軸(axis)がどの方向を示しているかについては、以下の記事を参考にしてみてください。

第二引数のnについてですが、こちらは複雑ですので後述します。基本的には、隣り合う要素間の差をとったものを要素とする配列を返すので、差分のとられた次元における配列の要素数は元の配列から1減少したものとなります。

このため、のちほど登場するcumsumを結果の配列に適用しても、元の配列を再現することはできません(初項がわからないので当たり前ではありますが)。

引数nについて

引数nというのは、n階微分の差分をとるということを意味しています。以下は、微分方程式などでよく使う差分法の公式となっています。

\frac{\partial u}{\partial x} \to \frac{u(x+\Delta x) - u(x)}{\Delta x} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \\ \frac{\partial^2 u}{\partial x^2} \to \frac{u(x+\Delta x) - 2u(x) + u(x-\Delta x)}{(\Delta x)^2} \ \ \ \ \ \ \ \\ \\ \frac{\partial^3 u}{\partial x^3} \to \frac{u(x+2\Delta x) - 3u(x+\Delta x) + 3u(x) - u(x-\Delta x)}{(\Delta x) ^3} \\

n=3にあたる3階微分までの公式を挙げてみました。この分子にあたる部分を自動的に計算してくれるのがnp.diff関数です。これを NumPyの配列上で行うと、以下のような処理を施してくれることになります。

配列xの差分を配列yに記録するとします。

n = 1 のとき\\ y[i]=x[i+1]-x[i] \\ n = 2 のとき\\ y[i]=x[i+2]-2x[i+1]+ x[i] \\ n = 3 のとき\\ y[i]=x[i+3]-3x[i+2]+3x[i+1]-x[i] \\ n = 4 のとき\\ y[i] = x[i+4] -4x[i+3] + 6x[i+2] - 4x[i+1] + x[i] \\

このあたりまで来ると勘がよい方は気づくかもしれませんが、これらの係数は二項係数と一致しています。ここに関して説明を始めると一気に長くなる上、少し本記事の趣旨から外れてしまうため、割愛します。詳しく知りたい方は、以下のページを参考にしてみてください。

インデックスが元の公式からズレていますので、注意してください。
では、実際のコードでnの値を変えてみましょう。

In [1]: import numpy as np

In [2]: a = np.array([1, 2, 4, 1, 6, 8, 3]) # 適当に値を並べた配列を用意する

In [3]: np.diff(a, n= 1) # まずはn=1から。
Out[3]: array([ 1,  2, -3,  5,  2, -5])

In [4]: np.diff(a, n=2) # 次はn=2。
Out[4]: array([ 1, -5,  8, -3, -7])

In [5]: np.diff(a, n=3) # 次はn=3
Out[5]: array([ -6,  13, -11,  -4])

In [6]: np.diff(a, n=4) # 最後にn=4
Out[6]: array([ 19, -24,   7])

axisを変えて差をとる方向を変更する

ついでに、他の引数であるaxisの値も変化させていきましょう。第4,5引数は省略します。

In [8]: b = np.random.randint(10, size = (5,5))

In [9]: b
Out[9]:
array([[4, 9, 9, 6, 3],
       [7, 6, 0, 3, 6],
       [8, 2, 2, 1, 2],
       [6, 5, 5, 8, 5],
       [6, 2, 2, 2, 9]])

In [10]: np.diff(b, axis=-1) # axis=1と同じ意味。
Out[10]:
array([[ 5,  0, -3, -3],
       [-1, -6,  3,  3],
       [-6,  0, -1,  1],
       [-1,  0,  3, -3],
       [-4,  0,  0,  7]])

In [11]: np.diff(b, axis=0) # 次は行方向
Out[11]:
array([[ 3, -3, -9, -3,  3],
       [ 1, -4,  2, -2, -4],
       [-2,  3,  3,  7,  3],
       [ 0, -3, -3, -6,  4]])

In [12]: np.diff(b, axis=1, n = 2) # n=2にしてみる。
Out[12]:
array([[-5, -3,  0],
       [-5,  9,  0],
       [ 6, -1,  2],
       [ 1,  3, -6],
       [ 4,  0,  7]])

np.cumsum, np.ndarray.cumsum

次は、np.diff関数と対になる機能を持つnp.cumsum関数について解説します。np.cumsum関数は配列内の要素を足し合わせていったものを順次配列に記録していていくものです。近似的な積分を行なっているイメージですね。

まずは、np.cumsumnp.ndarray.cumsumのAPIドキュメントをそれぞれみていきましょう。内容としては全く一緒ですし、行われる処理も全く同じです。

numpy.cumsum(a, axis=None, dtype=None, out=None)

params:

パラメータ名 概要
a array_like
配列に相当するもの
足し合わせをとりたい配列を指定します。
dtype dtype (省略可能)初期値None
返される配列の要素のデータ型、および計算処理を行うデータ型を指定します。初期値のNoneのままだとaの要素のデータ型が反映されます。例外として、aのデータ型がint8,int16,int32だった場合、計算の精度の高いint64で配列が返されます。
out ndarray (省略可能)初期値None
計算結果を入れていく配列を指定します。

returns:

指定された軸方向に足し合わされていった値を要素とする配列が返されます。
outが指定されているときはその結果が格納されたoutが返されます。

np.ndarray.cumsum(axis=None, dtype=None, out=None)

params:

パラメータ名 概要
dtype dtype (省略可能)初期値None
返される配列の要素のデータ型、および計算処理を行うデータ型を指定します。初期値のNoneのままだとaの要素のデータ型が反映されます。例外として、aのデータ型がint8,int16,int32だった場合、計算の精度の高いint64で配列が返されます。
out ndarray (省略可能)初期値None
計算結果を入れていく配列を指定します。

returns:

指定された軸方向に足し合わされていった値を要素とする配列が返されます。
outが指定されているときはその結果が格納されたoutが返されます。

cumsumはdtypeで計算処理を行うデータ型と出力される値のデータ型の両方を同時に指定することができます。また、axisを指定することで足し合わせを行なっていく方向を指定することができます。

では、実際のコードで使い方をみていきましょう。

In [1]: import numpy as np

In [2]: a = np.randint(10, size=20) # 0~9までの乱数を20個生成。
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-2-a51548630d0c> in <module>()
----> 1 a = np.randint(10, size=20) # 0~9までの乱数を20個生成。

AttributeError: module 'numpy' has no attribute 'randint'

In [3]: a = np.random.randint(10, size=20) # 0~9までの乱数を20個生成。

In [4]: a
Out[4]: array([6, 4, 1, 4, 8, 7, 4, 9, 1, 6, 7, 1, 2, 6, 0, 8, 5, 1, 1, 4])

In [5]: np.cumsum(a) # まずは普通に足し合わせてみる。
Out[5]:
array([ 6, 10, 11, 15, 23, 30, 34, 43, 44, 50, 57, 58, 60, 66, 66, 74, 79,
       80, 81, 85])

In [6]: a.cumsum() # この形でもok
Out[6]:
array([ 6, 10, 11, 15, 23, 30, 34, 43, 44, 50, 57, 58, 60, 66, 66, 74, 79,
       80, 81, 85])

次に、dtypeを指定していきましょう。

In [7]: np.cumsum(a, dtype='float32') # dtypeを"float32"に指定する。
Out[7]:
array([  6.,  10.,  11.,  15.,  23.,  30.,  34.,  43.,  44.,  50.,  57.,
        58.,  60.,  66.,  66.,  74.,  79.,  80.,  81.,  85.], dtype=float32)

In [8]: a.cumsum(dtype='float32')
Out[8]:
array([  6.,  10.,  11.,  15.,  23.,  30.,  34.,  43.,  44.,  50.,  57.,
        58.,  60.,  66.,  66.,  74.,  79.,  80.,  81.,  85.], dtype=float32)

In [9]: b = np.random.rand(3,4)*10

元の配列がintの64ビット未満の場合は、dtypeが特に指定されないとint64になって処理が行われます。

In [18]: c = np.random.randint(10, size = 10, dtype='int8')

In [19]: c
Out[19]: array([9, 4, 8, 1, 5, 7, 0, 0, 7, 0], dtype=int8)

In [20]: c.cumsum() # 特にdtypeを指定しないと
Out[20]: array([ 9, 13, 21, 22, 27, 34, 34, 34, 41, 41])

In [21]: c.cumsum().dtype # int64になる。
Out[21]: dtype('int64')

In [22]: d = c.cumsum(dtype='int8') # dtypeをしっかり指定すればこのような変化は起こらない。  

In [23]: d.dtype
Out[23]: dtype('int8')

次にaxisを指定していきましょう。

In [10]: b
Out[10]:
array([[ 1.17551018,  5.45865113,  1.7240211 ,  3.81506854],
       [ 3.33030478,  6.14776842,  5.71525196,  1.57526629],
       [ 4.44130255,  0.42665233,  5.81375779,  7.21976963]])

In [11]: np.cumsum(b) # axisに何も指定しないと
Out[11]:
array([  1.17551018,   6.63416131,   8.35818241,  12.17325095,
        15.50355573,  21.65132415,  27.3665761 ,  28.94184239,
        33.38314494,  33.80979727,  39.62355507,  46.8433247 ])

In [12]: np.cumsum(b, axis=1) # こうすると列方向に足し合わせてくれる。
Out[12]:
array([[  1.17551018,   6.63416131,   8.35818241,  12.17325095],
       [  3.33030478,   9.4780732 ,  15.19332516,  16.76859144],
       [  4.44130255,   4.86795488,  10.68171268,  17.90148231]])

In [14]: b.cumsum(axis=1)
Out[14]:
array([[  1.17551018,   6.63416131,   8.35818241,  12.17325095],
       [  3.33030478,   9.4780732 ,  15.19332516,  16.76859144],
       [  4.44130255,   4.86795488,  10.68171268,  17.90148231]])

In [15]: np.cumsum(b, axis=0) # 行方向に足し合わせる。
Out[15]:
array([[  1.17551018,   5.45865113,   1.7240211 ,   3.81506854],
       [  4.50581496,  11.60641955,   7.43927306,   5.39033482],
       [  8.94711751,  12.03307188,  13.25303085,  12.61010445]])

In [16]: b.cumsum(axis=0)
Out[16]:
array([[  1.17551018,   5.45865113,   1.7240211 ,   3.81506854],
       [  4.50581496,  11.60641955,   7.43927306,   5.39033482],
       [  8.94711751,  12.03307188,  13.25303085,  12.61010445]])

参考