祭囃子は遠く、

理系大学生のハッピーエヴリディを書いていきます。

ぷよぷよとランダムフィリング

こんにちは、すしくんです。

突然ですがみなさんこのツイート知っていますか?

統計物理専攻としてはなかなか興味を惹かれるものがあって、これをこうしたらどうなるんだろう・・・とか色々気になったものです。
しかし、当時はプログラミングが今よりも下手くそだったので、自分で解析することが叶いませんでした。最近は少しマシになったので自分でプログラムを組めたのでこの記事に至ります。


目次

はじめに

タイトルにもある通り今回は先行研究に従い、まずはランダムフィリングで色々条件を変えていったときに連鎖の分布がどう変化するかを見ていきます。具体的にはアスペクト比を保ったままスケールを変えていったときの平均連鎖数や、お邪魔ぷよを一定の割合で入れたときの変化です。
ランダムフィリングとは文字の通りですが、フィールドにぷよをランダムに敷き詰めていくことです。したがって一番の違いは現実のぷよぷよでは不可能な状況も発生するところです。(例えばフィールド真ん中から連鎖が始まる、12段3列目にぷよがあっても良いなど)
今回はとこぷよや対戦面を含めた盤面への主張はしないので、この場合こうなるんだとかそのくらいに思ってください。


プログラムの流れ

  1. あるスケールでのフィールドに4色のぷよをランダムに配置
  2. 連結を判定し、4連結以上でぷよを消す
  3. 連鎖数を測定
  4. 最初に戻る

拡張されたフィールドにおける連鎖数分布

まずは一応最初にやられている結果をフォローできるか確認しましょう。
また、記号についてですが13行6列のフィールドをL=1としています。26行12列はL=2、39行18列はL=3…と最初のフィールドからの倍率を表すものになっています。

L=1(先行研究)の結果です

f:id:ta_ichi:20180815142208p:plain

ピークの2−3連鎖がほぼ等しく、あとはlinearぽく落ちるとことかも似てるように見えますね。まぁこれでよしとしましょう

次にフィールドを等倍していったときの分布を見てみましょう。

f:id:ta_ichi:20180815142823p:plain

重なってて見えにくいですが、どのサイズもピーク左側が急峻で右側が緩やかな分布になっています。
分布の形自体はこの範囲内で収束しているように見えるため、どのサイズでも左右非対称で有限の歪度がある分布であることは間違いなさそうですね。

分布の形についてはわかったので次は量的なところをみていきます。
この比較だとLをどんどん大きくしていくと分布がどんどん右に動いていくだけで面白くありません。そこで各サイズでの最大連鎖数でスケールされた連鎖数で比較していきましょう。

f:id:ta_ichi:20180815143647j:plain

ピークに注目すると大きなサイズほど0に近づいていっているように見えますね。ここからわかることは、フィールドの最大連鎖数の増え具合よりも、ランダムフィリングしていったときの平均連鎖数の増え具合の方が緩やかということです。
最大連鎖数は\frac{13*6*L^2}{4}と計算していますが、せいぜいL^2のオーダーでしか増えません。一方平均連鎖数はこれを下回っていることが予測されます。

確認のためスケールされた連鎖数の期待値をみてみましょう

f:id:ta_ichi:20180815144838p:plain

緩やかにL→∞に向かって0に収束していきそうな感じが見えてますね。

片対数で見てみましょう(これでlinearならば指数関数的)

f:id:ta_ichi:20180815145006p:plain

うーん、綺麗な下凸でしょうか
両対数はどうでしょう(これでlinearならばべき関数的)

f:id:ta_ichi:20180815145134p:plain

これも微妙ですね、L=3以上でわりかし線形にも見えますが大まかには上凸でしょうか?

L=3以上でフィッティングかけてみましょう

f:id:ta_ichi:20180815145534p:plain

小さいサイズはものも見事に外れますね・・・
しかし大事なのはこのべき関数がどうなっているかです。

f:id:ta_ichi:20180815145821p:plain

ax^bという形でフィッティングしているので、−1.4程度で落ちていることがわかりました。
あんまり確かめるまでもなかった気はしますが、スケールサイズLに対して平均連鎖数はL^{0.6}程度のべき関数で抑えられそうということは一応やらないとわからなかったですね。
ただ、統計警察は気づいていると思いますが、これはフィッティング(回帰)させるモデルが多分ダメです。上の画像真ん中あたりのReduced chisquareに注目するとめちゃくちゃ小さいことがわかります、この値はだいたい1程度のオーダーだとよくフィッティングできてるなあという気持ちになるんですがこれは少し疑ったほうがいいですね。(多分)

ーおまけー
上の図から考えると、Lを大きくしていったときに早く収束するモードと遅く収束するモードがあると考えるのが自然でしょう。そこで補正項を入れてax^b+cx^dという形でフィッティング
f:id:ta_ichi:20180815150541p:plainf:id:ta_ichi:20180815150547p:plain

一応合わせるだけなら可能ですが、本当に合わせただけですね。もっと真面目にモデルを考えないと多分ダメです。

一定の割合でお邪魔ぷよを含むフィールドにおける連鎖数分布

今度はフィールドを固定してお邪魔ぷよをある一定の割合pで混ぜたときの分布をみてみましょう。
期待としては、ほどほどのお邪魔ぷよは多連結を切って連鎖になる方向に向かわせるはずなので、割合に対して平均連鎖数は上凸になる、と思ってました。(ました)

f:id:ta_ichi:20180815161207p:plain

0.05とか0.2とかがお邪魔ぷよの割合です、ダメそうですね、期待はずれです、普通に連鎖が減るだけです。

割合と平均連鎖数の図もみてみましょう。
f:id:ta_ichi:20180815161406p:plain

対数でみてみましょう
片対数
f:id:ta_ichi:20180815161423p:plain
両対数
f:id:ta_ichi:20180815161430p:plain

単調に減少していってますね・・・。

一応フィッティング、よくわからないので2次関数で合わせます。(どうせ適当な冪函数だし)

f:id:ta_ichi:20180815161903p:plain
f:id:ta_ichi:20180815161912p:plain

まぁ微妙ですねw


ただ、そもそもこのヒストグラムの出し方はもしかしたらよくないかもしれません。
おそらくもっとちゃんとやるならば、レプリカ法の精神に乗っ取って

  • ある割合でお邪魔を配置する
  • その隙間をぷよで埋める
  • 連鎖数を測定
  • 2、3を何度も繰り返してヒストグラムを取る
  • 同じ割合でまたお邪魔配置、2へ
  • こうしてできた「あるお邪魔配置でのヒストグラム」のアンサンブル平均をとる

とやるのが良さそうです。(めんどいからやらない)
今の解析でも割と大きなサンプルサイズではあると思うのでそこまで結果は変わらない気がします。

まとめと次の話

簡単にまとめると

  • ランダムフィリングにおいて最大連鎖数がL^2で増えるのに対し平均連鎖数はそれよりも緩やかにしか増えない
  • お邪魔ぷよは単純に平均連鎖数を下げる

といったところですね。

さて、次の進め方ですがアスペクト比を変えるとか、フィールドの倍率変更+お邪魔とかあると思いますが面白くなさそうなのでもっと違うことをします。
というか本来こっちがやりたかったことなのですが(専門なので)

「あるハミルトニアン(コスト関数)にしたがって統計力学的にフィールドを作る」

ということをやりたいと思います。(統計的にはq-state potts modelと等価)
どういうことかというと、例えばあるハミルトニアン(エネルギーと思ってください)を考えます。

\mathcal{H}=J\sum_{(i,j)_{NN}}^{}\delta_{p_i,p_j}

\delta_{p_i,p_j}はクロネッガーのデルタ、p_i,p_jは最隣接のぷよぷよです。このとき、Jが負だとするとこのエネルギーを最小化する解は「全てのぷよが同じ色」となることはすぐにわかると思います。正だと逆に連結が一切ないフィールドが解です。(piなどには色に対応する数字が入ってると考えてください)

このようにフィールドの状態はある程度この関数でコントロールすることができます。期待的には、フィールド一面に連結している状態と一切連結していないフィールドの中間を見るとそこそこ連鎖してくれそう。
さらに「温度」の概念も導入することで、温度揺らぎからくる効果により「盤面一色がエネルギー最小であるが何個かのぷよは違う色になる」ことができます。
これも連鎖の助けになりそうです。

このことから「ハミルトニアンの式そのもの」と「温度」をうまく組み合すことができればランダムフィリングで得られた平均連鎖数を越すことができるかもしれません。
連鎖数が多い方がエネルギーが低いという式を指定すれば連鎖数は増えるのでは?という意見もあると思います、しかしそれは結果が自明であり面白くないのでやりません。(そうしてできた形を見るのは面白いかもしれません)

また、今回やったことは「任意のハミルトニアンにおけるT→∞の性質」を見たと言い換えることもできます。このあたりも次の記事の導入で話したいと思います。

某での某某某

僕は怒っています。
就活のコーディングテストでよく使われる某某サイトの入力がゴミウンチ以下だからです。
普段競プロでしかほとんどプログラミングしなかったり、あんまり文字列扱わなかったりすると変な入力が来たらウンチを漏らしてしまうものです。

そこで某の注意点を少しだけまとめます。
主にc++想定です。Pythonも触れるかも。(あんまり触れませんでした)

言語が途中で変えられない

見ての通りです。しかしオンラインジャッジということもありTLE(時間オーバー)はあります。
問題をみて言語変えられないのは僕的に辛い。
某社のテストの時は、pythonを選んでしまいめちゃくちゃ最適化したものの、もう原理的に時間内に終わらないやろってレベルで無理でした。

結論
C、C++は正義、速度はパワー
でも難しい処理を普段Pythonに任せていると辛いものがある。

入力がcharだしクソ(僕が受けた範囲では)

これ僕のような人間からしたら面食らいます。(競プロの標準入力は優しい)
数字をどうこうしようってアルゴリズムを組めと言われたので、入力が数字だと思い込んで実装した後に気付く・・・ということもありました。
Pythonであればテキトーにやるだけで変換できるので特に問題はありません。
しかし、改行のエスケープ文字とかが文末に入っててややこしいことになったりします。(入力に含まれるのでそのままforとかで回すと色々やばい)
正直こんなことで頭を使いたくないのでcharで受け取ってintにキャストするやつを置いておきます。(テスト中に焦りながら調べた)

#include <stdlib.h> /* strtol */
#include <limits.h> /* INT_MAX, INT_MAX */

int main(){
char *endp;
char *str;
long r;

r = strtol(input, &endptr, 10);

if (input == endptr) {
   //文字列の処理;
}

if (r > INT_MAX || r < INT_MIN) {
   //オーバーフローの処理;
}

int val = (int)r; //intにキャスト
}

inputが標準入力で最終的にはvalが得たい値です。
atoiは??と思うかもしれませんが、あれは文字列をかますと0を返すはずなので入力が0の場合と分解できません。
後エラー時の挙動が環境で違うとかで危ないらしい。

少なくとも上のソースで僕は通したことがあるので大丈夫なはず。

感想を書かせてくる

これもわりとありがちなんですよね。
書かすのはいいんですよ、書かすのは。
ただ実装でカツカツのレベルの問題にぶち当たった時がクソ

「やべっ、時間足りない〜〜〜」

って時にあっ、感想(実装の注意点とか苦労した点とか書かせてくる)書かなきゃとなるわけです。
まぁこれに関しては初見でもどうにかできるとは思いますが、最初から計画を立てて置いたほうが良さそうですね。
事前説明にも多分書いてあるけど

終わりに

とまぁここまで書きましたが僕はもうこれといって頑張る必要もないのでみなさん頑張ってください。

データ解析諸々まとめ(R編)

例にもよって自分用、Rのコードでまとめていきます。

多分同じ内容でPython編も書きます。


以下目次

・はじめに

この記事ではパッケージの導入から基本操作をそれぞれまとめていく。
まだまだ触ったことがない手法が多いので使ったらどんどん追記していくスタイルにしておく。

あと、まだ使ってないけど使いたいやつは空でも章だけ作ってます。



~表記上の決まりごと~
train:訓練データ
test:テストデータ
y:目的変数


~環境(2018/7/11時点)~
・MacbookPrp HighSierra :  Rstudio ver.1.1.453 , R ver.3.5
・Windows10 : Rstudio ver.1.1.453 , R ver.3.5



あんまり解析手法の理論を理解していないのでヤバいこと書いちゃうかも>

・便利パッケージの導入と使い方

まずは補助的なパッケージの導入について

基本的に導入に詰まったところはRとRstudio のアップデートで解決した。

・R自体の更新方法
Rが3.5.0へメジャーアップデートしたので簡単アップデート - 清水誠メモ
・Rstudio の更新方法
HelpからCheck for Updates

・dplyr

フィルタリングとかグルーピングとかデータフレームの操作が早いやつ。
C++で書かれてるかららしい。

~導入手順~

  1. 各種ソフトアップデート確認
  2. R(Rstudio)のコンソールでinstall.packages("tidyverse")
  3. library(tidyverse)で確認

~基本操作~
ここに書くと膨大になってしまうので参考にした記事のリンクをはる。
dplyrを使いこなす!基礎編

行ごとに操作するのは簡単だが、列をある条件で動かしたりするのは厳しいっぽい。(多分、誰か知ってたら教えて)
例えば定数のカラムを消したい場合は

train_var <- train%>%dplyr::summarise_all(funs(var))#分散をそれぞれのカラムで計算し、代入
dropcol<-c()#落とすカラム名を入れていくやつ
for(i in 1:length(train)){
  if(dtrain_var[1,i]==0){
    dropcol<-append(dropcol,names(train)[i])#分散が0ならカラム名を追加
  }
}
train<-train[,-dropcol]#落とすやつ以外を訓練データに回す。

あと、%>%の記法もこいつらの中に入ってるぽい?
これは右に左のやつを引数として渡すという意味

scaled_data<-data%>%scale()

上のようにdplyr::を使わなくても良い(例はあるデータを正規化するコマンド)
可読性は高いか低いかわからないけどユーザーは割と見る。

・caret

ハイパーパラメータのグリッドサーチしてくれたりする。
対応してないパッケージもある(LightGBMとか・・・)

~導入手順~

  1. 各種ソフトアップデート確認
  2. R(Rstudio)のコンソールでinstall.packages("caret")
  3. library(caret)で確認

~基本操作~
これも膨大になるので参考記事を貼って主要部分のみ説明する。
Rによる機械学習:caretパッケージの使い方 | Logics of Blue

今の僕が詰まったのはtuneGrid
環境的な問題かもしれないけど、グリッドするパラメータを全指定する必要がある。

tuneGrid = expand.grid(mtry=1:20,splitrule="variance",min.node.size=5)

上はまた後ほど紹介するランダムフォレストのパッケージrangerのもの。(同様の内容を下でまた書きそう)
mtryは1から20の間でグリッドサーチ、splitrule、min.node.sizeは固定(最適化しない)
中身に適当なこと書いてエラーを見ると何を指定すべきか見れる。

~例文~
一応rangerでの例を挙げておきます。

modelRanger <- train(
  formula = y ~., 
  data = train, 
  method = "ranger", 
  tuneLength = 4,
  preProcess = c('center', 'scale'),
  num.trees = 500,
  tuneGrid = expand.grid(mtry=1:20,splitrule="variance",min.node.size=5),
  trControl = trainControl(method = "cv"),
  importance = 'impurity'
)
パラメータ名 デフォルト値 説明 推奨
method - 解析パッケージ指定 -
tuneLength - 各パラメータを何通りずつ試すか。 -
preProcess - 解析データの加工、この場合は正規化 -
tuneGrid - グリッドサーチするパラメータと範囲をあらわに指定 -
trControl - どのように解析するか、この場合はクロスバリデーション -

・doParallel

caretを使う場合にはセットで読み込んでおいたほうが良さそうな並列化のパッケージ。
現段階ではどの範囲で並列してくれるのかわかってません。(caret以外で)

~導入手順~

  1. 各種ソフトアップデート確認
  2. R(Rstudio)のコンソールでinstall.packages("doParallel")
  3. library(caret)で確認

~基本操作~
caretを使う前に以下を実行する

cl <- makePSOCKcluster(detectCores())#何コアで並列するか
registerDoParallel(cl)#並列するコアの登録?

detectCores()で使っているパソコンの最大コア数を取得できる。
上の例は実行環境での最大並列

また、原因はよくわかっていないが同じコードでもMac(4コア)では動かずWindows(8コア)では動いたことがあった。
初期化がどうこうとかよくわからないエラーが出るが、スペックの問題と思われる。(検索してもあんまりわからなかった)

・解析用パッケージの導入と基本操作

解析を実行するパッケージの導入方法と基本操作をまとめる。
この場合も各種アップデートは必須。あとは64bitじゃないパソコンはドブに捨てましょう。

・線形回帰

実はあんまりやったことないので使ったら書きます。

・主成分解析(PCA)

よくPCA(Principle Component Analysis)って略されるやつですね。
分散が大きい特徴量であるほどデータの本質を捉えているはず!という仮定から生まれた手法です。
直感的にもまぁそんな気がするなぁという。
そこで特徴量の分散が最大になる軸でデータを再構成し直す解析手法です。この辺の説明はネットに無限に転がっているので検索しよう。

一番参考にしたサイト
Practical Guide to Principal Component Analysis (PCA) in R & Python


~導入手順~
デフォルトで入ってます。

~基本操作~

result.pca <- prcomp(train,scale=T)

ほとんど上のコマンドで尽きています。
trainは解析したいデータで、scaleは正規化するかしないか。

~詰まったところ~
定数列があった場合スケールができないのでエラーが出る。
→上に書いた処理(分散が0の列を落とす)などを挟んでどうにかする。

~実行後~

#old.op <- options(max.print=999999)#特徴量がバカ多い場合は表示が溢れるのでこれをやる
sink("component.txt")#標準出力先を切り替える
round(result.pca$rotation, 3)#主成分を出力
sink("summary.txt")
summary(result.pca)#主成分が何パーセントデータを説明しているかがわかる。
sink()

f:id:ta_ichi:20180712010401p:plain

これのCumulative Proportionがそう

~その他応用~
成分解析→回帰に使う場合は目的変数を外して主成分解析、9割以上(一般的に何パーがいいのかは知らない)の主成分を拾ってきて回帰させる。

train <- data.frame(result.pca$x)
train<-train[,1:m]#mまでの主成分を拾ってくる(結果を見て変える)

これで適当な回帰にぶっこむ

・決定木系

決定木系の解析がkaggleとかだと一番使われてそうと思っている。
(決定木をスッキリ説明する語彙力がない)
理論の気持ちとかも充実させていきたいから勉強したらその辺も追記したい。

以下にわかりやすい図とかが乗っている気がします。
KaggleのHousePrices問題を決定木系のアンサンブルで解く - mirandora.commirandora.com

あと、決定木系は線形分離問題に弱いらしい。
そのためかXGBには線形用っぽい関数があります。

・決定木

実は使ったことがない、使ったら書きます。

・ランダムフォレスト(randomForest、ranger)

アンサンブル学習の一つ。弱学習機を〜〜って説明されます。
簡単に説明すると、多数決みたいなもんです。アンサンブル学習は全部多数決(アンサンブル平均をとる)ですけどね。
この場合は多数決取る人たち(決定木)は変わりません。変わる奴がXGBとかLGBMとかですね。

~導入手順(randomForest)~

  1. 各種ソフトアップデート確認
  2. R(Rstudio)のコンソールでinstall.packages("randomForest")
  3. library(randomForest)で確認

基本的に上のやつでいいとは思いますが、今はもっと高速な奴があります。
rangerって奴です。データが小さい場合は上のものでいいですが、特徴量が多い場合やデータが膨大な場合はこれが早いです。
あとcaretにも対応しているので使わない手はない(?)

速度比較
[1508.04409] ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R


~導入手順(ranger)~

  1. 各種ソフトアップデート確認
  2. R(Rstudio)のコンソールでinstall.packages("ranger")
  3. library(ranger)で確認

説明書(デフォルト値などはここで確認、しかしnum.treesのデフォルト値がよくわからない)
https://cran.r-project.org/web/packages/ranger/ranger.pdf


~基本操作~
・最小構成

ranger=ranger(y ~. , data=train.data)

分類も回帰もひとまずはこれで実行可能。分類の場合一応fcator(y)をかけといた方がいいかも?

・そこそこじっくり使うくらいの構成(回帰を想定した数値)

ranger = ranger(
    formula = y ~ .,
    data = train,
    num.trees = 100,
    mtry = 5,
    splitrule="variance",min.node.size=5,
    write.forest = TRUE,
    importance = 'impurity'
  )

これらのパラメータは重要そうなので説明します。(説明できないところはわかったら追記)

パラメータ名 デフォルト値 説明 推奨
num.trees - 決定木の数 -
mtry 特徴量のルート個数 1つの決定木がいくつの変数を持つか デフォルトの√n個らしい
splitrule "variance" - -
min.node.size 分類は1回帰は5 - -
write.forest True オブジェクトの保存をするかしないか Trueがいいと思います。
importance - 変数の重要度をどう表すか? -

赤くしたパラメータがcaretでグリッドサーチできる(しなければならない)ものです。

・caretでのグリッドサーチ

modelRanger <- train(
  formula = y ~., 
  data = train, 
  method = "ranger", 
  tuneLength = 4,
  preProcess = c('center', 'scale'),
  num.trees = 500,
  tuneGrid = expand.grid(mtry=1:20,splitrule="variance",min.node.size=5),
  trControl = trainControl(method = "cv"),
  importance = 'impurity'
)

上で同じコードを例に出していますが、なるべく章ごとに閉じた作りにもしたいのでもう一回貼ります。

パラメータ名 デフォルト値 説明 推奨
method - 解析パッケージ指定、今回はranger -
tuneLength - 各パラメータを何通りずつ試すか。 -
preProcess - 解析データの加工、この場合は正規化 -
tuneGrid - グリッドサーチするパラメータと範囲をあらわに指定 -
trControl - どのように解析するか、この場合はクロスバリデーション -

ランダムフォレストは結構使うので充実させたいね。

・Xgboost

まだ使ってないので、詳しくは使ったら。
導入と軽い調べ物だけしたので

~導入手順(Xgboost)~

  1. 各種ソフトアップデート確認
  2. R(Rstudio)のコンソールでinstall.packages("xgboost")
  3. library(xgboost)で確認

説明書
https://cran.r-project.org/web/packages/xgboost/xgboost.pdf

参考(パラメータなどについて)
xgboost のパラメータ - puyokwの日記

・LightGBM

XGBに対抗してマイクロソフトが作ったとか言われてるやつ。
基本はXGBと一緒らしい。そのライト版かつ正確な(と主張している)解析手法。
問題設定にもよりそうだけどkaggleでXGB と共にめちゃくちゃ流行ってるみたい。

残念ながらcaretには対応してないっぽい。

~導入手順(LightGBM)~
結構めんどい。
LightGBM/R-package at master · Microsoft/LightGBM · GitHub
一応ここに全部書いてあるので確認してください。あとMacでやったのでMac想定です。
gitとかもう入ってるよ〜って人はそこは飛ばして次に行ってください。

1.各種アップデート(R,Rstudio,Homebrewなど)
2.gitコマンドをターミナルで使えるようにする。

brew update
brew install git

3.Cmakeなどを使えるようにする。
gccがインストールされていないと多分動きません。

brew install gcc

この時バージョンも確認しておきましょう。僕は脳筋なのでgcc-7,gcc-8と打っていってコマンドが認識されたバージョンとしました。

brew install cmake

思ったよりインストールが長かった。待とう
4.公式のコマンドを入力

git clone --recursive https://github.com/Microsoft/LightGBM
cd LightGBM/R-package
# export CXX=g++-7 CC=gcc-7  # for macOS (replace 7 with version of gcc installed on your machine)
R CMD INSTALL --build . --no-multiarch

これでエラーでなければ多分大丈夫
5.Rのコンソールでlibrary(lightgbm)を確認

各種パラメータ
LightGBM/Parameters.rst at master · Microsoft/LightGBM · GitHub


~基本操作~
正直日本語でかつRで使ってる人全然いない(記事を書いてくれていない)ので、細かいところまでわかってないです。

・最小構成(最小もよくわかってないからとりあえず空)

empty

追記するので、それまではこの記事みてください、下の方にRで簡単なirisの分類をしてます。
LightGBM ハンズオン - もう一つのGradient Boostingライブラリ

・僕が今使ってるやつ(こっちは回帰)

#訓練用と評価用に分割
set.seed(123)
inTrain <- createDataPartition(y, p=.9, list = F)#caretの関数、均等な分割をしてくれるらしい
tr <- train[inTrain,]
va <- train[-inTrain,]
tr_ta <- y[inTrain]
va_ta <- y[-inTrain]
#lgb用のデータセット作成
lgb.train = lgb.Dataset(data.matrix(tr), label = tr_ta)
lgb.valid = lgb.Dataset(data.matrix(va), label = va_ta)

#ここから解析
params.lgb = list(
  objective = "regression"
  , metric = "rmse"
  , min_data_in_leaf = 5
  , min_sum_hessian_in_leaf = 100
  , feature_fraction = .5
  , bagging_fraction = .5
  , bagging_freq = 4
)

lgb.model <- lgb.train(
  params = params.lgb
  , data = lgb.train
  , valids = list(val = lgb.valid)
  , learning_rate = 0.003
  , num_leaves = 180
  , num_threads = 2
  , max_depth = 8
  , reg_alpha = .3
  , reg_lambda = .1
  , min_child_weight = 33
  , zero_as_missing = T
  , nrounds = 5000
  , early_stopping_rounds = 600
  , eval_freq = 50
)

実行すると環境によるかもしれませんがWarningがアホほど出てきますが、無視してオーケーです。
ちなみにこんなやつ
[LightGBM] [Warning] No further splits with positive gain, best gain: -inf

ハイパーパラメータの推奨値や意味はある程度まとまったら追記します。

パラメータ名 デフォルト値 説明 推奨
- - - -
- - - -
- - - -
- - - -
- - - -

SVM

使ったら書きます。

getch・getcheに関するめも(Mac/Unix用) 追記4/12(キーボード入力まわりについて)

とあるプログラムをCで書いてたら

#include <conio.h>
getch()

の部分でconio.hがないよと怒られる。

調べるとmac(というかunixには実装されてない模様)

他の手立てもあるらしいが、あんまり難しいことを考えたくないため、どうにかできないか検索

nonberanai.blog.fc2.com

僕「なーんだ、もうあるじゃんw簡単簡単w」

と言っていたのも束の間、includeフォルダにgconio.hを移動できない

su使って移動させればええやろ!と思ったが、それでもPermission deniedされる。

(ちなみにsuを有効化するのも1手順必要だった
Mac でルートユーザを有効にする方法やルートパスワードを変更する方法 - Apple サポート

おかしいなと思いつつ調べるとusr以下は最近のアップデートでルートユーザーでもアクセスできなくなっている模様

rootlessとかいうやつで。

やばくないですか?(やばくはない)

というわけで解除しました

qiita.com

これの通りにやればOKです。

ただ、rootlessはマルウェアからの防衛的な意味もあるようなので

作業が終わった同じ手順でcsrutil disableの部分を

csrutil enable

というコマンドにして有効化しておきましょう。

終わり。


追記4/12
・kbhitについて

#include <stdio.h>
#include <conio.h>

int main(){
	int key;
	int flag = 1;

	while(flag){
		if(_kbhit()){
			処理1
			flag = 0;
		}
	}
}

キーボード入力がヒットする→キーボード入力があると処理をするような関数

これもMacにはデフォルトでないみたいです、上のヘッダファイルにもありませんでした。

そこでここ参照

tricky-code.net

解決

・矢印キーの入力について
特殊入力はなんか色々特殊らしい、安直に調べたやつではうまく認識できなかった

Console Operation for RT-Series

結局ここで解決した。以下例文(unix用)インデントおかしいけど許して

#include <stdio.h>
#include <gconio.h>
#include <termios.h>
#include <unistd.h>
#include <fcntl.h>


int kbhit(void)
{
    ~略~
    }

    return 0;
}

int main(){

	switch(getche())
   {
    case 0x1b:
    switch(getche())
      {
			case 0x5b:
			switch(getche()){
			        case 0x41: printf("↑\n"); break;
		                case 0x42: printf("↓\n"); break;
		                case 0x44: printf("←\n"); break;
		                case 0x43: printf("→\n"); break;
			}
      }
    break;
   }

return 0;
}

ボルツマンマシン

この前のこれ

f:id:ta_ichi:20180223193050g:plain

の記事を書いていこうと思います。

はじめに

動画のこいつは何者なのか、とかいう質問が出てくると思いますが最初に答えます。

人類に限りなく近いゴリラです。

本題に入りましょう、今回のモデルは情報分野だと「ボルツマンマシン」とか「ホップフィールドネットワーク」とか言われてるやつです。

物理科の民は「え、ただのイジングモデルじゃん・・・。」ってなるやつです。

この辺の説明は後でやります。

流れとしては、説明書いてコードはって結果はって終わり!ってなると思います。

各モデルの説明

上で出てきたモデルの説明をしていきます。

物理の話から情報分野に下る流れでいきます。


イジングモデル(Ising-model)

統計物理とかやった人はわかると思いますが、磁性体のモデルの中で一番(?)簡単なやつです。

まず、棒磁石を思い浮かべてください。おそらく大勢の人は赤と青に塗られ、それぞれにN、Sと書いてあるやつを思い浮かべると思います。

こいつらを半分に切ったらどうなるでしょう?

答えは「長さが半分になった棒磁石が二つできる」です。知ってる人も知らない人もいると思いますが・・・。

では、それぞれもう半分にしたらどうなるでしょう?

もう簡単だと思いますが「さらに半分になった棒磁石が四つできる」ということになります。

これを繰り返していくと、これ以上分割できない「磁石の最小単位」にたどり着くことがわかると思います。

これが「スピン」です。

このスピンは矢印に例えられ、矢印がみんなで同じ方向を向いていると全体としてN極だとかS極だとかが生まれるわけです。

イジングモデルというやつもこれをモデル化したもので、矢印は上か下かのどちらかしか向きません。

全部上に向いていればその方向がN極で、逆はS極ということになります。

f:id:ta_ichi:20180223195528p:plain

大体こんなイメージです。全部上向いてますね

この時矢印(スピン)がどこを向きたいかと言うものを指定しなければいけません。

それがハミルトニアン(エネルギー関数)です。今回の場合

f:id:ta_ichi:20180223202644p:plain

こう書けます。本来は量子力学波動関数の重なりがどうとか議論して積分したりして得られるものですが、直感的にも良さそうなことがすぐにわかります。

またSはベクトルで書いていますがイジングモデルにおいては±1しかとりません(上か下か)、あとsummationはあらゆるijの組で取ることにしています。

例えばJijが全部+1の場合を考えましょう。しかし、どういう状態が「答え」でしょうか?

現実世界というのは基本的に万物が全体としてエネルギー安定な方向へ向かっています。高いところから物を落とすと下に落ちるのもそうです。

何が言いたいかというと、このエネルギー関数は現実では常に最小化される状態(向き)になるということです。(厳密には自由エネルギーが最小化)

全てがJij=+1となる時は簡単で全てが上、または下に揃えばエネルギーが最小になることがわかると思います。

これはスピンが何個ある時でも同じです。(N>2) 磁石は全体として同じ方向を向いているので、こういったJijによって振る舞いが決まっていそう、ということも納得できると思います。

例ではJijが全て等しい時を考えましたが、Jijが色々な値を取るとそれぞれの矢印がもはやどういった向きをとれば良いのかわからないと思います。実際この場合は解析的に解くことはできません。


・ホップフィールドネットワーク(Hopfield-network)/ボルツマンマシン(Boltzmann-machine)

まずはホップフィールド ネットワークについてです。これはもうイジングモデルそのものです「スピンをJijと言う結合定数で結んでいる」と考えるとなんとなくネットワーク感も出てくるのではないでしょうか。具体的には下図のようになっています。

f:id:ta_ichi:20180324005006p:plain

Jijの値が有限=繋がっている(繋がってないところは0)と考えるとイジングモデルはこのようになり、情報分野ではホップフィールド ネットワークと呼ばれます。

また、CNNなどは順方向に情報が伝わっていき、ニューロンが次々と発火していくような向きを持った構造になっています。しかしこのモデルは向きを持たない無向グラフの一種です。

次にボルツマンマシンについてですが、ホップフィールド とほとんど等価です。何が違うのかと言うと、「温度」の概念が取り込まれていることです。

温度を取り込む利点として、記憶を想起させる時に便利とかはあると思いますが、他にはあるんでしょうかねえ〜〜勉強してないのでよくわかりません。


・結局何が違うのか

結局一緒なんじゃん、となると思います。はっきり言うと

結局一緒

です。本質的には何も変わりません。

じゃあ情報と物理で何が違うかと言うと「問題設定」です。

物理では、物質構造や現象から「Jijを予め決めてスピンの向きを決める」のに対して

情報では基本的に「スピンの向きからJijの組みを推定する」と言う問題を解きます。

逆問題の関係になっているわけですね。

例えば今回の話だと「2値画像(1ピクセル=1スピン)からJijを推定(学習)」し、それを今度は物理の問題として「与えられたJijの組みに対してのスピンの向きを決める(学習画像を想起する)」と言う二段階になっています。

学習・アルゴリズム

ここまででイメージがついたと思うので(思うので)実装に入ります。

正直実装にあたって触れていない重要なことがあるんですが、1つの記事に書ききれないのでいくつか天下り的に紹介します。気になったら調べてください。


学習方法

ホップフィールド における二値画像の学習方法は古くから知られていて「Hebb則」ってやつを使います。

f:id:ta_ichi:20180324012032p:plain

Jij:i番目とj番目のピクセルの結合定数
N:総ピクセル
p:覚えさせたいp番目の画像
S:2値をとる学習させたい画像のピクセル

こんな感じです。今回の設定だと与えた画像のセットから決定論的に決まるので推定とかではないですが、立派な学習です。

また、Jiiなど添え字が同じやつは0にします。(自分自身への結合は0ということ)


想起方法

メトロポリスヘイスティングス法」を使います。

これは簡単に言うと「乱数で1つのピクセルの向きを決めて、それがエネルギーを下げる方向ならば採用する」というのを繰り返すモンテカルロ法の1種です。ボルツマンマシンの文脈では少し間違ったことを言っていますが、だいたいあってるのであとは調べてください。

もしかしたらこれ単体の記事を書くかもしれません。

実装

ではこれらの概念を使って実際に実装します。

import numpy as np
import cv2

def local_E(spin,i,j,img,img2,img_size):#ローカルな1スピンまわりのエネルギー計算

    E = 0
    for k in range(img_size):
        for l in range(img_size):

            N = img_size**2
            
            #2値(0 or 255)の画像からいちいち直接Jijを計算する。この際0は-1となり、255は1となるように調節
            #予めJijの配列を作成しても良いが、添え字を考えるのがめちゃくちゃ面倒
            #画像が二枚なので手打ちですが、画像自体をリスト化して、for回すのが賢いと思います
            tmp = 2*img[i][j]/255-1 #1枚目画像
            tmp2 = 2*img[k][l]/255-1
            tmp3 = 2*img2[i][j]/255-1 #2枚目画像
            tmp4 = 2*img2[k][l]/255-1
            jij = (tmp*tmp2+tmp3*tmp4)/N #Hebb則からJijを算出
            if i==k and j==l:
                jij = 0.0 #自分自身への結合は0

            E = -1*jij*spin[i][j]*spin[k][l] + E #p番目の画像の全ピクセルに渡ってsummationをとる

    return E

#hyper parameter
img_size = 82 #画像の大きさ(正方形を想定)
border = 105 #2値にする際のボーダーライン、お好みで調整
N = img_size**2 #総ピクセル数
seeds = 11 #乱数シード
sweep = 10 #画像の更新回数(sweep*N回更新)
temp = 0.06 #温度、小さな値だと綺麗に元画像が思い出せるが、計算コストがかかる

#1枚目画像読み込み
img = cv2.imread("pipo.png")
img = cv2.resize(img, (img_size, img_size)) #一応指定したサイズと形にリサイズ
img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) #グレースケール
ret, img = cv2.threshold(img, border, 255, cv2.THRESH_BINARY_INV) #2値化

#二枚目読み込み、上と同様
img2 = cv2.imread("ash.jpg")
img2 = cv2.resize(img2, (img_size, img_size))
img2 = cv2.cvtColor(img2, cv2.COLOR_RGB2GRAY)
ret, img2 = cv2.threshold(img2, border, 255, cv2.THRESH_BINARY_INV)

#2値画像をgazoフォルダに吐き出し(フォルダなどは適当に変えてください)
cv2.imwrite('./gazo/pipo.png', img)
cv2.imwrite('./gazo/ash.png', img2)

#スピン配列(画像を思い出させる用)を作成
spin = np.zeros((img_size,img_size))
tmp_image = np.zeros((img_size,img_size),dtype = int)

#+1-1を適当に配置
for i in range(img_size):
    spin[i] = np.random.choice([-1,1],img_size)

#どのスピン(ピクセル)を更新するかを決めるための配列
xy = np.zeros((img_size),dtype = int)
for i in range(img_size):
    xy[i] = i

np.random.seed(seeds)

for M_i in range(sweep):

    ran = np.random.rand(img_size**2)
    
    k = 0 #画像の吐き出し用

    for i in range(img_size):
        for j in range(img_size):

            #x座標とy座標をランダムにチョイス(普通にやる場合は順番に全部更新する)
            x = np.random.choice(xy,1)[0]
            y = np.random.choice(xy,1)[0]

            E_old = local_E(spin,x,y,img,img2,img_size) #1つのスピンが感じるエネルギーを計算
            spin[x][y] = - spin[x][y] #向きを反転させてみる
            E_new = local_E(spin,x,y,img,img2,img_size) #更新後の1スピン周りのエネルギーを計算
            dE = E_new - E_old #更新前後のエネルギー差を計算
            #本来更新することで系全体のエネルギーが下がるか上がるかを計算しなくてはならない
            #しかし、差分をとるため、注目しているスピンまわりだけ計算すれば十分

            if ran[k] > np.exp(-dE/temp): #乱数が更新採択確率を上回ればreject
                spin[x][y] = - spin[x][y] #元に戻す

            if k%50 == 0: #50回更新するごとに画像を吐き出し(gif用)
                for i in range(img_size):
                    for j in range(img_size):
                        tmp_image[i][j] = int((spin[i][j] * 255 + 255)/2) #スピンを0,255の2値画像に変換

                new_image_path = './gazo/'+str(M_i) +'_'+str(k)+ '.png';
               	cv2.imwrite(new_image_path, tmp_image) #書き込み

            k += 1

ざっとこんなもんです

色々書き足りていないところはありますが、上のプログラムをいい感じに書き直せばもっとたくさんの画像を記憶できるでしょう。

ただ、画像同士がほぼ直交していないとたくさんの画像は記憶できないようです。

結果

上のプログラムだと二枚記憶させているので、実は初期配置によっては

f:id:ta_ichi:20180324113347g:plain

違う人が出てきます。

まとめ

一番簡単な画像記憶をやりましたが、どうでしょうか。

多分説明が足らないところがあると思います。モチベがあればここに付け足すか別の記事で補完したいですね(やらない)

画像記憶以外でもボルツマンマシンですぐできそうなこと見つけたら実装してみようと思います。

藤江杯いってきた話

金沢レジャーランドで開催された第五回藤江杯に行ってまいりました。

詳細はこんな感じ

地方で開かれていますが、毎回参加人数が50人を超える人気のある大会です。

参加者のほとんどが遠征者という稀有な大会なので、色々な地方のぷよらーをつまみ食いできるいいチャンスです。

僕は今回(も)ぷよらーの方に車を出していただいたので、ワイワイガヤガヤしながら大阪から金沢へーーー



まず、到着後もりもり寿司という「質の高い寿司がたまたま回っているだけ」のお寿司屋さんにいきました。

ついつい3000円使ってしまうものの僕は旅先での予算は常に無限大、財布の紐が緩いどころか底が抜けている状態です。

爽やかな顔で会計を済まし、ACぷよ通のある施設へ出陣します。



最初に初心者大会が行われるわけですが、僕の友人であり犯罪者のぽかりすくんが予選で無双していました。決勝トナメですぐ負けた

僕は出ないことになっていたため一般大会の方に全てをかけていましたが、すんなり負けてしまい真顔で音ゲーコーナーに

f:id:ta_ichi:20180226012538p:plain

しかし、先週普通に落ちて凹んでいた中伝に合格することができ、一気にウキウキモードへ

大会は結局僕の友人であるキウくんが優勝しました。おめでとう



その後は食事会があり、新潟勢を迎え予約していたゆめのゆの大部屋宿泊会場に向かいました。

部屋はこんな感じ

f:id:ta_ichi:20180226010931p:plain

一人一人のスペースも確保されている上に何と言っても「安い

入館料や部屋代込みで一人当たり3700円程でした。遠征者のみなさん、おすすめです。

「仮眠スペースなら2000円で済むじゃん!」という方々もいらっしゃると思います。

確かに少人数での利用ならその方がいいかも知れませんが、10人集まるならばこのプランがいいと思います。

それはやはり「修学旅行感」に集約されるでしょう。

こんな大部屋に男だけで泊まったらやることは一つ

恋バナ

ですよね??



はい、嘘をつきました。僕はそういったものとは無縁の世界の住人なのでぷよぷよ対戦をしながらツイッターのオタクたちについて話し合っていました。

しかし、お話ししながらゲームしながら広々とスペースを使えるのはアドバンテージですし、テレビもあるので交渉次第ではうまく使えるかも知れません。

ゲーマーの皆さんは是非



翌日は早めに起床し、ジムで大胸筋を鍛えた後、港の食堂で食事をとりました。

これについても色々書きたいですが割愛

最後は忍者寺という全く忍者に関係ない寺に行き、施設内案内ツアーに参加したのですが、キウくんが係員に私語を注意されたところがピークでした。


帰りは車の中でドヒャドヒャ笑いながら運転していたらいつの間にか大阪にいました。

今回の藤江も楽しかったですし、早く就活を終わらせて次の大会にも笑顔で出場したいですね。

では

LSTMっていうオモチャを見付けました

はじめに

これはぼくのインターンメモです。

LSTM

標語的にはリカレントニューラルネットワーク(RNN)の一種でLong short-term memoryの略だとか言われてますね。

この手法の何が他と違うかというと、ある決められた(人の手で入れる)長さのデータを保持しつつ(たまに忘却しつつ)

学習していくというところみたいですね、すごそう。下のリンクの記事は読んだ事ないですが貼ります

qiita.com


何に強いかというと時系列データだとか、文章(前後の繋がりがあるため)でその辺メインで使われてるらしい。

qiita.com
deepage.net


画像認識もいけるみたいだけど、わざわざやる意味はいまいち分かりません。

動かしてみた

株価予想とか腐る程出てくるので他のでやりたい、文章どうこうはあんまり興味がない

んで、今回は


素数予想」


やります。

CNNで素数判定!みたいな話は結構出てくるんですが、時系列データとして生成するってのはあんまり出てこなかったんで。

まぁ多分やった結果面白くなかったから誰も言ってないないとか、数学的に不可能と証明されているとかそんな感じでしょう。

とりあえず素数をDLしましょう

素数2357: 素数一覧ダウンロード

10000個のやつでやりました、計算重いからね。

毎回思うんですが、こういうcsvファイルってプログラムに読み込ませる以外の用途思いつかないのに、なんで余計なヘッダー付いてるんでしょう。

エクセルで編集はあまりいい思い出が無いので、エディタで消しましょう。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math
from keras.models import Sequential
from keras.models import load_model
from keras.layers import Dense
from keras.layers import LSTM
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error


#データ読み込み
dataframe = pd.read_table('primelist10000.txt')

#変数の抽出ためのデータフレーム作成
df = pd.DataFrame()

#目的変数
df['prime'] = dataframe['Number']#csvのカラム名をぼくが勝手に変えました
max_prime = max(df.prime)#データ範囲での最大値を取得
df.prime = df.prime/max_prime#絶対値0~1でスケール、LSTMはこれが推奨されてるらしい

#説明変数+上と同じ操作
df['index'] = dataframe.index
max_index = max(df.index)
df.index = df.index/max_index

#ここでnumpy-arrayにする。型は適当に設定
df = df.values
df = df.astype('float32')

#データ分割
train_size = int(len(df) * 0.9)#データの9割で学習、1割でテスト
train, test = df[0:train_size,:], df[train_size:len(df),:]

#LSTMが受け取れる形にしましょう、[データ行数][変数の個数][ルックバックの個数]みたいな形でしか受け取れない
def make_dataset(dataset, look_back=1):
    dataX, dataY = [], []
    for i in range(len(dataset)-look_back-1):
        xset = []
        for j in range(dataset.shape[1]):
            a = dataset[i:(i+look_back), j]
            xset.append(a)
        dataY.append(dataset[i + look_back, 0])
        dataX.append(xset)
    return np.array(dataX), np.array(dataY)


look_back = 10#上で出てきましたが、保持しておくデータ数で今回は適当に10とおいた
trainX, trainY = make_dataset(train, look_back)#訓練データリシェープ
testX, testY = make_dataset(test, look_back)#テストデータリシェープ

Hidden_Layer = 300#隠れ層の数、多いほどコスト上がるけど多いほどいいとかなんとか
Epochs = 1000#学習回数
Batch_Size = 10#バッチサイズ(まとまった単位で最適化する感じかな?よくわかってなさそう)

#モデル作成
model = Sequential()
model.add(LSTM(Hidden_Layer, input_shape=(testX.shape[1], look_back)))#隠れ層の指定とインプットデータの指定、ルックバックの指定
model.add(Dense(1))#出力の数
model.compile(loss='mean_squared_error', optimizer='adam')#損失関数と最適化法の指定
model.fit(trainX, trainY, epochs=Epochs, batch_size=Batch_Size)#学習開始
model.save("prime_LSTM2.h5")#学習したモデル保存
#model = load_model("prime_LSTM.h5")#読み込む時はこっち

#予測値
trainPredict = model.predict(trainX)
testPredict = model.predict(testX)

#リスケール
trainPredict = trainPredict * max_prime
trainY = trainY * max_prime
testPredict = testPredict * max_prime
testY = testY * max_prime


# 誤差率評価
err_rate_train = (abs(trainPredict-trainY)/trainY)
err_rate_test = (abs(testPredict-testY)/testY)

#以下予測と実際の素数のプロット
test_len = len(testY)
x_pl = np.arange(test_len)
plt.scatter(x_pl,testY,label="actual")
plt.scatter(x_pl,testPredict,label="predict",color="r")
plt.legend()
plt.xlim([0,50])
plt.savefig("ac_pre.png")

#予測された数値を整数として吐き出し
np.savetxt('pre_train.csv',trainPredict.astype("int"),delimiter=',')
np.savetxt('pre_test.csv',(testPredict.astype("int")),delimiter=',')
np.savetxt('test.csv',testY.astype("int"),delimiter=',')

結果

f:id:ta_ichi:20180215231609p:plain

うーん、微妙そうだけど素数が書く軌跡?みたいなのは若干拾えてそう。わからん!

というか予測した値が実際素数かどうか判定しなきゃ意味がないですね!このプロットはなんなんだ!(一応平行移動して心眼で見る事は可能)

そのプログラムは面倒なのでまた今度回しておきます・・・。

まとめ

よくわからん。

というか、int型での推定(離散的な推定)ってそれ用に何かあるのかな

色々改善できそうだけどもうやらなさそうです。

おしまい