クランプ関数(0〜255に値を丸める)のやり方

またまたcodewarsの他人の解答より。 本題はRGBの変換だったけど,その前処理で0以下の数値を0に,255以上の数値を255にクランプする処理がある。

僕は富豪的に下記で。

    lst = map(lambda x: 0 if x<0 else x, lst)
    lst = map(lambda x: 255 if x>255 else x, lst)

でも,一行で書くやり方があって,他の人はみんなこれを使ってた。よくやる処理で当たり前なのかもしれない。

clamp = lambda x: min(255, max(x, 0))

文字列からの数値抜き出しとソート時のインデックス取得

CodeWarsより。 わかりやすいように少し改題。また,CodeWarsはPython2で動いているので下記のコードは通らない。

問題
入力:数値を含む文字列の配列が与えられる。ただし,各文字列には1~9の数字が1つだけ入っており,文字列間に含まれる数値の重複はない。
   例) sentence = ['hoge8hoge', 'aa2bed', '5foobar']
出力:文字列に含まれる数値順にソートせよ。
   例) sentence' = ['aa2bed', '5foobar', 'hoge8hoge']

練習になった点は,(1)文字列から数値部を抜き出すには?,(2)ソートした時に,そのインデックスは取得できる? (2)に関して,numpyだとargsortで簡単にできる。sortした時にそのargも返してくれれば良いけど,それは無いので,両方欲しいときはargsortをもとのnumpy配列にインデックスとして渡す。 matlabとかはmax演算の時にargmaxもくれるけど,その違いに似ている。でも,numpyもCとかに比べれば相当楽だ。Pythonは便利だなぁ。

文字列の抜き出し = re.subを使う。味噌はsubは抽出はできずに代替(substitute)するので,数値以外の部分を削る。

import re

string = 'hoge3'
re.sub(re.compile('\D'), '', string)

# 配列に適用したいならmap
map(lambda x: re.sub(re.compile('\D'), '', x), sentence)
imoprt numpy as np

d = [1, 3, 5, 2]
sorted_index = np.argsort(d)  # [0, 3, 1, 2]
sorted = np.array(d)[sorted_index]  # もしソートしたリストも欲しいなら

配列の中に奇数回出現する1つの整数値

codewarsというので遊んでいたら,他人の良い回答を見つけた。

問題
 整数の配列が与えられる。その中に1つだけ奇数回出現する整数値が存在する。それを求めよ。
 例)[1,3,5,2,6,8,4,3,1,5,2,6,8,3] -> 3

色々とやり方はあると思うけれど,XORを使う方法がなるほどと思った。 対象の整数以外は全て偶数回出現する,というのを利用して,XORとればそれは0になる。 reduceはpython3からfunctoolsに入ったみたい。map/fileterは標準なのに変な感じする。

from operator import xor
from functools import reduce
seq = [1,3,5,2,6,8,4,3,1,5,2,6,8,3]
return reduce(xor, seq)

並列計算の効果検証用プログラム

並列化の試行やJITの効果テストなどのために,Forをたくさん回す計算例が必要になる。 そのような例としてモンテカルロ法を置いておく。

ToDo 並列計算のテスト用としては,モンテカルロ法のような単純に並列化される(データ共有が無い)ものだけでなく,適時データ共有が必要となる例を別途置いておく。

モンテカルロ法による円周率計算

  • モンテカルロ法の例題でよく見る計算。[0,1]の範囲の乱数を生成して,その点が半径1の円の中に入るか否か?(点が円の中か否かは,点の距離dがsqrt(x2+y2)で計算できるから,円周率piを知らなくても判定できる)をランダムに試行して,N回サンプリングすると,円の中に入る比率が1:pi/4になることから推定する(円の面積がpirrなのは知っている)
# monte_carlo.py
import random

def calc_pi(iter):
    cnt = 0
    for _ in range(iter):
        x, y = random.random(), random.random()
        if x*x + y*y < 1:
            cnt += 1
    return (cnt*4.0)/iter
  • 全く同じ書き方でC言語でも。
#include <stdio.h>
#include <stdlib.h>

double calc_pi(int iter)
{
  int cnt=0;
  double x, y;
  
  for (int i=0; i<iter; i++) {
    x = (double)rand()/RAND_MAX;
    y = (double)rand()/RAND_MAX;
    if ((x*x + y*y) < 1)
      cnt++;
  }
  return cnt*4.0/iter;
}

int main(int argc, char** argv) {
  int iter;
  double pi;
  
  iter = atoi(argv[1]);
  pi = calc_pi(iter);
  printf("pi:%f", pi);

  return 0;
}
  • ちなみに,生のPythonは20倍くらい遅い。
# Python
$ time python3 monte_carlo.py 10000000 # 10M
real 0m3.838s                          # 0.4us/iter

# C(Clang, -O3)
$ time ./a.out 10000000                # 10M
real 0m0.214s                          # 0.02us/iter

モンテカルロ法1 (モンテカルロ積分)

Rによるモンテカルロ法入門のメモ,その1。

イントロ

  • 統計的推論では最適化と積分の計算,特に実用上は数値計算が必要になる。
  • 数値積分の手法としてモンテカルロ積分がある。シンプソン法などの区分法は多次元の場合に区分を細かくしようとすると計算量が増加してしまうので適用が困難である。そのような多次元の場合に対してモンテカルロ積分は有用。

モンテカルロ積分

  • モンテカルロ積分の考え方は非常に単純で,下記の期待値計算の積分を対象の確率分布からのサンプリングによる平均値で近似する方法。
  • 例えば,f(x)が確率密度関数で,xがf(x)に従う時のh(x)の期待値を計算したいときに,

 E_f(h(x)) = \int_{x}f(x)dx

 \frac{1}{m}\sum_{i=1}^{m}h(x_j)

で近似する。

  • 同様に考えて,定積分も対象の範囲に対する一様乱数を生成して,その時のh(x)の値を求めて平均取れば,その区間の平均値となるので,積分区間を掛けて上げれば積分値が得られる。(b-a倍するのを忘れないように)
import numpy as np
from scipy import integrate
import math

f = lambda x: math.sin(x)/x  # 被積分関数

def mc_integrate(f, a, b, n=10000):
    def scaling(x):
        '''
        np.random.randで生成する[0,1)の一様乱数を[a,b)に変換する関数
        '''
        return (b-a)*x + a

    xs = np.random.rand(n)
    xs = np.vectorize(scaling)(xs)
    return (b-a) * np.vectorize(f)(xs).sum()/n


def main():
    v = mc_integrate(f, 0, 0.0, 2.0*math.pi)
    
    # Scipyを使って検算
    v_val = integrate.quad(f, 0.0, 2.0*math.pi)[0]
    print('Scipy, v={:.5}, MC, v={:.5}'.format(v_val, v))

if __name__ == '__main__':
    main()    

パフォーマンス計測 in Python

timeモジュール(time.perf_counter)

  • 使いどころは特定のコードブロックに対する計測。精度はマイクロ秒程度。精度を求める場合はこれ。(逆にcProfile/line_profileなどのプロファイリングは基本的にクロック程度の精度,1/100秒程度しかなく,ボトルネックの発見が目的と考えてよいのかもしれない。)
  • 使い方は開始位置と終了位置でtime.time()ならエポック秒(UNIX時間)を,time.perf_counter()ならパフォーマンスカウンタの値を取得して差分として計算する。
  • Pythonのドキュメントに依ると,time.time()システムによっては1秒程度しか精度が無い場合もある。しかし,UNIXにおいてはgettimeofday()を使って時刻を取得するようなので,マイクロ秒単位程度の精度があるようだ。 # Web上にはtime.time()が1/60~1/100秒程度の精度の記述が見られるが,Python3のドキュメントにはその記述は見つけられなかった。
  • より明確にパフォーマンスカウンタ値を返す関数としてtime.perf_counter()がある。time.time()と同じような気もするけど,使い方は同じであるし,無難にこちらを使っておく。 # 実際に,time.get_clock_info('time')とtime.get_clock_info('perf_counter')のresolutionはどちらも1e-09。
import time
start_time = time.perf_counter()
# some heavy calculation
# ...
end_time =  time.perf_counter()
print('Elapsed time:{:.7}'.format(end_time - start_time))

cProfile

  • 使いどころは関数ごとの処理時間やコール回数を計測したい場合。精度は1/100秒程度。
  • 使い方はスクリプト実行時に,cProfileをインポートして実行するだけ。
  • 表示の意味をまとめておく。
    • ncalls: 呼び出し回数。再帰の場合が別途表示されていて,"再帰含む場合"/”再帰を含まない場合”の2つが出力される。
    • tottime: その関数で消費された時間。その関数が呼び出す関数の消費時間は含まない。
    • cumtime: その関数で消費された時間。その関数が呼び出す関数の消費時間は含む。
  • '-s'オプションでソートして出力することがほとんど。
    • time: tottimeでソート
    • cumulative: cumtimeでソート
    • calls: ncallsでソート
# tak.py (Takeuchi Function) is tested as an example.
hoge@foo:~/$ python3 -m cProfile -s time tak.py 

         12604868 function calls (8 primitive calls) in 3.457 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
12604861/1    3.457    0.000    3.457    3.457 tak.py:3(tak)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        1    0.000    0.000    3.457    3.457 tak.py:15(main)
        1    0.000    0.000    3.457    3.457 tak.py:1(<module>)
        2    0.000    0.000    0.000    0.000 {built-in method time.perf_counter}
        1    0.000    0.000    3.457    3.457 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

line_profiler

  • 使いどころは行単位のプロファイルを取りたい場合。cProfileであたりを関数単位で付けて,計測したい関数をline_profileで指定(デコレータ)する。
  • いくつかの使い方があるが,コードを修正してしまうので面倒ではあるが単純なのは関数に@profileデコレータを付けて,kerfnprof.py(line_profileと一緒にインストールされる)を使う。
# hoge.py
@profile
def f():
    # heavy calculation    
    ...
    ...
    ...
    return 
  • pipでインストールして,かつline_profilerをgit cloneしていて例があるが,kernprof.py自体はローカルにあるので不要。
$ pip3 install line_profiler --user # ~/.local/lib/python3.6/site-packages/にインストールされる。
$ python3 ~/.local/lib/python3.6/site-packages/kernprof.py -l -v hoge.py

訓練・検証・テストの3つにデータを分割する理由

MLの教科書とか記事を読んでいると,データを訓練,検証, テストの3つに分けている場合が多い。 訓練と検証の2つじゃだめなの?とついつい勘違いしてしまうのでメモ。

テストはなんのため?

  • ここの勘違いが全てなんだと思う。
  • テストの目的は実運用時の性能を推定すること。勘違いしがちなのは,モデル選択にテストデータでの評価結果を使うこと。
  • 繰り返すが,テストは性能を推定するだけ,モデル選択するためじゃない。

テストデータと検証データの使い方

  • MLのモデルはハイパーパラメタがある場合がほとんど。
  • よくある勘違いは訓練データで種々のモデル(モデル自体やハイパーパラメタが違うモデル)を訓練し,テストデータで性能を測定し,最も性能が良いものを採用してしまう方法。これはあくまでも,そのテストデータに対して最も良いモデルを選んだだけであって,その値を持って実運用時の性能を推定することはできない。
  • くどいけれど,あくまでも,テストはテストであって,選択基準にしてはならない。
  • 適切な使い方は,訓練データでパラメタを学習して,検証データでモデル選択をして,テストデータで性能を測定する。
  • だから,名前のままなんだ。訓練データで訓練して,検証データで検証してみて,テストデータでテストする。
  • そうか,検証とテストという言葉が似すぎているのか!検証は自分側に裁量がまだあるが,テストはもはや裁量権を失っているようなイメージなのかな。モデルを検証しましたので,テストお願いします,みたいな感じなのかな。

補足(クロスバリデーション)

  • データを分割すると,貴重なデータ数が目減りして見える。
  • より有効的に訓練と検証をする方法がクロスバリデーション。
  • データを訓練データとテストデータに分けて,訓練データに対して交差検証をする方法。これはデータをうまく使いつつ,良いモデルを選択できる。