任意関数の最小二乗法は EXCEL でも Pythonでも可能です。
精度、速度、記載の単純さで scipy.optimize.curve_fit の圧勝です。初期値によっては一発で求まらないこともありますが、 EXCELソルバーに比べると安定かつ高速で解を求めてくれます。変分推論は難しいですね。収束させるのにあたりのつけ方が必要になるようです。で、結局は最小二乗法に戻ってしまったり。慣れが必要なのでしょう。
Pandas の DataFrame でよく使う便利な乗算。
間違えやすいので、( ..)φメモメモ
**************************************
import numpy as np
import pandas as pd
df1=pd.DataFrame(np.ones(9).reshape(3, 3),
columns=['a','b','c'],
index=[1,2,3])
df1
| a | b | c | |
|---|---|---|---|
| 1 | 1.0 | 1.0 | 1.0 |
| 2 | 1.0 | 1.0 | 1.0 |
| 3 | 1.0 | 1.0 | 1.0 |
df2=pd.DataFrame(np.arange(9).reshape(3, 3),
columns=['a','b','c'],
index=[1,2,3])
df2
| a | b | c | |
|---|---|---|---|
| 1 | 0 | 1 | 2 |
| 2 | 3 | 4 | 5 |
| 3 | 6 | 7 | 8 |
#アダマール積
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | 0.0 | 1.0 | 2.0 |
| 2 | 3.0 | 4.0 | 5.0 |
| 3 | 6.0 | 7.0 | 8.0 |
index, column の名前、型が一致していないとダメ
df2.columns = ['d','e','f']
df1*df2
| a | b | c | d | e | f | |
|---|---|---|---|---|---|---|
| 1 | NaN | NaN | NaN | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN | NaN | NaN | NaN |
df2.columns = ['a','b','c']
df2.index = [4,5,6]
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN |
| 4 | NaN | NaN | NaN |
| 5 | NaN | NaN | NaN |
| 6 | NaN | NaN | NaN |
#indexを文字列として入力
df2.index = ['1','2','3']
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN |
| 1 | NaN | NaN | NaN |
| 2 | NaN | NaN | NaN |
| 3 | NaN | NaN | NaN |
#indexを数値として入力
df2.index = [1,2,3]
df1*df2
| a | b | c | |
|---|---|---|---|
| 1 | 0.0 | 1.0 | 2.0 |
| 2 | 3.0 | 4.0 | 5.0 |
| 3 | 6.0 | 7.0 | 8.0 |
dataframe * 行、列
#df1の各行に、df2の2行目をかける
df1*df2.iloc[1]
| a | b | c | |
|---|---|---|---|
| 1 | 3.0 | 4.0 | 5.0 |
| 2 | 3.0 | 4.0 | 5.0 |
| 3 | 3.0 | 4.0 | 5.0 |
# df1の各列に、df2のc列をかける・・・よりスマートな方法がありそうですが。
(df1.T*df2['c']).T
| a | b | c | |
|---|---|---|---|
| 1 | 2.0 | 2.0 | 2.0 |
| 2 | 5.0 | 5.0 | 5.0 |
| 3 | 8.0 | 8.0 | 8.0 |
df1*df2.iloc[1,2]
| a | b | c | |
|---|---|---|---|
| 1 | 5.0 | 5.0 | 5.0 |
| 2 | 5.0 | 5.0 | 5.0 |
| 3 | 5.0 | 5.0 | 5.0 |
ようやくMCMC を概ね理解し動かせるようになったので、地下水モデル10章の不確実性に関する部分へ少しだけ立ち戻ることに。
https://phreeqc.blogspot.com/2019/07/4.html
MCMC関連の文献を探してみると、簡単な計算をされている事例がありました。第2著者以降は日本の方ですね。
Julien Boulange, Hirozumi Watanabe, Shinpei Akai (2017) A Markov Chain Monte Carlo technique for parameter estimation and inference in pesticide fate and transport modeling
尤度を計算する前に採択されたパラメータでシミュを回す必要があります。大きなモデルを扱う実務向きではないですね。
実務で扱う地下水モデルで不確実性を検討するには、PEST一択なのでしょうか？確かに、市販のパッケージソフトには組み込まれており、使い勝手も良かった印象はあります。加えて、逆解析を古くから扱っているせいか、Bayes も良く出てきます。海外では当たり前なのかな？
もう少し、文献等を読んでみましょう。
Ohki et al.(2020)
Landslide detection in mountainous forest areas using polarimetry and interferometric coherence
https://earth-planets-space.springeropen.com/articles/10.1186/s40623-020-01191-5
ALOS-2 のデータ （Lバンド、空間解像度約６m)を用いた、土砂移動箇所検出の検討事例です。北海道胆振東部地震（2018年9月）と九州北部豪雨（2017年7月）を対象としています。
提案手法は change detection の範疇ですが、画像データを利用していません。波のデータのみです。そして最後はお決りの機械学習です（といっても決定木です）。
特徴
結果
留意点
最近の国交省さんの動向を見ていますと、今後、土砂災害分野における各種衛星データの利用は拡大して行くようですね。関連マニュアルも複数出ていますし、ALOS-3へ向けてのシンポジウムでも以下のような報告がされています。
https://www.pco-prime.com/alos-3_alos-4_sympo/pdf/1505.pdf
機械学習を用いた change detection の精度は高いといえない状況ですが、いずれ向上し、災害時の人手不足問題も（change detection に関しては）回避できるでしょう。
検出アルゴリズムの構築以外は特に難しい作業はないので、備えておきましょう。
EXCELをよく使っていた頃、データに多項式をフィッティングする場面で何次式を使うかは感覚的でした。
低次だと合いませんし、高次だとノイズまで合わせてしまいます（いわゆる過学習）。ちょうどよいところを探す方法は、、、古くからありましたね。
LASSO
https://scikit-learn.org/stable/modules/linear_model.html#lasso
\begin{align*}
\min_{w}{ \frac{1}{2n_{\text{samples}}} ||y - X w||_2 ^ 2 + \alpha ||w||_1}\end{align*}
この損失関数はモデルの適合性を表す誤差項と、複雑さを表すペナルティー項で構成されています。αは求めたいパラメータwを制御する上位のパラメータ（ハイパーパラメータ）。αを大きくするほど大きなL1ノルムに対してペナルティーがかかる仕組みです。するとL1ノルムは小さくならざるを得なくなり、スパースな解が導かれる（が適合性は低くなりやすい）という仕組み。逆にαを小さくするほどL1ノルムを大きく取れるようになり、適合性が向上しやすい（が複雑で、過学習に陥りやすい）。αをどのように調整・採用するかがキモです。
その判断基準は以下の通り。
・情報量基準（AIC、BICなど）
・cross‐validation
スパースを仮定した場合、非ゼロのwを増やすと組み込まれるXも高次まで多く含まれて適合性は良くなる（が複雑になる）といった表現が理想です。が、L0ノルムを扱うと組み合わせ（計算コスト）が膨大になるとのこと。そこでα＊L1ノルムで代用しているそうです。
ちなみに、L1、L2 は機械学習でも多用されていましたね。回帰だと以下の通りです。
・L0正則化・・・非ゼロノルムの個数
・L1正則化・・・L1ノルム利用、LASSO回帰
・L2正則化・・・L2ノルム利用、Ridge回帰
・L1＋L2正則化・・・Elastic Net
scikit-learn を利用すれば、これらのモデリングは容易。ありがたいですね。
==========================================
ラグランジュの未定乗数法を適用した場合の図や説明も web 上には多くあります。そちらの方が直感的で理解し易いでしょう。
関連内容はコチラ↓
https://phreeqc.blogspot.com/2019/06/blog-post_19.html
https://phreeqc.blogspot.com/2020/02/blog-post.html
SAR データを利用した Change Detection を整理。
以前利用していたDocker は －vオプションを使用しなければ動きました。機能が追加されているようです。
https://phreeqc.blogspot.com/2019/04/sar-change-detection-on-google-earth.html
GEE コミュニティのチュートリアルにもありました。Part１は動きましたが、Part２の後半はダメでした。
https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2
いずれも、GEE＋ Sentinel-1。解析エンジンもデータも海外製です。データは政府事業の Tellus 提供より揃っています。
https://phreeqc.blogspot.com/2019/02/google-earth-engine.html
国内の複数個所で Change Detection を試行しましたが、動作も問題ありません。きちんと結果は出てきます。
扱っているモノがモノなので、まだ精度面では物足りません。が、手軽に加工できるところはありがたいですね。将来的にはアルゴリズムが向上し、精度面でも耐え得るモノになるでしょう。
Tellus でも ALOS2 のデータが公開され始めました。が、未だ直近1年分くらいは公開されていません。過去分も、２～3年前くらいまで。直近のデータがないのは辛い。
しかも、開発環境を割り当ててもらえないと、自由に加工できない状況は継続中。https://phreeqc.blogspot.com/2020/09/tellusar.html
まだまだ海外製には追いつけない感が漂っています。頑張って欲しいですね。
MCMCを実装するライブラリはいくつかあるようです。
PyMC3:Theano
PyMC4:TensorFlow
Pyro:PyTorch
NumPyro:JAX
手元の環境で動いたのはPyroのみ。環境構築に関してはシビアなようです。
web上のサンプルを見てみますと、書き方はほぼPyMC3と同じでした。動かすのは簡単（traceplot がないのは残念。seaborn で書きました）。
が、写経し終わると理解できていない箇所が露わに。
y_model = pyro.sample('y_model', dist.Normal(a*x + b, sd), obs=y)
ここの obs の計算内容がわかりません。ここがキモなのはわかっているのですが、マニュアルに詳細が書かれていませんし、ソースも追えません。尤度×事前確率にどのようにかかわっているのでしょう？
おそらく、理解している方々が必要として使われているので、学び方が逆なのでしょうね。
ま、他の部分は理解できていることが分かりました。メモしておきましょう。
**************************************
20201123 MCMC メモ
データ駆動科学
取得データD[x,y]
予測分布モデル （例）物理モデル+ノイズ
物理モデル y=ax+b
ノイズ
予測分布モデルから θ[a:変形係数、b：切片、sd:標準偏差]の事後確率 p(θ|D)を求める。
マルコフ連鎖モンテカルロ法（MCMC法）
他手法との違い