---
robots: noindex, nofollow
lang: ja
dir: ltr
breaks: true
slideOptions:
theme: white
transition: slide
---
<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