Python3 – 対数

Pythonで対数を出すには、math.logを使います。import mathで使えるようになります。

import math
print(math.log(2, 10))

0.30102999566398114

これは、10を何乗したら2になるか?です。log102です。
試しに、10を0.30102999566398114乗してみます。

print(pow(10, 0.30102999566398114))

1.9999999999999998

100の0.5乗は、100の平方根なので10です。

print(pow(100, 0.5))

10.0

27の0.3333333333333333333333333333333333333333333乗は、大体27の立方根なので、3に近いはずです。

print(pow(27, 0.3333333333333333333333333333333333333333333))

3.0

おージャスト3.0になった。

print(pow(27, 0.333333333333333))

2.9999999999999964

print(pow(27, 0.3333333333333333))

3.0

少数15桁までだと2.99…になりますが、16桁にするとジャスト3.0になる。

10を底とする対数を常用対数といいます。

対数関数と指数関数は逆関数の関係にあります。逆関数はy=xに対して対称です。指数関数と対数関数のグラフを書いてみます。
y = axの指数関数と、y = logaxの対数関数を書きます。

numpyを使いたいので、math.logではなくnp.logを使います。aを2として、np.log2を使います。例えばnp.log2(4)とすると2になります。2を底とした4の対数を出しています。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0.1, 3, 0.1)
y = x
plt.plot(x, y, label='y=x')
y = 2 ** x
plt.plot(x, y, label='y=2**x')
y = np.log2(x)
plt.plot(x, y, label='y=log2 x')
plt.grid(True)
plt.legend()
plt.show()

対称なのか全然これじゃ分からない。

対数関数は、
対数が1より小さいとき、
底が1より大きければ、真数は底より小さい。
底が1より小さければ、真数は底より大きくなる。

対数関数を、底の大きさのパターンで2つ書いてみる。

import math
import numpy as np
import matplotlib.pyplot as plt

#底は1ではない正の数
#底が1より大きい場合
a = 2
#底が1より小さい場合
b = 0.5
#xは真数
x_list = np.arange(0.1, 3, 0.1).tolist()
#yは対数
y_a = []
y_b = []
for x in x_list:
    y_a.append(math.log(x, a))
    y_b.append(math.log(x, b))
plt.plot(x_list, y_a)
plt.plot(x_list, y_b)
plt.grid(True)
plt.show()

2100の桁数を出すには、常用対数を使います。
log102100 = 100log102 = 100 * 0.30102999566398114 = 30.102999566398114
よって31桁です。10の2乗は100なので3桁です。10の3乗は1000で4桁です。2.1乗は100より大きく、1000より小さいので3桁です。30.10乗は31桁です。

Python3 – xのn乗のグラフ(matplotlibのsubplotとアニメーション)

GIFアニメーション

matplotlibのアニメーション作成は2つ種類があって、ArtistAnimationとFuncAnimationとがある。

参考:matplotlib でアニメーションを作る

ArtistAnimation は、あらかじめ全てのグラフを配列の形で用意しておき、それを1枚ずつ流すというものである。FuncAnimationは予め完成したグラフを渡すのではなく、アニメーションの1フレームごとに関数を実行する。データが巨大すぎる場合や潜在的に無限に続く場合に便利。

ArtistAnimation

ArtistAnimationを使って、xの2乗~10乗のグラフを配列に入れて、アニメーションをつくってみる。

サンプルコード

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def xpown(s, e):
    graphs = []
    fig = plt.figure()
    for n in range(s, e + 1):
        x = np.arange(-100, 100, 0.1)
        y = x ** n
        g = plt.plot(x, y, label='y = x ** {}'.format(n))
        graphs.append(g)
    plt.legend()
    ani = animation.ArtistAnimation(fig, graphs, interval=500)
    ani.save('pow1.gif', writer='imagemagick')

if __name__ == '__main__':
    xpown(2, 10)

結果

これだと、10乗のグラフのyが大きすぎて、他がみんな直線みたいになってしまった。

FuncAnimation

FuncAnimationを使ってやってみる。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

x = np.arange(-100, 100, 0.1)
n = 1
max = 10
fig = plt.figure()

def plot_pow(data):
    plt.cla()
    global n
    n = n + 1 if n < max else 2
    y = x ** n
    im = plt.plot(x, y, label='y = x ** {}'.format(n))
    plt.legend()

if __name__ == '__main__':
    ani = animation.FuncAnimation(fig, plot_pow, interval=400, frames=9)
    ani.save("pow2.gif", writer="imagemagick")

結果

指数が偶数ならy軸に対称のグラフになり、奇数なら原点に対象のグラフになります。原点に対象っていうのか。

Subplot

matplotlibで複数のグラフを表示したいときにsubplotが使えます。
参考:[Python]Matplotlibで複数のグラフを描画する方法

サンプルコード

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-100, 100, 0.1)
fig = plt.figure()
ax = []
for i in range(9):
    y = x ** (i ++ 2)
    ax.append(plt.subplot2grid((3, 3), (i // 3, i % 3)))
    ax[i].set_title('y = x ** {}'.format(i + 2))
    ax[i].plot(x, y)
plt.show()

結果

Python3 – 素数

エラトステネスの篩

エラトステネスの篩

結果
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]

参考:素数と仲良しになる3つのプログラム

フェルマーの小定理

フェルマーの小定理で素数か判定できる。

Pythonを使って高速素数判定をしてみる」に詳しく書いてあった。フェルマーテストというので出せるけど、素数じゃないものも出しちゃう可能性があって、より高性能な方法がある。それもこのサイトに書いてある。

ベルマン方程式

ベルマン方程式は、動的計画法(動的な最適化問題)の最適性の必要条件を表す方程式らしい。必要条件は、再帰を使って部分を解くことで全体が解ける状態にあることと、メモ化を使うことです。最適化問題とは、集合内のすべての数値を、ある関数にいれたときに最小あるいは最大となる集合内の元を探すことを目的とする問題のことです。

動的計画法において、問題の目的を数式で表したものを目的関数といいます。

xは状態です。aは行動です。tは時間です。βは割引率です。F(x, a)は、状態xのときに行動aをとることで得られる報酬が出力される関数です。この式は、報酬を最大化させることが目的になっています。割引率は0 < β < 1で、0であればどの時間tに得られた報酬も同じ重さで扱います。1に近ければ未来のtで得られた報酬程重要視します。t100のときの報酬は、βの100乗が掛けられるからです。 この式は再帰可能な形態に変形できます。まず下記のようにします。

最初のx0のときの報酬だけ左に出しました。βも一つ左に出すことで、右の式はt0がt1になっただけで、あとは上記の式とまったく同じです。よって下記のように書けます。

これを下記のように簡単に書いたのが、ベルマン方程式だそうです。

なぜ、上の式のx1が、下ではT(x, a)に変わっているかというと、T(x, a)は、状態xのときに、行動aをとった場合に次に遷移する状態を返す関数として定義しているからです。変数の添え字を消したかったのでこれを使っています。

Python3 – 動的計画法(フィボナッチ数列)

動的計画法は、分割統治法メモ化を合わせた方法のことらしい。分割統治法は、問題を細分化して、細かい部分を順に解いていくことで全体を解明するようなことの総称らしい。

分割統治法は、コード的には下記のようになり、再帰することになる。

function conquer(x)
  if xは簡単にconquerできる then
    return conquerの結果
  else
    (x1, x2, ...) := divide(x)     // 複数の小さな副問題に分割
    (y1, y2, ...) := (conquer(x1), conquer(x2), ...)  // 再帰
    return synthesize(y1, y2, ...)  // 各副問題の解を合成(最大値を選ぶ、等)

メモ化とは、再帰の際に同じ副問題を複数回解いてしまうことによる効率低下を防ぐために、副問題の結果を保存して再利用することらしい。キャッシュみたいな感じ。

下記のフィボナッチ数列を表示するプログラムは、動的計画法の具体的な例らしい。

int fib(unsigned int n) {
    int memo[1000] = {0, 1}, i;
    for (i = 2; i <= n; i++) {
        memo[i] = memo[i - 1] + memo[i - 2];
    }
    return memo[n];
}

Python3だと下記のような感じ。

def fib(n):
    memo = [0] * n
    memo[0:2] = [0, 1]
    for i in range(2, n):
        memo[i] = memo[i - 1] + memo[i - 2]
    return memo[0:n]
print(fib(100))

ちなみに、ネットでフィボナッチ数列のPython版しらべたら下記がでてきたけど、これはまさしく動的計画法じゃない例だと思った。毎回メモらずに計算しているようだ。

def fib2(n, a=1, b=0):
    return b if n < 1 else fib2(n - 1, a + b, a)
print([fib2(i) for i in range(100)])

せっかくだから時間図ってみる。

import time

def fib(n):
    memo = [0] * n
    memo[0:2] = [0, 1]
    for i in range(2, n):
        memo[i] = memo[i - 1] + memo[i - 2]
    return memo[0:n]
start = time.time()
print(fib(900))
print(time.time() - start)

def fib2(n, a=1, b=0):
    return b if n < 1 else fib2(n - 1, a + b, a)
start = time.time()
print([fib2(i) for i in range(900)])
print(time.time() - start)

結果

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,...
0.001961946487426758
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,...
0.0782785415649414

おー40倍速いってことか。さすが動的計画法だ。

下記のような方法もある。

#include <stdbool.h>

int memo[1000] = {0, 1};
bool in_memo[1000] = {true, true};

int fib(unsigned int n) {
    if (!in_memo[n]) {
        memo[n] = fib(n - 1) + fib(n - 2);
        in_memo[n] = true;
    }
    return memo[n];
}

Python – 勾配法

勾配は、すべての変数に対する偏微分をベクトルにしたものです。勾配法は、勾配を使って関数が最小になるパラメタを探す方法のことです。偏微分はある点における、各パラメタに対する微分ですので、傾きです。傾きが分かるとパラメタをどう動かすとマイナスになるかが分かります。勾配はすべてのパラメタに対する傾きが入ってるので、これを使って全部のパラメタを徐々に関数の出力結果がマイナスになるように調整します。ただし、関数のグラフが凸凹している場合、ある部分においての凹の一番下についてしまった場合、出られなくなったりします。残念なことです。でも大体なんとかなるようです。ただしグラフの形状によっては、実際残念な結果に終わったり、最小値を発見するのにものすごく時間がかかったりするそうなのでケースバイケースで勾配法とは別の方法も試す必要があるそうです。

勾配によって、各パラメタの増減すべき方向が分かりますが、どの程度増減させるかは人間がいい感じの数字を設定する必要があります。この一回当たりにパラメタを増減させる数値は、学習率で決めます。偏微分結果に学習率を掛けたものが、増減数値になります。また、どこまでいけば本当の正解なのかは勾配法では分かりませんので、勾配法を何回繰り返すかも人間が決める必要があります。

勾配をPythonで書くと下記になります。

import numpy as np
 
def func(x):
    return x[0]**2 + x[1]**2
 
def gradient(f, x):
    h = 1e-4
    grad = np.zeros_like(x)
    for i in range(x.size):
        tmp = x[i]
        x[i] = tmp + h
        y1 = f(x)
        x[i] = tmp - h
        y2 = f(x)
        grad[i] = (y1 - y2) / (2 * h)
        x[i] = tmp
    return grad
        
g = gradient(func, np.array([3.0, 4.0]))        
print(g)

勾配法をPythonで書くと下記になります。

def gradient_descent(f, x, lr=0.01, num=100):
    for i in range(num):
        x -= lr * gradient(f, x)
    return x

コードサンプル

import numpy as np

def func(x):
    return x[0]**2 + x[1]**2

def gradient(f, x):
    h = 1e-4
    grad = np.zeros_like(x)
    for i in range(x.size):
        tmp = x[i]
        x[i] = tmp + h
        y1 = f(x)
        x[i] = tmp - h
        y2 = f(x)
        grad[i] = (y1 - y2) / (2 * h)
        x[i] = tmp
    return grad
        
def gradient_descent(f, x, lr=0.01, num=100):
    for i in range(num):
        x -= lr * gradient(f, x)
    return x

x2 = gradient_descent(func, np.array([3.0, 4.0]), 0.1, 100)
print(x2)

python – 損失関数

ディープラーニングで使う損失関数の、二乗和誤差と交差エントロピー誤差。

二乗和誤差

二乗和誤差は、差の二乗を全部足して2で割る。

def nijowagosa(y, t):
    return 0.5 * np.sum((y - t)**2)

y = np.array([0.2, 0.03, 0.8])
t = np.array([0, 0, 1])

print(nijowagosa(y, t))

交差エントロピー誤差


$$E = -\displaystyle \sum_{k}t_k \log y_k$$

logは底がeの自然対数を表すそうです。np.logでできます。ディープラーニングで使うときは、yが出力で、tが正解ラベルです。正解ラベルは正解を1とし、間違いを0とすると、正解ラベル以外のものは全部0になるので無視できて、-log yになります。間違いのやつ無視していいのかなと思いましたが、間違いのやつを高確率と判断した場合、ソフトマックス関数では正解の確立が低くなるので問題ないのかなと思いました。

def cross_error(y, t):
    delta = 1e-7
    return -np.sum(t * np.log(y + delta))

y = np.array([0.2, 0.03, 0.8])
t = np.array([0, 0, 1])
print(cross_error(y, t))

deltaは、np.log(0)が-infを出すことに対する対策です。超小さい数を加えています。

Python – 偏微分・勾配

偏微分

偏微分は、変数が2つ以上あるときに1つ以外を固定にして、固定じゃない1つに対して微分をすることだそうです。Pythonでやると簡単にできます。

下記の式で、x0が3、x1が4のときの偏微分


$$f(x_0, x_1)=x_0^2+x_1^2$$

コード

import numpy as np

def numerical_diff(f, x):
    h = 1e-4
    return (f(x+h) - f(x-h)) / (2*h)

def func_x0(x0):
    return x0**2 + 4**2

def func_x1(x1):
    return 3**2 + x1**2

d0 = numerical_diff(func_x0, 3) 
d1 = numerical_diff(func_x1, 4)
print(d0)
print(d1)

勾配

勾配は、すべての変数に対する偏微分をベクトルにしたものだそうです。

import numpy as np

def func(x):
    return x[0]**2 + x[1]**2

def gradient(f, x):
    h = 1e-4
    grad = np.zeros_like(x)
    for i in range(x.size):
        tmp = x[i]
        x[i] = tmp + h
        y1 = f(x)
        x[i] = tmp - h
        y2 = f(x)
        grad[i] = (y1 - y2) / (2 * h)
        x[i] = tmp
    return grad
        
g = gradient(func, np.array([3.0, 4.0]))        
print(g)

Python – Matplotlibで3次元グラフを書く

下記の式をグラフ化してみる。


$$f(x_0, x_1)=x_0^2+x_1^2$$

コード

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def func(x0, x1):
    return x0**2 + x1**2

x0 = np.arange(-3, 3, 0.25)
x1 = np.arange(-3, 3, 0.25)
X0, X1 = np.meshgrid(x0, x1)
Y = func(X0, X1)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X0,X1,Y)
plt.show()

参考:matplotlibで3Dグラフを描画する

Python – 数値微分

微分はある瞬間の変化量です。公式は下記です。超絶的に0に近い、というかもはや0という謎の極小値でもって計算することで瞬間変化量を出します。解析的にやるといろんな公式が編み出されているので真の微分が出せますが、コンピュータで公式通りにやろうとすると、謎の極小値を本当にすごい小さい値にすることになるので、誤差がでます。しかしその誤差が非常に小さいので多くの場合問題ありません。


$$\frac{df(x)}{dx} = \lim_{h \to 0}\frac{f(x+h) – f(x)}{h}$$

Pythonだと下記になります。

def numerical_diff(f, x):
    h = 1e-4
    return (f(x+h) - f(x)) / h

これよりも誤差を小さくする方法は下記です。

def numerical_diff(f, x):
    h = 1e-4
    return (f(x+h) - f(x-h)) / (2*h)

サンプル

import numpy as np
import matplotlib.pyplot as plt

def numerical_diff(f, x):
    h = 1e-4
    return (f(x+h) - f(x-h)) / (2*h)

def sample_func(x):
    return x**2

x = np.arange(-7.0, 7.0, 0.1)
x2 = np.arange(-1.0, 7.0, 0.1)
k = 2

y = sample_func(x)
d = numerical_diff(sample_func, k)
y2 = d * x2 + sample_func(k) - d*k

plt.plot(x, y)
plt.plot(x2, y2)
plt.show()

ソフトマックス関数

合計が1になるように、うまく数字を調整してくれる関数らしい。確率を確認するのに便利。


$$y_k = \frac{\exp(a_k)}{\displaystyle \sum_{i=1}^n \exp(a_i)}$$

aがn個あるときの、k番目のyを取得する。分母は全てのaの指数関数の和。分子はk番目のaの指数関数。

pythonで書くと下記になります。

a = np.array([0, 2.3, 3.5])
exp = np.exp(a)
sum_exp = np.sum(exp)
y = exp / sum_exp
print(y)

指数関数を使っているので数字がものすごく大きくなる可能性があり、それによってオーバフローで変な数値が返ってくる場合があるため、aの最大値をそれぞれのaの値から引きます。これで全部0かマイナスになります。

a = np.array([0, 2.3, 3.5])
c = np.max(a)
dexp = np.exp(a - c)
sum_exp = np.sum(exp)
y = exp / sum_exp
print(y)

ソフトマックス関数は、aのそれぞれの値の大きさの順番は変えません。なので、一番大きいaを取得したい場合は、ソフトマックス関数は不要です。

ディープラーニングの出力層で使う活性化関数は、回帰問題だと恒等関数、2クラス分類問題はシグモイド関数、多クラス分類問題はソフトマックス関数を使うのが一般的だそうです。シグモイドは0か1かを表すし、ソフトマックスは、複数の選択肢の確からしさを数値で表すからだと思います。回帰は数値を予測するので確率にする必要がないので恒等関数なんだと思います。

Python – 行列の掛け算

Numpyを使います。np.dot(A, B)で簡単にできます。
例えば、A([1, 2, 3], [4, 5, 6])と、B([7], [8], [9])を掛けます。

import numpy as np

A = np.array([[1, 2, 3], [4, 5, 6]])
B = np.array([[7], [8], [9]])
print(np.dot(A, B))

下記が出力されます。

[[ 50]
 [122]]


$$\begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{pmatrix}
\begin{pmatrix} 7 \\ 8 \\ 9 \end{pmatrix}$$