理科系の勉強日記

Linux/Ubuntu/Mac/Emacs/Computer vision/Robotics

B spline曲線と曲面

B spline曲線について勉強したかったんですけどずっと放置していたので、ついに実装して遊んでみました。 basic functionとknot vectorが曲線にどう関係してくるのか実装するまでわかっていませんでした。


S(t) = \sum_{i=0}^n{N_{i,p}(t)P_i}

Basic function


N_{i,0}(u) =  \left\{ \begin{array}{ll}
    1 & (u_i \leq u < u_{i+1}) \\
    0 & (otherwise)
  \end{array} \right. \\
N_{i,p}(u) = \frac{u-u_i}{u_{i+p}-u_i}N_{i, p-1}(u) + \frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u)

f:id:kenbell1988:20190304093622p:plainf:id:kenbell1988:20190304093638p:plain
B spline curve

f:id:kenbell1988:20190304093646p:plainf:id:kenbell1988:20190304093700p:plain
Knot vector, basic function

ついでにB spline曲面も実装しました。u方向の制御点m個とv方向の制御点n個をつかって、control netなる m \times nの行列 Pをつかって以下のように表現できます[1]。実装はgithubのjupyter notebook参照。


S(u, v) = \sum_{i=0}^n\sum_{j=0}^m{N_{i,p}(u)N_{j,q}(v)P_{i,j}}

f:id:kenbell1988:20190304093926p:plain
B spline surface

github.com

参考

[1]https://graphics.stanford.edu/courses/cs348a-09-fall/Handouts/surfaces8.pdf

ガウス過程回帰を実装して遊ぶ

ガウス過程回帰を実装して遊んだ

ガウス過程の本が発売されるということで復習として簡単にメモしようとおもったけど、はてなTeX記法にうまく馴染めなかった1ので 数式は載せずに結果だけ簡単に貼り付けて日記としておきます。

結果

ガウス過程回帰、ハイパーパラメタをうまく決めないとゴミみたな結果が得られるようです。 適当な回数イテレーションを回してハイパーパラメタを決めましょう。やり方はいろいろあると思いますが、MCMCで決めるのがいいかなと考えています。 勾配法はうまくいったためしがないので。。。

f:id:kenbell1988:20190211143121p:plainf:id:kenbell1988:20190211143126p:plainf:id:kenbell1988:20190211143205p:plain
ガウス過程回帰の結果.(左)入力、(中)適当なパラメタ、(右)よさそうなパラメタ


  1. インラインで数式がうまく反映されない。中身を変えずにブロックにすると正しく表示される。やる気失うタイプのトラブル。

続きを読む

Jupyter notebookのthemeを変えたら出力セルが見切れる件

最高なjupyter notebookをさらに最高にするためのjupyter notebook themesをインストールしたが、
どうも表示の調子がよくない。端が切れている。
f:id:kenbell1988:20180929172403p:plain
OS:Mac OS High Sierra
ブラウザ:Google Chrome v64
(Safariでも同様)

このままでは辛いので、もとのcssを直接編集して対応することにした。
~/.jupyter/custom/custom.cssにあるだろうと予想をつけて、
適当にoutputで検索しているとそれらしいものを発見した。
山勘でpaddingを増やしたところうまくいってしまった。

div.output_subarea {
 overflow-x: auto;
 padding: 1.8em !important;  /* もともとコレ padding: 0.8em !important;*/
 -webkit-box-flex: 1;
 -moz-box-flex: 1;
 box-flex: 1;
 flex: 1;
}

f:id:kenbell1988:20180929173018p:plain
1.8が妥当かはまったく検討なし。

以上。

Juliaで主成分分析

周りの人と読んでいる本(ベイズ推論による機械学習)の影響でJuliaをさわってみました。
まだあまりむずかしいことができないので、操作に慣れるために主成分分析をしてみました。
(主成分分析を一発で実行する関数もあるかもしれませんが、練習なので教科書通りやってみます)

using PyCall, PyPlot
using Distributions
# データの作成
n = 300
μ  = [1.0, 3.0]
Σ  = [10 3
        3 2]
X = rand(MvNormal( μ, Σ), n);

MvNormal関数で2次元のガウス分布に従うデータを作成しました。 全角の文字も変数にできるのでギリシャ文字を変数にできて楽しい。

# データから平均を計算して、データから引く
mu = [mean(X[1, :])
            mean(X[2, :])]

X_m = copy(X)
X_m[1, :] = X[1, :] - mu[1];
X_m[2, :] = X[2, :] - mu[2];

2x300のデータから2x1の形のベクトルをうまく引き算する方法がわからなかったのでダサい方法で。。

# 共分散行列の計算
Σ = 1.0 / n *  X_m * X_m'

2×2 Array{Float64,2}: 10.5301 3.24329 3.24329 2.1232

# 固有値分解
λ, S = eig(Σ)
println(λ)
println(S)

# 線形変換
y = S' * X_m;

[1.01741, 11.6359] [0.322704 -0.9465; -0.9465 -0.322704]

# 可視化用の第一成分、第二成分
p_ax1 = [mu[1] mu[1] + S[1,1]*4
                mu[2] mu[2] + S[2,1] *4]
p_ax2 = [mu[1] mu[1] + S[1,2]*4
                mu[2] mu[2] + S[2,2]*4 ]

2×2 Array{Float64,2}: 1.15528 -2.63072 3.10447 1.81366

figure("result")
clf()
subplot(121)
title("PCA result")
plot(X[1,:], X[2, :], "+")
plot(p_ax1[1, :], p_ax1[2, :])
plot(p_ax2[1, :], p_ax2[2, :])

axis("equal")

subplot(122)
title("Linear Transformed")
plot(y[1,:], y[2, :], "+")
axis("equal")

show()

プロットもPyPlotを使えば(pythonになじみがあれば)あまり迷うことなくできました。 次は画像を使って見ようかな、データ用意するの面倒だけど。

f:id:kenbell1988:20180512142948p:plain

セルオートマトンによる渋滞シミュレーション

はじめに

年末年始に渋滞学という本を読んだ。
車の渋滞だけではなく、緊急時の避難や蟻の行列、通信についても書かれていて勉強になった。

渋滞学 (新潮選書)

渋滞学 (新潮選書)

渋滞とセルオートマトン

交通をモデル化して解析することを交通流解析と呼ぶ。交通流モデルは連続モデルと離散モデルに分けられ、連続モデルでは最適速度過程がよく用いられる。離散モデルでは、セルオートマトンがよく用いられており、「渋滞学」でもこれについて説明されていた。

続きを読む

シェルスクリプトで可変長データを読み取る

はじめに

こんなデータに出くわした。data.txtとする。

1, hoge, foo, bar, piyo, [ID: 1; a; b; c; d;ID: 2; a; b; c; d;ID: 3; a; b; c; d;ID: 4; a; b; c; d;ID: 5; a; b; c; d;]
2, hoge, foo, bar, piyo, [ID: 1; a; b; c; d;]
3, hoge, foo, bar, piyo, [ID: 1; a; b; c; d;ID: 2; a; b; c; d;ID: 3; a; b; c; d;]
4, hoge, foo, bar, piyo, [ID: 1; a; b; c; d;ID: 2; a; b; c; d;]

[]で囲まれた部分が可変長のデータになっている。可変長部分から'b'だけを取り出して、縦1列に並べたい。

1   b
2   b
3   b
4   b
5   b
1   b
1   b
2   b
3   b
1   b
2   b

今回一回限りのデータ処理だったので、久々にシェルで遊ぶことにした。以下は考えた順番通りのメモ。もっといい方法はあると思うが。

実践

下ごしらえとしてawkで[ ]の部分だけを取り出す。[ ]の中と外でセパレータが違うのでありがたい。
簡単のため、最初の1行だけで処理を考える。

head -1 data.txt | awk -F ',' '{print $6}'
[ID: 1; a; b; c; d;ID: 2; a; b; c; d;ID: 3; a; b; c; d;ID: 4; a; b; c; d;ID: 5; a; b; c; d;]

awkでforを回すことも考えた。しかし今回は先に試した別の方法がうまく機能した。

head -1 awk -F, '{print $6}' | tr "ID:" "\nID:" | awk -F ';' 'NR>1{print $3}'

'ID'を検索して、'ID'が見つかるたびに改行する。すると

 [
ID 1; a; b; c; d;
ID 2; a; b; c; d;
ID 3; a; b; c; d;
ID 4; a; b; c; d;
ID 5; a; b; c; d;]

と複数の行であらわされた!あとは各行について、awkで'b'を取り出すだけだ。1行目は'['が邪魔なのでNR>1で回避しておく。

最後に、元データの各行についてこれを実行すればいいから、head -1 の部分をwhile readに変更して1行ずつ読み取るようにした。

cat data.txt | while read LINE;
do
    echo $LINE | awk -F, '{print $6}' | tr "ID:" "\nID:" | awk -F ';' 'NR>1{print $3}' | cat -n;
done
1   b
2   b
3   b
4   b
5   b
1   b
1   b
2   b
3   b
1   b
2   b

ついでに

一応関数化してみる。'cat -' で標準入力を受けられるものにした。

read_multi_frame(){
    cat - |
	awk -F, '{print $6}' |
	tr "ID:" "\nID:" |
	awk -F ';' 'NR>1{print $3}'
}
cat data.txt | while read LINE;
do
    echo $LINE | read_multi_frame | cat -n;
done

以上。

Git bashでSolarized Color

f:id:kenbell1988:20170115175537j:plain
Solarizedのdarkが大好き。
Solarized - Ethan Schoonover


自分が使うPCのターミナル(とEmacs)をすべてsolarized darkにすることで環境の差を小さくし、
会社にいながら家にいるような気持ちでリラックスしてPCに向き合える。*1

WindowsのPCにはGit bash(mintty)をインストールして、ここでシェルを書いたりコマンドを実行したりしている。
.bashrcに以下2行を書いて、対応するファイルを以下のように作ってホームディレクトリに転がしておけば、それでOK。

*1:家にいながら会社にいるような気持ちにもなる。

続きを読む