デジタル・デザイン・ラボラトリーな日々

アラフィフプログラマーが数学と物理と英語を基礎からやり直す。https://qiita.com/yaju

自然対数の底(ネイピア数) e を理解してみる

はじめに

前回、対数を理解してみました。 yaju3d.hatenablog.jp

今回は機械学習を学ぶ上で出てくる自然対数の底(ネイピア数 e)とは何かを理解していきます。
間違える人は結構多いですが、e は「自然対数」ではありません。「ネイピア数」あるいは「自然対数の底」と呼ばれる定数です。

自然対数の底とは

対数の底には、10進法を使う常用対数とネイピア数と呼ばれる数学定数を使う自然対数があります。
記号として通常は e が用いられ、その値は、e = 2.71828 18284 59045 23536 02874 71352\cdots と続く無理数超越数となっています。 とても不思議な値ですよね。
不思議な値といえば円周率 \pi = 3.14159265359\cdots も同様に無理数超越数となっています。

qiita.com

最後に

今回はQiitaの方に書いてしまったので、リンクするだけにします。

対数logを理解してみる

はじめに

機械学習を学んでいると対数logがでてくる。基礎的なことから対数を理解してみたい。
指数はイメージし易いが、対数は分かりにくいと思われている。指数と対数はペアの関係にあり、かけ算とわり算のように逆関係にある。 先ずは、指数の大きさを視覚的にイメージするために、アメリカ・ワシントンにある航空宇宙博物館で公開されていた9分半の映画「パワーズ・オブ・テン」を紹介する。

9分半の映画

パワーズ・オブ・テンです、10の冪(10^{n})の違いを視覚でご覧ください。


Powers of Ten with Japanese translation

対数とは

例えば、2を3回かけ算すると 2 \times 2 \times 2 = 8 になります。これを「2を3乗したら8になる」と言い、以下のように書きます。

  2 ^ 3 = 8

このとき、2 の右上に乗っている 3 のことを「指数」と言います。指数は「1つの数を何回掛けるか」を表しています。
一方、「◯を何乗すれば△になるか」を表す数のことを「対数」と言います。 例えば「2を何乗すれば8になるか」を表す数は以下のように表記され、これを「2を底とする8の対数」と言います。

  log_{2}8 = 3

2 を底をする 8 の対数は 3 ということになります。
もう少し身近な値としてみます。1000 は、10 の 3乗 で「10を底とする1000の対数」は以下のようになります。

  log_{10}1000 = 3

計算の簡略化

対数には大きな桁数でも簡単に計算できるようになるというメリットがある。
それには指数法則と対数法則を知っておく必要がある。

指数法則①

a^{p} \times a^{q} = a^{(p+q)}

3^{4} \times 3^{5} について考える。3^{4} \times 3^{5} を指数を使わない形で表現すると(3 \times 3 \times 3  \times 3)(3 \times 3 \times 3  \times 3 \times 3) となる。 3 を繰り返しかけ算する回数は、  4 + 5 = 9 回である。
よって、3^{4} \times 3^{5} = 3^{(4+5)} となる。

指数法則②

(a^{p})^{q} = a^{(p \times q)}

(5^{3})^{4} について考える。(5^{3})^{4} は、5^{3} を4回繰り返しかけ合わせるので、(5^{3})^{4} = 5^{3} \times 5^{3} \times 5^{3} \times 5^{3} である。指数の部分の合計は、3乗が4回、つまり (3 \times 4) によって計算できる。
よって、 (5^{3})^{4} = 5^{3 \times 4} となる。

指数法則③

(a \times b)^{p} = a^{p} \times b^{p}

(3 \times 7)^{5} について考える。(3 \times 7)^{5} は、(3 \times 7) を5回繰り返しかけ合わせるので、 (3 \times 7)^{5} = (3 \times 7) \times (3 \times 7) \times (3 \times 7) \times (3 \times 7) \times (3 \times 7)
 = (3 \times 3 \times 3 \times 3  \times 3) \times (7 \times 7 \times 7 \times 7  \times 7) 、つまり = 3^{5} \times 7^{5} によって計算できる。
よって、 (3 \times 7)^{5} = 3^{5} \times 7^{5} となる。

対数法則①

かけ算を足し算に変換

log_{a}(M \times N) = log_{a}M + log_{a}N

log_{2}(8 \times 4) = log_{2}32 = 5 = 3 + 2 = log_{2}8 + log_{2}4
注目すべきは、M \times N というかけ算が、log_{a}M + log_{a}N という足し算に変換されているということである。このことが対数を利用してかけ算を足し算へ簡略化して計算する際に生かされる。

対数法則②

わり算を引き算に変換

log_{a}(M \div N) = log_{a}M - log_{a}N

log_{2}(8 \div 4) = log_{2}2 = 1 = 3 - 2 = log_{2}8 - log_{2}4
今回は、M \div N というわり算が、log_{a}M - log_{a}N という引き算の形に変換される。

対数法則③

累乗を簡単なかけ算に変換

log_{a}M^{k} = k \times log_{a}M

log_{2}8^{3} = log_{2}512 = 9 = 3 \times 3 = 3log_{2}8
例えば、2^{29} のような大きな計算も、この対数法則③を使えば一瞬にして近似値を求められる。

底の変換公式

底も自由に変えられるようになります。
\displaystyle \log_{a}b = \frac{(\log_{c}b)}{(\log_{c}a)} = 対数

\displaystyle \log_{2}8 = \frac{(\log_{10}8)}{(\log_{10}2)}
\displaystyle = \frac{0.903089987}{0.301029995} = 3

C# 数学10 「常用対数、ネイピア数e-基本1、自然対数-基本1」e log ln exp

対数の起源

対数を考え出したのはジョン・ネイピアです。
対数(logarithm)の名前の由来は、logos (比、神の言葉)とギリシャ語のarithmos (数) を合わせて logarithms(ロガリズム) という造語でネイピアが考案しました。

ジョン・ネイピア

1550年、宗教戦争が続くスコットランドにネイピアは生まれました。時代はヨーロッパ列強諸国が覇権を争う大航海時代のまっただ中。ネイピアはマーキストン城の城主として官の仕事の他、領民のために農業土木、軍事技術など多くの発明を行うエンジニアとして活躍しました。(中略) 遠洋航海に必要な天文学(船の測位)を支える数学が球面三角法です。数学者でも天文学者でもないネイピアは現実の切実な問題解決を目標としていました。それが天文学的計算の克服です。三角関数の計算の中に現れる大きい数の計算は天文学者を苦しめました。大航海時代は計算との闘いでもあったのです。天文学者は直面する天文学的計算を克服する手立てを見つけることができませんでした。彼らの計算を助けるために、ネイピアはついに新しい計算法を見つけ出す決心をします。時にネイピア44歳、1594年。その20年後の1614年、ついに人類は青天の霹靂として「対数」を手にします。第5回 ジョン・ネイピア対数誕生物語

また、日々私たちが使っている小数点(.)を使う表記はネイピアの発明です。小数点はネイピアが対数を生み出す過程で考え出した副産物だったのです。第6回 ジョン・ネイピア小数点「・」誕生物語
※小数の考え方はネイピア(1550〜1617)とほぼ同時期のシモン・ステヴィン(1548-1620) ですが19.178と表す小数を「19⓪1①7②8」のように表記していた。ステヴィンの本は1585年に出版、これはネイピアが対数表の計算を開始する1594年の9年前のことです。ネイピアはステヴィンの小数は使いづらいと判断し、もっと機能的で使いやすい小数の表記法を求め続けたのです。

発明を行うエンジニア

ネイピアはマーキストン城の城主であり、困った人を助けるために自分の才能を使う、優秀なエンジニアだったのです。
自分の領地の収穫を増やすために肥料や揚水機の研究をしたり、スペインの.侵攻を恐れて潜水艦や戦車などの兵器も数多く開発しています。これも領地内の人民を安心させるためだったのでしょう。
また、対数の計算を効率化するためにネイピアの骨と呼ばれる、かけ算や割り算などを簡単に行うための計算器具を発明しています。これは1617年(ネイピアが亡くなる年)に論文「小さな棒による計算術」としてエディンバラで発表されました。

対数表

対数表には、指定した数を底とする対数(例 底=2)、常用対数(底=10)、自然対数(底e=2.718281828…)があります。試験などで底を省略した対数表が出た場合は常用対数となります。

ネイピアが作成した対数表は不思議な底(0.9999999)だったため、使いやすい底を10にした常用対数表をネイピアの意思をついだブリッグスが作成しました。1617年に1000まで計算して出版し、1624年には1から20000までと90000から100000まで、小数点以下14桁まで計算した対数表を出版しました。1628年にブリッグスの対数表で抜けていた20,000~90,000の対数表をオランダのアドリアン・ブラックが完成させました。ただ、ブラックの対数表は10桁のものでしたが実際に利用するには十分な桁数でした。対数の発見

qiita.com

対数をとるメリット

対数のメリットは、極めてケタ数の大きい数の面倒なかけ算や割り算を、易しい足し算や引き算に直してくれることにあります。但し、求めた値はあくまで近似値であり、正確な値ではありません。
電卓や計算機のなかった昔の時代において非常に便利な計算方法でした。ピエール・シモン・ラプラスからは「天文学者の寿命を2倍にした(天文学者が一生にこなせる計算量を2倍にした)」と絶賛されるようになりました。

対数表を使った計算

131 \times 219 \times 563 \times 608 を計算する場合

10進数を底とする対数で表すと、log_{10}(131 \times 219 \times 563 \times 608) となる。
log_{10}( (1.31 \times 10^{2}) \times (2.19 \times 10^{2}) \times (5.63 \times 10^{2}) \times (6.08 \times 10^{2}) )
= log_{10}((1.31 \times 2.19 \times 5.63 \times 6.08) \times 10^{8})
ここから、かけ算が足し算になります。
= log_{10}1.31 + log_{10}2.19 + log_{10}5.63 + log_{10}6.08 + log_{10}10^{8}
= log_{10}1.31 + log_{10}2.19 + log_{10}5.63 + log_{10}6.08 + 8

ここで、  log_{10}1.31 、log_{10}2.19 、log_{10}5.63 、log_{10}6.08 の値を対数表から読み取る。
log_{10}(131 \times 219 \times 563 \times 608) \fallingdotseq 0.1173 + 0.3404 + 0.7505 + 0.7839 + 8 = 1.9921 + 8  =  0.9921 + 9 ここで、0.9921 に近い値を対数表から読み取ると、0.9921 \fallingdotseq log_{10}9.82 となる。
log_{10}(131 \times 219 \times 563 \times 608) \fallingdotseq log_{10}9.82 + 9 = log_{10}9.82 + log_{10}10^{9} = log_{10}(9.82 \times 10^{9})
こうして得られた log_{10}(131 \times 219 \times 563 \times 608) \fallingdotseq log_{10}(9.82 \times 10^{9}) の両辺を見比べると
131 \times 219 \times 563 \times 608 \fallingdotseq 9.82 \times 10^{9} = 9820000000 (実際の値は、 9820359456)

2^{29} を計算する場合

対数法則③を使います。
2^{29}は、10を底とする対数で表すと log_{10}2^{29} = 29 \times log_{10}2 となる。
log_{10}2 の値は、常用対数表から 0.3010 である。よって、log_{10}2^{29} = 29 \times 0.3010 = 8.7290 = 0.7290 + 8
常用対数の値が、0.7290 に近い真数の値を表から読み取ると、0.7290 \fallingdotseq log_{10}5.36 となる。
log_{10}2^{29} \fallingdotseq log_{10}5.36 + 8 \times log_{10}10 = log_{10}5.36 + log_{10}10^{8} = log_{10}(5.36 + 10^{8})
よって、2^{29}は、約 5.36 \times 10^{8} = 536000000 と計算できる。(実際の値は、536870912)

基本情報技術者試験の問題

ソートの計算量

計算量は一般に、O記法(オー記法、もしくはビッグオー記法と呼ぶ)を使って表します。O記法では計算量を、O(n)O(n^{2}) のような形で記述します。
性能が良い : O(1) \lt O(\log n) \lt O(n) \lt O(n\log n) \lt O(n^{2}) \lt O(n^{3}) : 性能が悪い

  • O(1)
    データ数nにかかわらず一定なので、nが2倍、3倍、10倍となっても計算量は等しいままです。
  • O(\log n)
    計算量を求める場合は底を2として考えます。
    nが2倍、3倍と増えても、計算量は1、2と増える程度なのでデータ量が多い場合に有用です。
  • O(n)
    比例なので、nが2倍、3倍と増えると、計算量も同様に2倍、3倍となります。
  • O(n \log n)
    先述のO(\log n) にnの計算量がかけられたものです。
    例えば、nが2倍になると計算量は約2.5倍、nが10倍になると計算量は約20倍、nが100倍になると計算量は約300倍になります。
  • O(n^{2})
    n個の要素を用いた総当りの組み合わせを求める場合がこのパターンに当てはまります。
    二次関数なので、nが2倍、3倍と増えると、計算量は4倍、9倍と増えます。

O記法にはいくつかのルールがある。

  • 中身は一番大きな規模だけ残す(最大次数の項以外を除く) 例 2+n+2n^{2}→2n^{2}
  • 係数は除く 例 2n^{2}→n^{2}
  • 係数は1にする

例 1 ~ Nまでの整数の総和を求めよ

  • 普通に計算する → N-1回足し算 → O(N-1) → O(N)  
  • 公式 \displaystyle S=\frac{(n+1)n}{2} を使う → 足し算と掛け算と割り算 → O(3) → O(1)

その中で対数が使われるのが下記となります。

O記法 概要 使用例
O(logN) <対数時間>
処理をひとつ行うたびに処理対象を何割か減らせるようなアルゴリズム。データ量が増えても計算時間がほとんど増えない。多くの場合、これ以上の改善はする必要はない。
ソート済み配列の二分探索
O(NlogN) <準線形、線形対数時間>
ちょっと重いO(N)程度。マージソートのように二分木でデータを分割し、かつそれらをリニアサーチするようなアルゴリズムの計算量。二分木のオーダー(logN)×リニアサーチのオーダー(N)をかけあわせてNlogNになる。
クイックソート、マージソート、ヒープソート

二分探索の計算量

ソートされている前提で、中央の値より小さい(左側)か大きい(右側)が分かるので、次からは半分は探索する必要はなくなる。
これを2の指数で表した場合、指数の数が1つずつ減ることになるので対数にすると、O(\log n) になる。
f:id:Yaju3D:20180122012235p:plain

クイックソートの計算量

分割統治法に基づくアルゴリズムです。
f:id:Yaju3D:20180122013630p:plain
毎回半分ずつにデータが分かれたと仮定すると、平均計算量は n個 \times \log_{2} n 段 → O(n\log n) となる。
f:id:Yaju3D:20180122213121p:plain

n進数での最大値

14けたの16進数の最大値は、10進数で表すと何けたか。ここで、 log_{10}2 = 0.301 とする。

log_{10}(16^{14})
= 14 \times log_{10}(16)
= 14 \times log_{10}(2^{4})
= 14 \times 4 \times log_{10}(2)
= 14 \times 4 \times 0.301 = 16.856

繰り上げて答えは17桁になります。

身近な問題

地震の大きさや音階や星の等数や化石の年代測定や複利計算にも対数が使用します。他にも電気通信や音の単位に用いられるデシベル、音の大きさを表すホンでも使用されます。

地震のエネルギーの大きさ

地震のエネルギーの大きさを表す言葉として「マグネチュード」と言います。

マグニチュードと震度の違いは?
「マグニチュード」は、地震そのものの大きさ(規模)を表すものさしです。一方「震度」は、ある大きさの地震が起きた時のわたしたちが生活している場所での揺れの強さのことを表します。

マグネチュードを M、地震が発生するエネルギーを E (単位モジュール)とすると、logE =  4.8 + 1.5M という関係式があります。
この式より、E = 10^{4.8 + 1.5M} = 10^{4.8}10^{1.5M} となります。
したがって、 M が1増えるとエネルギーは10^{1.5}倍(およそ32倍)、 M が2増えるとエネルギーは10^{3}倍(およそ1000倍) になります。また、M が0.2増えるとエネルギーは10^{0.3}倍(およそ2倍) になります。

マグネチュードが2違うと1000倍で1違うと32倍って差が分かりにくいですが、x \times x = 1000 とすると x^{2} = 1000 = x = \sqrt{1000} = 31.6227  \fallingdotseq 32 となるわけです。

音階の違い

ピタゴラスの定理などで有名な古代ギリシャの数学者ピタゴラスが音階「ドレミファソラシド」を作りました。ピタゴラスが発見した音階は弦の長さで決められたものでしたが、我々が通常聞いている音階は振動数の比で決められています。

平均律音階は、1オクターブ高い音は振動数が2倍であり、その1オクターブを12音階に均等に分けたものである。
f:id:Yaju3D:20180128121348p:plain

音階 ド# レ# ファ ファ# ソ# ラ#
周波数比 2^{0}=1 2^{\frac{1}{12}} 2^{\frac{2}{12}} 2^{\frac{3}{12}} 2^{\frac{4}{12}} 2^{\frac{5}{12}} 2^{\frac{6}{12}} 2^{\frac{7}{12}} 2^{\frac{8}{12}} 2^{\frac{9}{12}} 2^{\frac{10}{12}} 2^{\frac{11}{12}} 2^{\frac{12}{12}}=2
周波数 fr fr^{1} fr^{2} fr^{3} fr^{4} fr^{5} fr^{6} fr^{7} fr^{8} fr^{9} fr^{10} fr^{11} fr^{12}
近似値 262 277 294 311 330 349 370 392 415 440 466 494 523

※周波数 f=262

r^{12} = 2 で、これを満たす r を求めると 1.06^{12} \fallingdotseq 2 となります。半音上がることに振動数を 1.06倍されている。
\displaystyle r = \sqrt[12]{2} = 2^{\frac{1}{12}} = 1.0594631 \fallingdotseq 1.06

つまり、音階も対数変換(底1.06)ととらえることができる。\log_{1.06}(周波数)=音階

星の等数

歴史上で星の明るさをランクづけしたのは、古代ギリシアの天文学者ヒッパルコスです。彼は、肉眼で見えるもっとも明るい星を1等星、かろうじて見える星を6等星として、1等星~6等星までの6段階に分けました。月日が経ち、1856年にイギリスの天文学者のノーマン・ロバート・ポグソンは、ジョン・ハーシェルの研究結果を元にし1 等星の明るさは 6 等星の 100 倍であり、1 等級ごとの明るさの違いは 100^{\frac{1}{5}} 倍であると定義しました。

f:id:Yaju3D:20180128204643p:plain

\displaystyle r = \sqrt[5]{100} = 100^{\frac{1}{5}} = 2.51188643151 \fallingdotseq 2.51

つまり、星の等数は対数(底2.51)となります。

化石の年代測定

化石の年代測定ってどうやっているのか不思議でした。
これは「放射性炭素年代測定法」というのが1974年にアメリカのリビーという人が発見した方法で、炭素14の性質を利用しています。
炭素14は放射性炭素ともいわれ、電子を放出して窒素14に変わる。この崩壊によって炭素14の数が半分になるまでの期間を半減期といい5730年となっている。
生きている生物(動物、植物)はこの炭素14を体内に取り込むので、体内の炭素14の割合は大気中の割合と同じとなる。また生物が死ぬと炭素14の供給がなくなり崩壊だけが続くので発見した資料の炭素14の割合を調べることで、その資料が死んでからの年数が推定できる。

ある木管に炭素14の割合を調べたら、75%に減っていた。
この場合、炭素14が1年で r 倍に現象するとして、この木簡が x 年前のものだとすると、
r^{x}=0.75 また r^{5730}=0.5

x\log r=\log 0.75 、② 5730\log r=\log 0.5
① と ② より
\displaystyle x=\frac{\log0.75}{\log r} = \frac{5730}{\log 0.5} \times \log 0.75
\displaystyle =\frac{5730(\log3 - 2\log2)}{- \log 2}
 = 5730 \times 0.4150 = 2378 年前

複利計算

サラリーローンで、1万円を1日1割の複利で2ヶ月借りる場合、
<複利計算>元利合計  = 元金 \times (1 + 利率)^{期間の数} の式となります。
1 \times (1+0.1)^{60} \fallingdotseq 305 なんと、約305万円にもなります。

これを対数で計算してみます。電卓がない状態で1.1を60回かけるんなんてうんざりですが、その場合は対数表を使えばいい。
\log_{10}x = \log_{10}1.1^{60}
対数の法則を使う
\log_{10}x = 60\log_{10}1.1
対数表では\log_{10}1.1 = 0.0414 となります。
\log_{10}x = 60 \times 0.414 = 2.484
このうち、小数部の 0.484 を対数表で見ると、3.05 です。
\log_{10}x = 2.484 = 2 + 0.484
= \log_{10}10^{2} + \log_{10}3.05
= \log_{10}100 \times 3.05
= \log_{10}305
 x = 305

自然対数

自然対数は、2.71828182845\cdots という無理数を底にした対数となります。一方、常用対数は底を10にした対数となります。
常用対数は電卓やコンピューターの登場により使用されることがなくなりましたが、自然対数はコンピューターで使用する上で標準と底となっています。これは次記事に書きます。 yaju3d.hatenablog.jp

機械学習での役割

  • かけ算を足し算または割り算を引き算にすることで計算を楽にさせる。
    \displaystyle \prod_{n=1}^{N} f(n) = \sum_{n=1}^{N}\log f(n)\prodは、かけ算を繰り返す記号
  • 機械学習では非常に大きな数を扱います。プログラムでこのような数を扱うと、オーバーフローを起こしてしまうときがあります。そのような数は、対数 をとって扱うことで、オーバーフローを防ぐことができます。例えば、100000000 = 10^{8}0.000000001 = 10^{-8} となります。
  • 機械学習では確率を使いますが、同時確率の場合には 1 以下の数の掛け算の連続になります。多い場合にはアンダーフローを起こしてしまうときがあります。そのような数は、対数 をとって扱うことで、アンダーフローを防ぐことができます。
    例としてサイコロを2回投げた場合、1回目に1の目が出て2回目に2の目が出る確率は、下記の計算となる。
    \displaystyle \frac{1}{6} \times \frac{1}{6} = \frac{1}{36}
  • 対数は単調増加関数であるため、ある関数 f(x) があって f(x) を最小にする値を求めた場合、対数をとった \log f(x) でも最小の値は同じになります。また、最大値を求める場合でも同じです。
    f(x)=(x-1)^{2}+2は、x=1の時、最小値をとる
    \log f(x)=\log ((x-1)^{2}+2) も、x=1の時、最小値をとる

単調増加関数

単調増加な関数は、任意な値  x1 と x2 があったとき 、 x1 \lt x2 ならば  f(x1) \lt f(x2) となるような関数  f(x) であるいうことです。
例えば、  x1=3、x2 = 5 x1 \lt x2 が成り立ちます。この関係性は対数に直しても  \log x1 \lt \log x2 が成り立ちます。

「サザエさんのじゃんけん データ分析」の2017年の結果

はじめに

明けましておめでとうございます。

昨年も東芝に暗雲が立ち込め、東証2部に降格して日経平均から外れ、ついには18年3月末にサザエさんからのスポンサー契約を降板する方向で調整とあいなりました。
さて、サザエさんのじゃんけんを長年記録しデータをWebスクレイピングさせて頂いた「サザエさんジャンケン学」が2017年6月25日を以て終了されていました。データを取得した際に7月以降がなかったので、あれっと思って見たらそういうことでした。長い間お疲れ様でしたm(_ _)m
www.asahi-net.or.jp

さてさて、2017年のサザエさんのじゃんけん結果はどうなったでしょう。 ちなみに、2016年のサザエさんのじゃんけん結果は、27勝11敗12分(勝率0.711)でした。

人工知能による予測化を断念

本来はDeepLearningを使用して予想をする予定でしたが、まだ勉強中のままです。あと少しだと思うので今年中にはなんとかなると思います。
昨年、R言語からPythonに切り替えました。
qiita.com

Python使うに際にはいつも Jyupter notebook だったので、この為に久しぶりに Visual Studio Code で Python を実行しようとしてバージョンアップ(1.19.1)したらタスク実行の設定がいろいろ変わって動かなくて、四苦八苦しながらなんとか動かせるようにしました。下記の記事を修正してあります。 yaju3d.hatenablog.jp

次の手の予測アルゴリズム

  • チョキが多いので、グー > チョキ > パーの優先順位とする
  • 前回と違う手を出すので、上記の優先順位で勝手を選ぶ
  • 二手前と一手前が違う手なら、残りの手を出すので勝手を選ぶ
  • 三手の中に同手がある場合、 残りの手を出すので勝手を選ぶ
  • 二手前と一手前が同じ手なら、勝手を出すので負手を選ぶ

2017年の勝敗結果

年月 サザエさんの手 予想の手 勝敗結果
01月01日 休み
01月08日 チョキ グー 勝ち
01月15日 チョキ グー 勝ち
01月22日 パー パー 引き分け
01月29日 グー パー 勝ち
02月05日 グー グー 引き分け
02月12日 パー グー 負け
02月19日 チョキ グー 勝ち
02月26日 チョキ パー 負け
03月05日 グー パー 勝ち
03月12日 パー チョキ 勝ち
03月19日 休み
03月26日 チョキ グー 勝ち
04月02日 チョキ パー 負け
04月09日 グー パー 勝ち
04月16日 パー チョキ 勝ち
04月23日 パー グー 負け
04月30日 チョキ グー 勝ち
05月07日 グー パー 勝ち
05月14日 パー チョキ 勝ち
05月21日 チョキ グー 勝ち
05月28日 グー パー 勝ち
06月04日 チョキ チョキ 引き分け
06月11日 パー チョキ 勝ち
06月18日 グー パー 勝ち
06月25日 グー グー 引き分け
07月02日 チョキ グー 勝ち
07月09日 パー チョキ 勝ち
07月16日 グー パー 勝ち
07月23日 パー グー 負け
07月30日 チョキ グー 勝ち
08月06日 パー パー 引き分け
08月13日 グー パー 勝ち
08月20日 チョキ グー 勝ち
08月27日 グー チョキ 負け
09月03日 パー チョキ 勝ち
09月10日 パー グー 負け
09月17日 チョキ グー 勝ち
09月24日 グー パー 勝ち
10月01日 チョキ チョキ 引き分け
10月15日 パー チョキ 勝ち
10月22日 グー パー 勝ち
10月29日 休み
11月05日 チョキ グー 勝ち
11月12日 チョキ チョキ 引き分け
11月19日 グー チョキ 負け
11月26日 パー チョキ 勝ち
12月03日 パー グー 負け
12月10日 チョキ グー 勝ち
12月17日 グー パー 勝ち
12月24日 パー チョキ 勝ち

結果は、32勝9敗7分(勝率0.780)となりました。

ちなみに、サザエさんじゃんけん研究所 公式ウェブサイトサザエさんの手の予想と勝負結果(2017年)が29勝8敗11分(勝率0.783)でした。

今回は勝数では上回ったのですが、勝率では僅差で負けました。 勝率の計算は、「勝ち / (勝ち + 負け)」で行っているのですが、負けの1つが響いたわけです。

【2017/01/09追記】
2017年冬版 サザエさんじゃんけん白書によるとクール(四半期)の初回(1月、4月、7月、10月の初回)はチョキが出やすいとのことです。もし、これを取り入れた場合、04月02日(負け)と10月01日(引き分け)なので、34勝8敗6分(勝率0.809)になるかな。

データ分析 研究所公式
2013 24勝13敗12分(勝率0.649) 25勝9敗17分(勝率0.735)
2014 30勝10敗11分(勝率0.750) 30勝9敗12分(勝率0.769)
2015 32勝9敗9分(勝率0.780) 33勝9敗8分(勝率0.785)
2016 27勝11敗12分(勝率0.711) 22勝13敗15分(勝率0.628)
2017 32勝9敗7分(勝率0.780) 29勝8敗11分(勝率0.783)

スライド

2013年に静岡Developers勉強会で機械学習を学び、2014年1月にネタとしてSlideShareに公開しました。

最後に

今年中にはTensorFlowを使って人工知能のディープラーニングで予測手を作りたいと思います。まだ組むだけの理解度(n-gramモデル)が足りないためですが、あと少しで理解度が進むはずです。

【2018/01/07 追記】
r_stdさんが、機械学習を用いて検証してくれました。2016年に統計検定準1級に合格している方で幾つかの手法を使っています。その中で、naive bayesで33勝9敗6分と好成績を残しています。タイトルが①なので次回も期待しています。
r-std.hatenablog.com

ディープラーニング(深層学習)を理解してみる(勾配降下法:計算確認)

はじめに

前回の続きです。 yaju3d.hatenablog.jp

前回の計算が本当に合っているのか、Pythonを使って実証してみたいと思います。

プログラム

重みの更新

①現在の重みで推測値を求める

import numpy as np
a = np.array([10, 20])
b = np.array([[1,3,5],[3,1,7]])
print(np.dot(a,b))

[ 70, 50, 190]

全体の誤差を求める。
二乗誤差は、mean_squared_error の名前です。

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

y = np.array([70, 50, 190])
t = np.array([190, 330, 660])

print(mean_squared_error(y, t))

156850.0

②勾配\Delta Eを計算する
勾配は、gradient の名前です。転置は x.T とTを付けるだけなんですね。

def gradient(y, t, x):
    return np.dot((y - t),x.T)
    
t = np.array([190, 330, 660])
y = np.array([70, 50, 190])
x = np.array([[1,3,5],[3,1,7]])

print(gradient(y, t, x))

[-3310 -3930]

③重みを更新する
②で求めた勾配\Delta E を用いて、現在の重みw を更新します。学習係数η は 0.02 としています。

def weight(w, lr, gd):
    return w - lr * gd

t = np.array([190, 330, 660])
y = np.array([70, 50, 190])
x = np.array([[1,3,5],[3,1,7]])
w = np.array([10, 20])

print(weight(w, 0.02, gradient(y, t, x)))

[ 76.2 98.6]

2回目の①現在の重みで推測値を求める

a = weight(w, 0.02, gradient(y, t, x))
b = np.array([[1,3,5],[3,1,7]])
print(np.dot(a,b))

[ 372. 327.2 1071.2]

全体の誤差を求める。
二乗誤差は、mean_squared_error の名前です。

y = np.dot(a,b)
t = np.array([190, 330, 660])

print(mean_squared_error(y, t))

101108.64

というのを繰り返す。

import numpy as np

# 二乗誤差を求める
def mean_squared_error(y, t):
    return 0.5 * np.sum((y - t)**2)

# 勾配を求める
def gradient(y, t, x):
    return np.dot((y - t),x.T)

# 重みを求める
def weight(w, lr, gd):
    return w - lr * gd

# 教師データ
t = np.array([190, 330, 660])
# 重み
w = np.array([10, 20])
# 入力データ(個数)
x = np.array([[1,3,5],[3,1,7]])
y = None
for i in range(100):
    if y is not None:
        # 勾配を求める
        grad = gradient(y, t, x)
        # 重みを求める
        w = weight(w, 0.02, grad)
    # 推測値を求める
    y = np.dot(w, x)
    # 誤差を求める
    l = mean_squared_error(y, t)
    if (i + 1) % 10 == 0:
        print("step=%3d, a1=%6.2f, a2=%6.2f, loss=%.2f" % (i + 1, w[0], w[1], l))
print("Estimated: a1=%6.2f, a2=%6.2f" % (w[0], w[1]))
step= 10, a1= 78.84, a2= 48.86, loss=4531.39
step= 20, a1= 89.41, a2= 32.85, loss=564.82
step= 30, a1= 94.92, a2= 27.91, loss=264.21
step= 40, a1= 97.30, a2= 26.05, loss=217.63
step= 50, a1= 98.28, a2= 25.30, loss=209.89
step= 60, a1= 98.68, a2= 25.00, loss=208.59
step= 70, a1= 98.84, a2= 24.88, loss=208.38
step= 80, a1= 98.91, a2= 24.83, loss=208.34
step= 90, a1= 98.94, a2= 24.81, loss=208.33
step=100, a1= 98.95, a2= 24.80, loss=208.33
Estimated: a1= 98.95, a2= 24.80

実行結果
リンゴ 98.95 ≒ 99、ミカン 24.80 ≒ 25 → 99 x 2 + 25 x 4 ≒ 297

以前、TensorFlowで組んだ際の結果とほぼ同じになったのでこれでいいかな。
リンゴ 98.81 ≒ 99、ミカン 24.90 ≒ 25 → 99 x 2 + 25 x 4 ≒ 297
yaju3d.hatenablog.jp

上記サイトでTensorFlowで組んだのに近付けるために、train_x と train_y を縦ベクトル に変更してみました。

f:id:Yaju3D:20160419005112p:plain f:id:Yaju3D:20160419002003p:plain

import numpy as np

# 二乗誤差を求める
def mean_squared_error(y, t):
    return 0.5 * np.sum((y - t)**2)

# 勾配を求める
def gradient(y, t, x):
    return np.dot((y.T - t.T), x)

# 重みを求める
def weight(w, lr, gd):
    return w - lr * gd

train_x = np.array([[1., 3.], [3., 1.], [5., 7.]])
train_y = np.array([190, 330, 660]).reshape(3, 1)
y = None
for i in range(100):
    if y is not None:
        # 勾配を求める
        grad = gradient(y, train_y, train_x)
        # 重みを求める(横ベクトルを一次元に変換)
        w = weight(w, 0.02, grad).reshape(-1,)
    else:
        # 初期重み
        w = np.array([10, 20])

    # 推測値を求める
    y = np.dot(w, train_x.T)
    # 誤差を求める
    l = mean_squared_error(y, train_y)
    if (i + 1) % 10 == 0:
        print("step=%3d, a1=%6.2f, a2=%6.2f, loss=%.2f" % (i + 1, w[0], w[1], l))
print("Estimated: a1=%6.2f, a2=%6.2f" % (w[0], w[1]))

ディープラーニング(深層学習)を理解してみる(勾配降下法:計算方法)

はじめに

前回の続きです。 yaju3d.hatenablog.jp

幾つかの人工知能関連の本やWebサイトを見ても、数式やプログラムのソースリストは記載されていても、数学が苦手な自分が理解できるようになるまでの説明が無い、そんな中でも下記3つの本(Kindle)がまだ理解できそうな感じで参考になりそうである。

    

基礎的な知識から、やっと実際の計算方法に入っていきます。

勾配降下法

関数 z=x^2 + y^2 について、その最小値を与える xy の値を勾配降下法で求めてみます。
ちなみに、正解は(x, y)=(0, 0)です。
f:id:Yaju3D:20171028224142p:plain

最初に勾配を求めておきましょう。
前回記事の偏微分で説明したように、関数 z=x^{2} + y^{2}偏微分すると指数を係数にした \displaystyle\frac{\partial z}{\partial x}=2x\displaystyle\frac{\partial z}{\partial y}=2y となります。

①勾配式 \displaystyle\left(\frac{\partial z}{\partial x}, \frac{\partial z}{\partial y} \right) = (2x, 2y)

それでは、ステップを追って計算を進めます。

1.初期設定

初期位置と学習係数 \eta を適当に与えます。
今回は初期位置を(3.00,2.00)、学習係数 \eta = 0.1 とします。

No 変位ベ クトル 関数値
i x_i y_i \partial z / \partial x \partial z / \partial y \Delta x \Delta y z
0 3.00 2.00

2.変位ベクトルを算出

現在位置 (x_i, y_i) に対して、勾配式から算出し、勾配降下法の基本式から変位ベクトル \Delta x = (\Delta x_i, \Delta y_i)を求めます。
前回記事の勾配降下法に適用で、2変数関数 z=f(x,y)の勾配降下法の基本式は次のように表しました。
基本式 \displaystyle({\Delta x}, {\Delta y}) = -\eta\nabla f(x,y)

これに①勾配式を当てはめたのが次の式となります。「\cdot」は掛け算の意味です。
②変位ベクトル (\Delta x_i, \Delta y_i) = - \eta(2x, 2y) = (- \eta\cdot2x, - \eta\cdot2y)

No 変位ベ クトル 関数値
i x_i y_i \partial z / \partial x \partial z / \partial y \Delta x \Delta y z
0 3.00 2.00 6.00 4.00 -0.60 -0.40 13.00
各計算結果の求め方

勾配
\partial z / \partial x = 2 \cdot x_0 = 2 \times 3.00 = 6.00
\partial z / \partial y = 2 \cdot y_0 = 2 \times 2.00 = 4.00
変位ベクトル
\Delta x = - \eta \cdot \partial z / \partial x = -0.1 \times 6.00 = -0.60
\Delta y = - \eta \cdot \partial z / \partial y = -0.1 \times 4.00 = -0.40
関数値
z = x_0^2 + y_0^2 = 3.00^2 + 2.00^2 = 9.00 + 4.00 = 13.00

3.位置を更新

勾配降下法に従って、現在位置 (x_i, y_i) から移動先 (x_{i+1}, y_{i+1}) の点を次の式から求めます。
②移動先 (x_{i+1}, y_{i+1}) = (x_i, y_i) + (\Delta x_i, \Delta y_i)

No 変位ベ クトル 関数値
i x_i y_i \partial z / \partial x \partial z / \partial y \Delta x \Delta y z
0 3.00 2.00 6.00 4.00 -0.60 -0.40 13.00
1 2.40 1.60

移動先
x_1 = x_0 + \Delta x_0 = 3.00 + (-0.60) = 2.40
y_1 = y_0 + \Delta y_0 = 2.00 + (-0.40) = 1.60

4.2と3の繰り返し

2と3の繰り返し(6~27は省略)、30回繰り返したときの座標 (x_{30}, y_{30}) の値です。
正解の (x, y) = (0, 0) と一致します。

No 変位ベ クトル 関数値
i x_i y_i \partial z / \partial x \partial z / \partial y \Delta x \Delta y z
0 3.00 2.00 6.00 4.00 -0.60 -0.40 13.00
1 2.40 1.60 4.80 3.20 -0.48 -0.32 8.32
2 1.92 1.28 3.84 2.56 -0.38 -0.26 5.32
3 1.54 1.02 3.07 2.05 -0.31 -0.20 3.41
4 1.23 0.82 2.46 1.64 -0.25 -0.16 2.18
5 0.96 0.66 1.97 1.31 -0.20 -0.13 1.40
28 0.01 0.00 0.01 0.01 0.00 0.00 0.00
29 0.00 0.00 0.01 0.01 0.00 0.00 0.00
30 0.00 0.00 0.01 0.00 0.00 0.00 0.00

Excel計算で小数誤差により微妙に値が違います。

バイアスについて

バイアス(bias)とは、一般に真値からの偏り、つまり系統的な誤差を指す。

切片(せっぺん)

ニューラルネットワークのパラメーターは重みとバイアスがセットとなります。
バイアスがイメージしやすいものとして、回帰直線があります。
f:id:Yaju3D:20171105230836p:plain
回帰直線は次のような1次式で表現されます。
回帰方程式 y = ax + b
f:id:Yaju3D:20171106001201p:plain
a を回帰係数(傾き)、bを切片の呼びます。切片が無いと必ず原点を通すことになってしまいます。
傾きが求まったところで切片で上下位置の調整をします。この切片がバイアスのことになります。

閾値(しきいち)

f:id:Yaju3D:20170528144340p:plainf:id:Yaju3D:20170528154127p:plain
左図が本物の神経細胞(ニューロン) で、右図が形式ニューロンです。
簡単に説明すると、入力が2つあり各入力に対して重みが掛け算され、その値が閾値を超えれば出力は「1」、そうでなければ出力は「0」となります。 たとえば、入力が(1,0)、重みが(0.5, 0.7)だとすると、1×0.5 + 0×0.7 = 0.5 を計算して閾値と比較します。
閾値未満だと出力無し(発火なし)、閾値を超えると出力有り(発火あり)となります。

出力信号無し (y=0):w_1x_1+w_2x_2 \lt \theta
出力信号有り (y=1):w_1x_1+w_2x_2 \geqq \theta

ここで \thetaニューロン固有の閾値です。 \theta を左に移行させた発火の式は次のように表現することができます。
発火の式 y=a(w_1x_1+w_2x_2 - \theta)
a は活性化関数となります。

活性化関数にシグモイド曲線を使った場合の閾値 \theta は生物的にはニューロンの個性を表現する値です。
\theta が大きければ興奮しにくく(すなわち鈍感)、小さければ興奮しやすい(すなわち敏感)という感受性を表します。
f:id:Yaju3D:20171106013953p:plain

発火の式 y=a(w_1x_1+w_2x_2 - \theta) にて、\theta だけマイナス記号が付いているのは数学的に美しくありません。美しさが欠けることは数学が嫌うところです。また、マイナスは計算ミスを誘発しやすいという欠点を持ちます。そこで、-\thetab と置き換えたのが次の式となります。
y=a(w_1x_1+w_2x_2 + b)
こうすれば式として美しく、計算ミスも起こりにくくなります。
この b をバイアス(bias)と呼びます。

バイアスの式表現

数式では b とするよりは w_0 として、次の式にします。
y=a(w_0 + w_1x_1+w_2x_2)

一般的に書き直すと、重みをwではなく \theta として、\theta_0 をバイアスとした場合、入力値 x とすると\theta_0x で次元が違うと扱いにくいので、最初の要素に 1 をセットする。
f:id:Yaju3D:20180103202510p:plain

それを x_0=1 と定義して、 x の最初の要素に x_0 を置くほうがより数学上では綺麗となる。
f:id:Yaju3D:20180103210343p:plain

\theta を転置したものと x を掛けたものを計算すると次の式になります。
\theta^Tx = \theta_0x_0 + \theta_1x_1+ \theta_2x_2 + \cdot + \theta_nx_n
更にこれを書き直すと 、f_\theta(x) = \theta^Tx と簡易的な表現の式となる。

やる夫で学ぶ機械学習 - 多項式回帰と重回帰 - · けんごのお屋敷

リンゴとミカン

リンゴとミカンを例題に勾配降下法を計算していきます。上記と違うのは、訓練データがあることです。 f:id:Yaju3D:20160417230234p:plain
f:id:Yaju3D:20160419005112p:plain f:id:Yaju3D:20160419002003p:plain

以前の記事を参考にします。
yaju3d.hatenablog.jp

前提

f:id:Yaju3D:20171029165625p:plain
1層の全結合ニューラルネットワークを用いて、勾配降下法による重みの更新例を示します。
入力層ユニット数は2、出力層はユニット数が1のシンプルなネットワークです。損失関数は二乗誤差を用います。

No 入力 データ 教師データ
i x_1 x_2 t
1 1 3 190
2 3 1 330
3 5 7 660

上表はトレーニングデータセットです。サンプル数は3で、1番目のサンプルを見ると、x_1=1x_2=3が与えられたときの教師データ t=190 (正解)になっています。

勾配\Delta Eの計算式

出力層の出力値 y は、入力層の出力値 x を用いて、次のような一次多項式で表すことができます。
a_1x_1+a_2x_2+b=y
パラメーター a_1,a_2,b

これにトレーニングデータセットを当てはめると、次のような連結方程式が出来上がります。
a_1+3a_2+b=y_1
3a_1+a_2+b=y_2
5a_1+7a_2+b=y_3

この連結方程式は行列とベクトルを用いて次のように表すことが出来ます。

\begin{pmatrix}a_1 & a_2 & b \end{pmatrix} 
\begin{bmatrix}
1 & 3 & 5 \\ 
3 & 1 & 7 \\ 
1 & 1 & 1 \\ 
\end{bmatrix} = \begin{pmatrix}y_1 & y_2 & y_3 \end{pmatrix}  
    w     X       Y
   重み  入力データ 出力データ(推測値)

今回はバイアスbは 0 とするので除外しました。そうしないリンゴとミカンの金額(重み)を求めたいのにバイアスを含む3つ値が求まってしまいます。

\begin{pmatrix}a_1 & a_2 \end{pmatrix} 
\begin{bmatrix}
1 & 3 & 5 \\ 
3 & 1 & 7 \\ 
\end{bmatrix} = \begin{pmatrix}y_1 & y_2 & y_3 \end{pmatrix}  
    w     X       Y
   重み  入力データ 出力データ(推測値)

誤差を含めた式と行列を表現すると下記のようなります。
y=a_1x_1 + a_2x_2 + 誤差 e

\begin{bmatrix}
190 \\ 
330 \\ 
660
\end{bmatrix}=
\begin{bmatrix}
1 & 3 \\ 
3 & 1 \\ 
5 & 7
\end{bmatrix}
\begin{bmatrix}
a_1\\ 
a_2
\end{bmatrix} +
\begin{bmatrix}
e_1\\ 
e_2\\ 
e_3
\end{bmatrix} 

入力データを X、重み(パラメーター)を w、出力データ(推測値) Y としています。wX を用いると、出力データ Y は次の式で表すことができます。
wX=Y
また、教師データは、 t を用いて次のように表します。
t=(190 \quad 330 \quad 660)

損失関数には二乗誤差を使用します。誤差は次のような式で表すことができます。
E=\displaystyle \frac{1}{2} (Y - t)^2
\quad=\displaystyle \frac{1}{2} (wX - t)^2
勾配\Delta Eは、誤差 E を重み w微分すると、勾配\Delta Eを次のような式になります。
\Delta E = \displaystyle \frac{\partial E}{\partial w} = (Y - t)X^T
X^TX の転置行列、(Y - t):誤差信号 \delta

今回、(Y - t)を誤差信号と呼び、記号\delta (デルタ)で表します。
勾配\Delta Eは、誤差信号 \delta と、入力データX から求めることができます。入力データはすでにわかっている値なので、誤差信号の値を求めれば、勾配\Delta Eを求めることができ、そして勾配\Delta Eがわかれば、式 w \leftarrow w-\eta\Delta E で重みを w を更新することができます。\eta:学習係数

初期値の設定

初めに、重みw と学習係数\etaに適当な初期値を設定します。ここでは学習係数\etaは 0.02 とし、重みwの初期値を10円と20円から、バイアスは 0 にして次のように設定します。
w=(a_1 \quad a_2) = (10 \quad 20)

重みの更新

①~③の手順で重みの更新を行います。

①現在の重みで推測値を求める

入力データXと現在の重みwから、推測値Yを求めます。
推測値

wX=\begin{pmatrix}10 & 20\end{pmatrix} 
\begin{bmatrix}
1 & 3 & 5 \\ 
3 & 1 & 7 \\ 
\end{bmatrix}
\quad\quad=(1 \times 10 + 3 \times 20 \quad 3 \times 10 + 1 \times 20 \quad 5 \times 10 + 7 \times 20)  
\quad\quad=(70 \quad 50 \quad 190)  
\quad\quad=Y  

ここで、全体の誤差をいったん計算してみます。損失関数には二乗誤差を使用するので、全体の誤差Eは次のように計算します。

現在の値
教師データ t=(190 \quad 330 \quad 660)
推測値   Y=( 70 \quad  50 \quad 190)
全体の誤差
E=\displaystyle \sum_{n=1}^{3} \frac{1}{2}(y_n - t_n)^2
\quad=\displaystyle \frac{1}{2}(y_1 - t_1)^2 + \frac{1}{2}(y_2 - t_2)^2 + \frac{1}{2}(y_3 - t_3)^2
\quad=\displaystyle \frac{1}{2}(70 - 190)^2 + \frac{1}{2}(50 - 330)^2 + \frac{1}{2}(190 - 660)^2
\quad=\displaystyle \frac{1}{2}(-120)^2 + \frac{1}{2}(-280)^2 + \frac{1}{2}(-470)^2
\quad=\displaystyle 7200 + 39200 + 110450
\quad=\displaystyle 156850

②勾配\Delta Eを計算する

現在の重みwに対する勾配\Delta Eを計算します。
入力データを転置行列X^Tにするのは、行列の掛け算の定義で「行列Aの列数(横の個数)と、行列Bの行数(縦の個数)が等しくないと掛け算できない」ためです。
行列A(Y-t) 1行3列 と 行列B(X^T) 3行2列 にして掛け算させます。行列の積

現在の値
教師データ t=(190 \quad 330 \quad 660)
推測値   Y=( 70 \quad  50 \quad 190)

入力データ X^T = \begin{bmatrix}
1 & 3 \\ 
3 & 1 \\ 
5 & 7 \\ 
\end{bmatrix}  

勾配
\Delta E=(Y-t)X^T
\quad\quad=(70-190 \quad 50-330 \quad 190-660)\begin{bmatrix}
1 & 3 \\ 
3 & 1 \\ 
5 & 7 \\ 
\end{bmatrix}
\quad\quad=(-120 \quad -280 \quad -470)\begin{bmatrix}
1 & 3 \\ 
3 & 1 \\ 
5 & 7 \\ 
\end{bmatrix}
\quad\quad=(-120 \times 1 + -280 \times 3 + -470 \times 5 \quad -120 \times 3 + -280 \times 1 + -470 \times 7)
\quad\quad=(-3310 \quad -3930)
③重みを更新する

②で求めた勾配\Delta Eを用いて、現在の重みwを更新します。学習係数\etaは 0.02 としています。

現在の値
現在の重み w=(10 \quad 20)
勾配   \Delta E=(-3310 \quad -3930)
更新後の重み
w \leftarrow w-\eta\Delta E
\quad=(10 \quad 20) - 0.02 \times (-3310 \quad -3930)
\quad=(10 \quad 20) - (-66.2 \quad -78.6)
\quad=(10 - (-66.2) \quad 20 - (-78.6))
\quad=(76.2 \quad 98.6)

これで重みwが更新されました。ここで、更新後の重みwを使って推測値を求め、全体の誤差を再び計算してみます。

現在の値
教師データ t=(190 \quad 330 \quad 660)

推測値

wX=\begin{pmatrix}76.2 & 98.6 \end{pmatrix} 
\begin{bmatrix}
1 & 3 & 5 \\ 
3 & 1 & 7 \\ 
\end{bmatrix}
\quad\quad=(1 \times 76.2 + 3 \times 98.6  \quad 3 \times 76.2 + 1 \times 98.6 \quad 5 \times 76.2 + 7 \times 98.6  
\quad\quad=(372 \quad 327.2 \quad 1071.2)  
\quad\quad=Y  

全体の誤差
E=\displaystyle \sum_{n=1}^{3} \frac{1}{2}(y_n - t_n)^2
\quad=\displaystyle \frac{1}{2}(y_1 - t_1)^2 + \frac{1}{2}(y_2 - t_2)^2 + \frac{1}{2}(y_3 - t_3)^2
\quad=\displaystyle \frac{1}{2}(372 - 190)^2 + \frac{1}{2}(327.2 - 330)^2 + \frac{1}{2}(1071.2 - 660)^2
\quad=\displaystyle \frac{1}{2}(182)^2 + \frac{1}{2}(-2.8)^2 + \frac{1}{2}(411.2)^2
\quad=\displaystyle 16562 + 3.92 + 84542.72
\quad=\displaystyle 101108.64

全体の誤差は 101108.64 となり、重みの誤差前の 156850 より小さくなっていることが分かります。

更新の繰り返し

以上の①~③の計算で、全サンプル(今回は2サンプル)の更新が1回終わりました。これで1エポックを終了したことになります。
更新後の重みwを用いて、①~③をもう一度繰り返せば、2エポックが終了となります。

最後に

計算の繰り返しはコンピューターの得意なところです。
次回はこの計算が本当に合っているのか、Pythonを使って実証してみたいと思います。

ディープラーニング(深層学習)を理解してみる(勾配降下法:ベクトル、内積、微分、偏微分)

はじめに

前回の続きです。 yaju3d.hatenablog.jp

幾つかの人工知能関連の本やWebサイトを見ても、数式やプログラムのソースリストは記載されていても、数学が苦手な自分が理解できるようになるまでの説明が無い、そんな中でも下記3つの本(Kindle)がまだ理解できそうな感じで参考になりそうである。

    

計算方法に入る前にベクトル等の基礎知識を先に説明していきます。

ベクトル

ベクトルとは大きさと向きを持つ量と定義されます。
f:id:Yaju3D:20171001213101p:plain

ベクトルの矢を座標平面上に置くことで座標のように表現できます。矢の始点を原点に置き、矢の終点の座標でそのベクトルを表します。これをベクトルの成分表示といいます。
成分表示されたベクトル a で下図のように表現します。また3D座標上でもz軸方向を追加して同じように表現します。
f:id:Yaju3D:20171001225146p:plain

下図のように矢印がたくさんある図を見たことがありませんか?
各データは、どこかの一点のデータを基準にして、大きさと向きを持つベクトルとなります。
f:id:Yaju3D:20171001213359p:plain

内積

ベクトルにも足し算と引き算があるように掛け算も…といいたいところですが、微妙に違うのが内積です。内積だとベクトル同士を掛け算っぽくするのに結果はベクトルではなくスカラという1つの値となります。内積以外にも外積がありますが、ここでは説明しません。

過去に3Dの勉強用に書いた記事がありますので、興味があれば読んでみて見ください。

なぜ内積の計算をそうするのかは次のサイトの説明が分かりやすいと思います。
www.yukisako.xyz

内積というのは2つのベクトルの“近さ”を表しているという程度に理解しておけばいいです。
2つのベクトルの近さを表すのにcos\thetaが使われます。

f:id:Yaju3D:20171001231022p:plain

高校までの数学だと「内積の値がcos\thetaを使って表される」と教わるが、逆にcos\thetaの値が内積によって定義されるのだ」と考えるようにするといいでしょう。

cos\thetaの値は2つのベクトルの向きが反対なら「-1」、同じであれば「1」となる。 f:id:Yaju3D:20171002005747p:plain

(ア)2つのベクトルが反対向きのときに内積は最小値
(イ)2つのベクトルが平行でないとき、内積は平行の場合の中間の値
(ウ)2つのベクトルが同じ向きのときに内積は最大値

そして、(ア)の内積が最小値というのが、勾配降下法の基本原理となります。

手書き数字認識でどの数字か分類できるのは、より似ているかで判断しています。
方向が似ているベクトルを「似ている」と判断するなら、2つのベクトルが似ていると内積は大きくなるのです。
f:id:Yaju3D:20171002005805p:plain

内積の計算

平面(2次元)
f:id:Yaju3D:20171002005819p:plain

立方空間(3次元)
f:id:Yaju3D:20171002005835p:plain

内積は2次元や3次元に限らず、何次元のベクトルについても計算することができます。現実世界では4次元や10次元ベクトルで表せるものは存在しないですが、数学上は計算できるところが数学のおもしろいところでもあります。

この内積が何次元(多次元)であっても計算できるということが重要なポイントです。
リンゴとミカンの値段を求める際は2変数です、それにブドウを付け加えると3変数になり、レモンも付け加えると4変数になります。変数が多くなることを多変数と言います。この多変数というのは機械学習における多次元と同じ意味を持っています。

参照

微分

微分とは、ある関数の各点における傾き(変化の割合)のことです。
そして勾配降下法においては、傾きの最小値を求めることが目的です。 f:id:Yaju3D:20171007211313p:plain

最小値

微分において、傾きが水平になる 0 が最小値といいたいところですが、下図右のように最大値においても水平になり 0 となります。
f:id:Yaju3D:20171008132102p:plain

このことは、勾配降下法においては厄介な性質になります。
f:id:Yaju3D:20170917223314p:plain

導関数

関数 y=x^{2}微分して得られた導関数は指数を係数にした {y}'=2x となります。
微分の定義式では極限の記号\displaystyle \lim_{h \to 0}を使って、x の変化量 h を限りなく 0 に近づけます。

\displaystyle \lim_{h \to 0}\frac{f(x+h)-f(x)}{h}={f}'(x)

※現代の数学では、0 で割る(分母が0)ことは禁止事項とされているので、極限まで 0 に近づければいいという屁理屈です。

ちなみに、'の名称は「プライム」で、\limの名称は「リミット」と読み、極限(limit)を取る記号です。

https://sci-pursuit.com/math/differential-1.htmlsci-pursuit.com

偏微分

関数 y=x^{2} は変数が1つでしたが、リンゴとミカンの値段を求める場合は2変数による関数z=f(x,y)となります。このように変数が2つ以上の関数を多変数関数といいます。
多変数関数でも微分法を適用できます。ただし変数が複数あるので、どの変数について微分するかを明示しなければなりません。この意味で、ある特定の変数について微分することを偏微分(partial derivative)といいます。
例えば、2変数x,yから成り立つ関数を考えてみましょう。変数xだけに着目してyは定数と考える微分を「xについての偏微分」と呼び、次の記号で表します。

\displaystyle\frac{\partial z}{\partial x}=\frac{\partial f(x,y)}{\partial x}=\lim_{x \to 0}\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}

yについての偏微分も同様です。

\displaystyle\frac{\partial z}{\partial y}=\frac{\partial f(x,y)}{\partial y}=\lim_{y \to 0}\frac{f(x,y+\Delta y)-f(x,y)}{\Delta y}

ちなみに、\partialの名称は「デル」、\Deltaの名称は「デルタ」です。

関数 z=x^{2} + y^{2}偏微分すると指数を係数にした \displaystyle\frac{\partial z}{\partial x}=2x\displaystyle\frac{\partial z}{\partial y}=2y となります。

remedics.air-nifty.com

www.youtube.com

勾配降下法に適用

上述で説明してきたベクトルと内積偏微分が勾配降下法にどのように使われるのか?

勾配降下法は、少しずつ下りながら、場所ごとに最も急な勾配を探すことになります。
P0から最も勾配の急な点P1の位置を求める。その位置P1から更に最も勾配の急な点P2の位置を求める。このように繰り返すことで最も早く最小点にたどり着くようになります。 f:id:Yaju3D:20170924211620p:plain

関数z=f(x,y)において、xをΔxだけ、yをΔyだけ変化させたとき、関数z=f(x,y)の変化は下記式となります。
\displaystyle \Delta z = f(x + {\Delta x},y + {\Delta y}) - f(x,y)
これを近似公式から次の関係が成立します。
\displaystyle \Delta z = \frac{\partial f(x,y)}{\partial x}{\Delta x} + \frac{\partial f(x,y)}{\partial y}{\Delta y}

f:id:Yaju3D:20171010235143p:plain

そして、これは2つのベクトルの内積を表す式になっているのです。 f:id:Yaju3D:20171010234406p:plain

ベクトルが反対向きになっている時に、内積が最小値となっている。 f:id:Yaju3D:20171011001024p:plain f:id:Yaju3D:20171011001230p:plain

ベクトルが反対向きは、-1となり、勾配降下法の基本式は次のように表せます。(\etaは正の小さな定数)
\displaystyle({\Delta x}, {\Delta y})=-\eta\left(\frac{\partial f(x,y)}{\partial x}, \frac{\partial f(x,y)}{\partial y} \right)

ちなみに、\etaはイータと読むギリシャ文字です。ローマ字のiに対応します。
なぜループカウンタ変数のほとんどに “i”が使用されるのか? - Qiita

重要なのは、内積は何次元(何変数)のベクトルについても計算することができます。

\displaystyle({\Delta x_1}, {\Delta x_2}, \cdots, {\Delta x_n})=-\eta\left(\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \cdots, \frac{\partial f}{\partial x_n} \right)

式が長いので数学のベクトル解析で用いられるハミルトン演算子\nablaのナブラを使って省略します。
\displaystyle \nabla f = \left(\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \cdots, \frac{\partial f}{\partial x_n} \right)

2変数関数 z=f(x,y)の勾配降下法の基本式は次のように表せます。
\displaystyle({\Delta x}, {\Delta y}) = -\eta\nabla f(x,y)

3変数関数 z=f(x,y,z)の勾配降下法の基本式は次のように表せます。
\displaystyle({\Delta x}, {\Delta y}, {\Delta z}) = -\eta\nabla f(x,y,z)

次回は計算方法について説明していきます。

参照

ディープラーニング(深層学習)を理解してみる(勾配降下法:最急降下法と確率的勾配降下法)

はじめに

前回の続きです。 yaju3d.hatenablog.jp

勾配降下法をどうして使うのかは理解できたのですが、その計算方法がまだ理解が足りてなくて、微分の本とかを読んでいました。
数学は苦手なんですが、理解はしたい。おまじないとかそういうルールだからとかで済ましたくないんですよね。

勾配降下法の他にも最急降下法とか確率的勾配降下法の用語がありますので、この違いも理解したい。
今回は下記の本(Kindle)を参考にしている。

  

勾配降下法

f:id:Yaju3D:20170917211137p:plain
勾配には、傾斜の程度・斜面という意味があります。
上図では1変数による関数y=f(x)でグラフは放物線となります。

リンゴとミカンの値段を求める場合は2変数による関数z=f(x,y)となり下図のようにワイングラスの底のような形になります。ちなみに3変数以上は図で表わすのは困難です。
f:id:Yaju3D:20170924173848p:plain
下図は関数のグラフの一部を拡大し、斜面に見立てた図となります。その斜面上のある点Pにピンポン玉を置き、そっと手を放すと玉は最も急な斜面を選んで転がり始めます。
f:id:Yaju3D:20170924181321p:plain

この操作を何回も繰り返せば、ピンポン玉は最短な経路をたどってグラフの底、すなわち関数の最小値にたどり着くはずです。この玉の動き(誤差関数という坂道、つまり「勾配」を下へ下へと「降下」していく)を真似たのが勾配降下法です。
f:id:Yaju3D:20170924183956p:plain

正確には最小値を探す場合を勾配降下法(gradient descent mothod)と言い、最大値を探す場合を勾配上昇法(gradient ascent mothod)と言います。descentは、SQLで降順にソートさせる時に略称(desc)で、order by column descと使いますね。

最短で下るには、少しずつ下りながら、場所ごとに最も急な勾配を探すことになります。
P_0から最も勾配の急な点P_1の位置を求める。その位置P_1から更に最も勾配の急な点P_2の位置を求める。このように繰り返すことで最も早く最小点にたどり着くようになります。
f:id:Yaju3D:20170924211620p:plain

数値解析の分野では勾配降下法を最急降下法と呼びますが、勾配降下法の中にもいくつかの方法が存在します。

最急降下法とは

最急降下法は学習データのすべての誤差の合計を取ってからパラメーターを更新します。学習データが多いと計算コストがとても大きくなってしまいます。また、学習データが増えるたびに全ての学習データで再学習が必要となってしまいます。

確率的勾配降下法とは

確率的勾配降下法は学習データをシャッフルした上で学習データの中からランダムに1つを取り出して誤差を計算し、パラメーターを更新をします。勾配降下法ほどの精度は無いが増えた分だけの学習データのみで再学習する(重みベクトルの初期値は前回の学習結果を流用)ため再学習の計算量が圧倒的に低くなります。
運が良ければ最急降下法よりも早く最適解にたどり着けますが、運が悪ければいつまで経っても適した答えを導けません。

ミニバッチ確率的勾配降下法とは

ミニバッチ確率的勾配降下法最急降下法確率的勾配降下法の間を取ったような形となります。 最急降下法では時間がかかりすぎ、確率的勾配降下法では一つ一つのデータにかなり揺さぶられることになるので、学習データの中からランダムにいくつかのデータを取り出して誤差を計算、パラメーターを更新をします。このときの一回に取り出すデータの数をバッチサイズと呼びます。

最急降下法確率的勾配降下法の違い

パラメーター更新の式

最急降下法のパラメーター更新の式
\displaystyle \theta_j := \theta_j - \eta\sum_{i=1}^n(f_{\theta}({x}^{(i)}) - y^{(i)}){x_j}^{(i)}

確率的勾配降下法のパラメーター更新の式
\displaystyle  \theta_j := \theta_j - \eta(f_{\theta} ({x}^{(k)}) - y^{(k)}){x_j}^{(k)}
(式中のkは、パラメーター更新毎にランダムに選ばれたインデックス)

大きな違いとして、確率的勾配降下法ではシグマ\sum (1~nまで合計)が取れています。
その分、計算コストは少なくなっていることが分かりますよね。

最小値へのたどり着き方

最急降下法がほぼ直線的に最小値にたどり着くのに対し、確率的勾配降下法はランダムなのでくねくねしながら最小値にたどり着きます。
f:id:Yaju3D:20170925010546p:plain

局所解に陥る可能性

誤差関数が素直なカーブのみならいいですが、下図のようにいったん下がってまた上がるみたいなうねうねした場合があります。本当は④が最小値なのに別の小さな谷(これを局所解といいます)に捕まってしまうことがあります。
f:id:Yaju3D:20170917223314p:plain

確率的勾配降下法はこのうねうねのいろいろなところから勾配を下ろうとするため、最急降下法よりも局所解に陥る可能性が小さくなります。

定数η(イータ)は学習係数と呼ばれます。このηは人が移動する際の「歩幅」と見立てられます。このηで決められた値に従って次に移動する点が決められるからです。その歩幅が大きいと最小値に達しても、それを飛び越えてしまう危険があります(下図左)。歩幅が小さいと、極小値で停留してしまう危険があります(下図右)。
f:id:Yaju3D:20170925010945p:plain

最後に

言葉で説明だけだと分かりにくいので、実際にPythonでプログラムされているサイトを見つけました。 sinhrks.hatenablog.com

これを、組み直ししてみました。 qiita.com

次回は、勾配降下法を実際に計算してみます。

参照

スポンサーリンク