Juliaで散布図行列を表示する

Julia v0.6.2 + StatPlots.jlで散布図行列を表示する。

using CSV
d = CSV.read("data-attendance-1.txt")

using StatPlots
@df d corrplot([:A :Score :Y])

png("corrplot")

以下の画像が得られる。

f:id:YuK_Ota:20180510003420p:plain

相関を楕円で出す奴が好きとか色々あると思うので、結局自分で作ることになるんですかね。

すべてのものを直立にする

Juliaでヒストグラムを描く。

環境はJulia v0.6.2
入力データはCSVとする。

using CSV
d = CSV.read("data-attendance-1.txt")
# Scoreというラベルのデータをプロットする。
score = convert(Array{Int64}, d[:Score])

using StatsBase
# nbinsで棒の数、第三引数でstart:step:maxでエッジの指定ができる。
h = fit(Histgram, c)

using Plots
gr()
bar(h.edges[1][2:end], h.weights)

概形みたいだけとかならExcelとかのほうが良さそう。

StanとJuliaでベイズ統計モデリング

確率的プログラミング言語を学ぶ機運を感じたので勉強している。
以下の本を元に勉強しているが、写経するのも面白くないのでJuliaでやっていく。
環境はUbuntu16.04、Julia0.6.2。

インストール

Stan.jlはCmdStanを呼び出すため、CmdStanとStan.jlの両方のインストールが必要である。
Stan.jlのインストール時にいい感じにやってくれんかなと思ったけどそんなことはないし、やられたらそれはそれでうざい。

CmdStanのインストール

aptとかにもないので公式ページから落としてくる必要がある。
Stan - CmdStan

私は~/bin以下に解凍した。
次に.zshrcに以下を追記。

export CMDSTAN_HOME=/home/youraccount/bin/cmdstan-2.17.1

いつも通りsource .zshrcして完了。

Stan.jlのインストール

CmdStanが設定されていればPkg.addするだけである。

Pkg.add("Stan.jl")

Juliaのlm関数で推定

上記書籍でも一旦、Rのlmを使用して回帰させているので、こちらも一旦それに習う。
GLM.jlを使用し同じことを行う。

using CSV
using GLM
using DataFrames

data = CSV.read("data-salary.txt")
profit = lm(@formula(Y ~ X), data)
new_x = DataFrame(X=23:60)
                                                                                      
@show predict(profit, new_x) 

ちなみに予測区間を出すにはpredict(profit, new_x, :confint)とやれば良さそうにREADMEには書いてあるが、
DataFrameからMatrixの変換はStatsModelsで実装されており、StatsModels側の実装がSymbolを受け取るようなI/F
ではないため、できない。READMEにはかいてあるのに。
issueをみると1年くらい前に治そうとしている形跡はある。
私も直してやろうかと思ったらそもそもテストが落ちているので、まずそこを直そうとなり、これがそもそもこの仕様やめね?という話に飛び火してつらみを出してる。
巷でJulia最高と言っている人間はこういう辛い面を全く言わないので多分違うJuliaを使っている気がする。

Stanで実装、Juliaからの実行方法

上記書籍のmodule4-5.stanを写経する。
下からダウンロードしても良い。
RStanBook/model4-5.stan at master · MatsuuraKentaro/RStanBook · GitHub

後はStan.jlからモデルを呼び出す。

using CSV
using Stan
d = CSV.read("data-salary.txt")
data = [Dict("N"=>size(d, 1), "X"=>d[:X], "Y"=>d[:Y])]
f = open("model4-5.stan", "r")
model = readstring(f)
close(f)
stanmodel = Stanmodel(model = model, random=Stan.Random(1234))
rc, sim = stan(stanmodel, data)

実行結果であるsimはMambaというMCMCのライブラリで扱われるデータ型になる。

実行結果を見る。

using Mamba
describe(sim)

最後にグラフに表示する。
MambaはGadflyでグラフ表示することを意識して作っているが、
Plots(GRBackend)が好き(というかそれ以外であまり上手く行ったことがない)なため、これで描画する。
今回もGadflyで描画しようとしてクラッシュしたので逃げたとも言う。

Traceを表示する。
といっても表示要素が多いためどれを表示したいかを選ぶ。

sim.names

表示させたいものが決まったらグラフに出す。

using Plots
gr()

nrows, nvars, nchains = size(sim.value)
i = 9 # bを表示
Plots.plot(1:nrows, sim.value[:, i, :], title=sim.name[i])

bが21付近をうろうろしているグラフが得られ、1年ごとに21万も昇給するわけねぇだろという気持ちになったので終わります。

Mastodonインスタンス立ち上げ後にやること

なんとなくDockerでMastodonインスタンスを作れるが、まだやるべきことがあります。
デイリー、ウイークリー、マンスリーでやるべきタスクがあるのでこれをcronに登録する必要があります。
ユーザーがいる場合、サービスが止まらないか気にする必要があるので予めやっておくほうが良いです。

mastodon:dailyのcron設定

mastodon:dailyというコマンドがここでやれと言われているコマンドをまとめてやってくれています。
documentation/Maintenance-Tasks.md at master · tootsuite/documentation · GitHub
dailyで実行するものなのでcronに登録しましょう。

@daily cd /your/gitcloned/mastodon && /usr/local/bin/docker-compose run --rm web rake mastodon:daily | logger -t mastodon -p local0.info

mastodon:media:remove_remoteのcron設定

1週間以上前の添付ファイルのキャッシュを消します。

@weekly cd /your/gitcloned/mastodon && /usr/local/bin/docker-compose run --rm web rake mastodon:media:remove_remote 2>&1 | logger -t mastodon -p local0.info

docker-composeへのパスは絶対パスで指定しておくか、cron実行時のPATHを確認したうえで指定してください。
cron実行時とsshのときのPATHは異なる場合が有ります。

Let's Encryptの証明書の更新のcron設定

Let's Encryptの証明書の期限は90日です。
公式で60日での更新を奨めています。
FAQ - Let's Encrypt - Free SSL/TLS Certificates
monthlyで更新しましょう。
以下は最初にwebrootでやった場合です。standaloneの場合は一度nginxを止める必要があると思います。

@monthly /your/gitcloned/letsencryption/certbot-auto renew && /bin/systemctl reload nginx

一度手元でdry-runで実際に動作するのか試しておいたほうが良いです。

certbot-auto renew --dry-run

私の環境だとはnginxのconfが雑だったため動きませんでした。
大体下のリンクからパクったconfですが、rootの場所を設定していませんでした。
documentation/Production-guide.md at master · tootsuite/documentation · GitHub
nginxのインストール->証明書取得->mastodonインストールの順で作業してるとなんとなくうまく行ってしまうポイントなのでアレ。


ちなみに立ち上げたインスタンスです。
mastodoll.net

JuliaImagesで画像の平行移動をする

Julia v0.5.1 + JuliaImages/ImageTransformations.jl v0.2.2で画像の平行移動を行う。

using Images 
using CoordinateTransformations

img = Gray{Float64}[                                                                                                                                                          
 0.0 0.0 0.0 0.0 0.0;
 0.0 0.5 0.5 0.5 0.0;
 0.0 0.5 1.0 0.5 0.0; 
 0.0 0.5 0.5 0.5 0.0;
 0.0 0.0 0.0 0.0 0.0;]

tfm = recenter(Translation(2,0), center(img_pyramid))
translated_img = warp(img, tfm)

このとき、以下のように表示され、なんか移動してなくね?となる。

 3:7×1:5 OffsetArray{Gray{Float64},2}:
 Gray{Float64}(0.0)  Gray{Float64}(0.0)  Gray{Float64}(0.0)  Gray{Float64}(0.0)  Gray{Float64}(0.0)
 Gray{Float64}(0.0)  Gray{Float64}(0.5)  Gray{Float64}(0.5)  Gray{Float64}(0.5)  Gray{Float64}(0.0)
 Gray{Float64}(0.0)  Gray{Float64}(0.5)  Gray{Float64}(1.0)  Gray{Float64}(0.5)  Gray{Float64}(0.0)
 Gray{Float64}(0.0)  Gray{Float64}(0.5)  Gray{Float64}(0.5)  Gray{Float64}(0.5)  Gray{Float64}(0.0)
 Gray{Float64}(0.0)  Gray{Float64}(0.0)  Gray{Float64}(0.0)  Gray{Float64}(0.0)  Gray{Float64}(0.0)

だが、この行列はOffsetArrayで表現されており、これでちゃんと平行移動している。
1行目の型情報にオフセットの情報が書いてあり、この値域の情報だけが下に表示されている。
ちなみにtranslated_img[1, 1]でアクセスするとBoundsErrorになるのでindicesなりで座標の値域をとってからアクセスする必要がある。


私がOffsetArrayに気付かず、テストコードやサンプルも回転移動のものしかないため少し詰まったのでメモ。
関係ないが、スーパーpre記法がjuliaに対応していた。
はてなブログ最高。

Ubuntu16.04でKinnect v1を動かす

本日ThinkPadが届き、これまで放置してしまっていたKinnectでViolin演奏中の姿勢の確認する奴作るかということでubuntu16.04でKinnect v1の動作を確認した。
とく難しいことをせず動いてしまったので、動いたよという報告だけになる。

参考にしたのは以下のスライドである。参考というかこのままである。
https://staff.aist.go.jp/kanezaki.asako/pdf/SSII2016_AsakoKanezaki_tutorial.pdf

ROSインストール

以下のマニュアルに従ってインストールする。*1
kinetic/Installation/Ubuntu - ROS Wiki

freenectインストール

Kinectを扱うためのコンポーネントをインストールする。。

sudo apt-get install ros-kinetic-freenect-launch

これをインストールすることでfreenect(ドライバ)もついでにインストールされるので予めfreenectを入れないと〜ということに悩まされない。
入るパッケージはros-kinetic-libfreenectのバージョンが0.5.1とubuntuリポジトリ経由で入るlibfreenectの0.5.2より古い。
同じドライバが違うパッケージ名で競合してるの関わりたくない感じがしたので見なかったことにした。

動作

rvizでなんか写ってることを確認する。
点群だけ表示してもは???って感じなのでrvizでみたほうが良い。
ありがちなのだが、PointCloud2とか表示に加えないと左のメニューにも表示されない。
この手の、最初設定したらもう知らん系にちょいちょい悩む。

roslaunch freenect_launch freenect.launch
rosrun rviz rviz

おわりに

前日OpenNI入れろとかNiTEを闇から入手しろみたいなブログばかり見まくって死ぬか思った。
インターネット怖い。
手順的には産総研のスライドと同じなのでこんな記事いらないのだが、16.04でもうごいたで!ありがとや!という気持ちを伝えたかっただけである。

*1:そっからかよと思うかもしれないが、僕のThinkPadは今日届いたのだ。

Golearnとplotで散布図を書く

Go入門中です。
Goleanで読み込んだデータをgonum/plotでプロットします。

package main

import (
	"fmt"
	"image/color"
	"math/rand"

	"github.com/sjwhitworth/golearn/base"

	"github.com/gonum/plot"
	"github.com/gonum/plot/plotter"
	"github.com/gonum/plot/vg"
)

func main() {
	rawData, err := base.ParseCSVToInstances("datasets/iris.csv", false)
	if err != nil {
		panic(err)
	}

	_, rows := rawData.Size()

	spal_length_attrs := rawData.AllAttributes()[0].(*base.FloatAttribute)
	spal_length_spec, _ := rawData.GetAttribute(spal_length_attrs)

	spal_width_attrs := rawData.AllAttributes()[1].(*base.FloatAttribute)
	spal_width_spec, _ := rawData.GetAttribute(spal_width_attrs)

	species_attrs := rawData.AllAttributes()[4].(*base.CategoricalAttribute)
	category_values := species_attrs.GetValues()
	category_num := len(category_values)

	fmt.Println(category_values)
	fmt.Println(category_num)

	class_attrs := rawData.AllClassAttributes()[0].(*base.CategoricalAttribute)
	class_spec, _ := rawData.GetAttribute(class_attrs)

	category_vs_raw_map := make(map[string]plotter.XYs)
	fmt.Println(category_vs_raw_map)

	// plot
	for i := 0; i < rows; i++ {
		// category
		category_bytes := rawData.Get(class_spec, i)
		category_name := class_attrs.GetStringFromSysVal(category_bytes)
		fmt.Println(category_name)
		_, isExist := category_vs_raw_map[category_name]
		if !isExist {
			fmt.Println("Add ")
			fmt.Println(category_name)
			category_vs_raw_map[category_name] = make(plotter.XYs, 0)
		}

		spal_length_bytes := rawData.Get(spal_length_spec, i)
		spal_length := spal_length_attrs.GetFloatFromSysVal(spal_length_bytes)

		spal_width_bytes := rawData.Get(spal_width_spec, i)
		spal_width := spal_width_attrs.GetFloatFromSysVal(spal_width_bytes)

		append_xy := plotter.XYs{{spal_length, spal_width}}
		pts := category_vs_raw_map[category_name]
		category_vs_raw_map[category_name] = append(pts, append_xy[0])

	}

	p, err := plot.New()
	if err != nil {
		panic(err)
	}

	p.Title.Text = "Iris Spal Length vs Width"
	p.X.Label.Text = "Length"
	p.Y.Label.Text = "Width"
	p.Add(plotter.NewGrid())

	for name, xys := range category_vs_raw_map {
		fmt.Println(name)
		fmt.Println(xys)
		plot_data, _ := plotter.NewScatter(xys)
		r := uint8(rand.Float32() * 255)
		g := uint8(rand.Float32() * 255)
		b := uint8(rand.Float32() * 255)
		plot_data.GlyphStyle.Color = color.RGBA{R: r, G: g, B: b, A: 255}
		p.Add(plot_data)
		p.Legend.Add(name, plot_data)
	}

	// save pict
	p.Save(4*vg.Inch, 4*vg.Inch, "iris_spal.png")
}


これで次のグラフがプロットされます。
f:id:YuK_Ota:20160427232438p:plain


Go、チュートリアル中なので微妙なところがあれば教えて欲しいです。