やり残し事項 2019

今年は研究職を半分担当しました。

内業が主体となり、時間的にも体力的にも楽になりました。報酬の対象となる勤務時間だけを見ると確実に短縮されています。働き方改革と銘打って時短を推奨されている御時世ですが、素直な後輩たちの動向を見ていると将来への危機感を抱かざるを得ません。

人のことを言えないのですが、やり残し事項は2年連続ほぼそのまま。
かわりではないですが、3Dレーザー関連を扱えるようになりました。python を利用したデータ分析、機械学習関連のテクニックもさらに身に付けることができました。


以下、やり残し事項です。昨年とほぼ同内容ですが、一応書いておきましょう。

道具
帯磁率計
ガンマ線測定器（携帯型）

技術
動的解析（耐震、液状化）
斜面設計

コード・ソフト
優先度低：DtransuのCUDA化・・・100万で購入可。
優先度低：PSIｎSAR・・・700万で購入可
優先度中：地表流＋地下水・・・GSFLOW、HydroGeoSphere
優先度中：GMS・・・J-SHISの取り込みまではできました。いつか続きを。
優先度中：OpenMI・・・PHREEQCとDtransuの連成等


中期目標3年目はもがき中。あと2年。間に合うかな。

的中率、見逃し率

砂防分野では「的中率」や「見逃し率」などの用語をよく聞きます。

独自に進化したのか、機械学習分野と同じ定義でも異なる用語が用いられています。毎回、変換に戸惑います。

先日、以下の論文に砂防分野での定義の記載がありました。書き留めておきましょう。
堤ほか「ストリームチューブによる地形分割を基にした表層崩壊解析手法」2019, 砂防学会誌, 72, 2, 343, 3-13

TP:A
FN:B
FP:C
TN:D

的中率=(A+D)/（A+B+C+D）
捕捉率=A/(A+B)
空振り率=C/(C+D)
見逃し率=B/(A+B)

正答率：accuracy rate = 的中率
再現率:recallのうち、感度:sensitivity = 捕捉率
　　　　　　　　　特異度：pecificity = 空振り率
適合率：precisionの利用で、1-陰性的中率=見逃し率

imbalanced data のままだと、accuracy は機能しません。個々の評価指標を押さえておく必要があるのですが、重視すべきもののみ砂防分野で再定義されているということでしょう。

***********************************
20200304追記
「的中率」よりは「適中率」をよく見ますね（私は後者を使っています）。
また、上記の文献と気象庁では、一部の用語の定義が異なっています。ややこしや。
https://www.data.jma.go.jp/fcd/yoho/kensho/explanation.html

		予 　　報
		降水あり	降水なし	計
実況	降水あり	A	B	N1 = (A+B)
	降水なし	C	D	N2 = (C+D)
	計	M1 = (A+C)	M2 = (B+D)	N = (A+B+C+D)

• 適中率（降水の有無の適中率）（％）＝（A＋D）÷N×100
• 「降水あり」予報の適中率（％）　　＝A÷M1×100　（一致率に同じ）
• 「降水なし」予報の適中率（％）　　＝D÷M2×100
• 見逃し率（％）　　　　　　　　　　＝B÷N×100 
• 空振り率（％）　　　　　　　　　　＝C÷N×100
• 捕捉率（％）　　　　　　　　　　　＝A÷N1×100
• 一致率（％）　　　　　　　　　　　＝A÷M1×100　（「降水あり」予報の適中率に同じ）

機械学習モデルの解釈可能性

機械学習モデルの解釈可能性に関する、素晴らしいスライドと kaggle マイクロコースです。スライドの１、マイクロコースの２が昨日の特徴量評価にあたります。

スライド
https://speakerdeck.com/line_developers/machine-learning-and-interpretability

1. どの特徴量が重要か
　Feature Importance
2. 各特徴量が予測にどう影響するか
　Partial Dependence
　Surrogate Model
3. ある予測結果に対して特徴量がどう寄与するか
　LIME, SHAP
　Grad-CAM, Grad-CAM++

kaggle
https://www.kaggle.com/learn/machine-learning-explainability

1. Use Cases for Model Insights
　Why and when do you need insights?
2. Permutation Importance
　What features does your model think are important?
3.Partial Plots
　How does each feature affect your predictions?
4.SHAP Values
　Understand individual predictions
5.Advanced Uses of SHAP Values
　Aggregate SHAP values for even more detailed model insights


両者に出てくるSHAP、発想は面白いですね。単純な数式でモデルの推定結果を模擬し、各特徴量にかかる係数で定量的に結果への寄与度を表現しています。これは決定木ベース以外の説明できなかったモデルに対する解釈性の要望を反映したものでしょう。

が、正解を模擬しないのはなぜでしょう。あくまでモデルの解釈性を加えることが目的だからでしょうか。それとも理解できていないのでしょうか？
精度を犠牲にしつつ解釈性を加えた簡易モデルの追加。そこまでして「ブラックボックス」を避ける必要はないと思います。最後に判断するのはヒトなので。

このあたり、機械学習分野の今後の動向が注目されます。


********************************
20200104追記
技術評論社「Kaggleで勝つデータ分析の技術」
6.2特徴量選択および特徴量の重要度
6.2.2 特徴量の重要度を用いる方法
・モデルから出力される重要度：RF、GBDT
・他の方法：Permutation Importance、null importance、boruta、特徴量の機械的大量生成

特徴量選択

機械学習の特徴量選択について。

機械学習にかける前に有効な特徴量を選択するのですが、その手法は以下のように区分されています。https://www.oreilly.co.jp/books/9784873117980/

1. filter method、単変量統計: 特徴量個々の分離性・相関係数・分散などを利用。
2. wrapper method、反復特徴量選択: 再帰的特徴量削減(Recursive Feature Elimination, RFE)、Boruta など。feature_importance 利用。計算コスト高。

決定木ベースの手法による feature_importance を利用すれば、モデル作成と重要度の評価が同時に実施できます。３つ目の手法としてモデルベース特徴量選択と呼ぶそうです。

分類問題では、決定木ベースのモデルに依存した選択手法が主体でしょう。ライブラリも提供されており、環境は整っています。

逆に、決定木ベース以外のモデルでは feature_importance 属性を利用できませんので、RFE はダメ、モデルベース特徴量選択も当然ダメです（XGBoost や LightGBM 等のモデルが分類問題に対し高精度に機能することを示してきた中で、あえて他の手法に積極的に移行する動機はないのですが）。

決定木ベース以外では、モデルに依存しない指標を利用するようですね。例えば、Permutation importance（組もうかと思ったら、wrapper がありました）。

これは素人でも容易に思いつく手法なのですが、高コスト。完走させるために軽い手法を選択しても結果が良くないと、その中での指標の低下量、重要度を比較するだけになります（木を見て森を見ず状態）。

いずれにしても、何らかの評価・解釈を行って重要度の低い特徴量を破棄し、モデル全体の計算コスト低下、予測精度向上に努める必要があります。

特徴量の評価・選択は前処理に続き必要な作業だと、あらためて認識させられました。

********************************
20200104追記
技術評論社「Kaggleで勝つデータ分析の技術」でも同じ3区分でした。
6.2特徴量選択および特徴量の重要度
6.2.1 単変量統計を用いる方法
6.2.2 特徴量の重要度を用いる方法
6.2.3 反復して探索する方法

手法に関しては以下の記載がありました。
「なお、scikit-learn のドキュメントにはモデルの選び方として図4.8がありますが、データが十分にあるテーブルデータの回帰や分類のタスクでは、GBDT でうまくいくことが多く、このように細かに条件で分岐させる必要はあまりないでしょう。」p230

SVM

延長戦突入。

今日はサポートベクターマシン（support vector machine, SVM）。SVM はデータ数が少ない場合に用いられる手法で、昔よく利用されていました。が、今回は大量のデータを処理する必要があります。力技でこなそうと、scikit-Learn で組んでみました。が、やはり遅い。

SVCではなく、SVRだとどうか？と試しましたが、終わりません。

GPU版がないかな？と思い RAPIDS の中を覗いてみると、cuML の中に SVC がありました。が、fit()を使うと GTX1060 3GB のメモリに載りません（想定内）。scikit-Learn の cross validation に放り投げてみるも連携しません。

デモデータだと計算速度が40倍も違ったので、GPU は利用したいところです。

RTX8000 48GB なら動くでしょう。汎用性はなくなりますが、このマシンのみで動かすことを前提に、組んでみましょうか。

********************************
20200104追記
結局、他のライブラリとの連携の関係で、scikit-Learn を選択。年末から調整を続けていますが、最終モデルは40時間以上経過した今でもまだ完走しません。遅すぎ。特徴量を減らさないとダメということはわかりました。
技術評論社「Kaggleで勝つデータ分析の技術」では以下の記載がありました。SVM の使われない理由が明記された書き物を見たのは初めてです。

本章では、テーブルデータを扱う分析コンペにおいて使われるモデルとして以下を紹介します。・・・・・・ なお、教師あり学習のモデルとして良く紹介されるモデルのうち、サポートベクターマシンは精度や計算速度が見劣りするため、あまり使われることはありません。

********************************
20200109追記
O'REILLY「Pythonで始める機械学習」にも、以下のような記述がありました。

• サンプル数が大きくなるとうまく機能しない（10万サンプルぐらいまではうまく機能する）
• 10万サンプルぐらいになると、実行時間やメモリ使用量の面で難しくなる。
• 注意深い前処理とパラメータ調整が必要。→決定木ベースが良く用いられる理由の一つ。
• 検証が難しい。予測された理由を理解・説明することが困難。

 Ｃ、gamma の影響についてはｐ98に図解あり。

不飽和透水試験

不飽和透水試験機のデモに立ち会いました。

試験機自体は古くからあるのですが、小型化、高速化したと宣伝されてた機種です。小型化はまだしも、機器の改良だけで高速化？という点が気になっていました。

実際に見せていただくと、測定自体は確かに速い。が、不飽和透水試験で時間がかかるのは測定前の浸透。今回は肝心のプレ浸透が済まされており、試験だけの営業でした。
何時間浸透させたのかと聞いてみると、ほとんどしていないとのこと。つまり、測定が速い=高速化と宣伝されているようでした。うーん。

ま、プレ浸透を長くすれば飽和透水係数を得られるわけでもありません。が、雑に扱うのは論外。短時間のプレ浸透で飽和とみなせる理屈はないでしょう。
https://phreeqc.blogspot.com/2019/09/blog-post_17.html

残念ながら私の勘違いだったのですが、コンセプトは良いと思われます。これからも開発が続くことでしょう。一つぐらいは手元に置いておきたい機材でした。

正弦波

正弦波
=================================

振幅u0、角速度ω のとき、

u(t)=u0sin(ωt)

これはx=0での変位。
x=x位置での変位は、ちょっとしたアイデアで求まる。

x=0からx=xまで波が到達するのにかかる時間はx/c秒。
(x/c秒前の原点の状態を求めれば良い）

u(t,x)=u0sin(ω(t-x/c))

T=2πr/rω=2π/ω
λ=c/f
f=1/T
k=2π/λ　より

u(t,x)=u0sin(ωt-2πx/Tfλ)
　     =u0sin(ωt-2πx/λ)
　     =u0sin(ωt-kx)

RIVeR

PIV、STIV で河川の表面流速を比較しました。

以前、Python で組んでいたPIV＋STIV。これに動画をよませるだよませるだけです。
https://phreeqc.blogspot.com/2019/10/piv-stiv.html
https://phreeqc.blogspot.com/2019/10/stiv.html

今回は流速を実際に計測していますので、答え合わせができます。
PIV は静かな流れだと流速を掴みづらいようです。パラメータを調整しましたが、目立つ浮遊物がない限り、実際よりも低速と判定されます（というよりも認識しません）。
STIVは安定しています。浮遊物がなくても安定して流速を把握できます。
両者とも正しい答えを出すのですが、安定性で STIV 優位でしょうか。結果的には土研さんや国交省さんと同じ選択となりました。
http://www.jsece.or.jp/event/conf/abstract/2014/pdf/P2-56.pdf

表面流速が得られても、流量を出すにはひと手間必要です。複数のラインで平均流速を出したのち、断面形状と位置をあわせてから、それぞれの流速が受け持つ範囲を決めて流量を出す必要があります。このあたり、1から組むのは面倒。
ソフトがないかな？と探してみると、ありました。またしてもUSGSが噛んでいるのでしょうか。

RIVeR
https://www.sciencedirect.com/science/article/pii/S0098300417307045?via%3Dihub#!
http://riverdischarge.blogspot.com/

使ってみたところ、優秀。計算はもちろん、静止画への変換から幾何補正まで一通りの機能が備わっています。Unshake と言って、画像のブレを補正してくれる機能まであるので、スマホ撮影でも対応できそうです。しかも、PIV 版と STIV 版の両方が備わっています。
残念ながら、PIV で一部の流速を得られないのは同じ。これはソフトに起因するものではなく、手法に起因するものでしょう。一方、STIV では安定した結果を得られます。良いですね。

PIV で指定横断に直行方向の流速が得られる点もgood。Python では u,v 成分から指定横断位置での直交方向の流速を出す必要がありました。これを自動で算出してくれます。

残る課題は高さの把握。
洪水時等に幾何補正を為せるだけの高さ情報が必要になります。これを効率的に得るためにはどのようにすればよいかアイデアが必要があるでしょう。この点では、水路のような形状だと楽なのですが。

ひとまず、理屈を理解し、ソフトを準備しました。
残る課題をクリアし、本番に備えましょう。

ネットワークカメラ

ここ２か月で、３台のネットワークカメラを設置。

カメラの世界、進んでいますね。

http のみならず様々なプロトコルに対応しています。容易にデータを入手・転載できるほか、警報メールなども出せるようになっています。当然、発報にはトリガーが必要で、音の異常、動きの検出などを判定するためのアルゴリズムが入っています。精度は別として、技術レベルとしては高くないですからね。
https://phreeqc.blogspot.com/2019/10/blog-post_27.html

ハード面では、夜も見えるように近赤外領域を撮れるカメラも普及しているようです。併設するライトまで「近赤外LED」です。これは、自動車にも使われていますので、カメラの世界では浸透しているのでしょう。
高解像度化や省電力化も進んでいるようで、ライトを合わせて数Wの消費で済むように工夫されているものもありました。

知らない間に、ハードやソフトが変化しています。
さらに未来の状況は容易に想像できます。それを考えると、ここで脱落するわけには参りません。
発展させる立場にないことはやや残念ですが、開発者に敬意をもちつつ、積極的に利用させていただきましょう。

Fourier Transform in Python

先週末から触っていたフーリエ変換。結局、Python を使いました。

有限複素フーリエ係数Ck。Nで割る前の値が出ます。

Ck=np.fft.rfft(data)/N

周波数。ここまではいつも通り。

freq=np.fft.rfftfreq(N,d=1./fs)

躓いたのがランニングスペクトル（スペクトログラムという方が正解でしょうか）。scipyでの計算結果が EXCEL の計算結果と一致せず、最初は何の値が吐き出されているのかわかりませんでした。

from scipy import signal

freqs, length, Sx = signal.spectrogram(data, 

                                       fs=fs, 

                                       nperseg=segment,

                                       noverlap=segment-1,

                                   #raw-data

                                       window=('tukey', 0.), 

                                   # Remove linear trend along axis from data.

                                        detrend = False,

                                  # power spectral density(V**2/Hz)

                                  # or power spectrum (V**2) 

　　　　　　　　　　    scaling='spectrum', 

　　　　　　　　　   #mode='magnitude'
                                    )

ｗindowを削り、傾き補正を削り、スペクトルを指定すると、結果が一致。デフォルトではこれらが指定されていたため、明記しないと外れませんでした。

得られる値は Power。継続時間 T=ndt は乗じられていません。

この T や N が乗じられているか否かは式を見ないとわからないのですが、明記されていない場合があります（今回もそうでした）。

プロに聞けば、それぞれの流儀があるそうです。地震分野でも「私は乗じない」と言われる方がいらっしゃるそうで、決まりはないとのこと。ま、Nの場合は逆算する際に気を付けておけば良いだけですので、あまり気にする必要はないのでしょう。

EXCELでのミスも見つかりました。Power の最初と最後は2倍しない、など。数式を見直すと、確かにそのようになっていました。

初心者は手を動かさないと正解にたどり着けません。EXCEL で整理し、Python で答え合わせをしました。これで次は大丈夫でしょう。

Fourier Transform in EXCEL

この週末、EXCEL でのフーリエ変換を試していました。

これまで FFT は Python を使用しており、仕事でも Python を使うつもりでいたのですが、今回は EXCEL。理由は「後輩君が EXCEL で始めたから」。

私は空間、彼は時間を扱っており、別の作業です。が、DFT 部分は共通。彼が EXCEL で作り始めた（が正解にたどり着けていない）ので、今回はそれをフォローできる体制を整えておくことにしました。

EXCELでは、以下の3つの方法が考えられます。
1．分析ツールの「フーリエ解析」を使用する。
2．関数を使用する。
3．VBAで組む。

まず、分析ツールのフーリエ解析。試したところ、一番手軽でした。
当然、後輩君はこれを使用しています。が、 FFT なのでデータ数は2の累乗に限定されます。たとえパディングでかわしたとしても、採用されている計算式が明記されていないため、結果をデータ数で割るべきかどうかわかりません。故に振幅も出せません。
（VBA での結果と比較すると、データ数で割っていない値が表示されていました。）

関数では複素数を扱うのがやや面倒です。四則演算のみでもIM○○といった関数を使う必要に迫られます。さらに、k次の答えを求めるためには、k列必要です。これでは非現実的。

VBA では複素有限フーリエ変換を組もうとしました。が、一旦躓きました。VBAでも複素数を扱いづらい。WorksheetFunction でもエラーが出るので原因究明しかけました。が、結論としてはその時間をとらずに有限フーリエ変換に変更しました。今回はデータ数が700以下と少ないので、FFTでない方が効率的なのかもしれません。
ここまで来れば簡単。k、mの2重ループで有限フーリエ係数を出して、それから振幅、フーリエ振幅、Power、複素振幅を計算します。

計算後、スペクトルを図化する際にどの振幅を使うか迷いました。
調べてみると、フーリエ振幅ですね。今まで、スペクトル表示の際はただの振幅を利用しているものと思っていましたが、基本はフーリエ振幅を使っているようです（ランニングスペクトルではPowerが一般的）。慣習のようですね。分野にもよるのでしょう。


EXCEL では VBA が楽でした。分析ツールの出す答えが何を示すのかも分かったので、今後はコチラも使えます。
彼が週明けに自力で正解にたどり着いていたなら、それに越したことはありません。どうなっているでしょうか。
ま、頭の整理ができましたので、手を動かして良かったと思います。

*******************
20191217
結局、私は Python を利用しました。これが楽でした。

TREND-POINT

地上レーザーのデータを可視化しています。

サーフェス化して断面を切っていますが、同時に航空LPのデータも切っています。大体、同じ形状になります。両者の精度、思ったより良いですね。

レーザーのデータから地表を抽出するには、木や草を除去する必要があります。また、座標の変換も必要な場合があるでしょう。これ、ReCap ではできません。

プロに聞いたら、ライセンスに空きのあるソフトがありました。
TREND-POINT
https://const.fukuicompu.co.jp/products/trendpoint/

Leica のデータも読めました。
が、樹木等を除去するのはかなり困難。プロもトライ＆エラーで解決するそうです。航測会社さんの処理はとても速いのですが、きっと独自のノウハウがあるのでしょう。
これはあきらめました。

座標変換は簡単。
図面から公共座標を拾って点群に与えるだけ。それを計算させると独立座標から変換してくれます。この機能だけでも使えそうです。

プロに聞いたところ、点群を平面図に変換できるソフトはないようです。が、いずれできるでしょうね。フルサーフェス化したり、ソリッド化したり。

新たなハードが出ると新たなデータが生まれ、それを効率よく処理するソフトが整備されます。昔からそうですが、近年は浸透が速く、多様化しているように思えます。
頑張って、ついて行きましょう。