それでは毛玉諸君、これにて失敬

日々の精進を備忘録的に綴ります。

ベルマンフォード法とダイクストラ法の概念を完全に理解する

概要

社内の勉強会でダイクストラ法(dijkstra)とベルマンフォード法(bellman-ford)の紹介をしました。その際にアルゴリズムの証明を色々調べていましたが、解説されている記事があまり多くないことに気づきました。 「なぜそうなるの?」という部分は個人的に大事な要素だと思うので、その辺も踏まえてこれらのアルゴリズムを紹介したいと思います。

筆者はどんな人?

プログラミング歴:1年ちょい(Pythonのみ) Atcoder茶(よわよわ) 2020年4月からデータ分析委託の企業に勤めているが、なぜかアルゴリズムの勉強ばっかりやっている。

このアルゴリズムはいつ使うの?

重み付きグラフの単一始点最短経路問題を解く際に使います。

それぞれの用語ですが、 重み付き:辺を通るのに必要なコスト 単一始点:スタート地点が一か所 最短経路:スタート地点から目的地まで最小コストで向かう経路 です。

下のグラフを例にすると、 image.png スタート地点$A$から全ての辺への最短経路を求めることが目的になります。 ちなみに$A→E$の最短経路は $(A→D→B→F→E)$ です。これをアルゴリズムを使用して求めていきます。

アルゴリズムの前に、グラフの大事な性質

1.最短路の部分道は最短路

最短路の部分道はどこをとっても最短路になります。

<証明> image.png 仮に上のグラフで$A→E$の最短路を $(A→B→C→D→E cost:5)$ とします。

このうち$(A→B→C cost:3)$ はこの最短路の部分道ですが、この区間にはB'という頂点を利用すれば コストを1少なく移動できます。 するとこの区間の最短路は $(A→B'→C cost:2)$ となるため、$A→E$を $(A→B'→C→D→E cost:4)$ で移動できてしまい、仮定と矛盾します。 最短路の部分道は必ず最短路になります。

2.経路緩和性

始点に関連する頂点から順番に最短路を求めていくと、最終的にすべての最短経路が求められます。 これは結局、始点から徐々に求めた最短路がその先に到達する頂点の部分道となるためです。

ベルマンフォード法とは

全頂点間の最短路を求めます。ダイクストラ法と大きく異なる点は、負の重みが許容されることです。

最短路を求めるプロセスは、

  1. 全ての辺に注目し、その辺を通った先の頂点のコストを小さくできれば値を更新
  2.  1.のプロセスを(頂点の数-1)回行う です。

下の図だと、

image.png

1.のプロセスで、$A$を出発し、$C$にコスト$5$で到達できるので、$C:5$とメモします。 次に、$A→B$は$-1$で到達できるので$B:-1$、$B→C$にはコスト$-1+2=1$で到達できます。 この時、$B$を経由するとより小さなコストでCに到達できることが分かったので、$C:1$と更新します。 この作業を全ての辺に対して行います。これが1回目。 あとは(頂点の数-1)回繰り返せば、$A$から出発して全ての頂点への最短路を求めることが出来ます。

なぜ(頂点の数-1)回行うの?

始点から目的地までに通る頂点の最大の個数は、始点を除く頂点、すなわち(頂点の数-1)個です。 前述したプロセスでは、1回目のループで始点から移動できる頂点)の最短路が必ず求められます。 2回目のループでは(1回目で決まったルートから移動できる頂点)の最短路が必ず求められます。 ・ ・ ・ 結局このループは、始点をスタートして順番に最短路を求める作業になるので、経路緩和性により最短路を求めることができます。

負の重みによって困る場合

下の図を見てください。 image.png 丸で囲んだ部分はサイクルになっていますが、ここを1周するとコストが-2下がります。 つまりこのサイクルを使えば無限にコストを下げることが可能です。 これを負閉路といいます。

ベルマンフォード法では、負閉路の検出も可能です。 先ほど、(頂点の数-1)回のループで最短路を求めることを確認しましたが、再度全ての頂点のコストを計算し、コストを更新できる頂点がある場合は負閉路が存在するということになります。

いいからソースコード見せろよ

隣接リスト形式(競プロとかでよく利用する、隣接行列というのもある)の入力に対応し、ベルマンフォード法により最短路を検出します。 コードの下に負閉路のある入力例を示します(上の図のグラフ)。負閉路を通る頂点のコストを$-inf$に置き換えています。

クリックしたら見れるよ

def shortest_path(s,n,es):
    #s→iの最短距離
    # s:始点, n:頂点数, es[i]: [辺の始点,辺の終点,辺のコスト]
    
    #重みの初期化。ここに最短距離をメモしていく
    d = [float("inf")] * n
    
    #スタート地点は0
    d[s] = 0

    #最短路を計算
    for _ in range(n-1):
        # es: 辺iについて [from,to,cost]
        for p,q,r in es:
            #矢印の根が探索済、かつコストができるルートがあれば更新
            if d[p] != float("inf") and d[q] > d[p] + r:
                d[q] = d[p] + r
    
    #負閉路をチェック
    for _ in range(n-1):
        # e: 辺iについて [from,to,cost]
        for p,q,r in es:
            #更新される点は負閉路の影響を受けるので-infを入れる
            if d[q] > d[p] + r:
                d[q] = -float("inf")
    return d
################
#入力
n,w = map(int,input().split()) #n:頂点数 w:辺の数

es = [] #es[i]: [辺の始点,辺の終点,辺のコスト]
for _ in range(w):
    x,y,z = map(int,input().split())
    es.append([x,y,z])
    # es.append([y,x,z])

print(shortest_path(0,n,es))

'''
入力
7 12
0 1 2 
0 2 -5
0 3 5
1 2 3
1 4 1
1 5 3
2 4 6
3 2 -2
3 5 4
4 6 -5
5 4 -2
6 2 -3
'''
>>[0, 2, -inf, 5, -inf, 5, -inf]

計算量は$O(V2+E)$(頂点の数:$V$、辺の数:$E$)です。

ダイクストラ法とは

非負の重みを持つグラフで使用できます。

最短路を求めるプロセスはヨビノリさんの動画 )が死ぬほど分かりやすいです。 (ヨビノリさんの守備範囲広すぎる…)

一応こちらでもご説明します。 image.png

頂点のコストを(初期値:$S = 0$, その他$ = ∞$)とします。

  1. 一番コストが小さい頂点を始点とします。最初は$S$が始点です。
  2. 始点のコストを確定します(図では緑になる)。この頂点のコストはこれ以降更新されません。
  3. 始点から移動可能な頂点に出発点のコスト+辺重みを頂点にメモします。中段列の頂点に$7、3、1$とメモされました。コストはまだ更新されるかもしれないので未確定です。

  4. すべての頂点が確定するまで1~3を繰り返します。メモされた頂点の中で1のコストを持つ頂点が一番小さいので確定し、その頂点から移動可能な頂点をメモします。もし、より小さなコストで移動できるのならコストを更新します。

なぜ最短路が求められる?

未確定の頂点のうち最小のコストを持つ頂点は、それ以上小さなコストに更新することができないためです。 非負の重みしか存在しないので、別ルートを通るにしても0以上のコストが発生します。そのため未確定のうち、最小のコストは最短路と確定することが出来るのです。

で、ソースコードは?

クリックしたら見れるよ

import heapq

def dijkstra(s, n, es):
    #始点sから各頂点への最短距離
    #n:頂点数, cost[u][v] : 辺uvのコスト(存在しないときはinf)
    d = [float("inf")] * n
    #探索リスト
    S = [s]
    #優先度付きキューを使うと計算量が少なくなるよ
    heapq.heapify(S)

    d[s] = 0
    
    while S:
        #最小コストの頂点を取り出し
        v = heapq.heappop(S)
        # print("-----------v", v)
        for u, w in es[v]:
            # print(u, w)
            #より小さなコストで到達可能であれば更新し、リストに追加
            if d[u] > d[v]+w:
                d[u] = d[v]+w
                heapq.heappush(S, u)
        # print("d", d)
    return d

n,w = map(int,input().split()) #n:頂点数 w:辺の数

cost = [[float("inf") for i in range(n)] for i in range(n)] 
es = [[] for _ in range(n)]

for i in range(w):
    x, y, z = map(int, input().split())
    es[x].append([y, z])
    es[y].append([x, z]) #有向グラフの場合はこの行を消す
# print(es)
print(dijkstra(0, n, es))

オリジナルの計算量はベルマンフォード法と同じ$O(V2)$ですが、優先度付きキューを使用すれば$O((E+V)logV)$に削減できます。

参考文献とかサイト

アルゴリズムイントロダクション 第3版 第2巻: 高度な設計と解析手法・高度なデータ構造・グラフアルゴリズム (世界標準MIT教科書) (T.コルメン,R.リベスト,C.シュタイン,C.ライザーソン,浅野哲夫,岩野和生,梅尾博司,山下雅史,和田幸一) かなり重い本ですが、グラフ理論の証明がとても丁寧でした。

グラフ理論⑤(ダイクストラアルゴリズム) (予備校のノリで学ぶ「大学の数学・物理」) https://www.youtube.com/watch?v=X1AsMlJdiok こちらの動画を見るだけで、ダイクストラ法は完全に理解できます!

・蟻本 python 単一最短経路法2(ダイクストラ法) 競技プログラミング (じゅっぴーさんのブログ、コードはこちらを参考にさせて頂きました) https://juppy.hatenablog.com/entry/2018/11/01/%E8%9F%BB%E6%9C%AC_python%E5%8D%98%E4%B8%80%E6%9C%80%E7%9F%AD%E7%B5%8C%E8%B7%AF%E6%B3%952%EF%BC%88%E3%83%80%E3%82%A4%E3%82%AF%E3%82%B9%E3%83%88%E3%83%A9%E6%B3%95%EF%BC%89%E7%AB%B6%E6%8A%80%E3%83%97

おまけ

記事中のグラフはnetworkxで作図しました。 気になる人向けにコードを置いておきます。

クリックしたら見れるよ

#ダイクストラ法のグラフ
import matplotlib.pyplot as plt
import networkx as nx

#有向グラフのセット
G = nx.DiGraph()

#辺を追加していく
# G.add_edges_from([("A", "B", {"weight":2)])
# nx.add_path(G, ["A", "C"])
# G.edges["A", "B"]["weight"] = 2
G.add_weighted_edges_from([("A", "B", 5), ("A", "C", 7), ("A", "D", 1), ("B", "C", 3), ("B", "F", 3), ("D", "B", 2)])
G.add_weighted_edges_from([("C", "E", 3), ("D", "F", 8), ("F", "E", 2), ("B", "E", 6)])

#図中で頂点の座標を指定しないとランダムに配置される
pos1 = {}
pos1["A"] = (0, 0)
pos1["B"] = (1, 0)
pos1["C"] = (1, 2)
pos1["D"] = (1, -2)
pos1["E"] = (2, 1)
pos1["F"] = (2, -1)

#そのまま描写するとweightの文字が邪魔なので消す
edge_labels = {(i, j): w['weight'] for i, j, w in G.edges(data=True)}

#重みを書き込む
nx.draw_networkx_edge_labels(G, pos1, edge_labels=edge_labels, font_size=18)
#作図
nx.draw_networkx(G, pos1, node_size=1000, label_pos=100, node_color="r", font_size=18)
#図の保存
plt.savefig("dijkstra.png")
plt.show()

#networkxにもダイクストラ法が実装されている。計算量はO(V^2+E)なので遅い
pred, dist = nx.dijkstra_predecessor_and_distance(G, "A")
print(pred, dist)

>>{'A': [], 'B': ['D'], 'C': ['B'], 'D': ['A'], 'F': ['B'], 'E': ['F']} 
>>{'A': 0, 'D': 1, 'B': 3, 'C': 6, 'F': 6, 'E': 8}

Numpy100本ノックで出会った「こりゃあ、便利だなぁ~」と思った関数一覧

https://fontmeme.com/permalink/200609/09f69cd77b1ab658e128cee41158569a.png

はじめに

Qiitaで全く同じ内容の記事を載せたら100以上のいいねを貰って嬉しかった記事です。

承認欲求が満たされるってこんなに幸せなことなのね…(2週間くらい幸せ期間が続きました)

概要

会社の研修中にNumpy100本ノックに取り組みさせられましたが、そもそも関数を知らないと倒せない問題が数多くあり、知識不足を痛感しました。 その、僕が初めましての関数の中でも「便利だなぁ~」と思ったものを備忘録を兼ねてご紹介するのがこの記事になります。

筆者プロフィール

プログラミング歴:約1年(Pythonのみ) Numpyはこの本で勉強しました。

初めに

初めての関数に出会ったら、ググるよりも先にhelpを見た方が圧倒的に理解が早いです。 helpは

np.info([関数名])
ex) np.info(np.add)

jupyter notebookであれば

[関数名]?
ex) np.add?

で表示されます。 関数の説明のみならず、引数の種類や使用例まで載っており、情報量がとても多いです。 英語なので慣れていないと苦痛に感じるかも知れませんが、是非参考にしていただきたいです。

「便利だなぁ~」な関数

np.flip

配列を反転します。

a = np.arange(10)
np.flip(a)
-> array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

ベクトルであればスライス([::-1])で事足りますが、行列の場合などに重宝しそうです。 行列の場合、axisを指定すると任意の方向にひっくり返すことができます。

a = np.arange(9).reshape(3,3)
np.flip(a)
-> array([[8, 7, 6],
          [5, 4, 3],
          [2, 1, 0]])

np.flip(a, axis=0)
->array([[6, 7, 8],
         [3, 4, 5],
         [0, 1, 2]])

np.flip(a, axis=1)
-> array([[2, 1, 0],
          [5, 4, 3],
          [8, 7, 6]])

np.eye

単位行列を生成します。

np.eye(3)
-> array([[1., 0., 0., 0.],
          [0., 1., 0., 0.],
          [0., 0., 1., 0.],
          [0., 0., 0., 1.]])

np.diag

対角成分を抽出します。

a = np.random.randint(0,10, (3,3))
print(a)
-> [[5 2 8]
    [2 7 5]
    [5 1 0]]
p.diag(a)
-> array([5, 7, 0])

np.tile

配列を敷き詰めます。

a = np.array([[0, 1], [1, 0]])
np.tile(a, (2,2))
-> array([[0, 1, 0, 1],
          [1, 0, 1, 0],
          [0, 1, 0, 1],
          [1, 0, 1, 0]])

np.bincount

配列中の数値(非負、int型)をカウントし、その値のindexに格納します。

a = np.random.randint(0, 10, 10)
print(a)
-> array([8 6 0 8 4 6 2 5 2 1])
np.bincount(a)
-> array([1, 1, 2, 0, 1, 1, 2, 0, 2], dtype=int64)

np.repeat

要素を指定した数だけ繰り返します。

np.repeat(3, 4)
-> array([3, 3, 3, 3])

np.roll

配列を指定した数だけ右にシフトします。

a = np.arange(10)
np.roll(a, 2)
-> array([8, 9, 0, 1, 2, 3, 4, 5, 6, 7])

np.flatten

元の配列を1次元配列に直します。

a = np.arange(9).reshape(3,3)
print(a)
-> [[0 1 2]
    [3 4 5]
    [6 7 8]]
b = a.flatten()
print(b)
-> [0 1 2 3 4 5 6 7 8]

np.nonzero

0以外の値が含まれているindexを教えてくれます。

a = np.array([1,0,0,1])
b = np.array([1,1,1,1])
print(np.nonzero(a))
-> (array([0, 3], dtype=int64),)
print(np.nonzero(b))
-> (array([0, 1, 2, 3], dtype=int64),)

np.copysign

第1引数の符号を第2引数と同じ符号に変換します。

a = np.arange(10)
b = np.repeat([1, -1], 5)
np.copysign(a, b)
-> array([ 0.,  1.,  2.,  3.,  4., -5., -6., -7., -8., -9.])

np.intersect1d

2配列の共通する要素を取り出します。

a = np.arange(10)
b = np.arange(5, 15)
np.intersect1d(a, b)
-> array([5, 6, 7, 8, 9])

最後に

こういった関数は暗記するものではありませんが、そもそも知らないと使うことができません。 ご紹介した関数達が皆様の頭の片隅にでも残っていただけるととても嬉しいです!

正規分布の確率密度関数を理解してみる

あらすじ

f(x)=\frac{1}{\sqrt{2\piσ^2}}e^{-\frac{(x-μ)^2}{2σ^2}}

正規分布は代表的な分布の一つであり、統計の分野で頻出する分布です。 しかし正規分布確率密度関数は↑式のようにとてもややこしく、多くの初学者を葬り去ってきました。 かく言う僕も統計の勉強を始めて、初めてこの式を見た時は頭痛と眩暈を覚えた記憶があります。

本記事では、↑式の係数部分($\frac{1}{\sqrt{2\piσ2}}$)と指数部分($-\frac{(x-μ)2}{2σ2}$)に分けて、なぜこのような形になったのかを記述していきます。

指数部分

世の中の多くの事象は平均値を取る確率が最も大きく、平均値から離れるにつれその値を取る確率は小さくなります。

これを簡単に表す式が、

f(x)=e^{-x^2}

です。 グラフで表すと、

import numpy as np
import matplotlib.pyplot as plt

def normal_dist(x, ave = 0, disp=1):
    return np.exp((-(x-ave)**2)/(2*disp**2))

x = np.linspace(-3, 3)
y = normal_dist(x)

plt.plot(x,y)
plt.grid(axis="both")
plt.show()

image.png こんな感じで山なりの形になります。

この式では一通りのグラフしか書けないので汎用性に劣ります。 そこで、 ①グラフの左右移動 ②グラフ幅の変更 の機能を追加します。

グラフの左右移動

$x$の部分を$x-μ$に変更することで、極大値を取る$x$の位置をずらすことが出来ます。

f(x)=e^{-(x-μ)^2}

こんな式になります。 $μ$の値を変えていき、グラフの変化を見てみましょう。

color = ["b", "g", "r", "c", "m"]
for i, col in enumerate(color):
    y = normal_dist(x, i)
    plt.plot(x, y, color=col)
    
plt.show()

image.png $μ$が大きくなるにつれ、グラフが右にずれていきました。

グラフ幅の変更

指数部分に$\frac{1}{2σ2}$を掛けると幅の変更が可能になります。

f(x)=e^{-\frac{x^2}{2σ^2}}
x = np.linspace(-10, 10)
color = ["b", "g", "r", "c", "m"]
for i, col in enumerate(color):
    y = normal_dist(x, 0, i+1)
    plt.plot(x, y, color=col)
    
plt.show()

image.png

幅を変えることに成功しました。 σを二乗にしたのは正負どちらも扱えるようにするため、$2$を掛けているのは計算しやすくするためです。

係数部分

確率密度関数なので面積の合計が$1$でなければいけません。 そこで都合の良い値を関数に掛けてやります。 都合の良い係数を$c$と置き、

\int_{-\infty}^{\infty}ce^{-\frac{(x-μ)^2}{2σ^2}}dx=1

を計算し、$c$を求めます。

一見やっかいな計算に見えますが、

\int_{-\infty}^{\infty}e^{-ax^2}dx=\sqrt{\frac{\pi}{a}}

このガウス積分の公式を用いて$a=\frac{1}{2σ2}$とすれば

c=\frac{1}{\sqrt{2\piσ^2}}

となり、簡単に計算が可能です。 解は一番上の式の係数になりました。

結局この式は、$f(x)=e^{-x2}$に汎用性を持たせ、係数を調整することで確率密度関数を定義するという式なのでした。 めでたしめでたし。

参考にしたサイト

https://to-kei.net/distribution/normal-distribution/density-function-derivation/ 正規分布の密度関数を意味的に理解する

https://mathtrain.jp/gauss ガウス積分の公式の2通りの証明

https://bellcurve.jp/statistics/course/7797.html 14-1. 正規分布