numpy v1.6以下でもrandom.choiceの重み付け

## GAE/py でnumpy

Google Cloud Platform(GCP)のAppEngine standard environment with pythonでは
様々な制約がありますが、numpyを使うことができます :+1:

* [組み込みサードパーティ ライブラリ ](https://cloud.google.com/appengine/docs/standard/python/tools/built-in-libraries-27?hl=ja)

ただし、versionは1.6.1に固定されています :scream: 
(2017/8/31現在の最新は1.13.1？)

そのため、様々な機能が使えなかったりしますが、今回はnp.random.choiceを使えなくて困りました。

## numpy.random.choice

1.7以上で追加されたrandom抽出の便利機能です。
[numpy.random.choice](https://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.random.choice.html)

### 基本機能

基本的には配列or配列長を受け取って、乱数抽出で配列要素orインデックスを返すメソッドです。
第二引数(またはsize)にランダム抽出する数を指定します。
サイズは未指定の時は1で、返り値が配列でなく要素になります。

```py
# 配列長を指定してindexを返す
>>> np.random.choice(5)
3
# サイズを指定すれば配列でindexを返してくれる
>>> np.random.choice(5, 3)
array([3, 2, 4])

# 配列を指定して要素を返す
>>> np.random.choice(['alpha', 'beta', 'gamma'])
'alpha'
# サイズを指定すれば配列要素を取り出して返してくれる
>>> np.random.choice(['alpha', 'beta', 'gamma'], 2)
array(['alpha', 'beta'],
      dtype='|S5')
```

実際にはsizeや配列には多次元もサポートされています。
これも1.6に持っていくのは面倒そう。。

しかし標準機能だけならnp.random.randintとかを駆使して再現できそうですね。
randintは配列の範囲とsizeを指定できるので、choiceのindex選択の部分は共通と言えるでしょう。

### オプション機能: 重複制御
サイズ指定して、ランダムに抽出した内容は、デフォルトで重複が許可されています。

```py
# デフォルトだと重複許可なので、同じ値を選んでしまうことも。
>>> np.random.choice(5, 3)
array([0, 0, 0])
>>> np.random.choice(['alpha', 'beta', 'gamma'], 2)
array(['gamma', 'gamma'],
      dtype='|S5')

# replace=Falseで重複を禁止できます。
>>> np.random.choice(5, 3, replace=False)
array([0, 1, 4])
>>> np.random.choice(['alpha', 'beta', 'gamma'], 2, replace=False)
array(['gamma', 'alpha'],
      dtype='|S5')
```

サイズが選択対象の配列長よりも大きいとエラーになります。

今回はこの重複禁止は1.6以下のnp.randomだと実装されていないように見えたので困りました。
しかもrandintを重複しないようにloopさせると運が悪いとだいぶ時間がかかりそう。
(擬似乱数わかってないのでそんなことないのかもしれませんが。。)

なので今回はpythonのrandomモジュールのsampleを使ってしまいます。
random.sampleは配列から重複なし・サイズ指定ありで抽出ができます。
乱数生成方法が別かもしれないので厳密に同じではないと思いますが、、、

### オプション機能: 重み付け

pという名前で抽出対象と同じ長さの配列を渡すことで重み付けをすることができます。
pはprobabilityの略で、合計が1になる必要があるようです。
デフォルトでは均等に分布します。

```py
# 重み付けをすることで偏りをつけられます。
>>> np.random.choice(5, 3, p=[0.1, 0, 0.3, 0.6, 0])
array([3, 2, 3])
>>> np.random.choice(['alpha', 'beta', 'gamma'], 2, p=[0.1, 0.1, 0.8])
array(['gamma', 'gamma'],
      dtype='|S5')

# この記事では対応できてませんが、replaceと併用もできます
>>> np.random.choice(5, 3, p=[0.1, 0, 0.3, 0.6, 0], replace=False)
array([2, 3, 0])
>>> np.random.choice(['alpha', 'beta', 'gamma'], 2, p=[0.1, 0.1, 0.8], replace=False)
array(['gamma', 'alpha'],
      dtype='|S5')
```

この重み付けはstackoverflowで質問している人がいたので、そのスレッドへの回答を拝借しました。

[How do I “randomly” select numbers with a specified bias toward a particular number](https://stackoverflow.com/questions/25507558/how-do-i-randomly-select-numbers-with-a-specified-bias-toward-a-particular-num)

理解が不十分ですが、pを累積配列に変換して、そこに一様分布の乱数を配置して実際の乱数の値にしているようですね。

## 実装例

だいぶダサい条件分岐ですが、実装してみるとこんなイメージになるかと思います。
重複なしの抽出でpythonのrandomを使ってしまっているのでそこだけパフォーマンス劣化や乱数の偏りとかあるかもしれません。
もっといい実装があったら教えて欲しいです。（実は1.6系でも相当のapiがあったりしそう）

一番いいのはnumpy1.7以上を使うことですね。

```py
def numpy_choice(a, size=1, replace=True, p=None):
    # 1.6だとchoiceがないので再現する. replaceは重複許可. pは重みづけ
    # 整数ならrange配列、そうでないなら配列化する.
    values = np.arange(a) if isinstance(a, int) else np.asarray(a)
    if p:
        # TODO: pの長さとaの長さを検証する必要があります。
        # あとreplaceとの併用は今は対応してません。。(たまたま不要だったので)
        choiced = weighted_choice(values, p, size)
    else:
        length = len(values)
        if replace or size > length:
            # 重複ありならrandintで対応
            idx = np.random.randint(0, length, size)
        else:
            # 重複なしはpythonのrandom.sampleを使ってしまう
            idx = random.sample(np.arange(length), size)
        choiced = values[idx]
    if size == 1 and len(choiced) == 1:
        # size 1の時は要素を返却する
        return choiced[0]
    return choiced


def weighted_choice(values, p, size=1):
    # 重みありのchoice. stackoverflowのままです。
    values = np.asarray(values)

    cdf = np.cumsum(np.asarray(p))
    cdf /= cdf[-1]

    uniform_samples = np.random.sample(size)
    idx = cdf.searchsorted(uniform_samples, side='right')
    sample = values[idx]

    return sample
```

### 課題

- randomとnp.randomを併用していてよくない(randomだけの方がいいかも)
- pの検証ができてない
- pとreplaceの併用の実装ができていない