理科系の備忘録

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

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:家にいながら会社にいるような気持ちにもなる。

続きを読む

jediによるpythonコード補完 with Emacs on Windows

f:id:kenbell1988:20170122114754g:plain

背景

年始なのでEmacsの設定を見直すことにした。大学時代に作った環境をずっと使っていたが、あれからpackage.elによるパッケージ管理が主流となり、自分の環境が時代遅れのものとなっていた。(auto-completeよりはcompany-modeがいいとか)

pythonのコーディング環境はipythonで満足していたが、Emacsでも関数の補完やヘルプの表示がしたくなったので、パッケージをインストールした。インストールするパッケージはjedi。

jediとは

github.com

pythonの静的解析ライブラリで、コードの補完などに利用できる。色々なエディタで使うことができ、当然Emacsでも利用できる。

続きを読む

線分と平面の交点を求める

線分の両端点(a, b)と平面の法線ベクトル(nv)と平面上の任意の点(p)から、線分と平面の交点を求めるC++のプログラム。

線分abと平面の交点が線分abを内分する点となることから交点の座標を計算する。

#include <iostream>
#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <string>

using namespace std;

cv::Mat Intersect3D(
		 cv::Mat& nv, // 平面の法線ベクトル cv::Mat(3,1)
		 cv::Mat& p,  // 平面上の任意の1点  cv::Mat(3,1)
		 cv::Mat& a,  // 線分の端点         cv::Mat(3,1)
		 cv::Mat& b   // 線分の端点         cv::Mat(3,1)
		 ){
  
  // 返り値
  cv::Mat ret = (cv::Mat_<double>(3,1) << 0,0,0);

  // ベクトルの準備
  cv::Mat pa, pb; // p->a, p->bのベクトル
  pa = a - p;
  pb = b - p;

  // 内積計算
  double dot_a = pa.dot(nv); // paベクトルと法線ベクトルの内積
  double dot_b = pb.dot(nv); // pbベクトルと法線ベクトルの内積
 
  // 内積が0の場合交点がなく、直線が平面に含まれる
  double MIN_DOT_TH = 0.000001;
  if( fabs(dot_a) < MIN_DOT_TH ){dot_a = 0;}
  if( fabs(dot_b) < MIN_DOT_TH ){dot_b = 0;}
  
  // 交点なし
  if( dot_a==0 && dot_b==0){
    return ret; // [0,0,0]'で返す
  }

  // 交点なし
  if( dot_a*dot_b > 0 ){
    return ret; // [0,0,0]'で返す
  }

  // 交点あり
  cv::Mat ab = b - a;
  double ratio = fabs(dot_a) / (fabs(dot_a) + fabs(dot_b) );
  ret = a + ab * (ratio);

  return ret; // 計算結果を返す
}



int main(int argc, char *argv[]){

  cv::Mat a  = (cv::Mat_<double>(3,1) << 0,0, 0); // 端点1
  cv::Mat b  = (cv::Mat_<double>(3,1) << 2,2,-2); // 端点2
  cv::Mat p  = (cv::Mat_<double>(3,1) << 0,0,-1); // 平面上の点
  cv::Mat nv = (cv::Mat_<double>(3,1) << 0,0, 1); // 法線ベクトル
  cv::Mat c = Intersect3D(nv, p, a, b); // 交点

  cout << "result: " << c << endl;
  
  return 0;
}

LabelingクラスをOpenCVで利用

学生時代にも一度お世話になったラベリングクラスhttp://imura-lab.org/products/labeling/を久々に使いました。
OpenCVで使える関数を作ったのでメモ。

2値化が大津の手法なので、琵琶湖の画像にしてみました。
素敵な写真はこちらから使わせていただきました。綺麗ですね。
https://www.itoen.co.jp/itoen-motherlake/photocontest.html



続きを読む