KaggleのTitanicデータに対してsklearnの決定木を試してみる
(注)コンペ目的でない素人の備忘録です。参考になるかは不明ですがコメントは歓迎します。
Kaggleとは
KaggleはHPにあるとおり機械学習/データサイエンスに携わる人や企業などのコミュニティです。
正直始めたてなのでざっくりとしかわからないですが,Competitionという、企業や政府などがコンペ形式で課題を提示し、賞金と引き換えに最も制度の高い分析モデルを買い取るといったシステムが有名です。
勉強用としてのKaggle(Titanic)
機械学習は教科書で深層学習やパーセプトロンなどをざっくり読んだだけなのでコンペに参加することはしばらくないでしょうが,勉強用にはいくつか使えそうな点があります。
例えば,練習用としてTitanicというタイタニック号事件のデータが提供されています。 構造もシンプルで初学者にとっても分かりやすいです。
また,人によってはnotebookをわかりやすくコメント付きで投稿してくれていたりします。 www.kaggle.com
他にも,ブラウザ上からアクセスするnotebook形式ではデフォルトでsklearnやlightBGMなどのライブラリがインストールされてる状態で始まるので環境構築の心配が現状なさそうです。
なお,真面目に入門するなら以下の解説が懇切丁寧でよかったです。(英語)
私個人の目的としてはスプラトゥーン2の編成格差問題を自前のデータから解析してみようかなと思っています。(別記事にて記載予定)
メモ:Notebookの作成
初めてKaggleに登録して開くと何をしていいのか本当によくわからなくなりますね。 コンペを探してそれに参加するというのが本来でしょうが,今回はKaggleの環境とデータで練習をするという目的でNotebookを作ります。なお,PythonやPandasなどの知識はどこか別のところで勉強してください。
コンペから探して作成
上のSearchでtitanic
を検索すればtitanicのページに移動します。
右上の「New Notebook」から始めます。
なお,notebookのタブには他の人のnotebookがありさんこうになります。
新規作成→Notebookにあとからデータを追加
sklearnの分類木ライブラリを用いた分類
決定木の説明は教科書等にに投げます。(Qiitaなら以下のリンクなんかが読みやすいのでは?)
[入門]初心者の初心者による初心者のための決定木分析 - Qiita
今回はsklearnというライブラリを用いて決定木のうち,分類木 (classification tree)を用いて学習します。
import pandas as pd import numpy as np import warnings warnings.filterwarnings('ignore') train = pd.read_csv('../input/titanic/train.csv') test = pd.read_csv('../input/titanic/test.csv') sample_submission = pd.read_csv('../input/titanic/gender_submission.csv')
trainは訓練用,testは評価用のデータです。
データの前処理
まずはデータの前処理をします。 機械学習のデータ前処理備忘録 - Qiita
主に,以下のようなことをします。
- 訓練に必要なデータの抜き出し,加工
- 空白データの補間,削除
データ形式の変換と不要なデータの削除
# データの変換 ##maleを0に、femaleを1に変換 train["Sex"] = train["Sex"].map({"male":0,"female":1}) test["Sex"] = test["Sex"].map({"male":0,"female":1}) # EmbarkedのOne-Hotエンコーディング train = pd.get_dummies(train, columns=['Embarked']) test = pd.get_dummies(test, columns=['Embarked']) ## 不要な列の削除 train.drop(['PassengerId', 'Name', 'Cabin', 'Ticket'], axis=1, inplace=True) test.drop(['PassengerId', 'Name', 'Cabin', 'Ticket'], axis=1, inplace=True)
- One-Hotエンコーディング:[A,B,C]のような3つ以上の値のある要素をAの有無としてBoolean値に変換する。
欠損データの確認と対処
ところで,データにはところどころ欠損(NaN)があり,学習に悪影響を及ぼすので取り除いたり穴埋めをしたいです。 pandasの機能でどの量に欠損があるのか簡単に割り出せます。 例えばNanを確認して取り除く場合は以下のように書くことができます。
# NaN の存在確認 と 除去 print(train.isnull().sum()) train2 = train.dropna() test2 = test.dropna()
一応結果を示すと「年齢」の項が抜けているのがわかりますね。
Survived 0 Pclass 0 Sex 0 Age 177 SibSp 0 Parch 0 Fare 0 Embarked_C 0 Embarked_Q 0 Embarked_S 0 dtype: int64
他にも中央値や平均値,ランダム値などで埋める手法があるようです。 参考:平均値でNanを埋めるタイプの方
学習
さて,いよいよ学習のフェーズです。まずは与えられたデータを入出力に分けます。ここでの出力(=予測したい値)は生存したかどうかなのでSurvivedの列を抽出します。
from sklearn.model_selection import train_test_split from sklearn.metrics import accuracy_score # devide train input and output X_train = train2.drop(['Survived'], axis=1) # X_trainはtrainのSurvived列以外 y_train = train2['Survived'] # y_trainはtrainのSurvived列 # X_trainとY_trainをtrainとvalidに分割 train_x, valid_x, train_y, valid_y = train_test_split(X_train, y_train, test_size=0.20, random_state=0)
最後の実行では訓練用のデータを更に分けて4:1に分けて訓練の評価をします。 訓練用データ内でのフィッティング結果と実際の評価用データ内でのフィッティング結果を分けたいがためです。(過学習とかで検索してください。)
訓練用データを2つに分けるのをホールドアウト法,3つ以上に分けるのを交差検証法と言います。(引用)
これで訓練用の入力と出力を分けられたのでsklearnの分類木ライブラリの学習で学習します。
import sklearn.tree as tree # 分類木だからClassifier,Regressorもある clf = tree.DecisionTreeClassifier(max_depth=4) # データを用いて学習 model = clf.fit(train_x, train_y) # データを用いて予測 predicted = model.predict(valid_x) print(accuracy_score(predicted,valid_y))
なんと,これだけです。
結果は0.7902097902097902
になりました。約8割程度の評価です。これは訓練用データに8割程度の正答率で合致しているという解釈になります。
補足:交差検証法による評価
先程の結果はデータのどの部分を訓練用に使うかにより左右されそうです。
したがってよりきちんとやる場合は,交差検証法で訓練データと評価用データをうまくローテーションして平均値などでモデルを評価します。
# 3分割交差検証を指定し、インスタンス化 from sklearn.model_selection import KFold kf = KFold(n_splits=3, shuffle=True, random_state=0) # スコアとモデルを格納するリスト score_list = [] models = [] # 各分割ごとに評価 for fold_, (train_index, valid_index) in enumerate(kf.split(X_train, y_train)): print(f'fold{fold_ + 1} start') train_x = X_train.iloc[train_index] valid_x = X_train.iloc[valid_index] train_y = y_train.iloc[train_index] valid_y = y_train.iloc[valid_index] ## 分割データで学習・予測・評価 clf = tree.DecisionTreeClassifier(max_depth=4) model = clf.fit(train_x, train_y) # データを用いて予測,記録 predicted = model.predict(valid_x) score_list.append(accuracy_score(predicted,valid_y)) models.append(model) print(score_list, '平均score', round(np.mean(score_list), 3))
結果は以下の通り,データによってばらつきがあるということがわかりますね。
[0.8067226890756303, 0.7647058823529411, 0.819327731092437] 平均score 0.797
混合行列による予測結果の視覚化
予測結果の成否をよりわかりやすくするために混合行列を用いた評価ができます。
from sklearn.metrics import confusion_matrix #混同行列の作成 cmatrix = confusion_matrix(valid_y,predicted) #pandasで表の形に df = pd.DataFrame(cmatrix,index=["actual_died","actual_survived"],columns=["pred_died","pred_survived"]) print(df)
pred_died pred_survived actual_died 131 14 actual_survived 29 64
ざっくりというと,死んだという予測は18.1%
で間違っており,生きていたという予測は17.9%
で間違っていたということになります。
誤差のようなものですが,生きていたという予測のほうがかろうじて得意と言えそうです。
木構造の図示
学習した木を評価する手法にはいくつかあります。まずは木の構造を見たいと思うので木を図示します。
# 図示その1 import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(30, 30)) # whatever size you want tree.plot_tree(model, ax=ax) plt.show()
ちょっと見づらい感じになります。
より良く見せたいならば,graphvizやdtreevizなどを使うと良いと思います。
木構造のかっこいい視覚化
Kernelに以下のモジュールをインストールします。
!pip install dtreeviz pydotplus
# graphvizによる視覚化 import pydotplus as pdp file_name = "./tree_visualization.png" dot_data = tree.export_graphviz(model, # 決定木オブジェクトを一つ指定する out_file=None, # ファイルは介さずにGraphvizにdot言語データを渡すのでNone filled=True, # Trueにすると、分岐の際にどちらのノードに多く分類されたのか色で示してくれる rounded=True, # Trueにすると、ノードの角を丸く描画する。 feature_names=train_x.columns, # これを指定しないとチャート上で特徴量の名前が表示されない class_names=['Dead','Sruvived'], # これを指定しないとチャート上で分類名が表示されない special_characters=True # 特殊文字を扱えるようにする ) graph = pdp.graph_from_dot_data(dot_data) graph.write_png(file_name) # dtreevizによる視覚化 from dtreeviz.trees import dtreeviz viz = dtreeviz( model, train_x, train_y, target_name='alive', feature_names=train_x.columns, class_names=['Dead','Sruvived'] ) #viz.view() # 外窓で画像を開くのでKaggleのKernelではエラーがでる? viz.save(”dtreeviz.svg”)
木構造の解析:変数の重要性
では,解析の結果どの変数が生存に大きな影響を及ぼしたのでしょうか。 関数を叩くと一発ですが,以下のワードくらいは抑えておくとします。
- ジニ不純度: ノードによる分類のバラケ方(≒エントロピーみたいな)を数値化したもの→半々で最大。「1 - 各分類木のノードの分類割合の2乗和 」の値。CARTで使われている。→ sklearnでのデフォルト。
- 変数重要度:”ある特徴量で分割することでどれくらいジニ不純度を下げられるか。”を表す量。分割後のジニ不純度のサンプル数の分の和が,分割前のジニ不純度のサンプル数和よりどの程度小さくなったかを表す。
詳細な説明はこちらの記事に投げます。
treeオブジェクトでは変数重要度は.feature_importances_
というメンバに入っています。
# 重要度を表示 print(dict(zip(train_x.columns, model.feature_importances_))) # bar plot fig, ax = plt.subplots() plt.grid() ax.bar(train_x.columns,model.feature_importances_) fig.autofmt_xdate() # make space for and rotate the x-axis tick labels plt.show()
確かタイタニックでは女性と子供が優先的に救命ボートに乗せられたはずなのでこの洞察は合っていると言えるでしょう。
ただし重要度は,如何にクラスを上手く抽出できたかを表す無方向性の値なので,年齢が幼い方が良いのか女性の方が良いのかはわかりません。どうすればよいかちょっと調べてみます。(TODO)
テスト提出
とりあえず提出用のデータを作ります。
# テスト これを提出 ## test の中身 print(test.isnull().sum()) ## 中央値で埋める test = test.fillna(test.median()) ## 予測結果 test_predicted = model.predict(test) ## 提出用データ sample_submission['Survived'] = test_predicted sample_submission.to_csv('submission.csv',index=False)
このcsvをダウンロードしてSubmit Predictionの該当箇所にUploadします。
本当はNotebookでSubmitしたいんですが,Submitできそうなボタンがないような?(@2020年4月20日)なんでやろ?
スコアボードを見ると全く同じスコアの人々がいるので多分テストコードを走らせただけの人々がそれなりにいるのでしょう。 なお,正答率100%の人々もかなりいます。だって答えわかってるもんね。
最後に一応My Notebookを載せます。大してきちんと書いてないですけどね。
補足・TODO
気分がアガり次第書きます。スプラトゥーンの記事は書いたまま塩漬けになってるので今度ちゃんと書き起こします…
TODOs:
今回は最も単純な分類木を用いましたが,
補足:今後の参考
素晴らしい記事をありがとうございます。後ほど参考にさせて頂く予定です。
ボツデータ:欠損値をどのように埋めるか問題
欠損値をどのように埋めるかですが,3番めに重要と思われるAgeが抜けておりました。
絶対適当に埋めるべきではない気がするのですが,中央値,平均値,それとも削除して無視するのどれが良いのか比較しようと試みました。
結論から言うと結果は全く変わらず。なんか違いが出ると思ったんですがね…
# 欠損値除去 VS 補填 (なお,比較のためここでやっているが,本来はデータの前処理の段階でやるべき) ## train_yのうちtrain_xがNanを含む行をdropする train_y_dropna = train_y.drop(train_x[train_x['Age'].isnull()].index.tolist()) train_x_dropna = train_x.dropna() ## 埋める方 shallow copy train_x_fillmed = train_x train_x_fill0 = train_x train_x_fillmed["Age"] = train_x_fillmed["Age"].fillna(train_x["Age"].median()) train_x_fill0["Age"] = train_x_fill0["Age"].fillna(0) # それぞれFit model_dropna = clf.fit(train_x_dropna, train_y_dropna) model_fillmed = clf.fit(train_x_fillmed, train_y) model_fill0 = clf.fit(train_x_fill0, train_y) # 予測 ## validの欠損値補間はMedianで valid_x["Age"] = valid_x["Age"].fillna(valid_x["Age"].median()) pred_dropna = model_dropna.predict(valid_x) pred_fillmed = model_fillmed.predict(valid_x) pred_fill0 = model_fill0.predict(valid_x) ## 評価 print(accuracy_score(pred_dropna,valid_y), accuracy_score(pred_fillmed,valid_y), accuracy_score(pred_fill0,valid_y))
追記:ちゃんと変化しました。Medianで適当に埋めるほうが少し精度が下がりそうです。
[0.8148148148148148, 0.7811447811447811, 0.7542087542087542] 平均score 0.783
柏キャンパス-江戸川台駅 のグルメ攻略ガイド<Maniac>
どうも。本日柏とお別れしましたが,今後は終身栄誉柏民です。
さるネイティブカシワリアンの書いた記事に触発(挑発)されたので, 柏キャンパス-江戸川台のエリアにフォーカスしたグルメ攻略ガイドのManiac版を公開します。
TX民?ララポに行きなさい。
柏キャンパスのグルメがわかる記事
渋谷・本郷と暮らしてきた人々からすれば柏の葉は不毛の地と思われがちですが,よく見ると柏キャンパス周辺にはたくさん食事処があります。
以前簡潔に書いた柏キャンパス民のためのグルメガイド記事から約2年が経ちました。
その後,優秀なネイティブカシワリアン?氏に以下の記事にて補足(と煽り)をいただきました。
ということで,本記事の趣旨は上記の記事の店も含めよりManiacな 表通り(通称グルメストリート)にないお店を紹介することです。
和食
突然ですが,柏キャンパスには和食があまりありません。 チェーンを除くと基本中華の勢力が強いんですね。
鳥一や理彦などは必修ですので他の和食店を紹介します。 といっても実質2店舗くらいしかないんですけど。
わっ嘉
インドカレー屋の向かいの路地裏あたりにあるお店。昼はランチ,夜は飲み屋らしいです。
- ○:小皿の料理が美味しい
- ○(✖):店主がとてもフレンドリーでアットホームな雰囲気(苦手な人はいると思う)
- ✖:学生にはちと高い
- △:ランチしかやってない/日曜空いてない
”グルメ”としては結構完成度高いので自分へのご褒美としてどうぞ。夜も安く飲ましてくれるそうです。
12品の彩り定食が魅力の「わっ嘉」は、車いす料理人が江戸川台にオープンした和風ダイニング | リビングかしわWeb
美こま・お食事処
すごく評判が良さげなので行きたかったのですが,ついぞ行けなかったので名前だけの紹介にします。 レビューの感じもすごく良さそうなので期待を込めて入れています。
だんらん食堂
どうやらチェーンの類のようですがあまり知られていないのでここに記します。 自分でおかずを選んでご飯,豚汁と組み合わせるタイプの店です。
- ○:比較的手軽な値段(600円程度~)で健康的な定食がを自分で選んで食べれる。
- ○:250円くらいでご飯&豚汁を無限おかわりできる。うまい。
- ✖:江戸川台駅から西へチャリで10分。超遠い。
洋食
江戸川台駅の北のエリアには実はフレンチのビストロがあったりします。うっかり柏に来た友達などを接待してやる時に使うと良いでしょう。
フレンチ(コンツェルト・ブラッスリータケ)
駅から北に歩いて5分程度の位置に並んで2店舗があります。2店舗の違いはそんなに通ってないのでなんとも言えないですが個人的にはコンツェルトの方のメニューが印象に残っています。
- ○:駅から近い
- ○:1500円くらい出せばちゃんとしたコースが食べれる
- △:昼遅めになるとお得なコースがなくなってたりするので注意
- あくまで接待用??
イタリアン
イタリアンはいくつかあるようですが,行ったことのある店のみにします。
ボスコロッソ
駅の広場のすぐ裏にあるカフェっぽいイタリアンです。 1000円前後でランチが食べられます。
- ○:駅から近い
- ○:イタリアンなので基本はずれない
- △: 値段からして学生が通う店ではない気がする
ひるふぁ~む
みどり台郵便局の所にあるピザの美味しいお店。 ランチは1000円強くらいです。
- ○:ピザがうまい。サラダもうまい。おそらくパスタもうまいんだろう。
- ○:健康的な食事をした!という気分になる
- △:通うには高いなぁ…
イタリアンにはコスパ神サイゼリアがあるのでなかなかよっしゃ食べに行こうとはならないのが残念なところです。
中華
中華は激戦区です。イチオシは華記食府ですが,個人的には小満堂も好きです。 ホントはもっといろいろあるんですが紹介したい店だけにします。あと,ラーメンも同列です。
ラーメン食いたければ柏駅に行け。
ほんわ華屋
ヨークマートの近くにある品の良い中華屋という感じです。
- ○:胃に優しいタイプの中華
- △:1000円前後でドリンク付きランチ。やはりちょっと高い
親御さんなど連れてきてもいいでしょう。
ラーメンイヨマンテ
おじいちゃんが経営しているラーメン屋さん。ベーシックだけど結構好きな味です。
- ○:値段の割に量が多め。お腹いっぱいになれる。
- △:場所を知っている人は殆ど居ないだろう。近くに住んでいるならどうぞ。
さくらミニ
駅前北側にあるラーメン屋さん。意外と気づかないような位置にある。 まぜそばが人気のようだが,個人的には味噌ラーメンもイケる。
その他
その他の雑多なグルメ情報です。
ケーキ・パン屋
パン屋として必修なのはサフランですが,付近のちょっと路地を入った所にミレーというケーキ屋さんがあります。
ケーキのデザインも良く,値段も高くないので気になる教授にケーキを買う時に是非(?)
配達系
知っているか。柏キャンパスにUber eatsは来ない(2020年4月現在)。
配達系はピザか寿司,かわりものでパエリア程度しかないのでおうちで何か注文するには絶望的に適していない。ご注意。
スーパー
新人はわかりやすい駅前の”かどや”やヨークマートに行きがちですが,より食品の鮮度・コスパで優れているのはフーズマーケットセレクション西原店です。 位置はほぼ初石ですが,果物系の鮮度や安さ,暴力的な弁当のコスパ(~300円でそこそこ)を考えると行くべきはここ! 頑張って自転車を漕ぎましょう。
なお,近くの西原小学校では毎週日曜日の昼に有志でバドミントンをやっています。興味ある方はそちらもどうぞ。
16号沿い
キャンパスの裏から高校の脇の道を通って国道16号線に出れます。キャンパスから最も近い牛丼屋である”すき家”他よくあるチェーンやラーメン屋もあるので頑張って足を伸ばしてみましょう。 もつ煮太郎は必修な。
おわりに
柏キャンパス-江戸川台付近にはもっともっとお店があるのですが,ある程度書いて満足したのでこの辺で終わりにします。 みなさんもぜひお店を開拓して見てください。
【雑談】宝石をじゃらじゃらしたいのでAmazonUSで100カラット注文してみた話
この記事は私が疲れている時に勢いで書いたものになります。
読まれる方もどうか勢いをつけて読んでください。
追記:後に冷静になってから貴石達が届きました。石に詳しい人に是非助言をいただきたいです。
宝石じゃらじゃら欲
追記:この章は茶番です。結局AmazonUSで買いました。
宝石じゃらじゃらしたい。これは石に興味のある子供なら誰しも思ったことがあるはずである。
というかキラキラ輝く宝石を手に収めて眺めていたい欲は多分結構普遍的に人類の持っている欲だと思う。
昔は道端で透明な石英を拾って喜んでいたちびっこだった私はその後もビー玉を集めたりしたものだった。ラムネの瓶から抜いてくるやつがちょっとマットな質感があって好きだった。 鉱物に関する本も結構読んだ気がする(ただし小学生レベル)。
小さい頃はそれで我慢できたかもしれない。が,ちょっと大人になった今なんだか急に宝石をじゃらじゃらしたい欲が湧いてきた。
アクリルでもガラスでもない。「本物」をただじゃらじゃらしたいのだ。 はい,というわけで探してきました。安価にじゃらじゃらする方法を。
え?安価に??
サブスクリプションサービス
デザインアトリエカケラというサブスクリプションサービスがあることを最近になって知った。
コンセプトも夢があって,加工ミスや品質の問題により値段があまり付けられないような宝石を定期的にランダムに送り届けてくれるというものだ。
- 本来価値の高いものを訳ありということで安価に
- ある程度ランダム性をもって届けられる
という点において「訳あり高級カニ」系統に劣らないワクワク感を演出してくれる。
このランダムというのがおみくじ(ガチャ)感あってとても好きだが,私は今すぐじゃらじゃらしたいので今回は保留とした。
なお,こちらのショップはバラ売りも行っている。定期的に眺めたい感じはある。
ネットショッピング
こういう掘り出し物的なのはネットで買うのが一番である。どうやら「天然石」「ルース」「訳あり」などのワードでサーチすると良いらしい。
専門店で買う
普通に「天然石」「ルース」「訳あり」などのワードでサーチすると次のようなページに辿り着く。
まずページがずるい。きれいな形の宝石がじゃらっと並んでるだけでIQが下がってしまう。 しかも一個あたり結構安いのでホイホイ買えてしまうではないか。
ただし実際のサイズより写真ではかなり大きく見えていることに注意。
また,ぶっちゃけどの宝石をどの割合でとかを考えるのも面倒だ。 私は「じゃらじゃら」が注文したいんであって最適な石の選定には興味がなかった。たくさんじゃらるには結構値も張るので今回は保留とした。
ヤフオク・メルカリ・楽天etc
上記の専門店のような業者は一部フリマにも進出している。 友人曰く交渉次第で良いものを手に入れられるらしいが,私は宝石をじゃらじゃらしたいチンパンなので今回はパス。
しかし,メルカリって無限に時間潰せちゃうよね。
Amazonで買う
今や偽レビューや謎業者によって信頼性がガタ落ちしてはいるが腐ってもネットショッピング大手。Amazonにはこの世の半分くらいのものは揃っている。
たまたまそれっぽいのを見つけたので筆者は最終的にAmazonで買うことにした。
日本のAmazon
自分が探した中だと以下のようなランダムミックスが面白そうと思った。他にもいくつかあったがアクリルが混じっているなどの悪いレビューがあったりするので気をつけてみるべし。
5000円は少し買うのに勇気がいる値段だがサブスクリプションで毎月待たずとも一気にじゃらじゃらできる量が手に入る。 具体的にレビューの写真を見ていこう。
料金は5000円くらい,納期は2週間といったところか。
すでにありなのだが海外から輸入になるので仲介業者の存在により高くなっている場合がある。
英語で同様のショップを調べた所アメリカのAmazonにも同じ商品を見つけた。
アメリカのAmazon
海外のAmazonからものを買う手法は
があるが,目当てのものはAmazonGlobalで探せなかったので仕方なくアメリカのAmazonのアカウントを作成してみた。
これの良いところは仲介の店舗を介さないのでより安く買える点である。
https://www.amazon.com/gp/product/B006C4VWX2/ref=ppx_od_dt_b_asin_title_s00?ie=UTF8&psc=1
というか納期も日本より早いのでアカウント作る手間を惜しまなければアメリカのAmazonで買ったほうが良さげである。
なおコメントにもあるのだが,これはロットによって当たり外れが大きく,外れた人は低評価をつけている。 曰く,
- 欠けや傷が存在している
- 色や大きさに偏りがある
- 思ったより少ない
ということらしい。 しかしまぁ「本来価値のあるものバルクで扱ってみたい」という本趣旨を考えるにそこまで問題にはならなそうではある。
新年まだやっていないおみくじのつもりで買ってみることにする。ガチャっとな。
追記:届いたものについて
茶番終了です。 こちらが届いたものになります。アメシスト多めですが各色揃っているのがわかります。
内訳について
鑑定に関しては素人ですが多い順に並べていきます。
アメシスト(紫水晶):27個
一番多かったです。 手前から奥に向けて濃くなるように並べました。
ガーネット?:15個
次に多かったのが暗赤色の宝石です。多分ガーネットと思います。
水晶?:11個
石英(クォーツ)のうち特に透明度の高いものをクリスタルと言います。
他にも,トパーズなんかは透明のものもあるようです。 【ダイヤ】無色透明の宝石一覧【珍しい石】 - NAVER まとめ
以下の写真は少し自信がないので色ごとに分けています。
でもバルクで売ってる以上全部クォーツだと思います。
スモーキークォーツ(煙水晶):5個
水晶が放射線を浴びたり中にアルミニウムイオンなどが含まれている時にできると言われています。
シトリン(黄水晶):3つ + 何か1つ
紫水晶を加熱することで得られるのがシトリンらしいです。こちらは鉄イオンが関係しているとか。
残りの1つはなんでしょう?ペリドット,イエローベリル,トパーズあたりでしょうか?
翡翠?:2つ
濁り具合が翡翠を想起したので翡翠としましたがよくわかりません。 大粒で入っているのでそこまで貴重な石ではなさそうですが…
ということで袋に詰め替えて終了です。 じゃらじゃらは数回して満足してしまったので鑑定の方に時間を費やしました。
やはり福袋的な感じがして楽しいのですが正解が非常に気になる…石に詳しい人に聞きたいですね。
化学系の研究室なら紫外線ライトとかがあるんでしょうが…。
布教コーナー
宝石の国は良いぞ。
RICOH THETA_SCで多重露光・合成をしてPixel4みたいにクリアな星空を撮りたい
この記事はたまたま集まった技術&もの作りが好きな人たち Advent Calendar 2019の記事です。
はじめに
それはゴミのような論文を泣く泣く提出した日の夜でした。デスマを終えて見上げた空は信じられないくらい澄んでいてそれはもう最高に晴れやかな気分でした。
『こら撮らな!』と思いたって家からTHETAを引っ張り出してきて撮ったのが以下の写真です。
何度かパシャパシャ撮ってから帰りましたが,やっぱり後でビューワーで視るとまだ少し物足りない感じがします。
そこで,Google Pixel4でも使われている多重露光を行って合成してより精細な写真を撮るという手法を試してみたくなりました。
以下がPixel4で撮った写真の一例らしいです。カッコよすぎんか。私も撮りたい。
戦略
暗い夜空をきれいに取るためには露光時間を長くすればよいのですが,長すぎると星々が動いてしまいます。
Google Pixel4の戦略としては短めの露光を複数回撮って,それらを重ね合わせてノイズを除去しているようです。
ということで私も以下のような流れをとることにしました。
- 360°カメラで星空を複数回撮影
- 複数の画像を合成してノイズを落とす
- 星空の回転に合わせて画像を補正
諸々の処理のプログラムはPythonで書きます。
Step0: 星空を撮る
撮影に使っているのはRICOHの全天球カメラのTHETA_SCです。
RICOH 360度カメラ RICOH THETA SC (ブルー) 全天球カメラ 910743
- 発売日: 2016/10/28
- メディア: エレクトロニクス
扱いやすく個人的には360度カメラの入門にはもってこいと思います。
撮像素子はデジカメでは標準的な 1/2.3型
,レンズのF値は2.0
とまぁまぁの明るさなので15秒程度露光すれば都心の夜空は十分に綺麗に撮れます。
都市郊外で空を撮る感覚だと
露光20-30秒 ,ISO:100-160程度
が向いていそうです。
当日は小型の三脚に固定して撮影しました。
全天球カメラで撮った画像はEquirectalgular (正距円筒図法)という形式で保存されます。
本稿ではこの画像の中心を原点とします。
横方向:水平方向へは一周で2π,縦方向は緯度方向なので一周でπのスケールがあります。 したがってざっくり画像座標→正距円筒座標の変換を書くと以下のようになります。
hei,wid,_ = image.shape cx,cy = wid/2,hei/2 Equilecpoints = (np.array(polarpoints)-np.array([cx,cy]))/np.array([wid/2/pi,-hei/pi])
Step1:多重露光を合成する
ということで多重露光をした画像を合成する段階に入ります。 本来なら露光やISOを統一するべきなのですが,合成を思いついたのが撮影会後なので設定がバラバラな画像になってしまいました。
何も考えずにピクセル輝度の平均をとると以下のようにぼやけた画像になってしまいます。
そこでOpenCVにあるHDR合成(ハイダイナミックレンジ合成)のライブラリを使用してみることにしました。
HDR合成とは簡単に言うと露光条件の異なる画像を組み合わせて白飛び黒飛びを軽減した画像を撮像することです。
丁寧にまとめてくれている方がこちら(文献[1]): [CV]High-Dynamic-Range(HDR) Imagingについて - Qiita
今回は最も楽そうなMergeMertens
という手法を用いました。
Pythonで書くなら以下のような感じです。
import cv2
merge=cv2.createMergeMertens()
proc1 =merge.process([img1,img2,img3])
先ほどのようなモヤが生じる問題を回避できましたが,よく見ると星の動きが軌跡となっています。これはこれで乙ですがくっきりとした星空が見たいですね。
実際オリオン座の周辺に着目すると,星空が動いているのがよくわかります。
露光を繰り返した3分強の間にも星々はよく動いていることがわかります。
Step 2 星空の回転を補正する
先程のGIF画像の様に星空は一日で天球を一回転します。回転中心は北極星(Polaris)近辺の極北にあたります。
OpenCVの既存関数による位置ずれ補正
OpenCVにもAlignMTBというクラスがあり,calculateShiftという関数で位置ずれを補正できるそうですがプログラムの中身を見た感じこれは微小並進量のみの対応です。
また,一般の平面画像であれば射影変換によって回転ズレを補正できます。 ossyaritoori.hatenablog.com
しかし,Equirectangular画像においては回転操作は非線形変換になります。
Equirectangular画像における三次元回転
回転操作がどの様に影響するか,次の座標変換から導出することができます。
逆変換は以下の通りです。
星空はrが変化しない無限遠にあるので,変換後のEquirectangular画像のは天球上での三次元回転行列を用いて
と計算できます。
実際変形する際にはに対応する元画像の座標を逆変換から求め,ピクセル値をBilinear補間で求めます。
ロドリゲスの回転行列の計算
回転行列はOpenCVでサポートされているロドリゲスの回転行列を用います。
星空は原点と極北≒北極星を通る回転軸の反時計回りに回転します。1秒あたり1/240度なので微小のように思われますが,先程のGIFのように複数回露光をすると動きとして効いてきます。
変化を補正するために手動で極北の座標を決定した後,画像のExif情報から撮像の間隔を取得して回転変換をかけます。
# polarw :Equirectangularでの極北のtheta,phi nx = np.cos(polarw[1])*np.cos(polarw[0]) ny = np.cos(polarw[1])*np.sin(polarw[0]) nz = np.sin(polarw[1]) # omega = time_diff/3600/24*2*pi Rrod,_ = cv2.Rodrigues(omega*np.array([nx,ny,nz]))
このRrodは先程の式における回転行列Rの転置になります。
変換&重ね合わせの結果
ということで重ね合わせの結果です。
重ね合わせが成功して満足してしまいましたが,見えないものが見えるわけではないのでもっと暗い所に行って多重露光合成の威力を試したくなりました。 まずは休暇が欲しい…
STEP3 :自動で極北を抽出(途中)
SIFTやAKAZE特徴点でマッチングした画像は案の定,空と建物の境界にハマります。 星々をマッチングできれば回転を検出できますが,やはり空と建物・地面の分離を事前にやる必要がありそうです。
夜空部分の抽出
先程のマッチングでは信頼度の低いマッチング結果を除去していました。具体的には距離情報で縛りをかけています。
良いマッチングを得るためにユークリッド距離の比を元に選別するコードがよく書かれています。
for m,n in self.matchings: if m.distance < 0.75 * n.distance: #good.append()
この選別部分をなくした場合のマッチング点(AKAZE)はいくつか空の星々を捉えています。
したがって,最初の選別有りのマッチングにて建物と空の協会となるY座標を抽出し,(今回はY=619)Y座標がその閾値以下,すなわちカメラから見てより上方にある箇所を空の領域と仮定しました。 その結果絞り込まれたのが以下の17点の対応です。
3次元座標での回転推定
Equirectangular平面内ではお互いの変換は非線形になってしまうので,一度求めた特徴点の対応を三次元球面極座標に写像してから対応する回転行列を求めます。
三次元回転の最適化といえば金谷先生の教科書を参照するのが良いでしょう。
また,一応英語のテキストが無料で転がっていたりします。
https://igl.ethz.ch/projects/ARAP/svd_rot.pdf
を最小化する回転Rは対応点x_iとy_iを列方向に並べた行列X,Yを用いて表した
を特異値分解して得たから,
のように計算できます。なお,ここで推定したRはデータによっては厳密な回転行列にならない場合があります。それを補正するために同様の操作を繰り返します。 こちらは金谷先生の教科書に記述があるのでそちらを確認してください。
numpyでは以下のように書きます。np.linalgのsvdがVを転置の状態で出すことには注意です。
def estimate3dR(X,Y): S=np.dot(X,Y.T) U,W,V=np.linalg.svd(S) Rhat_ = np.dot(V.T,U.T) Rhat = np.linalg.multi_dot([V.T,np.diag([1,1,np.linalg.det(Rhat_)]),U.T]) # additional Step For Optimization U,W,V=np.linalg.svd(Rhat) Rhat_ = np.dot(V.T,U.T) Rhat = np.linalg.multi_dot([V.T,np.diag([1,1,np.linalg.det(Rhat_)]),U.T]) return Rhat
ランダムサンプルによるロバスト推定
ところで上記の手法でそのまま回転を計算すると,外れ値の影響を受けてしまいます。
この外れ値を除去するためにランダムサンプルによる手法が数多く提案されており,有名なRANSACなどはopenCVでの推定アルゴリズムにも採用されています。
二次文献ですが以下のブログがわかりやすいと思います。
RANSACの完全な実装はなかなかに面倒で,以下のブログを参考にしました。 クラスの継承を有効に使うコードを初めて書いた気がします。これはこれで長いので別記事で。
その結果,
- 手動設定の回転行列
array([[ 0.99989064, -0.00811519, -0.01236295], [ 0.00808042, 0.99996327, -0.00285925], [ 0.0123857 , 0.00275904, 0.99991949]])
- RANSACで推定した回転行列
array([[ 0.99994968, -0.00893819, -0.00455381], [ 0.00895058, 0.99995627, 0.00270865], [ 0.0045294 , -0.00274927, 0.99998596]])
はい。そこそこまともな値が抽出できてきました。 観測にノイズが混じっているのでこれは厳密には回転行列ではないので,補正が必要になります。(TODO:もう1つSubsectionが必要)
とりあえず,これをこのサイトにて分解して回転軸と回転角(rad)を抽出してみると,
[ -0.2625027, -0.4368636, 0.8603735 ], 0.0103961
これは,約143秒の露光という判定です…うーん,一分ほど短い? 北極星の位置も微妙にずれた検出をしているのでマッチングの精度が悪かったかRANSACの調整ミスですね。(TODO)
また,回転中心の位置推定もやはり手動には敵いません。これはなかなか手ごわいです。
とまぁこれ以上詳細に書くとめちゃ長くなるので別記事にします。 とりあえず,
- 星空の領域抽出
- 星空の運動の推定
がネックとなるため最初持っていたイメージよりかは簡単には自動化できないという感覚を得ました。 手動で北極星さえ指定すればExif情報で基本的に推定できるのでこの頑張りがそもそも必要かというのも少し怪しげですが気力があれば続きも記事にします。
TODO
- 空の領域とそうでない領域の分離:雑にDone
- 回転中心の自動抽出(特徴点マッチングを用いる):雑にDone
- 三次元回転推定まわりを整形:多分やる
- 公開ようにプログラムを整形:やる?
おまけプログラム's
PILでgif画像を保存する関数
from PIL import Image def saveasGif(images,filename='out.gif',duration=500,loop=0): imlist = [] for img in images: img_array = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) imlist.append(Image.fromarray(img_array)) imlist[0].save(filename, save_all=True, append_images=imlist[1:], duration=duration,loop=loop)
Exifデータから撮影間隔を取得する関数
from PIL.ExifTags import TAGS from datetime import datetime # Exif データを取得 def get_exif_of_image(file): """Get EXIF of an image if exists. 指定した画像のEXIFデータを取り出す関数 @return exif_table Exif データを格納した辞書 """ im = Image.open(file) try: exif = im._getexif() except AttributeError: return {} exif_table = {} for tag_id, value in exif.items(): tag = TAGS.get(tag_id, tag_id) exif_table[tag] = value return exif_table def get_diff_time(file1,file2): """Calculate capture interval of images from EXIF 指定した2つの画像ファイル間の撮影時間の間隔を計算する """ FMT = '%H:%M:%S' exif1=get_exif_of_image(file1)['DateTime'] exif2=get_exif_of_image(file2)['DateTime'] dtime = datetime.strptime(exif2[-8:], FMT) - datetime.strptime(exif1[-8:], FMT) return min(3600*24-dtime.seconds,dtime.seconds)
多くなりすぎたのでStep3に関しては別記事に移行する気がします。
参考文献等
[1] HDRなどに関して:[CV]High-Dynamic-Range(HDR) Imagingについて - Qiita
[2] 球面極座標に関する図表などの引用:360度パノラマ画像から平面画像への変換手法 - LASTMILE WORKS / DYNAMO TECH - R&D Project - Medium
[3] パノラマ画像から矩形画像を切り取るプログラム:GitHub - NitishMutha/equirectangular-toolbox: Handy tool for equirectangular images
他余裕でき次第追記
VSCodeでコマンドをショートカットから実行する設定方法 -docmuteで個別コンパイルができない問題の解決-
TexはずっとTexstudioで書いてきましたが,docmuteで分割コンパイルするときに困ったのでVScodeも併用します。
VScodeでコマンドをショートカットから実行する
ほぼ以下の記事の通りです。
keybindings.jsonを開く
まずは,設定を書き込むkeybindings.jsonを開きます。
コマンドパレット(ctrl+shift+p)を開き'key'と入力し「基本設定:キーボードショートカットを開く(JSON)」を選びます。
設定の書き方
開いたjsonファイルは[]
で囲まれた記述が連なっていると思います。
下のような表記を中括弧[]
の中に追記します。
{ "key": "F5", "command": "workbench.action.terminal.sendSequence", "args":{ "text": "echo ${file}\n" } }
これはworkbench.action.terminal.sendSequence
というコマンドを介してtext: に連なるターミナルのコマンドを実行しています。
ここでいうターミナルはVSCodeの下にあるターミナルのことで,実行結果や経過を確認することができます。
引数について
textの引数には以下のURLにあるような変数が使用できます。
例えば上の例の${file}
は現在VSCode上で編集しているファイルのフルパスとファイル名です。
Visual Studio Code Variables Reference
自分の設定
後半で解説する問題を解決するために自分が書いたファイルは以下のとおりです。
{ "key": "F5", "command": "workbench.action.terminal.sendSequence", "args":{ "text": "cd ${fileDirname}\n platex ${fileBasenameNoExtension}.tex \n dvipdfmx ${fileBasenameNoExtension}.dvi \n SumatraPDF ${fileBasenameNoExtension}.pdf \n" } }
「自分の今開いているファイルのディレクトリに移動」「platexとdvipdfmxでPDF作成」「SumatraPDFで表示」という順になっています。
docmuteパッケージの個別コンパイルができない問題
docmuteとは端的に言うとTexファイルの分割コンパイルのためのパッケージです。
以下のような構造で,mainがcahpter1,2をincludeないしinputしているときにmainでの全体コンパイルとchapter1単体のコンパイルをサポートしてくれます。
root/ ├ main.tex ├ chapter/ ├ chapter1.tex └ chapter2.tex
確認している問題
TexStudioのオートコンパイルの機能を使っている前提ですが, 他のファイルをコンパイルする前にmainをコンパイルした場合,chapter1のファイルを開いてコンパイルしても間違って全体がコンパイルされてしまう現象が起きてしまいます。
おそらく,main.texのコンパイルで生成されたファイルの依存関係を参照してしまうためと思われ,ターミナルで必要なコマンドを直打ちするか,一度中間生成ファイルをすべて取り払うことで想定した動作を実現できます。
VSCodeを使った解決
TexStudioで編集しているファイルをVSCode同時に開いて全体コンパイルと個別コンパイルでエディタを分けて対応する様にしてやろうというのが今回の解決策です。
この解決策は最適ではきっとないでしょうが,VSCodeとは今後長いこと付き合っていきそうなのでできそうなことはVSCodeにやってもらおうかなと思っています。
Barrierを用いてUbuntuとWindows間でマウス・キーボード・クリップボード共有
前回の続きです。これはかなり嬉しい。
Barrier with Ubuntu16.04
前回の記事ではBarrierでWindowsのPC同士でマウスとキーボードを共有しましたが,今度はWindowsとUbuntu間で試しました。 Ubuntu16は公式ではサポートしていないのでコミュニティのログを血眼になって探したところ,一つ前のバージョンで動作をさせることができました。
インストールする手順
インストールされるのはBarrier2.0です。現在の最新は2.1で2.2がbetaなのでギリギリ大丈夫だとは思います。
以下のコマンドを走らせることで,Barrierのソースを落としてきて自前の環境でコンパイルすることができます。
sudo apt-get -y install build-essential cmake qt5-default wget libxtst-dev libxinerama-dev libice-dev libxrandr-dev libavahi-compat-libdnssd-dev libcurl4-openssl-dev libssl-dev dh-make wget https://github.com/truatpasteurdotfr/barrier/archive/v2.0.0.tar.gz && tar xzvf v2.0.0.tar.gz cd barrier-2.0.0 && dpkg-buildpackage -us -uc | tee debian.log
barrierの実行ファイルはbarrier-2.0.0/build/bin
に作成されるのでディレクトリ移動して起動するかPathを設定しましょう。個人的にはディレクトリ移動して起動するスクリプトを書くほうが無難な気がします。
動作確認
さて,公式のissueなどではUbuntuマシンではクリップボード共有などはうまく行かないなどの報告が多いですが,偶々?私の環境ではきちんと動作させることができました。
欠点はWindows同士と比べて若干ラグを感じることです。タッチパッドのスクロールなんかは少し怪しいです。(ASUSのドライバがあんまり出来が良くないという説もあります。)
それでも無料のソフトを使ってマウス・キーボード・クリップボードをLinux-DOS間で共有できるというのは割とすごいことの様に思われます。 ちなみにLogicoolのFlowでも今の所対応できていません。
製品版のSynergyにも期待が持てますね。
とはいえ安心=課金
フリーソフトマニアが言うのもなんですが,安心は金で買うべきです。 Flow対応のキーボードやマウスをLogicoolから買えばもっと楽にファイル共有なども行えるでしょう。
というわけで広告を貼ります。これ別に私にお金はいるわけではないですけどね。
ロジクール ワイヤレスキーボード 無線 キーボード K380BK Bluetooth マルチOS:Windows Mac iOS Android Chrome OS対応
- 出版社/メーカー: Logicool(ロジクール)
- 発売日: 2015/11/05
- メディア: Personal Computers
- この商品を含むブログ (2件) を見る
- 出版社/メーカー: Logicool(ロジクール)
- 発売日: 2017/06/15
- メディア: Personal Computers
- この商品を含むブログを見る
matplotlibを用いた作図テンプレート(PDF図)
より良い記法をみつけ次第更新。こだわりポイントがある人はコメントで教えてください。
結論
折れ線プロット
PDFでplotをきれいに吐きたいならとりあえずこれ。
# setup import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages import numpy as np x = np.arange(0,10,0.5) y = np.sin(x) # plot template fig = plt.figure(figsize=(5,4),dpi=80) plt.plot(x,y,'-*',label='this is used for legend') plt.grid() plt.rcParams['font.family'] = 'Arial' plt.title('Title above') plt.ylabel('ylabel') plt.xlabel('Time [s]') plt.legend(loc='best', borderaxespad=0, fontsize=12) plt.show() # Save figure as PDF pp = PdfPages('pdf_test.pdf') pp.savefig(fig) pp.close()
写真とplotを重ねる
写真と重ねたplotはベクターにする意義が薄いがまぁ統一するならPDFか。
fig = plt.figure(figsize=(10,10),dpi=200) # Image Plot img = plt.imread(’Filename’) plt.imshow(img) plt.grid() # Scatter Plot plt.scatter(ConvertedX,ConvertedY,c=color_code) plt.colorbar() plt.xlabel('longitude') plt.ylabel('lattitude')
個々の機能について
Font変更
Times New Roman教の方はこちらから。 matplotlib のフォントはデフォルトで Bitstream Vera Sansですが,ArialやNew Romanなどがいい人は以下のように変更すればいいと思います。
自分用メモ: matplotlibを使ったグラフの調整 - Qiita
plt.rcParams['font.family'] = 'Arial' #plt.rcParams['font.family'] = 'Times New Roman'
figureのサイズなどをいじる
plt.figure()の中身をいじればfigureのレイアウトや解像度を指定可能。
fig = plt.figure(figsize=(5,4),dpi=80)
PDFにエクスポート
多用するなら関数化しても良いかも。なお,ただの画像ならplt.savefig('figure.png')
でOK。
from matplotlib.backends.backend_pdf import PdfPages def savefigPDF(fig,filename): pp = PdfPages(filename) pp.savefig(fig) pp.close()
subplot①
fig = plt.fgure(figsize=(10,10),dpi=200) ax=fig.add_plot(1,2,1) plt.plot() ...
subplot②
こちらの書き方の場合,plt.~~
で設定していた箇所が ax.set_~~
になる。図ごとに設定をいじれるが書き方が面倒と思ってしまう。
img = plt.imread("Map3forsave.png") fig, ax = plt.subplots(1, 2,figsize=(10,10),dpi=200) ax[0].imshow(img) ax[1].scatter(sigpos_img[:,0],sigpos_img[:,1],c=signal_freq) ax[0].set_xlabel() ax[0].set_ylabel() ax[0].set_title()
参考サイト
https://shimolab01.blogspot.com/2018/02/matplotlibpdf.htmlshimolab01.blogspot.com