はじめに

livedoor提供のお天気Webサービスに興味を持ったため、Djangoを利用したお天気検索ページを作成することにしました。(実装環境はPaizaCloud)

APIの概要：市町村ごとのID(以後、「市町村ID」と呼称)をリクエストパラメータとして引き渡すと、IDに紐づくお天気情報(JSONデータ)が返却されるAPI。
http://weather.livedoor.com/weather_hacks/webservice

以下がサンプルコード。

import requests

url = 'http://weather.livedoor.com/forecast/webservice/json/v1'
payload = {'city':'290010'} # 奈良県の市町村ID
tenki_data = requests.get(url, params=payload).json() 
print('地方：', tenki_data['location']['area'])
print('都道府県：', tenki_data['location']['prefecture'])
print('市町村：', tenki_data['location']['city']) 
print()
print(tenki_data['forecasts'][0]['dateLabel'] + ' の天気は ' + tenki_data['forecasts'][0]['telop'] + ' です。')
print(tenki_data['forecasts'][1]['dateLabel'] + ' の天気は ' + tenki_data['forecasts'][1]['telop'] + ' です。')

以下が実行結果。

地方： 近畿
都道府県： 奈良県
市町村： 奈良

今日 の天気は 曇り です。
明日 の天気は 晴れ です。

作成予定の画面

【1】地域情報一括更新画面
・市町村IDを管理しているxmlを取り込み、テーブルとして管理する為の画面。
xml：http://weather.livedoor.com/forecast/rss/primary_area.xml

・作成するテーブルは「都道府県テーブル(親)」と「市町村テーブル(子)」。「市町村テーブル(子)」で「市町村ID」を管理する。

・画面上の更新ボタンを押下することで以下の処理を行う。
　1. 「都道府県テーブル(親)」、「市町村テーブル(子)」の初期化。
　2. xmlファイル読み込み、解析。
　3. 解析したxmlファイルのデータを「都道府県テーブル(親)」と「市町村テーブル(子)」へ格納する。

【2】お天気検索画面
・画面上の「都道府県プルダウン」「市町村プルダウン」を選択し、確定ボタンが押下されたタイミングで天気情報を画面表示する。「都道府県プルダウン」には「都道府県テーブル(親)」のデータを、「市町村プルダウン」には「市町村テーブル(子)」のデータを設定する。

・「都道府県プルダウン」の選択値に応じて「市町村プルダウン」のデータを制御する。
(例：「都道府県プルダウン」が「奈良県」の場合、「市町村プルダウン」には「奈良、風屋」を選択肢として表示する。)

・《可能であれば》GoogleMapを表示し、選択された市町村の箇所にピン付けする。
(プルダウン選択後か確定ボタン押下後かは検討中。なお、フィジビリティー確認も不十分な為、見送る可能性あり。)

Next

まずは【1】地域情報一括更新画面を作成します。

はじめに

Mouse DictionaryをPDFで使うために、PDF.jsを用いるという方法がある。(@kroton様の参考記事)

しかし、私のようなめんどくさがりは、これすらも面倒くさいので以下のようなツールを作った。

• ワンラインでネット上のPDFをPDF.jsで開く
• ローカルに落としたpdfをわざわざ削除しなくてもよい

完成物

https://github.com/monado3/pdfjs-broker

開発環境

• Manjaro linux
• Python 3.7.3
• Docker 18.09.5-ce

実現方法

1. Dockerコンテナ内にpdfを落とす。
2. そのままDockerコンテナ内でPDF.jsが機能するPHPビルトインサーバを立てる。
3. ローカルのブラウザからDockerコンテナ内のpdfを閲覧する。
4. ダウンロードしたpdfはDockerコンテナが終了すれば(再起動時や明示的にコンテナを停止させた際)自動的に消える。

また、複数のpdfを同時閲覧することを前提とし、自動で動的にDockerコンテナとローカルのポートフォワーディングをするために、Pythonプログラムにコンテナ立ち上げを仲介させることにした。

課題点

DockerイメージのEntrypointとなっているシェルスクリプトが終了すると、例え-dオプションをつけてrunしてもコンテナはなぜか終了してしまう。
それにより、コンテナ内でPDFのサービングが開始したことをPython側にシンプルに伝える方法がなく、コンテナ内の必要処理が終了する前に、Pythonがブラウザを開いてしまう。

最後に

実装があまりスマートではないので、どしどしアドバイスやPR待ってます。

はじめに

Mouse DictionaryをPDFで使うために、PDF.jsを用いるという方法がある。(@kroton様の参考記事)

しかし、私のようなめんどくさがりは、これすらも面倒くさいので以下のようなツールを作った。

• ワンラインでネット上のPDFをPDF.jsで開く
• ローカルに落としたpdfをわざわざ削除しなくてもよい

完成物

https://github.com/monado3/pdfjs-broker

開発環境

• Manjaro linux
• Python 3.7.3
• Docker 18.09.5-ce

実現方法

1. Dockerコンテナ内にpdfを落とす。
2. そのままDockerコンテナ内でPDF.jsが機能するPHPビルトインサーバを立てる。
3. ローカルのブラウザからDockerコンテナ内のpdfを閲覧する。
4. ダウンロードしたpdfはDockerコンテナが終了すれば(再起動時や明示的にコンテナを停止させた際)自動的に消える。

また、複数のpdfを同時閲覧することを前提とし、自動で動的にDockerコンテナとローカルのポートフォワーディングをするために、Pythonプログラムにコンテナ立ち上げを仲介させることにした。

課題点

DockerイメージのEntrypointとなっているシェルスクリプトが終了すると、例え-dオプションをつけてrunしてもコンテナはなぜか終了してしまう。
それにより、コンテナ内でPDFのサービングが開始したことをPython側にシンプルに伝える方法がなく、コンテナ内の必要処理が終了する前に、Pythonがブラウザを開いてしまう。

最後に

実装があまりスマートではないので、どしどしアドバイスやPR待ってます。

TL; DR

日本語文書分類タスクを機械学習で解くとき、下記の文書ベクトル表現手法ごとの精度を比較しました。
Github: https://github.com/nekoumei/Comparison-DocClassification/tree/master/src
※最近のGithub、jupyter notebookのレンダリングがよく失敗するのでnbviewerも貼っておきます
https://nbviewer.jupyter.org/github/nekoumei/Comparison-DocClassification/tree/master/src/

• Bag of Words
• TF-IDF
• Word2Vecの平均値
• Doc2Vec
• SCDV
• SWEM

結論としては下記の3点です。

• 問題によるので銀の弾丸はない（あたりまえ）
• BoW, TF-IDFのような古典的手法も案外悪くない
• 個人的にはSWEM推し

なんて定性的な結論なんだ、、

問題設定

以下2つの問題について実験しました。

Livedoorニュースコーパスのカテゴリ分類

利用データ
https://www.rondhuit.com/download.html

9カテゴリの日本語ニュース記事です。1カテゴリ500~1000記事くらいあります。
日本語の文書分類についてググるとよく出てくるデータセットです。
各記事の文書から、どのカテゴリに属するかを解くマルチクラス分類問題として扱います。
記事ごとにファイルが分かれていてちょっとめんどくさいので1csvに変換するNotebookもつくりました。
https://nbviewer.jupyter.org/github/nekoumei/Comparison-DocClassification/blob/master/src/NewsParser.ipynb

Twitterツイートネガポジ分類

利用データ
http://bigdata.naist.jp/~ysuzuki/data/twitter/

詳細は上記ページに書いてありますが、特定のトピックに関するツイートのネガポジをアノテーションしたデータです。
ツイート本文はTwitter API規約のため公開されていないため、statusIDを用いて自分でAPIを叩いて取ってくる必要があります。
今回は、ルンバに関するツイートがネガティブか、ポジティブか予測する2クラス分類問題を考えます。

APIを使ってツイートを取得するNotebookも下記に公開しています。（API KEY, ツイート本文は載せていません）
https://nbviewer.jupyter.org/github/nekoumei/Comparison-DocClassification/blob/master/src/GetTweetText.ipynb
ルンバに関するツイートは2000件ほどありますが、現時点でapiで取ってこれたのは半分くらいでした。ネガポジも半々くらい。

比較する文書ベクトル表現手法

下記のとおりです。細かい説明はしないです。

• Bag of Words
• TF-IDF
• Bag of Words + Truncated SVD
• TF-IDF + Truncated SVD
  • BoWやTF-IDFはすごい高次元になるのでSVDで次元削減して使うこともあるのでやりました（kaggleでたまに見る）
• Word2Vecの平均値
  • 文書中に登場する単語ベクトルの平均値を文書のベクトルとする手法です。思考停止でよくやる（私が）
• SWEM -max
  • 文書中に登場する単語ベクトルのmaxを取る手法です。arXivにもあります。詳しくは下記のブログがとても参考になりました。
  • https://nykergoto.hatenablog.jp/entry/2019/02/24/%E6%96%87%E7%AB%A0%E3%81%AE%E5%9F%8B%E3%82%81%E8%BE%BC%E3%81%BF%E3%83%A2%E3%83%87%E3%83%AB%3A_Sparse_Composite_Document_Vectors_%E3%81%AE%E5%AE%9F%E8%A3%85
  • arXiv: https://arxiv.org/abs/1805.09843
• SCDV
  • 実装についてはSWEMで紹介した記事での実装を参照しました。本当にありがとうございます！
• SCDV -raw
  • 上記ブログ記事にも書かれていますが、SCDVの元実装には最後に小さいベクトルをゼロに押しつぶす(compress)処理が入っています。理由がよくわからない（スパースにして計算量を減らすためか？）ので上記記事同様この処理を除いた場合もやりました。
• Doc2Vec
• Doc2Vec - epochs=30
  • なぜかよく精度が悪いDoc2Vecですが、下記記事の解説によるとデフォルトのエポック数だと不十分らしいです。
  • そこでエポック数を30にした場合も追加しました。
  • https://buildersbox.corp-sansan.com/entry/2019/04/10/110000

また、上記手法のうちword2vecを用いる手法（w2v平均、SWEM、SCDV）については対象データセットを用いてword2vecを学習した場合と、学習済ベクトルを利用する場合の2パターン試しました。
学習済ベクトルは下記のfasttextでWikipedia日本語を学習したものを利用しています。ありがとうございます。
https://qiita.com/Hironsan/items/513b9f93752ecee9e670

前提

環境

Ubuntu18.04 LTS
Python3.7
メモリ32GB
下記記事で自作したPCです。
https://qiita.com/nekoumei/items/9304f355d66a247215f3

比較方法

1. 文書を一切前処理せずMeCab+NEologdで分かち書きする。
2. 上述の文書ベクトル表現手法でベクトルにする。
3. LightGBMで5-folds CVする。 (Stratified Kfold)
4. Out Of FoldのAccuracy5つずつと、モデルの学習+推論にかかった時間を比較する。

実験Notebook

ニュースカテゴリ分類: https://nbviewer.jupyter.org/github/nekoumei/Comparison-DocClassification/blob/master/src/Classification_News.ipynb
ルンバツイートネガポジ分類: https://nbviewer.jupyter.org/github/nekoumei/Comparison-DocClassification/blob/master/src/Classification_Roomba.ipynb
どちらもやってることはほぼ同じです。ちゃんとメソッドかクラスにすればよかった。

結果

Notebook
https://nbviewer.jupyter.org/github/nekoumei/Comparison-DocClassification/blob/master/src/Evaluate.ipynb

ニュース記事カテゴリ分類

Accuracyのboxplot+swarmplot

BoWが一番精度良い、という結果になりました。面白いですね。
カテゴリの分類って、特定の単語を含むかどうかで分類できそうなので確かにそれはそうな気もします。
BoW、TF-IDF以外ではcomporessしないSCDV（学習済ベクトル利用）が良さそうです。このへんは上述のGOTOさんの記事と同様ですね。
個人的には、w2vのmeanが比較的良くないこと、SWEMがシンプルな手法なのにめっちゃ良い感じなのが印象的です。
また、BoW、TF-IDFについては未知語へのロバストネスがちょっと懸念なのでこの結果を見て「BoW最強！！」ってわけではないと思っています。
ちなみにBoW、TF-IDFはだいたい80,000次元くらい、SCDVが18,000次元、SWEM他w2v系手法は300次元です。

Accuracy(mean)とモデル学習時間のscatterplot

学習にかかった時間と精度の関係性を可視化しました。
かかった時間にはベクトル化する時間は含めてません。純粋にfitしてpredictする時間のみです。ベクトル化時間を含めるとSCDVが更に増えます。
BoW, TF-IDFが次元の割に計算時間短いのは、疎なのでLightGBMがうまくやってくれたんでしょうね。
これを見るとSCDVの時間のかかり方が特徴的ですね。
学習時間に制約があるケースだとSWEMとかの方が良いかもしれません。（問題による）

ルンバツイートネガポジ分類

Accuracyのboxplot+swarmplot

データがちょっと少ないのもあってやや分散が大きいですね。
この問題ではw2v(学習済fastText)のmeanが一番よさそうです。
ニュース記事と違い、ツイートなので1文書あたりの単語数が少ないのでBoW等よりもベクトル表現の方がより文書を表現できてるってことでしょうか。
先程と比較して、SCDVはあまり良い結果が得られていません。GMMでのクラスタリングによる効果があまりワークしないお気持ちはなんとなく分かる気がします。
また、こちらでもSWEMは安定して良い精度が得られています。
Doc2Vecに関しては、エポック数を増やすことで大幅に精度が改善していますが、他の手法よりは劣っています。
Wikipediaとかの巨大データで学習済のDoc2Vecでやってみるとまた違った結果が得られるかもしれませんね。

Accuracy(mean)とモデル学習時間のscatterplot

あいかわらずSCDVが突出していますね。
こちらではcompress処理を入れるか否かで処理時間が大きく異なっています。やっぱり処理時間短縮のための処理だったんでしょうか。

おわりに

いろいろな手法を試してみましたが、個人的にはSWEMが良いなあと思いました。
シンプルな実装、低次元、低計算時間で十分な精度がいずれの問題でも得られています。
勿論、解く問題によって何が最適なのかは変わるので、都度検討する必要はあります。
あと、日本語の文書分類に使えるデータが全然なくてそこが一番詰まりました。

参考文献

Livedoorニュースコーパス: https://www.rondhuit.com/download.html
Twitter日本語評判分析データセット: http://bigdata.naist.jp/~ysuzuki/data/twitter/
SWEM論文: https://arxiv.org/abs/1805.09843
SWEM, SCDV実装: https://nykergoto.hatenablog.jp/entry/2019/02/24/%E6%96%87%E7%AB%A0%E3%81%AE%E5%9F%8B%E3%82%81%E8%BE%BC%E3%81%BF%E3%83%A2%E3%83%87%E3%83%AB%3A_Sparse_Composite_Document_Vectors_%E3%81%AE%E5%AE%9F%E8%A3%85
Doc2Vecについて： https://buildersbox.corp-sansan.com/entry/2019/04/10/110000
fastText学習済ベクトル: https://qiita.com/Hironsan/items/513b9f93752ecee9e670

はじめに

グリッドサーチを行うにあたり細かい制御を行う必要があり、パラメータの組み合わせ生成等が面倒だったため調べたときのメモ

環境

• Windows 10
• Python 3.6
• scikit-learn 0.20

やり方

グリッドサーチをやる場合、各パラメータの探索範囲をDictで定義し、それをGridSearchCVの引数に与えると思う。同じ形のDictをインプットとして、グリッドサーチを自力でやろうとした場合、Dictからモデルのset_paramsの引数に与えられる形式で、全パラメータの組み合わせを生成する必要があるが、この実装がかなり面倒くさい。

scikit-learnのGridSearchCVクラスのソース等をいろいろ調べたところ、ParameterGridクラスを用いると簡単にできることが分かった。

百聞は一見に如かず、とりあえずソースを見てみよう。まずは必要なものをインポートする。

import pandas as pd
import numpy as np

from sklearn.model_selection import ParameterGrid
from sklearn.model_selection import KFold

from sklearn.metrics import matthews_corrcoef
from sklearn.metrics import roc_auc_score

from sklearn.ensemble import ExtraTreesClassifier
from sklearn import metrics

次にモデルおよびGridSearh用のDictの定義。そして、ParameterGridを用いてDictに対する全パラメータの組み合わせのリストを生成。

   clf = ExtraTreesClassifier()
   param_grid = {
        "class_weight" : ["balanced", "balanced_subsample", None],
        "max_features" : ["sqrt", 0.3, 0.6, 0.9],
        "n_estimators" : [10]
       }

   param_dicts = list(ParameterGrid(param_grid))

そして、全パラメータの組み合わせのリスト毎にモデルの評価を実施することで、自前のグリッドサーチを実施。

   auc_all = []

   for i, param_dict in enumerate(param_dicts):
        auc_cvs = []

        for train_index, test_index in KFold(n_splits=10, random_state=42, shuffle=True).split(X):
            X_train, X_test = X[train_index], X[test_index]
            y_train, y_test = y[train_index], y[test_index]

            clf.set_params(**param_dict)
            clf.fit(X_train, y_train)
            y_predict = clf.predict(X_test)
            y_predict_score = clf.predict_proba(X_test)[:, 1]
            auc_cvs.append(metrics.roc_auc_score(y_test, y_predict_score))

        auc_all.append(np.array(auc_cvs).mean())

   print(auc_all)

decision threshold等を変化させて指標を最適化するなど、より細かいグリッドサーチを行う場合は、こんな感じグリッドサーチを実装すればでできますよ、という話でした。

記事の趣旨

理系元カノとの間に生じたコミュニケーションに関する課題に、簡単なアプリケーションのプロトタイプを開発して、対応していく過程を日記として残しています。第二弾です。

第一弾はこちらです。
https://qiita.com/tantan0327/items/8d02a75233c2c54db013

課題

淡白な文章やネガティブなLINEを相手に送ってしまい、相手を機嫌悪くさせてしまうこと。
・相手からの丁寧なメッセージに対して、「了解」だけみたいなそっけないメッセージを送ってしまう。
・ネガティブな文章を返す。（例：「最近忙しすぎて辛いな。」等）

解決策

相手が機嫌を悪くする文章になっていないかを判定するシステム。入力された文章について以下の判定を行うシステムを実装します。
①10文字以内になっていないか。(文字数チェック)
②ネガティブな文章になっていないか。(ネガポジ分析)

プロトタイプ

日本語でのネガポジ分析について、GitHubで調べ、今回はGitHubの以下のレポジトリが良さそうだったので、こちらを利用させていただきました。
https://github.com/liaoziyang/negapoji

様々な単語についてスコアを割り当てた辞書を用いて、文章に含まれる単語のスコアの合計値で判定しています。

結果

①文章が短かったら以下のようになります。

②10文字以上のネガティブな文書だったら以下のメッセージを利用します。

③10文字以上のポジティブな文書だったら以下のメッセージを利用します。

　プログラムを動かしてみて、スコアを数値で出力した方が、入力した文書に修正を加えた時にどう変化するかわかりやすくてUXが良くなると感じました。辞書のスコアの割り当て方についても改良できないか検討してみようと思います。

感想

　より丁寧で、相手のことを思った、メッセージであるかを判定するための追加条件も検討中です。具体的には「ありきたりの文章になっていないか」を判定するアルゴリズムをどう実装するか考えています。一旦、こちらはUIを実装する際に「自分にしか言えない言葉で話していますか?」と書いてある確認モーダルを表示するようにしようかと思っています。
　まだまだNLP初心者なので、まずは少しづつ既存のものを改良して、技術も磨いていければと思います。

背景

• 以下の記事に沿って Gmail API を叩く Python スクリプトを作成
  • https://news.mynavi.jp/article/zeropython-22/
• スクリプトを Lambda で動くようにして、実行したところ [Errno 30] Read-only file system が発生
  • ローカルで実行したときには起きなかった

原因

• 権限ファイル (client_id.json) の書き込みでエラーになっている模様
• エラーのスタックトレースを見ると、oauth2client ライブラリ内で発生しているようだった

"errorMessage": "[Errno 30] Read-only file system: 'client_id.json'",
  :
    [
      "/var/task/oauth2client/file.py",
      79,
      "locked_put",
      "f = open(self._filename, 'w')"
    ],
  :

• self._filename は Storage クラスのコンストラクタ引数
  • https://github.com/googleapis/oauth2client/blob/master/oauth2client/file.py#L33
  • 作成したスクリプトでは client_id.json のパスを渡しており、client_id.json はスクリプトと同じディレクトリに置いていた
• Lambda が書き込みできるのは /tmp 配下のファイルのみなので、Read-only エラーとなった

対応

• client_id.json を一旦 /tmp に退避させてから Storage クラスを作成する

before

    user_secret_file = 'client_id.json'
    store = Storage(user_secret_file)
    credentials = store.get()

after

import shutil
 :
    user_secret_file = 'client_id.json'
    user_secret_tmp_path = '/tmp/' + user_secret_file
    shutil.copy2(user_secret_file, user_secret_tmp_path)

    store = Storage(user_secret_tmp_path)
    credentials = store.get()

参考

• https://stackoverflow.com/questions/45608923/aws-lambda-errormessage-errno-30-read-only-file-system-drive-python-quic

背景

• 以下の記事に沿って Gmail API を叩く Python スクリプトを作成
  • https://news.mynavi.jp/article/zeropython-22/
• スクリプトを Lambda で動くようにして、実行したところ [Errno 30] Read-only file system が発生
  • ローカルで実行したときには起きなかった

原因

• 権限ファイル (client_id.json) の書き込みでエラーになっている模様
• エラーのスタックトレースを見ると、oauth2client ライブラリ内で発生しているようだった

"errorMessage": "[Errno 30] Read-only file system: 'client_id.json'",
  :
    [
      "/var/task/oauth2client/file.py",
      79,
      "locked_put",
      "f = open(self._filename, 'w')"
    ],
  :

• self._filename は Storage クラスのコンストラクタ引数
  • https://github.com/googleapis/oauth2client/blob/master/oauth2client/file.py#L33
  • 作成したスクリプトでは client_id.json のパスを渡しており、client_id.json はスクリプトと同じディレクトリに置いていた
• Lambda が書き込みできるのは /tmp 配下のファイルのみなので、Read-only エラーとなった

対応

• client_id.json を一旦 /tmp に退避させてから Storage クラスを作成する

before

    user_secret_file = 'client_id.json'
    store = Storage(user_secret_file)
    credentials = store.get()

after

import shutil
 :
    user_secret_file = 'client_id.json'
    user_secret_tmp_path = '/tmp/' + user_secret_file
    shutil.copy2(user_secret_file, user_secret_tmp_path)

    store = Storage(user_secret_tmp_path)
    credentials = store.get()

参考

• https://stackoverflow.com/questions/45608923/aws-lambda-errormessage-errno-30-read-only-file-system-drive-python-quic

タイトルはこれで適当につけています。

きっかけ

一回試してみたかったので、以下を読みながらLightGBMを試してみた。
https://www.codexa.net/lightgbm-beginner/

インストールするのがめんどくさかったので、環境はGoogle Colabにて。

座学

アンサンブル

アンサンブルの手法は大きく分けると「バギング」「ブースティング」「スタッキング」となる
参考：https://www.codexa.net/what-is-ensemble-learning/

• バギング
  • 並列的に学習
• ブースティング
  • 前の弱学習器の結果を次の学習データに反映
• スタッキング
  • 二段階に学習を積み上げる

※今回はアンサンブルの各手法を掘り下げるのが目的ではないので、ふんわりとさせておく。

LightGBMって？

勾配ブースティングは複数の弱学習器（LightGBMの場合は決定木）を一つにまとめるアンサンブル学習の「ブースティング」を用いた手法

仮に最初に訓練を行う決定木を1号、次に訓練を行う決定木を2号としましょう。まずは決定木（1号）でモデル訓練を行い推測結果を評価します。決定木（1号）の推測結果と実際の値の「誤差」を訓練データとして、決定木（2号）の訓練を行います。N号の決定木はN-1号の決定木の誤差（Residuals）を学習するわけです。

実践

せっかくなので、こっちのデータでやってみる
https://prob.space/c

ただ、codexaとコンペだとファイルの数が違うけど、いったん無視して続けてみる。

kmnist-test-imgs.npz
kmnist-train-imgs.npz
kmnist-train-labels.npz

ファイルアップロード

# filesモジュールでアップロード
from google.colab import files
uploaded = files.upload()

これでファイル選択できるようになるので、対象ファイルをアップロードする

ファイル読み込み

# 訓練データの読み込み
X_train = np.load('kmnist-train-imgs.npz')['arr_0']
y_train = np.load('kmnist-train-labels.npz')['arr_0']
# テストデータの読み込み
X_test = np.load('kmnist-test-imgs.npz')['arr_0']

# データのサイズ確認
print(X_train.shape)
print(X_test.shape)

アップロードしたファイルを読み込み、想定通りの件数になっているかを確認

データの確認

print(X_train.min())
print(X_train.max())

print(y_train[423])

plt.imshow(X_train[423], cmap = plt.cm.gray)
plt.show()

どういう中身かを確認

データ前処理

正規化

# 処理前データ確認
X_train[0,10:15,10:15]

# 訓練/テストデータの正規化
X_train = X_train / 255 
X_test = X_test/ 255

# 処理後データ確認
X_train[0,10:15,10:15]

リサイズ

# 処理前データ確認
print(X_train.shape)
print(X_test.shape)

X_train = X_train.reshape(X_train.shape[0], 784)
X_test = X_test.reshape(X_test.shape[0], 784)

# 処理後データ確認
print(X_train.shape)
print(X_test.shape)

LightGBM モデル訓練

# 訓練・テストデータの設定
train_data = lgb.Dataset(X_train, label=y_train)
eval_data = lgb.Dataset(X_test, label=y_test, reference= train_data)

ハイパーパラメータはとりあえず初期値としておく。

このタイミングでvalidation用のデータとラベルが必要なのだが、どちらも無いことに気が付く。
なので訓練データを分割して作成することにする。
件数的にはとりあえず3割ぐらいにしておく。

validation_size = 0.3
np.random.seed()
perm_idx = np.random.permutation(len(X_train))
_X_train = X_train[perm_idx]
_y_train = y_train[perm_idx]

# split train and validation
validation_num = int(len(y_train) * validation_size)

validation_imgs = _X_train[:validation_num]
validation_lbls = _y_train[:validation_num]

train_imgs = _X_train[validation_num:]
train_lbls = _y_train[validation_num:]

このデータを使って、改めてモデル設定

# 訓練・テストデータの設定
train_data = lgb.Dataset(train_imgs, label=train_lbls)
eval_data = lgb.Dataset(validation_imgs, label=validation_lbls, reference= train_data)

params = {
'task': 'train',
'boosting_type': 'gbdt',
'objective': 'multiclass',
'num_class': 10,
'verbose': 2,
}

boosting_typeのgbdtは「Gradient Boosting Decistion Tree」の略です。分類クラスは3つ以上ですので多項分類（multiclass）を目的（objective）に指定します。またターゲットクラスは10種類の平仮名ですのでnum_classは10となります。

学習

gbm = lgb.train(
params,
train_data,
valid_sets=eval_data,
num_boost_round=150,
verbose_eval=5,
)

サンプルでは100回訓練していたが、もちっと繰り返してもよさそうだったので150回にしてみる。

検証

y_pred = []
eval_predicts = gbm.predict(validation_imgs)
for x in eval_predicts:
    y_pred.append(np.argmax(x))

confusion_matrix(validation_lbls, y_pred)

accuracy_score(validation_lbls, y_pred)

0.9426111111111111

思ってた以上よりいい感じ。高すぎる気もするけど。

予測

_predicts = gbm.predict(X_test)
predicts

predicts = []
for x in _predicts:
    predicts.append(np.argmax(x))

ファイル作成

submit = pd.DataFrame(data={"ImageId": [], "Label": []})
submit.ImageId = list(range(1, _predicts.shape[0]+1))
submit.Label = predicts


submit.to_csv("submit.csv", index=False)
files.download('submit.csv')

これで予測した結果をDLできる。

提出してみた結果

0.862

こんなもんか。過学習しちゃっているような感じなのかなぁ。
もちっとパラメータをいじってみたり、XGBoostも試してみての差を見てみたりしたいけど、とりあえず連休はもうおしまいなのである。。

今回作ったものは以下
https://github.com/RyoNakamae/LightGBM

tips

502 Bad Gateway

HTTPのエラーコードの1つ．
400番台がクライアント側のエラー（File not found とか）であるのに対し，500番台はサーバ側のエラーを示す．
今回は，ゲートウェイやプロキシが不正なリクエストを受け取ったことを示している．

参考：HTTP502（Bad Gateway）の原因と対処法

SQLインジェクションの1つの解

管理者でのログインを試みる際，ユーザ名がadminで決め打ちだった場合，

admin'--

と入力すればパスワードが何であれログインに成功する．

調べた知識

Base64

英数字と「+」「-」「=」のみを使って，マルチバイト文字等をエンコードする方式．これによって，英数字+@しか使えない通信環境でもマルチバイト文字等を扱うことができる．
エンコード・デコードには以下のサイトが便利である．

オンラインBase64でデコーダ

ヴィジュネル暗号

換字式暗号の1つ．
アルファベットの対応表をもとに平文を暗号化する．
1つ1つ目で追って解読しても解けるが，スクリプトを書いて楽にしたい．

vigenere_cipher.py

import sys


class MyVigenere:
    def encode(text, key):
        ret = ''
        for i in range(len(text)):
            offset = (ord(text[i]) + ord(key[i]) - 2 * ord('a')) % 26
            ret += chr(ord('a') + offset)
        return ret

    def decode(text, key):
        ret = ''
        for i in range(len(text)):
            offset = (ord(text[i]) - ord(key[i])) % 26
            ret += chr(ord('a') + offset)
        return ret


if __name__ == '__main__':
    plaintext = sys.argv[1]
    key = sys.argv[2]
    flag = MyVigenere.decode(plaintext, key)
    print(flag)

John The Ripper

パスワードクラッキングツール．
以下，Ubuntuでの実行例を示す．

// ツールのインストール
# apt install john
// パスワードファイルとシャドウファイルを結合する
$ unshadow passwd shadow > tmp
// johnコマンドに渡す
$ john --show tmp

robots.txt

検索エンジンのクローラーのアクセスを制限するためのファイル．
見せたくないページはここに記述するらしい．

picoCTF

きなこ(@kinako_software)さんのブログ(こちら)を拝見して，cpawCTFの次はこれをやるといいよと書いてあったので取り組んでみる．
ヒントが書かれているので，もしわからなかったとしても調べればなんとかなりそう．

Logon

ユーザ名とパスワードを入力してログインする．
何を入力してもログインに成功するが，adminでログインしたい．
Cookieを見ると，ログイン成功後にadminという変数がFalseで格納されている事がわかる．
これをTrueに書き換えてリロードするとログイン成功．

Recovering From the Snap

ディスクイメージファイルが渡される．
fileで見てみると，FATであることがわかる．
Ubuntuで，以下のコマンドを実行してマウントする．

# mount -t vfat animals.dd /mnt/animal/

これで/mnt/animal/以下に4つの画像が現れた．
Linuxでは，jpgファイルはeogコマンドで開くことができる．
しかし，どの画像を見てもフラグに関係するようには思えない．
ヒントを見てみると，いくつかのファイルが削除されてしまっているらしい．
FTK Imagerというツールを使うと，ディスクイメージをマウントしたときに削除済みのファイルにもアクセスできるらしい．
MacやLinuxではコマンドラインツールしか利用できず使い方が難しかったので，Windows版をダウンロードした．
与えられたイメージを開くと，フラグが書かれている画像が削除されていることがわかる．
そのファイルを開けばフラグを獲得できる．

admin panel

与えられたpcapファイルをWiresharkで開く．
/loginにPOSTしている部分があるので，そのパケットを見ると，パスワードのところにフラグが書かれている．

hertz

換字暗号だったので，便利そうなツールを調べたところ，quipqiupというのを見つけた．
「Puzzle:」の部分に暗号化されたテキストを貼り付けて「Solve」を押すだけ．

今までチュートリアルのコンペに参加してみましたが、今回は初めてちゃんとしたコンペに参加してみます。今回はMercari Price Suggestion Challengeです。

前回の記事

https://qiita.com/HayatoYamaguchi/items/30e93a43d519f8c8aea0

Mercari Price Suggestion Challengeとは

なんか一年くらい前にすごい話題になったものらしく、カーネル（お手本の投稿みたいなやつ）もすごく多いのでこれをやってみた。データの特徴としては、特徴量はそんなに多くないが、とにかくデータが今までのと比べてものすごく多い。訓練データだけでも1482531件もあります、しかもその上にItem descriptionという特徴があり、これは言語なので自然言語処理をしないといけず、むちゃくちゃデータが重くて待つのが大変だった（）

前処理

今回はkaggleのカーネルを参考にやっていきました。
まずデータをsalesPrice（予測すべきもの）だけ落として訓練データとテストデータを連結させます。

path = "../input/"
train = pd.read_csv(path+"train.tsv", header=0, sep='\t')
test = pd.read_csv(path+"test.tsv", header=0, sep='\t')
train_df = train
test_df = test
drop_train_df = train_df.drop('price',axis=1)
allData = pd.concat([drop_train_df,test_df],axis=0)

特徴量は７つあるので一つずつ見ていきます

単純な処理の特徴量

brand_name
これは単純にブランドの名前です。ブランドの名前は大きく価格に影響するので、今回はこれの全てをLabelEncoderを用いて数値化していきます、あとでまとめてやります

category_name
Men/Tops/T-shirts、Women/Tops & Blouses/Blouseのように最大２つのスラッシュによって３つに分かれています。また、欠損値になっているものも多く存在するので今回はこれをスラッシュの最初と真ん中と最後の三つの要素に分解し、それぞれ新たな特徴量としていきます。

def split_cat(text):
    try: return text.split("/")
    except: return ("No Label", "No Label", "No Label")
allData['general_cat'], allData['subcat_1'], allData['subcat_2'] = zip(*allData['category_name'].apply(lambda x: split_cat(x)))
allData.head()

item_condition_id
商品の状態ごとのスコアです。１から４までありますがそれぞれがどういう意味なのかは分かりません笑

shipping
客と店どちらが送料を負担しているかというものです。０と１の値をとります。

ここまではやっていることはそこまで難しくありませんでしたが、次からが自然言語処理を行うのでやったことがない自分にとってはとても大変でした（汗）

自然言語処理を行う特徴量

item_description
その商品の説明が英語で書かれたものです。No description yetと書かれたものも多くあるほか、４つほど欠損データがテストにありました
これらは全て訓練データにあるので、データを落とします。

allData = allData.dropna(subset=['item_description'])

次に新しくis_descriptionという、descriptionがあるかどうかを判定する特徴量を作成します。

def if_description(row):
    if row == 'No description yet':
        a = 0
    else:
        a = 1
    return a

allData['is_description'] = allData["item_description"].apply(lambda row : if_description(row))

さて、ここからが自然言語処理の本番です。正直自分はまだまだ全然理解できていないですが、今回はTfidfVectorizerという手法を用いてdescriptionに存在する単語の出現回数を数えてそれをベクトル化して、Tfidfを算出します。この辺からむちゃくちゃ重くなりましたw
その後にSVDという手法を使います。これはPCAと類似したもので、次元削減をします。これにより、たくさんの単語をベクトル化した後、そのベクトルの次元削減をすることができます。
ちなみに時間が出てきますが、これはどれくらい時間がかかったかをプリントしてるだけで、特にそれ以上の意味はありません。
ちなみに自然言語処理はよくわからなかったので
https://aiacademy.jp
↑以下のサイトで無料でやってみました、今回やっているのとは内容は異なりましたがどんな感じでtokenizeするのかとかの感覚は分かりました。

from sklearn.feature_extraction.text import TfidfVectorizer

train_df = allData[allData['test_id'].isnull()]
test_df = allData[allData['train_id'].isnull()]

from sklearn.decomposition import TruncatedSVD
import time
start = time.time()
tfidf_vec = TfidfVectorizer(stop_words='english', ngram_range=(1,1))
full_tfidf = tfidf_vec.fit_transform(allData['item_description'].astype('U').values.tolist())
print("tfidfallData終わり")
train_tfidf = tfidf_vec.transform(train_df['item_description'].astype('U').values.tolist())
test_tfidf = tfidf_vec.transform(test_df['item_description'].astype('U').values.tolist())
print("tfidf終わり")
n_comp = 20
svd_obj = TruncatedSVD(n_components=n_comp, algorithm='arpack')
svd_obj.fit(full_tfidf)
print("svdのfit終わり")
train_svd = pd.DataFrame(svd_obj.transform(train_tfidf))
test_svd = pd.DataFrame(svd_obj.transform(test_tfidf))
print("svdのtransform終わり")
train_svd.columns = ['svd_item_'+str(i) for i in range(n_comp)]
test_svd.columns = ['svd_item_'+str(i) for i in range(n_comp)]
train_df = pd.concat([train_df, train_svd], axis=1)
test_df = pd.concat([test_df, test_svd], axis=1)
end = time.time()
print("time taken {}".format(end - start))

name
このネームも上のitem_descriptionと同じような処理をします。Tfidfを算出した後にSVDを算出します。
ちなみにどうして他の次元削減の方法やトピックモデリングを使わずに今回これを選択したかというと、いくつか試した中でこれが一番早く動いたからですw
他のだと１時間かかっても次元削減すら終わらないみたいなことが普通に起きました、、、

from sklearn.decomposition import TruncatedSVD
import time
start = time.time()
tfidf_vec = TfidfVectorizer(stop_words='english', ngram_range=(1,1))
full_tfidf = tfidf_vec.fit_transform(allData['name'].values.tolist())
print("tfidfallData終わり")
train_tfidf = tfidf_vec.transform(train_df['name'].values.astype('U').tolist())
test_tfidf = tfidf_vec.transform(test_df['name'].values.astype('U').tolist())
print("tfidf終わり")
n_comp = 20
svd_obj = TruncatedSVD(n_components=n_comp, algorithm='arpack')
svd_obj.fit(full_tfidf)
print("svdのfit終わり")
train_svd = pd.DataFrame(svd_obj.transform(train_tfidf))
test_svd = pd.DataFrame(svd_obj.transform(test_tfidf))
print("svdのtransform終わり")
train_svd.columns = ['svd_name_'+str(i) for i in range(n_comp)]
test_svd.columns = ['svd_name_'+str(i) for i in range(n_comp)]
train_df = pd.concat([train_df, train_svd], axis=1)
test_df = pd.concat([test_df, test_svd], axis=1)
end = time.time()
print("time taken {}".format(end - start))

その他の前処理

ここまでで難しい前処理は終わりです。数値化できる特徴を数値化して、数値化した特徴量以外の特徴量を消し、その上でデータを標準化します。
使わない特徴量をリストに入れて、それを削除します。

from sklearn.preprocessing import LabelEncoder
le = LabelEncoder()
allData['general_cat'] = le.fit_transform(list(allData['general_cat']))
allData['subcat_1'] = le.fit_transform(list(allData['subcat_1']))
allData['subcat_2'] = le.fit_transform(list(allData['subcat_2']))
allData['brand_name']=allData['brand_name'].fillna('missing').astype(str)
allData['brand_name'] = le.fit_transform(list(allData["brand_name"]))

do_not_use_for_training = ['test_id','train_id','name', 'category_name','item_description']
feature_names = [f for f in train_df.columns if f not in do_not_use_for_training]
X_train = train_df[feature_names]
X_test = test_df[feature_names]
X_test.head()

次に、まとめて標準化します。外れ値が多そうなので、正規化ではなく標準化を使用しました。ちなみにこの標準化の処理はRidgeモデルを使用した時には役に立ちましたが、今回はXGboostを使用したので全く必要ではありませんでした笑

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X_train)
columns = X_train.columns.values
X_train_scaled = pd.DataFrame(scaler.transform(X_train),columns=columns)
X_test_scaled = pd.DataFrame(scaler.transform(X_test),columns=columns)
X_train_scaled.head()

最後にPriceについてです。訓練データのPriceを見ると右に長い分布になっています。そのため、この形が綺麗な二項分布になるようにします。

y_train = np.log(train['price'].values + 1)

ここまでで前処理は終了です！
ここからモデルを設計していきます。

モデルの設計

本当はこのような大量で複雑なデータを使用するときはニューラルネットワークを使用するのが基本であると思われます。しかし、今回はあまりにもこのデータが重かったので単純に時間を取られるなと思い、他の方法によってモデルを作り上げました、他にもやり方調べるのがだるry

まず試したのがRidgeです、いくつかのカーネルでこれが使用されていたため使ってみました。これは30秒もあれば回せる代わりに、スコア予想が低くなりました。
次に試したのがXGboostingです。これは実際のデータを回そうとすると提出までに１時間かかりましたwこれをprediction用にlogなどを直して提出します。

import time
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_log_error
model = xgb.XGBRegressor(max_depth=15,colsample_bytree=0.5,eta=0.05,min_child_weight=20)
model.fit(X_train_scaler,y_train)
preds = model.predict(X_test_scaler)

test["price"] = np.expm1(preds)
test[["test_id", "price"]].to_csv("submission_ridge.csv", index = False)

ちなみにこのパラメータの設定は本当はグリッドサーチしないといけませんが、するとそれだけで多分10時間くらいパソコンを回さないといけないので、今回はカーネルにあったパラメータをそのままパクってきました。とりあえずこれで提出して、スコアは0.50315となり、上位40数パーセントくらいになりました。

最後に

多分もっとしっかりパラメータを調整したり、学習曲線書いたりすれば前処理で次元削減の数とかを調整したりすればもっといいスコアが出ると思いますが、一回回すだけでも１時間くらいかかってしまうのでここで今回は終わりにしますw
自然言語処理ってコードが重くて大変ですね、、、
今度は画像処理とかもやってみようと考えています。

ちなみに参考までに交差検証を行った時のそれぞれのスコアです
ちなみに評価指標はRMSLEです、よく分からない人はググってくださいw
特徴量20個２つでRidge
0.17040984745724053
0.16994236305495025
特徴量20個２つでXGboosting(30分かかる)
0.10012834069970207
0.12882562545649595
特徴量30個２つでRidge
0.16885835791622594
0.16846895351513075

if,else,elifによるデータの扱い

皆さんはお店のアルバイトをしている時にレジのバイトをしていると考えてください。それで130円のアイスを買いに来ているお客様がいるとします。その時もらったお金によって、「ちょうどお預かりいたします。」なのか、「お釣りのお返しです。」なのか「お金が足りません」と返すのかということをif文でそれぞれ処理をしてくれます。elifは条件２〜みたいなものです。

#もし支払った金額がちょうどだったら、、、ここは払ってもらっている金額をかく。
cost=130
#もらった金額が130円だった時の処理
if cost==130:
    print('ちょうどお預かりいたします。')
#もらった金額が130円より少ない時
elif cost < 130:
    print('お金が足りません。')
#もらった金額が130円より多い時
elif cost > 130:
    print('お釣りのお返しです。')
#それ以外（お金をもらってない時）
else:
    print('金額をおっしゃってください。')

このようにして条件文はできています。それでは少し応用のinput()を使ってみましょう。これは自分で値段を書いてその時に応じてかくというものです。
一回例を書いてみますね！！

while True:
    number=input('支払金額: ')
#ここで数字を文字列化していく
    cost=int(number)
#もらった金額が130円だった時の処理
    if cost==130:
        print('ちょうどお預かりいたします。')
#もらった金額が130円より少ない時
    elif cost < 130:
        print('お金が足りません。')
#もらった金額が130円より多い時
    elif cost > 130:
    print('お釣りのお返しです。')
#それ以外（お金をもらってない時）
    else:
        print('金額をおっしゃってください。')

このようにすれば自分が金額を指定してそれに応じた金額を出力することが可能です。ここで使っているのがwhile文になるwhile文はその条件があっている限り処理を繰り返していくというもの。例えば５以下の場合はずっと同じ処理を繰り返していくなどというもの。簡単な例を示していきます。

cost=2
while cost<11:
   print(cost*cost)
   cost+=1

4
9
16
25
36
49
64
81
100

ここでは数字の２から１０までの数字をcostに代入してそれを二乗するということです。

breakによるループ中止

やばい、渡すお釣りが無くなっちゃた！これはもうちょうど支払って欲しいなって思った時にこのように言いましょう！！！

while True:
    cost=input('支払金額: ')
    number=int(cost)
#もらった金額が130円だった時の処理
    if number==130:
        break
    print('ちょうどいただきたいです')

これだったらお釣り渡さなくてもいいですよね？？？

continueによる次のループの開始

continueは大雑把にいうとスルーみたいなものです！！ここでは釣りが出るとダメだよ！！みたいなことを言ってちょうどの時はもう処理を終了する。そしてお金が足りない時は何も言わないことにしよう

while True:
    cost=input('支払金額: ')
    number=int(cost)
    if number==130:
        break
    if number<130:
        continue
    print('ダメ！！！')

こんな感じでやってみればできると思います！！

もしこの表現ムカつくとかありましたら編集リクエストもしくはコメントにておしらせください！お願い申しあげます。GWは投稿休んでしまい申し訳ございませんでした

if,else,elifによるデータの扱い

皆さんはお店のアルバイトをしている時にレジのバイトをしていると考えてください。それで130円のアイスを買いに来ているお客様がいるとします。その時もらったお金によって、「ちょうどお預かりいたします。」なのか、「お釣りのお返しです。」なのか「お金が足りません」と返すのかということをif文でそれぞれ処理をしてくれます。elifは条件２〜みたいなものです。

#もし支払った金額がちょうどだったら、、、ここは払ってもらっている金額をかく。
cost=130
#もらった金額が130円だった時の処理
if cost==130:
    print('ちょうどお預かりいたします。')
#もらった金額が130円より少ない時
elif cost < 130:
    print('お金が足りません。')
#もらった金額が130円より多い時
elif cost > 130:
    print('お釣りのお返しです。')
#それ以外（お金をもらってない時）
else:
    print('金額をおっしゃってください。')

このようにして条件文はできています。それでは少し応用のinput()を使ってみましょう。これは自分で値段を書いてその時に応じてかくというものです。
一回例を書いてみますね！！

while True:
    number=input('支払金額: ')
#ここで数字を文字列化していく
    cost=int(number)
#もらった金額が130円だった時の処理
    if cost==130:
        print('ちょうどお預かりいたします。')
#もらった金額が130円より少ない時
    elif cost < 130:
        print('お金が足りません。')
#もらった金額が130円より多い時
    elif cost > 130:
    print('お釣りのお返しです。')
#それ以外（お金をもらってない時）
    else:
        print('金額をおっしゃってください。')

このようにすれば自分が金額を指定してそれに応じた金額を出力することが可能です。ここで使っているのがwhile文になるwhile文はその条件があっている限り処理を繰り返していくというもの。例えば５以下の場合はずっと同じ処理を繰り返していくなどというもの。簡単な例を示していきます。

cost=2
while cost<11:
   print(cost*cost)
   cost+=1

4
9
16
25
36
49
64
81
100

ここでは数字の２から１０までの数字をcostに代入してそれを二乗するということです。

breakによるループ中止

やばい、渡すお釣りが無くなっちゃた！これはもうちょうど支払って欲しいなって思った時にこのように言いましょう！！！

while True:
    cost=input('支払金額: ')
    number=int(cost)
#もらった金額が130円だった時の処理
    if number==130:
        break
    print('ちょうどいただきたいです')

これだったらお釣り渡さなくてもいいですよね？？？

continueによる次のループの開始

continueは大雑把にいうとスルーみたいなものです！！ここでは釣りが出るとダメだよ！！みたいなことを言ってちょうどの時はもう処理を終了する。そしてお金が足りない時は何も言わないことにしよう

while True:
    cost=input('支払金額: ')
    number=int(cost)
    if number==130:
        break
    if number<130:
        continue
    print('ダメ！！！')

こんな感じでやってみればできると思います！！

もしこの表現ムカつくとかありましたら編集リクエストもしくはコメントにておしらせください！お願い申しあげます。GWは投稿休んでしまい申し訳ございませんでした

はじめに

こんにちは。
僕はプログラミング歴3ヶ月の初心者なのですが、今回はそんな初心者の目から見た「覚えておくと便利なライブラリ」を紹介したいと思います！
すでにいくつか記事を書いてるので、よかったらそちらもついでに見てやってください。下のリンクの記事から繋げやすい内容になっています。
プログラミングしてみたい方はまずPythonのインストールを！(Mac向け)

ライブラリってなんだ？

プログラミング(この記事でいうとPython)初心者の方はまずライブラリという言葉を聞き慣れていないかもしれません。
ライブラリとはプログラミングをする上で欠かせないものであり、画像を編集したいときや、ファイルを整理するときなどに使えるいわば飛び道具のようなものです。すでに完成されているので僕たちは決められた書き方をするだけでその飛び道具が使えてしまうんです。
import 〜という決まり文句で呼び出すと使えるようになります。

ではそんな便利なライブラリをいくつか見ていきましょう！

Numpy

Numpyは数値計算をするときや、画像をピクセルごとに分解して配列にしたい時に使えるライブラリです。ディープラーニングで画像分類をする時などに必ず必要になってくる知識なので、AI関係のことに興味のある方は必ず覚えておきましょう！

こんな感じで使うのですが、今は全く分からなくて構いません。「へー、そんなものがあるのか」ぐらいの認識で結構です。
Numpyを使うにあたって学習（人によってはおさらい）しておくべきなのが行列です。高校の数学で習った方もいるとは思うのですが、指導要領の変更で外れている年もあるので、習ったことのない方はぜひ学習してみましょう。全然難しい内容じゃないですし、ググったら簡単に噛み砕いて教えてくれているサイトが山ほど出てきますよ！

Pandas

Pandasはデータ分析を行う上で必須のライブラリです。Pandasを使うとデータの読み込み(エクセルのデータなど）や不必要なデータの除外など、機械学習をする上では欠かせない作業が一気に簡単になります。機械学習に必要な作業の大部分がこのデータ分析になってくるので、覚えておくとめちゃくちゃ強いです！データサイエンティストになりたい方には特に必須ですよ。
Pandasの発展的な内容では株価分析などもできるようになるので、株に興味がある方にもおすすめです。プログラミングと組み合わせてトレードなんかできたらカッコよくないですか？（笑）

Matplotlib

MatplotlibはPythonの描画ライブラリで、線グラフや棒グラフなどデータを可視化したい時に使えるライブラリです。
例えばパワーポイントなどを使ってプレゼンをしないといけない場面でグラフを使いたい時、Matplotlibが大きな力を発揮します。データをグラフにしたい時、シンプルなやり方で緻密なグラフが描画できるのでとても便利ですよ！

これはMatplotlibの公式サイトから引用したグラフです。こんなグラフが数行のプログラミングで実現できるってすごくないですか！？
勉強し始めの頃はなんの意味もないグラフをかくだけでも構いません！「こういう使い方をするんだな」ぐらいの軽い気持ちで色々書いてみましょう。まずは自分で書くことが大事ですし、こういうグラフの形で目に見えるとシンプルに楽しいです。

Keras

Kerasは今までの3つとは少し毛色が違って、本格的に機械学習をする時に使うライブラリです。僕の今までの記事でも何度か出ている画像分類などで力を発揮します。
また、Kerasは公式ドキュメントに使い方が多く載っていたり、Kerasを使ったコードレシピもネット上にたくさん転がっているので何かを形にしたい時に即座に使いやすいライブラリです。
なんでこんなにサンプルが転がっているかというと、そのシンプルさが理由なんです。これまでにあったライブラリよりもはるかに少ない記述量で機械学習が簡単にできちゃう優れものであり、みんなが使うから必然的にネットに情報量も多いので勉強しやすいのがおすすめポイントです。
「エラーメッセージが出てきたけど何が書いてあるか分からない?」なんて時も同じエラーメッセージが出た人が掲載した解決法をググればなんの問題も無しです。安心して勉強しちゃいましょう！

僕が書いたこんな記事もありますので、ある程度勉強が進んだ方はこのやり方を真似て画像分類してみると面白いかもしれません。
VGG16で画像分類(2クラス)

その他のライブラリ・モジュール

上の4つを特におすすめしたかったので詳しく書きましたが、その他にも以下のようなライブラリがあるので気になった方はぜひググってみてください！今は使わないとしても結局使うことになるライブラリばかりだと思うので、覚えておいて損はないです。
・os(ファイル操作で使う）
・dlib(顔認証などで使える）
・OpenCV(画像の処理はこれさえあれば困りません！）
・sklearn
・sys
・glob

etc.....

こんな感じでたくさんあるのでぜひ調べてみてください！

参考にしたサイト

Numpy公式サイト
Matplotlib公式サイト

終わりに

いかがだったでしょうか？
この記事では書ききれない量の便利なライブラリがあるので、各自で調べてみるとさらに勉強の質が上がります！少しでも気になることはどんどんググって自分のものにしていきましょう！
また、この記事で書いたものはどれも僕の最初の記事で紹介したネット講座で勉強できるものばかりなのでそちらの記事も見ていただけたら嬉しいです！
プログラミングやディープラーニングの勉強、何から手をつければ良い！？

では、また次の記事で！

自分がRENATの開発者である。今回、ユーザの立場から、RENATの使い方を説明する 。
元々、RENATはネットワーク装置の検証の目的で作られた（参照資料）。今回、ソフトウェア・ルータを用い、renatの使い方を簡単に説明する。物理ルータや他のネットワーク装置へ対しての使い方は同じである。

仮想ルータはなんでもいいが、昔からBIRDを使おうと持っていたので、今回birdを使ってみる。

環境

• osX Mojave
• docker desktop Version 2.0.0.3
• renat 0.1.14

概要

以下の流れで説明する

• renatとは
• ネットワークトポロジー
• renatのインストール
• birdルータのインストール
• シナリオ作成・実行
• 実行結果の確認
• まとめ

RENATとは

詳細を省略するが、RENATはRobotframeworkを拡張したネットワーク検証用のフレームワークである。簡単なテキストのシナリオを用い、ラボで検証する内容をほぼ再現でき、実際の現場、サービスで既に多く使われている。

ネットワークトボロジー

作りたいトポロジーはとても簡単なもので、ソフトウェアルータ（bird）とrenatサーバが直結されるものである。

+-------+     +-------------+
| ルータ +-----+ renatサーバ  |
+-------+     +-------------+

birdルータのインストール

今回、cloudRouter(fedora/bird)を試す。以下の手順で一つのルータを立ち上げる。

root権限でcontainerを開始する

$ docker run --privileged --publish-all=true cloudrouter/bird-fedora /sbin/init

設定を変更するため、containerへloginする

$ docker exec -it 22f73220b705 /bin/bash --login

telnet-serverを追加し、rootでloginできるようにする

$ rm -f /etc/yum.repos.d/cloudrouter.repo
$ yum install -y telnet-server
$ sed -i '/pam_securetty.so/d' /etc/pam.d/remote
$ systemctl start telnet.socket

このcontainerのIP(eth0)を確認する。この場合は、ルータのIPが 172.17.0.2 である。

[root@22f73220b705 etc]# ip addr
1: lo: <LOOPBACK,UP,LOWER_UP> mtu 65536 qdisc noqueue state UNKNOWN group default qlen 1
    link/loopback 00:00:00:00:00:00 brd 00:00:00:00:00:00
    inet 127.0.0.1/8 scope host lo
       valid_lft forever preferred_lft forever
2: tunl0@NONE: <NOARP> mtu 1480 qdisc noop state DOWN group default qlen 1
    link/ipip 0.0.0.0 brd 0.0.0.0
3: ip6tnl0@NONE: <NOARP> mtu 1452 qdisc noop state DOWN group default qlen 1
    link/tunnel6 :: brd ::
34: eth0@if35: <BROADCAST,MULTICAST,UP,LOWER_UP> mtu 1500 qdisc noqueue state UP group default
    link/ether 02:42:ac:11:00:02 brd ff:ff:ff:ff:ff:ff
    inet 172.17.0.2/16 brd 172.17.255.255 scope global eth0
       valid_lft forever preferred_lft forever

rootアカウントのパスワードを system へ変更する

$ echo "cloudrouter" | passwd root --stdin

BIRDを開始する

$ bird

renatのインストール

以下のようにrenat環境を構築する。RENAT(github)のサイトでインストール方法を書いてあるが、ためすだけならば、dockerのイメージは最もは早い方法である。

dockerからインストール

RENAT(docker)の手順に従い、renatサーバを動かす（dockerhub.comから最新なイメージをダウンロードしてくる）

$ mkdir /opt/renat
$ docker run --rm -d --privileged -v /opt/renat:/opt/renat -p 80:80 -p 10022:22 --name renat bachng/renat:latest

ここで、renatサーバが動き始め、http port 80 で実行結果が確認できるようになる。

renat環境を調整

renatの設定ファイルは全て$RENAT_PATH/configのは以下で、以下ファイルである。なお、環境変数 $RENAT_PATHはデオフォルトで /home/work/robot/renatである。

• template.yaml: サポートしているplatformの情報（プロムプトなど）
• auth.yaml: 認証情報
• device.yaml: システムワイドの設備情報（主にIP情報）

元々、renatではJuniper,CiscoなどのNW装置をサポートしているが、birdのサポートは入っていない。以下のように、birdのサポートを追加する。

必要なのは、template.yamlファイルにBIRDのサポートを追加するのみである。

access-template:
    bird:
        access: telnet
        auth: plain-text
        profile: default
        prompt: ".*(>|#) "
        login-prompt: "login: "
        password-prompt: "Password: "
        init:
            - birdc
        finish:
            - quit
            - exit
    ssh-host:
...

ここで、ホストへのアクセスをtelnetとし、ホストにログイン成功したら、birdcを呼び、birdのCLIに入る。また、認証はplain-textで、認証情報のプロフィルはdefaultとする。

auth.yamlを編集し、BIRDルータの認証情報設定する。

auth:
    plain-text:
        default:
            user: root
            pass: cloudrouter
        samurai:
...

ここでは、BIRDルータのパスワードは前節でbird containerのrootパスワードと同じにする

device管理ファイルdevice.yamlに bird01 のルータを追加する。ここに追加されたdeviceはrenatサーバを共有したテストitemから利用できるようになる。

device:
    # Cloudrouter
    bird01:
        type: bird
        description: a cloudrouter BIRD version
        ip: 172.17.0.2
    # Testers
...

ここで IPアドレスは、birdコンテナのIPアドレスである。

シナリオ作成・実行

renatサーバへrobotユーザとしてログインする。デフォルトで /home/robot/work がワーキングフォルダーになる。

$ docker exec -it --user robot renat /bin/bash --login

まず、~/workは以下で、テストprojectを作成する。

[robot@9985510e9a87 work]$ cd /home/robot/work/
[robot@9985510e9a87 work]$ $RENAT_PATH/tools/project.sh renat-sample
created test project:  renat-sample
entered project folder:  renat-sample
use item.sh to create test case
[robot@9985510e9a87 work]$ cd ~/work/renat-sample
[robot@9985510e9a87 renat-sample]$

そして、batchモードでbird01というデバイスを含めたテストitemを作成する。

[robot@9985510e9a87 renat-sample]$ $RENAT_PATH/tools/item.sh -b -l -n bird01 test01


=== Created `test01` test item ===
Case scenario:     /home/robot/work/renat/tools/test01/main.robot
Case run file:     /home/robot/work/renat/tools/test01/run.sh
Local config file: /home/robot/work/renat/tools/test01/config/local.yaml
Tester config file:/home/robot/work/renat/tools/test01/config/
Check and change the `local.yaml` local config file if necessary
[robot@9985510e9a87 renat-sample]$

test01は以下にシナリオ・ファイル main.robot を以下のように編集する（既存のテンプレートは消しても良い）

 20
 21 *** Test Cases ***
 22 01. 基本的なコマンドの確認
 23     Router.Switch       bird01
 24     Router.Cmd          show status
 25     Router.Cmd          show memory
 26

まず、文法チェックを行う

[robot@9985510e9a87 test01]$ ./run.sh --dryrun
Current time:       Mon May  6 09:46:07 UTC 2019
Current RENAT path: /home/robot/work/renat

### Current folder is /home/robot/work/renat-sample/test01 ###
Run only once

Run: 001
Current local.yaml: /home/robot/work/renat-sample/test01/config/local.yaml
Loaded extra library `Hypervisor`
==============================================================================
test01 :: This is a sample test item
==============================================================================
01. 基本的なコマンドの確認                                            | PASS |
------------------------------------------------------------------------------
test01 :: This is a sample test item                                  | PASS |
1 critical test, 1 passed, 0 failed
1 test total, 1 passed, 0 failed
==============================================================================
Output:  /home/robot/work/renat-sample/test01/result/output.xml
Log:     /home/robot/work/renat-sample/test01/result/log.html
Report:  /home/robot/work/renat-sample/test01/result/report.html

全てPASSならば、次に実際シナリオを実行する（--dryrunなしで）

### Current folder is /home/robot/work/renat-sample/test01 ###
Run only once

Run: 001
Current local.yaml: /home/robot/work/renat-sample/test01/config/local.yaml
Loaded extra library `Hypervisor`
==============================================================================
test01 :: This is a sample test item
==============================================================================
RENAT Ver:: RENAT 0.1.14
------------------------------------------------------------------------------
README:
Write you readme file here


------------------------------------------------------------------------------
00. Lab Setup
------------------------------------------------------------------------------
01. 基本的なコマンドの確認                                            | PASS |
------------------------------------------------------------------------------
99. Lab Teardown
------------------------------------------------------------------------------
test01 :: This is a sample test item                                  | PASS |
1 critical test, 1 passed, 0 failed
1 test total, 1 passed, 0 failed
==============================================================================
Output:  /home/robot/work/renat-sample/test01/result/output.xml
Log:     /home/robot/work/renat-sample/test01/result/log.html
Report:  /home/robot/work/renat-sample/test01/result/report.html

実行結果の確認

前節でシナリオの実行が終わったら、dockerホストからのブラウザーでテストitemの結果を確認できる。renatの各ユーザのworkフォルダーが ~<username>へマッピングされる。今回、このurlで結果が確認できる http://127.0.0.1/~robot/renat-sample/test01/result/log.html

ここで、+のマークを選択すると、各ステップの詳細情報が見える。また、以下のurlで、書く装置の実行ログも確認できる
http://127.0.0.1/~robot/renat-sample/test01/result/bird01.log

まとめ

RENATを使い、birdルータの基本コマンドを確認テストitemの方法について簡単に説明した。シナリオを実行した例でしめしたように、ユーザがlogin/logoutやログ収集について一切気にせず、renatへ任せることができた。ユーザが、検証したい内容を注目し、シナリオを書き、テストを実行することによって、従来の検証プロセスをより効率的に実現できるだろう。

また、元々、birdをサポートしていないrenat対して、簡単にbirdのサポートを入れる方も説明した。今回、コンセプトを見せるため、テストの内容は基本のコマンドを確認するだけにしたが、renatを用い、他のNW装置、トラヒックジェネラータを含めた複雑なテスト・シナリオも作成、確認できる。また、今度別の記事を書きたいと思う。

ここでは簡単な使い方を説明したが、renatの詳細は以下のサイトへ。

参考資料

• RENAT(github)
• RENAT(dockerhub)
• RobotFrameworkについて（英語）
• RobotFrameworkのドキュメント（日本語）
• Janog#41で紹介したrenat
• How to configure CloudRouter containers as a Route Server

squre.py

import sys
import pygame
from pygame.locals import *

width, height = 800, 600

pygame.init()
screen = pygame.display.set_mode((width, height))
pygame.display.set_caption("square")
clock = pygame.time.Clock()

def main():
    while True:
        for event in pygame.event.get():
            if event.type == QUIT:
                pygame.quit()
                sys.exit()
        screen.fill((0, 0, 0))
        square = pygame.draw.rect(screen, (255, 0, 0), (width/2, height/2, 100, 100))
        pygame.display.update()
        clock.tick(30)

if __name__ == "__main__":
    main()

sys とは

sysはPythonのインタプリタや実行環境に関する情報を扱うためのライブラリです。
使用しているプラットフォームを調べる時や、
スクリプトの起動パラメータを取得する場合などに利用します。

インストール

標準ライブラリなので、はじめから入っています。
pipなどでインストールする必要はありません。

sys.version

使用しているPythonのバージョンを文字列型で。

> sys.version
'3.7.3 (default, Mar 27 2019, 23:40:30) \n[GCC 6.3.0 20170516]'

sys.platform

sys.platform

sys,platformはスクリプトを実行したOSを表示します。
私の場合はlinuxと表示されました。

例えば、カレントディレクトリ/home/user1に、test.pyとPython packageディレクトリtestがあるとします。

sys.argv

sys.argvはコマンドライン引数を扱います。
コマンドライン引数とは名前の通り、コマンドラインで指定した引数のことです。

Pythonでプログラムを作成し利用する場合、
条件に応じた値を引数で指定する場合があります。
その時にsys.argvが役に立ちます。

これのことです。

jovyan@bb26596c4f15:~$ python hoge.py koko kore

上記のプログラムではhoge.pyというPythonプログラムに対して、
kokoとkoreというコマンドライン引数が指定されています。
sys.argvを使用することで、これらの引数をPythonプログラム内で扱うことができます。

これからsys.argvを使用してコマンド引数の扱いを確認します。
はじめにhoge.pyというPythonファイルを下記のように作成します。

hoge.py

import sys

#↓ここでコマンドライン引数をargsという変数に格納
args = sys.argv

if len(args) == 3 :
    print(args) #コマンドライン引数を出力
    print(type(args)) #コマンドライン引数のデータ型を出力
    print(type(args[2])) #コマンドライン第一引数のデータ型
    print(len(args))
    print("-"*30) #区切ってるだけ。意味はないです。

    print("コマンドライン引数の第一引数は " + args[1] )
    print(float(args[2])*3.14) #floatの意味は下に書きます。

else:
    print("コマンドライン引数を２つ入力してください \n 第一引数はstr型、第二引数はfloat型です。")

そして実行してみると……

bash

jovyan@bb26596c4f15:~$ python hoge.py hello 10
['sys.py', 'hello', '10']
<class 'list'>
<class 'str'>
3
------------------------------
コマンドライン引数の第一引数は hello
31.400000000000002

出力の一行目から見ていきましょう。
argsはsys.argvによってコマンドライン引数を格納したものです。
リストの０番目はスクリプトのファイル名が入っています。
そして１番目以降はコマンドライン引数が順番に入っています。

出力の二行目はargsの、つまりsys.argvのデータ型を確認しています。
<class 'list'>と出力されているのでリスト型です。

出力の三行目はargs[2]、つまり２番目のコマンドライン引数のデータ型を確認しています。
２番目の引数は10ですが……<class 'str'>と出力されています。
数字を入れているので、整数型か浮動小数点数型になりそうな気がしますが、
コマンドライン引数は全て文字型になるようです。

出力の四行目はargsの長さを出力しています。
スクリプト名 + コマンドライン引数 なので3と出力されます。

sys.exit

このメソッドはSystemExit例外を発生させます。
この例外は特に処理をしない限り、スタックトレースを表示しません。
例えば、下記のようなpythonコードを書くと……

exit.py

import sys

for i in range(100):
  print(i)
  if i == 10:
    sys.exit()

jovyan@bb26596c4f15:~$ python exit.py 
0
1
2
3
4
5
6
7
8
9
10

このようにexit()到達時点でプログラムが終了します。

参考文献&参考サイト

今回の記事制作にあたり参考にさせていただいた記事や書籍です。
とても勉強になりました。ありがとうございました。

色々なPythonのライブラリを「知る」ための本
今回の記事を書いたきっかけ。今後もこの書籍を参考に記事を書きます。(Primeなら無料で読めます！)
LIFE WITH PYTHON 「ライブラリ：sys」
基本的な構成を参考にしました。とてもわかりやすいです。
python公式ドキュメント sys --- システムパラメータと関数
THE 本家。
Python学習講座　コマンドライン引数
sys.argvに関する記述が非常にわかりやすい。

bytes オブジェクトを生成できるらしい。

>>> b'\xe3\x81\x82'
b'\xe3\x81\x82'

>>> b'\xe3\x81\x82'.decode('UTF-8')
'あ'

>>> 'あ'.encode('UTF-8')
b'\xe3\x81\x82'


>>> b'A' == b'\x41'
True

>>> b'A' == b'A'
True

>>> type(b'A')
<class 'bytes'>

>>> type('A')
<class 'str'>

Original by Github issue

https://github.com/YumaInaura/YumaInaura/issues/1671

はじめに

• 「ガウス過程に基づく分類モデル」で行った計算を、R, Python, Juliaで実装して処理時間を比較。
• ついでにC#とC++でも実装して比較してみた。
• 事後分布をラプラス近似して、周辺尤度を最大化するパラメータを求める。
• 最適化手法は、勾配情報の不要なNelder-Meadアルゴリズムを採用。
• Windows 10 Pro, 16.0GBメモリ, i7-8550U CPU @ 1.80GHz
• The Julia Programming Language
• JuliaNLSolvers/Optim.jl
• SciPy
• Accord.NET Framework
• Math.NET Numerics
• Eigen: C++ template library for linear algebra

結果

評価関数の呼び出し回数を考慮すると、Juliaが最速で、その後にC++、Python、R、C#と続く。

言語環境	秒数	呼び出し回数	平均秒数	比率
Julia(1.1.0)	10.86	70	0.155	1
C++(VS2017)	8.51	47	0.181	1.17
Python(3.6.7/WSL)	13.06	59	0.221	1.43
R(3.5.3)	38.6	47	0.822	5.3
C#(VS2017)	47.6	47	1.013	6.5

Pythonで測定したところ、CPUの利用率が100%となり、24秒くらいかかった。
OMP_NUM_THREADS=1として使用するスレッド数を1に制限したら速くなった。
C#が遅すぎるが、Install-Package MathNet.Numerics.MKL.Win-x64として、
Intel Math Kernel Library (MKL)を使用すると劇的に速くなる。

言語環境	秒数	呼び出し回数	平均秒数	比率
C+++MKL	4.58	47	0.097	0.63
C#+MKL	6.59	47	0.140	0.90
R Open(3.5.2)+MKL	7.51	47	0.160	1.03
Python(3.7.3)+MKL	10.9	59	0.184	1.19

Microsoft R Openは、デフォルトでMKLを使用する。使用スレッド数を1に制限。
MinicondaでPythonをインストールするとnumpy+MKLを使用するため速くなる。
最速はC++/Eigen/MKL。
Julia+MKLは未評価。MKLを使用するためにはソースからコンパイルする必要あり。

処理の大部分は行列(767x767)のコレスキー分解であり、MKLの効果が大きい。

サンプルデータ

USPS handwritten digit data
http://www.gaussianprocess.org/gpml/data/

R

## データの入力

library("R.matlab")

mat_data <- readMat('usps_resampled.mat')

index3 = (mat_data$train.labels[3+1,] == 1)
index5 = (mat_data$train.labels[5+1,] == 1)

y_train = ifelse(index3,1,-1)[index3 | index5]
x_train = mat_data$train.patterns[,index3 | index5]

## カーネル関数とシグモイド関数

x_train_L2 = apply(x_train,2,function(a){colSums((x_train-a)^2)})
my_kernel_train = function(k_sigma,k_scale){
  exp(x_train_L2 / (-2 * k_scale^2)) * k_sigma^2
}

my_sigmoid = function(y,f){
  pnorm(y * f)
}
my_sigmoid_log_1st = function(y,f,prob){
  y * dnorm(f) / prob
}
my_sigmoid_log_2nd = function(y,f,prob){
  a = dnorm(f) / prob
  a * (a + y * f)
}

## 周辺尤度の最大化

my_evidence_optim = function(par){
  k_scale = exp(par[1])
  k_sigma = exp(par[2])
  K11 = my_kernel_train(k_sigma,k_scale)
  n = length(y_train)
  f = seq(from=0,to=1,length.out=n)
  while(TRUE){
    prob = my_sigmoid(y_train,f)
    V1 = my_sigmoid_log_1st(y_train,f,prob)
    V2 = my_sigmoid_log_2nd(y_train,f,prob)
    V2q = sqrt(pmax(V2,0))
    B = diag(n) + t(K11 * V2q) * V2q
    U = chol(B) # t(U) %*% U
    b = V2 * f + V1
    a = b - backsolve(U,forwardsolve(t(U),(K11 %*% b) * V2q)) * V2q
    f0 = f
    f = drop(K11 %*% a)
    if(max(abs(f - f0))<1e-5)break
  }
  -sum(a * f) / 2 + sum(log(prob)) - sum(log(diag(U)))
}
start = proc.time()
ret_optim = optim(c(2.85,2.35),my_evidence_optim,control=list(fnscale=-1))
proc.time() - start
ret_optim

Julia

Rとは異なり、whileでも変数のスコープが構成される。
sum()したら次元をドロップして欲しい。

using MAT
using Distributions
using Optim
using LinearAlgebra

## データの入力

matfile = matopen("usps_resampled.mat")

train_labels = read(matfile, "train_labels")
train_patterns = read(matfile, "train_patterns")

index3 = (train_labels[3+1,:] .== 1)
index5 = (train_labels[5+1,:] .== 1)

y_train = ifelse.(index3,1,-1)[index3 .| index5]
x_train = train_patterns[:,index3 .| index5]

## カーネル関数とシグモイド関数

x_train_L2 = mapslices(t -> dropdims(sum((x_train .- t) .^ 2, dims=1), dims=1), x_train, dims=1)
function my_kernel_train(k_sigma,k_scale)
  exp.(x_train_L2 / (-2 * k_scale^2)) * k_sigma^2
end

function my_sigmoid(y,f)
  cdf.(Normal(), y .* f)
end
function my_sigmoid_log_1st(y,f,prob)
  y .* pdf.(Normal(), f) ./ prob
end
function my_sigmoid_log_2nd(y,f,prob)
  a = pdf.(Normal(), f) ./ prob
  a .* (a .+ y .* f)
end

## 周辺尤度の最大化

function my_evidence_optim(par)
  k_scale = exp(par[1])
  k_sigma = exp(par[2])
  K11 = my_kernel_train(k_sigma,k_scale)
  n = length(y_train)
  f = Array(range(0,1,length=n))
  evidence = 0
  while true
    prob = my_sigmoid(y_train,f)
    V1 = my_sigmoid_log_1st(y_train,f,prob)
    V2 = my_sigmoid_log_2nd(y_train,f,prob)
    V2q = sqrt.(max.(V2,0))
    B = I + K11 .* V2q .* V2q'
    C = cholesky(Symmetric(B))
    b = V2 .* f + V1
    a = b - (C.U \ (C.L \ ((K11 * b) .* V2q))) .* V2q
    f0 = f
    f = K11 * a
    if maximum(abs.(f - f0)) < 1e-5
      evidence = -dot(a,f) / 2 + sum(log.(prob)) - log(det(C)) / 2
      break
    end
  end
  evidence
end

@time result = maximize(my_evidence_optim, [2.85,2.35], NelderMead())

Python

np.arrayとnp.matrixの二種類があり挙動が異なる。

import os
os.environ["OMP_NUM_THREADS"] = "1"

import time
import numpy as np
from scipy import io
from scipy.stats import norm
from scipy.linalg import cholesky
from scipy.linalg import solve_triangular
from scipy.optimize import minimize

## データの入力

matdata = io.loadmat('usps_resampled.mat')

train_labels = matdata["train_labels"]
train_patterns = matdata["train_patterns"]

index3 = (train_labels[3,:]==1)
index5 = (train_labels[5,:]==1)

y_train = np.where(index3,1,-1)[np.logical_or(index3,index5)]
x_train = train_patterns[:,np.logical_or(index3,index5)]

## カーネル関数とシグモイド関数

x_train_L2 = np.apply_along_axis(lambda x: np.sum(np.power(x_train-x.reshape(-1,1),2),0),0,x_train)
def my_kernel_train(k_sigma,k_scale):
  return np.exp(x_train_L2 / (-2 * k_scale * k_scale)) * (k_sigma * k_sigma)

def my_sigmoid(y,f):
  return norm.cdf(y * f)

def my_sigmoid_log_1st(y,f,prob):
  return y * norm.pdf(f) / prob

def my_sigmoid_log_2nd(y,f,prob):
  a = norm.pdf(f) / prob
  return a * (a + y * f)

## 周辺尤度の最大化

def my_evidence_optim(par):
  k_scale = np.exp(par[0])
  k_sigma = np.exp(par[1])
  K11 = my_kernel_train(k_sigma,k_scale)
  n = len(y_train)
  f = np.linspace(0,1,n)
  evidence = 0
  while True:
    prob = my_sigmoid(y_train,f)
    V1 = my_sigmoid_log_1st(y_train,f,prob)
    V2 = my_sigmoid_log_2nd(y_train,f,prob)
    V2q = np.sqrt(np.maximum(V2,0))
    B = np.identity(n) + K11 * V2q * V2q.reshape(-1,1)
    U = cholesky(B)
    b = V2 * f + V1
    a = b - solve_triangular(U,solve_triangular(U,np.dot(K11,b) * V2q,trans=1)) * V2q
    f0 = f
    f = np.dot(K11,a)
    if np.max(np.abs(f - f0)) < 1e-5:
      evidence = -np.dot(a,f) / 2 + np.sum(np.log(prob)) - np.sum(np.log(np.diag(U)))
      break
  return -evidence

start = time.time()
par = np.array([2.85,2.35])
ret = minimize(my_evidence_optim,par,method='Nelder-Mead')
print(time.time() - start)
print(ret)

C#

.matファイルの入力にはAccord.Mathを使用。
行列計算と最適化にはMathNet.Numericsを使用。
データの抽出と結合に少し手間がかかる。

class BinaryHandwrittenDigit
{
    (Vector<double>, Matrix<double>) Select(Int16[,] labels, Double[,] patterns)
    {
        int nrow = patterns.GetLength(0);
        int ncol = patterns.GetLength(1);
        int count_3 = 0;
        int count_5 = 0;
        for (int j = 0; j < ncol; j++)
        {
            if (labels[3, j] == 1) count_3++;
            else if (labels[5, j] == 1) count_5++;
        }
        Vector<double> y = new DenseVector(count_3 + count_5);
        Matrix<double> x = new DenseMatrix(nrow, count_3 + count_5);
        int i = 0;
        for (int j = 0; j < ncol; j++)
        {
            if (labels[3, j] == 1)
            {
                for (int k = 0; k < nrow; k++) x[k, i] = patterns[k, j];
                y[i++] = 1;
            }
            else if (labels[5, j] == 1)
            {
                for (int k = 0; k < nrow; k++) x[k, i] = patterns[k, j];
                y[i++] = -1;
            }
        }
        return (y, x);
    }
    Vector<double> y_train;
    Matrix<double> x_train;
    void ReadMat()
    {
        string path_to_file = "usps_resampled.mat";
        var reader = new MatReader(path_to_file);
        var train_labels = reader.Read<Int16[,]>("train_labels");
        var train_patterns = reader.Read<Double[,]>("train_patterns");
        (y_train, x_train) = Select(train_labels, train_patterns);
    }
    Matrix<double> x_train_L2;
    Matrix<double> MakeMatrixL2(Matrix<double> x)
    {
        int m = x.ColumnCount;
        var result = new DenseMatrix(m, m);
        for (int j = 0; j < m; j++)
        {
            var cj = x.Column(j);
            for (int i = 0; i < m; i++)
            {
                var s = x.Column(i) - cj;
                result[i, j] = s * s;
            }
        }
        return result;
    }
    Matrix<double> kernel_train(double k_sigma, double k_scale)
    {
        return (x_train_L2 / (-2 * k_scale * k_scale)).PointwiseExp() * (k_sigma * k_sigma);
    }
    Vector<double> sigmoid(Vector<double> y, Vector<double> f)
    {
        return y.PointwiseMultiply(f).Map((x) => Normal.CDF(0, 1, x));
    }
    Vector<double> sigmoid_log_1st(Vector<double> y, Vector<double> f, Vector<double> prob)
    {
        var df = f.Map((x) => Normal.PDF(0, 1, x));
        return y.PointwiseDivide(prob).PointwiseMultiply(df);
    }
    Vector<double> sigmoid_log_2nd(Vector<double> y, Vector<double> f, Vector<double> prob)
    {
        var df = f.Map((x) => Normal.PDF(0, 1, x));
        var a = df.PointwiseDivide(prob);
        return a.PointwiseMultiply(a + y.PointwiseMultiply(f));
    }
    Matrix<double> bXbPlusOne(Matrix<double> X, Vector<double> b)
    {
        int n = X.RowCount;
        var result = X.Clone();
        for (int j = 0; j < n; j++)
        {
            for (int i = 0; i < n; i++)
            {
                result[i, j] *= b[i] * b[j];
            }
            result[j, j] += 1;
        }
        return result;
    }
    double GetEvidence(Vector<double> par)
    {
        double k_scale = Math.Exp(par[0]);
        double k_sigma = Math.Exp(par[1]);
        var K11 = kernel_train(k_sigma, k_scale);
        int n = y_train.Count;
        Vector<double> f = new DenseVector(n);
        for (int i = 0; i < n; i++)
        {
            f[i] = 1.0 / (n - 1) * i;
        }
        double evidence = 0;
        while (true)
        {
            var prob = sigmoid(y_train, f);
            var V1 = sigmoid_log_1st(y_train, f, prob);
            var V2 = sigmoid_log_2nd(y_train, f, prob);
            var V2q = V2.PointwiseMaximum(0).PointwiseSqrt();
            var B = bXbPlusOne(K11, V2q);
            var b = V2.PointwiseMultiply(f) + V1;
            var c = K11.Multiply(b).PointwiseMultiply(V2q);
            var chol = B.Cholesky();
            var a = b - chol.Solve(c).PointwiseMultiply(V2q);
            var f0 = f;
            f = K11 * a;
            if ((f - f0).AbsoluteMaximum() < 1e-5)
            {
                evidence = -(a * f) / 2 + prob.PointwiseLog().Sum() - Math.Log(chol.Determinant) / 2;
                break;
            }
        }
        return -evidence;
    }
    public void MarginalLikelihoodMaximization()
    {
        ReadMat();
        x_train_L2 = MakeMatrixL2(x_train);
        var evidence_optim = ObjectiveFunction.Value(GetEvidence);
        Stopwatch sw = Stopwatch.StartNew();
        var ret = NelderMeadSimplex.Minimum(evidence_optim, new DenseVector(new[] { 2.85, 2.35 }));
        sw.Stop();
        Console.WriteLine("{0}sec", (double)(sw.ElapsedMilliseconds) / 1000);
        var par_m = ret.MinimizingPoint;
        Console.WriteLine("par={0},{1}", par_m[0], par_m[1]);
        Console.WriteLine("Value={0},Reason={1},Iterations={2}", ret.FunctionInfoAtMinimum.Value, ret.ReasonForExit, ret.Iterations);
    }
}
static void Main(string[] args)
{
    Control.MaxDegreeOfParallelism = 1;
    var sample = new BinaryHandwrittenDigit();
    sample.MarginalLikelihoodMaximization();
}

C++

データは、C#で出力したものを入力。省略。
最適化関数nmminは、Rのソース/src/appl/optim.cを流用して作成。省略。
Eigenでautoを使用すると型が変わってしまう。

MKLを使用する場合は、#define EIGEN_USE_BLASとして、
mkl_intel_lp64.lib,mkl_sequential.lib,mkl_core.libをリンクする。
Intel® Math Kernel Library Link Line Advisor

#include <stdio.h>
#include <direct.h>

#include <chrono>
#include <cmath>

//#define EIGEN_USE_BLAS
#include "Eigen/Dense"

using namespace std;
using namespace Eigen;

typedef double optimfn(int n, const double *x0, void *ex);

void nmmin(int n, double *Bvec, double *X, double *Fmin, optimfn fminfn,
    int *fail, double abstol, double intol, void *ex,
    double alpha, double bet, double gamm, int trace,
    int *fncount, int maxit);

double Normal_PDF(double value)
{
    const double pi = 3.1415926535897932384626433832795;
    return exp(-value * value / 2) / sqrt(2 * pi);
}
double Normal_CDF(double value)
{
    return erfc(-value / sqrt(2)) / 2;
}

class BinaryHandwrittenDigit {
private:
    int nrow;
    int ncol;
    MatrixXd x_train;
    ArrayXd y_train;
    MatrixXd x_train_L2;
public:
    void read_data();
    void make_matrix_L2() {
        x_train_L2.resize(ncol, ncol);
        for (int j = 0; j < ncol; j++) {
            for (int i = 0; i < ncol; i++) {
                double t = 0;
                for (int k = 0; k < nrow; k++) t += (x_train(k, i) - x_train(k, j)) * (x_train(k, i) - x_train(k, j));
                x_train_L2(i, j) = t;
            }
        }
    }
    MatrixXd kernel_train(double k_sigma, double k_scale) {
        return (x_train_L2 / (-2 * k_scale * k_scale)).array().exp() * (k_sigma * k_sigma);
    }
    ArrayXd sigmoid(const ArrayXd& y, const ArrayXd& f) {
        ArrayXd result(y.size());
        for (int j = 0; j < y.size(); j++) {
            result(j) = Normal_CDF(y(j) * f(j));
        }
        return result;
    }
    ArrayXd sigmoid_log_1st(const ArrayXd& y, const ArrayXd& f, const ArrayXd& prob) {
        ArrayXd result(y.size());
        for (int j = 0; j < y.size(); j++) {
            result(j) = y(j) * Normal_PDF(f(j)) / prob(j);
        }
        return result;
    }
    ArrayXd sigmoid_log_2nd(const ArrayXd& y, const ArrayXd& f, const ArrayXd& prob) {
        ArrayXd result(y.size());
        for (int j = 0; j < y.size(); j++) {
            double a = Normal_PDF(f(j)) / prob(j);
            result(j) = a * (a + y(j) * f(j));
        }
        return result;
    }
    MatrixXd bXbPlusOne(const MatrixXd& X, const ArrayXd& b)
    {
        int n = X.rows();
        MatrixXd result(X);
        for (int j = 0; j < n; j++)
        {
            for (int i = 0; i < n; i++)
            {
                result(i, j) *= b[i] * b[j];
            }
            result(j, j) += 1;
        }
        return result;
    }
    double get_evidence(const double* par) {
        double k_scale = exp(par[0]);
        double k_sigma = exp(par[1]);
        MatrixXd K11 = kernel_train(k_sigma, k_scale);
        int n = y_train.size();
        ArrayXd f = ArrayXd::LinSpaced(n, 0, 1);
        double evidence = 0;
        while (true) {
            ArrayXd prob = sigmoid(y_train, f);
            ArrayXd V1 = sigmoid_log_1st(y_train, f, prob);
            ArrayXd V2 = sigmoid_log_2nd(y_train, f, prob);
            ArrayXd V2q = V2.cwiseMax(0.0).cwiseSqrt();
            MatrixXd B = bXbPlusOne(K11, V2q);
            ArrayXd b = V2 * f + V1;
            ArrayXd c = (K11 * VectorXd(b)).array() * V2q;
            LLT<MatrixXd> chol(B);
            ArrayXd a = b - chol.solve(VectorXd(c)).array() * V2q;
            ArrayXd f0 = f;
            f = K11 * VectorXd(a);
            if ((f - f0).abs().maxCoeff() < 1e-5) {
                evidence = -(a * f).sum() / 2 + prob.log().sum() - log(chol.matrixL().determinant());
                break;
            }
        }
        return -evidence;
    }
public:
    static double get_evidence_optim(int n, const double *x, void *ex) {
        auto p = (BinaryHandwrittenDigit*)ex;
        return p->get_evidence(x);
    }
};

int main()
{
    auto sample = new BinaryHandwrittenDigit();
    sample->read_data();
    sample->make_matrix_L2();

    double ninf = -INFINITY;
    double eps = 2.220446e-16;

    ArrayXd xini(2); xini << 2.85, 2.35;
    ArrayXd xout(2);
    double fmin = 0.0;
    int fail = 0;
    double abstol = ninf;
    double intol = sqrt(eps);
    void *ex = sample;
    double alpha = 1.0;
    double bet = 0.5;
    double gamm = 2.0;
    int trace = 0;
    int fncount = 0;
    int maxit = 500;

    auto start = chrono::system_clock::now();
    nmmin(xini.size(), xini.data(), xout.data(), &fmin, BinaryHandwrittenDigit::get_evidence_optim, &fail, abstol, intol, ex, alpha, bet, gamm, trace, &fncount, maxit);
    auto end = chrono::system_clock::now();
    long elapsed = chrono::duration_cast<chrono::milliseconds>(end - start).count();
    printf("fmin=%f, fncount=%d, fail=%d\n", fmin, fncount, fail);
    for (int i = 0; i < xout.size(); i++) printf("%d, %f\n", i, xout[i]);
    printf("%d msec\n", elapsed);

    delete sample;
}