<style> .reveal, .reveal h1, .reveal h2, .reveal h3, .reveal h4, .reveal h5, .reveal h6 { font-family: "Source Sans Pro", "Helvetica Neue", Helvetica, Arial, "Microsoft JhengHei", Meiryo, sans-serif; } h1, h2, h3, h4, h5, h6 { text-transform: none !important; } .color-yellow{ color: yellow; } .alert { padding: 15px; margin-bottom: 20px; border: 1px solid transparent; border-radius: 4px; text-align: left; padding: 10px 0; } .alert-info { color: #31708f; background-color: #d9edf7; border-color: #bce8f1; } .alert-success { color: #3c763d; background-color: #dff0d8; border-color: #d6e9c6; } .alert-danger { color: #a94442; background-color: #f2dede; border-color: #ebccd1; } .reveal .slides span { text-align: left; display: inline-block; } p, li { font-size: 0.88em !important; } li>p { font-size: 1em !important; } </style> # JuliaTokai #8 (2020/11/14) ###### tags: `JuliaTokai` `prezentation` --- # MCMCで求める<br>Holy Numbers [堀川 由人, ほりたみゅ, @Hyrodium](https://hyrodium.github.io/profile) ---- ### おしながき * 自己紹介 * Holy Numbers * MCMC * モデリング * 計算結果 * まとめ等 --- ## 自己紹介 私の学生時代、知りたいですよね..？ * 高専時代 * 大学以降 ---- ### リンク機構と私 (高専時代) * 高専ロボコン(1年-3年) * リンク機構大好き少年      ---- ### リンク機構と私 (大学以降) * まだまだリンク機構は面白い * 3年前のストローで作った作品    最近の気持ち: **Juliaによる数値計算がアツい！！** --- ## Holy Numbers 本LTの題材の説明です * Strand Beest * Holy Numbers * やりたいこと！ ---- ### Strand Beest * オランダの芸術家Theo Jansen氏による[作品](https://www.youtube.com/watch?v=MYGJ9jrbpvg) * 「風を食べて動く生命体」  ---- ### Holy Numbers (1) * 「聖なる13の数」 * Strand Beestのリンク比のこと * Tシャツも販売されていた    ---- ### Holy Numbers (2) * リンク比マジ重要 * 少しでも変わると変な軌跡になる →滑らかに歩けない →**疑似直線**運動が必要     ---- ### やりたいこと！ * このリンク比(Holy Numbers)を自分で求めたい！ * あわよくば..もっと良いリンク比を見つけたい！   --- ## MCMC 最適化のための道具の紹介！ * MCMCってなに？ * MH法 * 焼き鈍し法 ---- ## MCMCってなに？ (1) * **M**arkov **C**hain **M**onte **C**arlo * Markov Chain (マルコフ連鎖) * 前の状態を記憶する * Monte Carlo (モンテカルロ) * 確率的に実行する * サンプリング手法の一つ * ベイズ推論の応用でもよく使われる * (らしいですが、まだ詳しくないです..) *(最近勉強し始めた分野なのでマサカリ歓迎です！！)* ---- ## MCMCってなに？ (2) * 与えられた確率分布$p(x)$を元にサンプリング * ヒストグラムが$p(x)$に一致する (すごい！)  ---- ## MH法 * Metropolis-Hastings法 (MCMC法の一つ) 1. 初期点$x_{0}$を決める 2. 点$x_{i}$をもとにランダムに$x'$を決定 3. 以下の漸化式で確率的に$x_{i+1}$を決定 ここで$\operatorname{rand}()$は$[0,1]$区間の一様乱数 $x_{i+1} = \begin{cases} x' & \left(\operatorname{rand}()<\frac{p(x')}{p(x)}\right) \\ x_{i} & (\text{otherwise}) \end{cases}$ * これだけで確率分布$p(x)$に従うサンプリングが可能 ---- ## 焼き鈍し法 (1) * MCMCって最適化問題にも使えるんですか? * 使えます！(焼き鈍し法) * 最適化問題 * エネルギー関数$x\mapsto E(x)$が最小になる$x$を求める * ハイパーパラメータ$\beta>0$を用意して以下で変換 * 分母は定数なので分母の計算は不要 $\displaystyle p(x) = \frac{\exp(-\beta E(x))}{\int_{x' \in D}\exp(-\beta E(x')) dx'}$ ---- ## 焼き鈍し法 (2) * ハイパーパラメータ$\beta>0$って？ * 逆温度と呼ばれる量 (統計力学からの輸入) * $\beta$を徐々に大きくすると.. * エネルギーの小さい所に存在する確率が上がる！ <iframe src="https://www.desmos.com/calculator/cvonetwtmk" width="600" height="400" frameborder="0" style="border:0" allowfullscreen></iframe> --- ## モデリング 「MCMCで求めるHoly Numbers」をどう実現するか？ * 未知変数の設定 * エネルギー関数の設計 * 最適化手法の検討 ---- ### 未知変数の設定 * リンク機構を決定する変数は13個 * $(w_{1},\dots,w_{13})$とおきます * 全体の定数倍は同じ軌跡(相似)になる * 原動節の長さ$w_{13}$を固定すれば残るは12個 * つまり$(w_{1},\dots,w_{12})$が未知変数   ---- ### エネルギー関数の設計 * 悪い軌跡には高いエネルギー:arrow_up: * 良い軌跡には低いエネルギー:arrow_down:     10　　　　　60　　　　　150　　　　　80 エネルギー関数の最小化問題に帰着される $E : \mathbb{R}^{12}\to \mathbb{R};(w_{1},\dots,w_{12})\mapsto E(w)$ ---- ### エネルギー関数の例 * 足先の位置ベクトル: $\displaystyle\newcommand\Pare[1]{\left(#1\right)} \newcommand\od[2]{\frac{d#1}{d#2}} \newcommand\pd[2]{\frac{\partial#1}{\partial#2}} \newcommand\sod[2]{\sfrac{d#1}{d#2}} \newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}} \newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2} \newcommand\bm[1]{\boldsymbol{#1}}\bm{c}(\theta; w)=(c_x(\theta;w), c_y(\theta;w))$ * リンク比: $w=(w_{1},\dots,w_{12})$ * 原動節の角度: $\theta\in[0,2\pi)=[0,\tau)$ * 位置エネルギー (重力の一般化) $\displaystyle\newcommand\Pare[1]{\left(#1\right)} \newcommand\od[2]{\frac{d#1}{d#2}} \newcommand\pd[2]{\frac{\partial#1}{\partial#2}} \newcommand\sod[2]{\sfrac{d#1}{d#2}} \newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}} \newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2} \newcommand\bm[1]{\boldsymbol{#1}} E_{\text{potential}}(w)=\int_{0}^{\tau} \Pare{c_y(\theta; w)-\min_{\theta'\in [0,\tau)}c_y(\theta';w)}^r d\theta$ * 面積一定 (風船) $\displaystyle E_{\text{area}}(w)=\Pare{\frac{1}{2}\int_{0}^{\tau} \Pare{\bm{c}(\theta; w)\times\od{\bm{c}(\theta; w)}{\theta}} d\theta-A_0}^2$ ---- ### 最適化の出番！ * 次元の呪い * 10分割ずつ総当りだと$10^{12}$パターン * 他の方法ダメなの？ (e.g. 最急降下法, Newton法) * 微分可能性: 無いかも (あるいは面倒) * 大域的性質: 局所最適なんて要らねえ(ドン！) * そこでMCMCですよ！ $\displaystyle p(w) = \frac{\exp(-\beta E(w))}{\int_{w' \in \mathbb{R}^{12}}\exp(-\beta E(w')) d\mu(w')}$ --- ## 計算結果 色々な面白い結果が揃いました！ * 使用言語 * デモ * 計算例 * 目的関数の設計について ---- ### 使用言語は..もちろんJulia！ * 速い！ * LLVMベースでJITコンパイル * 基本的にCやRustなどと同等の速度 * 書きやすい！ * 動的実行 * 多重ディスパッチ /Julia_(programming_language)-Logo.wine.png =400x) ---- ### デモ ここで計算します:computer: * 10,000回の反復計算が1秒くらい * グラフィックスは[Luxor.jl](https://juliagraphics.github.io/Luxor.jl/stable/)で出力 * (最適化の計算より画像出力の方が遅い..) ---- ### 計算例 (1/5: 位置エネルギー最小化) とりあえず位置エネルギーを小さくしてみるか:thinking_face: $\displaystyle\newcommand\Pare[1]{\left(#1\right)} \newcommand\od[2]{\frac{d#1}{d#2}} \newcommand\pd[2]{\frac{\partial#1}{\partial#2}} \newcommand\sod[2]{\sfrac{d#1}{d#2}} \newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}} \newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2} \newcommand\bm[1]{\boldsymbol{#1}} \begin{aligned}E_{\text{potential}}(w)&=\int_{0}^{\tau} \Pare{c_y(\theta; w)-\min_{\theta'\in [0,\tau)}c_y(\theta';w)}^r d\theta \\ E(w)&=E_{\text{potential}}(w)\end{aligned}$ ---- ### 計算例 (1/5: 位置エネルギー最小化)  潰れた:hushed:！？ ---- ### 計算例 (2/5: 面積一定) 軌跡は閉曲線だから面積が計算できるのでは:smirk: $\displaystyle\newcommand\Pare[1]{\left(#1\right)} \newcommand\od[2]{\frac{d#1}{d#2}} \newcommand\pd[2]{\frac{\partial#1}{\partial#2}} \newcommand\sod[2]{\sfrac{d#1}{d#2}} \newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}} \newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2} \newcommand\bm[1]{\boldsymbol{#1}} \begin{aligned}E_{\text{area}}(w)&=\Pare{\frac{1}{2}\int_{0}^{\tau} \Pare{\bm{c}(\theta; w)\times\od{\bm{c}(\theta; w)}{\theta}} d\theta-A_0}^2 \\ E(w)&=E_{\text{area}}(w) \end{aligned}$ ---- ### 計算例 (2/5: 面積一定)  潰れないが、これでは歩行できん:worried: ---- ### 計算例 (3/5: 組み合わせ) 組み合わせたらいけそう:bulb: $\displaystyle\newcommand\Pare[1]{\left(#1\right)} \newcommand\od[2]{\frac{d#1}{d#2}} \newcommand\pd[2]{\frac{\partial#1}{\partial#2}} \newcommand\sod[2]{\sfrac{d#1}{d#2}} \newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}} \newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2} \newcommand\bm[1]{\boldsymbol{#1}} \begin{aligned}E_{\text{area}}(w)&=\Pare{\frac{1}{2}\int_{0}^{\tau} \Pare{\bm{c}(\theta; w)\times\od{\bm{c}(\theta; w)}{\theta}} d\theta-A_0}^2 \\ E(w)&=E_{\text{potential}}(w)+E_{\text{area}}(w)\end{aligned}$ ---- ### 計算例 (3/5: 組み合わせ)  結構ええんちゃう！？膝がめり込んでるが..:anguished: ---- ### 計算例 (4/5: 面積の調整) 軌跡の面積が大きすぎたか..？:face_with_rolling_eyes: $\displaystyle\newcommand\Pare[1]{\left(#1\right)} \newcommand\od[2]{\frac{d#1}{d#2}} \newcommand\pd[2]{\frac{\partial#1}{\partial#2}} \newcommand\sod[2]{\sfrac{d#1}{d#2}} \newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}} \newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2} \newcommand\bm[1]{\boldsymbol{#1}} \begin{aligned}E_{\text{area}}(w)&=\Pare{\frac{1}{2}\int_{0}^{\tau} \Pare{\bm{c}(\theta; w)\times\od{\bm{c}(\theta; w)}{\theta}} d\theta-A_0}^2 \\ E(w)&=E_{\text{potential}}(w)+E_{\text{area}}(w)\end{aligned}$ ↓↓↓ $A_0$を半分くらいにしてみよう！ ---- ### 計算例 (4/5: 面積の調整)  かなりええんちゃう！？ちょっと細長いけれども..:grin: ---- ### 計算例 (5/5: 負の面積) 軌跡の面積って負の数にしてもいけるんやろか..:exclamation::question: $\displaystyle\newcommand\Pare[1]{\left(#1\right)} \newcommand\od[2]{\frac{d#1}{d#2}} \newcommand\pd[2]{\frac{\partial#1}{\partial#2}} \newcommand\sod[2]{\sfrac{d#1}{d#2}} \newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}} \newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2} \newcommand\bm[1]{\boldsymbol{#1}} \begin{aligned}E_{\text{area}}(w)&=\Pare{\frac{1}{2}\int_{0}^{\tau} \Pare{\bm{c}(\theta; w)\times\od{\bm{c}(\theta; w)}{\theta}} d\theta-A_0}^2 \\ E(w)&=E_{\text{potential}}(w)+E_{\text{area}}(w)\end{aligned}$ ↓↓↓ $A_0$の代わりに$-A_0$にしてみよう！ ---- ### 計算例 (5/5: 負の面積)   逆回転の軌跡ができた！:star-struck: ---- ### 計算結果の一旦まとめ * 疑似直線は意外と簡単に作れる * 目的関数(=エネルギー関数)次第で形状が大きく変化 * そもそも良い軌跡ってどんな形？ * 同値な問い: 最適な目的関数とは..？    ---- <!-- .slide: data-background="https://i.imgur.com/yEDcmq3.gif" --> ### 目的関数の設計 * 目的関数(=エネルギー関数)への要求 * 接地部の**軌跡は平坦**な方が良い:walking: * スムーズに歩くため * 接地部の**速度は一定**の方が良い:runner: * もっとスムーズに歩くため * **軌跡は高い**方が良い:arrow_down: * 重心が低い方が安定 * 軌跡の**上下幅は小さい**方が良い:battery: * 足を上げるのは無駄なエネルギー消費 * 軌跡の**上下幅は大きい**方が良い:mountain: * 不整地走破性のため * 各**リンクは短い**方が良い:balloon: * 軽量化したい * *全部を考慮するのは大変！！* <small>..(でも"それっぽい"ものは作れた)</small> --- ## まとめ等 長々とお付き合いありがとうございました。 <small>(10分に収まっていたでしょうか..？)</small> * Theo Jansen氏の計算との比較 * 色々実装した感想 * 参考文献など ---- ### Theo Jansen氏の計算との比較 * オリジナルのHoly Numbersは得られなかった * 目的関数が分からん.. * 新しい結果 * 逆向きの回転方向 * 軌跡の高さを(多少は)調整可能 * 軌跡の面積を(多少は)調整可能     ---- ### Theo Jansenマジ天才！！ * 1991年にはリンク比の計算を終えていた * [当時の計算環境](https://youtu.be/FFS-2axFo1Y?t=85)・インターネット環境..:floppy_disk: * GA(遺伝的アルゴリズム, ≠MCMC)を使ったらしい * リンク機構だけじゃない * 風力で歩く * 水を感知して制御 * 圧縮空気も動力源に  ---- ### 色々実装した感想 * リンク機構は楽しい！ * 単純な回転運動から複雑な軌跡を生成 * MCMCは面白い！ * しかし最適化での目的関数の設計は泥臭い.. * Juliaは最高！ * 超速い！書きやすい！気持ちいい！   /Julia_(programming_language)-Logo.wine.png =x220) ---- ### 参考文献など * Theo Jansen Japan * 公式サイト * https://theojansen.net/ * 計算統計Ⅱ * MCMC法とその周辺 * https://www.amazon.co.jp/dp/400730789X * 平面上の2円の交点の座標について * リンク機構の計算に使いました * https://hyrodium.github.io/pdf * HolyNumbersByMCMC * 今回の計算のリポジトリ * https://github.com/hyrodium/HolyNumbersByMCMC