地下水 と MCMC その2

Using Markov Chain Monte Carlo to quantify parameter uncertainty and its effect on predictions of a groundwater flow model - ScienceDirect 

数千回もの浸透流計算が必要となるパラメータ推定手法は、実務向きではありません。https://phreeqc.blogspot.com/2020/12/mcmc.html
この文献では、既存の知見により浸透流計算を省略しています。具体的には、尤度関数内の二乗誤差、誤差分散計算の代替えとして、涵養量と透水係数の比 R/K との相関性を利用します。

overcome this obstacle, 300 realizations of Hassan et al.’s(2002) model were studied and it was found that the terms needed in the likelihood function, namely the sum of squared errors and the error variance are highly correlated to the recharge–conductivity ratio, R/K.

この文献の前に、300 realizations の計算があるのでしょう。それに気づいて、途中で細かく読むのをやめてしまいました。R/K などの相関関係は他の現場で使えたとしても、やはり事前の計算による確認が必要なのでしょうね。

おっ！と思ったり、うーん、と思ったりした文献でした。

 

やり残し事項 2020

GEE コミュニティのチュートリアルを再度実施してみようと手を動かしました。
https://phreeqc.blogspot.com/2020/12/gee-sentinel-1.html

1か月でPart3まで増えていました。仮説検定を用いた判定がメイン。今回は任意の場所でも１～3まで走りました。が、結果はイマイチ。
これまで、IｎSAR、Change Detection ともに、良い結果を得られたことがありません。イマイチなんですよね。やはり、有償ソフトの出番でしょうか。この辺りは来年の短期目標にしましょう。

職種が変わったこともあり、地質関連の作業がなくなりました。帯磁率計などはやり残したままで終わりました。

あらためてやり残し事項を整理しておきましょう。

優先度高：粒子法コード のGPGPU対応
優先度高：SAR・画像の Change Detection
優先度低：PSIｎSAR・・・700万で購入可
優先度低：DtransuのGPU対応・・・100万で購入可
優先度低：地表流＋地下水・・・GSFLOW、HydroGeoSphere
優先度中：OpenMI・・・PHREEQCとDtransuの連成等

中期目標４年目は惨敗。新たな目標を考えましょう。

DataFrames.jl vs. pandas

 Python は私のようなノンプロが手を出しやすい簡素な言語です。

一方、ライブラリへの依存が高いため、それらの互換性や version に応じて複数の仮想環境の構築を迫られることがあります。将来的には、作成していたコードがライブラリのVer.UP や開発中止で修正を余儀なくされることもあるでしょう（前者は既に経験済みです）。
個人的には 精度や速度が必須となる科学技術計算 を Fortran 、データ整理や可視化を Python、といったような使い分けをしています。現状維持で良いともいえるのですが、そろそろ代替え言語を試しておくべきと思い始めました（リスク管理です）。

先日、ある方から Julia を勧められました。Python と同じくらい容易な書き方で、計算が速く、コンパイルができるとのこと。
それは良い。

ということで、試してみました。
少し古いのですが、LTS の v1.0.5 + win10 （i5のラップトップ）です。インストールはプレコンパイル版を入れるだけ。そのままでは Jupyter での利用は不可でしたが、何度か試行し、動くようになりました。
今回の比較は csv 読み込みからの DataFrame 操作。頻繁に使う機能です（Juria と Python の比較ではなく、DataFrames.jl と pandas の比較になりました）。
結果は以下の通り。

csv読み込み	df=DataFrame(CSV.File(open(read,	1170万行	84.8s
	df6=DataFrame(CSV.File(open(read,	6万行	0.50s
	df=pd.read_csv	1170万行	51.4s
	df6=pd.read_csv	6万行	0.14s
カラム名変更	rename!(df,		0.21s
	df.rename(columns= ,inplace=True)		0.02s
inner_join	dfj=innerjoin(df, df6	290万行	15.5s
	dfj=pd.merge(df, df6	290万行	5.8s
型変換	dfj.＿=Date.(dfj.＿,"yyyy/mm/dd")		86.7s
	dfj['__']=pd.to_datetime(dfj['__'])		1.4s
期間抽出	df0=dfj[dfj[:,"__"] .< _‗, :]	190万行	4.6s
	df0=dfj[dfj['＿’] < ‗_]	190万行	3.6s
集約	df0_st=stack(df0,	380万行	5.6s
	df0_pv=pd.pivot_table(df0,values=	190万行	1.4s

pandas の圧勝でしたね。強い。

Julia でもパッケージを使うので Python と同じような問題は潜在的に有しているのでしょう。
イタレーションでは Julia の方が速くなるそうですので、扱う作業内容により言語を選択する形になるでしょうか？（今回の結果を見る限り、当面の出番はないかもしれませんが）
https://julialang.org/benchmarks/ 

引き続きその動向に注目しましょう。

常微分方程式

振動などの物理現象を解く際には、常微分方程式の解法が利用されています。

が、既に学んだことすら記憶にありません。
冬休み中に再読。以下にメモ。学ぶべき時に励んでおけばよかった。

*******************************

未知関数 y(x) と x および導関数 y'(x),y"(x),・・・yn(x) からなる方程式を n階常微分方程式という。
F(x,y(x),y'(x),y"(x),・・・yn(x))=0
このとき、y(x) を常微分方程式の解という。（方程式を解く＝式を求める）

常微分方程式の分類
・線形・非線形
・定数係数・変数係数
・同次・非同次

登坂宣好「微分方程式の解法と応用」より常微分方程式の例
線形‐変数係数‐非同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+x^2y\left(x\right)=e^x\end{align*}
線形‐定数係数‐非同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+3y\left(x\right)=e^x\end{align*}
線形‐変数係数‐同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+x^2y\left(x\right)=0\end{align*}
線形‐定数係数‐同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+3y\left(x\right)=0\end{align*}
非線形‐非同次形
\begin{align*}y\left(x\right)\frac{dy}{dx}\left(x\right)+3y\left(x\right)=x^3\end{align*}

解法（最近、使用した解法のみ記載）
線形-1階-定数係数-非同次形
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=f\left(x\right)\end{align*}
直接積分形：両辺をxで積分（両辺にdxを乗ずる形にして積分)
\begin{align*}f\left(x\right)=x : y\left(x\right)=\int xdx=\frac{1}{2}x^2+c\end{align*}\begin{align*}f\left(x\right)=sinx：y\left(x\right)=\int sinxdx=-cosx+c\end{align*}
線形-1階-定数係数-同次形
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=ay\left(x\right)\end{align*}
指数関数解:
\begin{align*}\frac{d}{dx}e^{\lambda x}=\lambda e^{\lambda x}=ae^{\lambda x}\end{align*}\begin{align*}\left(\lambda-a\right)e^{\lambda x}=0\rightarrow\lambda=a\end{align*}
解として\begin{align*}e^{ax}\end{align*}を有する。定数倍した\begin{align*}Ce^{ax}\end{align*}も解となる。（一般解）

線形-n階-定数係数-同次形
\begin{align*}\frac{d^ny}{dx^n}\left(x\right)+a_1\frac{d^{n-1}y}{dx^{n-1}}\left(x\right)+\ldots+a_ny\left(x\right)=0\end{align*}左辺に\begin{align*}y=C_1y_1+C_2y_2+\ldots+C_ny_n\end{align*}を代入すると\begin{align*}\frac{d^n}{dx^n}\left(C_1y_1+C_2y_2+\ldots\right)+a_1\frac{d^{n-1}}{dx^{n-1}}\left(C_1y_1+C_2y_2+\ldots\right)+a_n\left(C_1y_1+C_2y_2+\ldots\right)\end{align*}\begin{align*}{=C}_1\left({\frac{d^n}{dx^n}y}_1+a_1\frac{d^{n-1}}{dx^{n-1}}y_1+\ldots+a_ny_1\right)+C_2\left({\frac{d^n}{dx^n}y}_2+a_1\frac{d^{n-1}}{dx^{n-1}}y_2+\ldots+a_ny_2\right)+\ldots\end{align*}\begin{align*}{=C}_1・0+C_2・0+\ldots C_n・0=0\end{align*}よって\begin{align*}y=C_1y_1+C_2y_2+\ldots+C_ny_n\end{align*}は一般解となる。
解の形が\begin{align*}y=Ce^{\lambda x}\end{align*}のとき、
\begin{align*}C\ \left(\lambda^n+a_1\lambda^{n-1}+\ldots+a_n\right)e^{\lambda x}=0\end{align*}
C＝0のとき、y＝0は自明な解。\begin{align*}\lambda^n+a_1\lambda^{n-1}+a_2\lambda^{n-2}+\ldots+a_n=0 ：特性方程式\end{align*}のとき、非自明な解\begin{align*}y=C_1e^{α_1ｘ}+C_2e^{α_2ｘ}+⋯+C_ne^{α_nｘ}\end{align*}が存在。

非線形（特別な形式）
変数分離型解法：両辺をg(y)で割ってから積分(g(y)≠0)
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=f\left(x\right)g(y\left(x\right))\end{align*}\begin{align*}\frac{1}{g\left(y\left(x\right)\right)}dy=f\left(x\right)dx\end{align*}\begin{align*}\int\frac{1}{g\left(y\left(x\right)\right)}dy=\int f\left(x\right)dx\end{align*}

Mann-Whitney's U test

Mann-Whitney's U test を実施。

scipy.stats.mannwhitneyu を利用していたのですが、ランキングのつけ方が書かれていません。Σ の計算に n(n+1)/2 が含まれていますので、１位、２位、２位、４位のような同順位のつけ方はまずい。そこは正しく計算してくれているのだろうと思いながら、初めて使うので確認。

ソースでは rankdata が使用されています。

ranked = rankdata(np.concatenate((x, y)))

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rankdata.html?highlight=rank#scipy.stats.rankdata

>>> from scipy.stats import rankdata
>>> rankdata([0, 2, 3, 2])
array([ 1. ,  2.5,  4. ,  2.5])
>>> rankdata([0, 2, 3, 2], method='min')
array([ 1,  2,  4,  2])
>>> rankdata([0, 2, 3, 2], method='max')
array([ 1,  3,  4,  3])
>>> rankdata([0, 2, 3, 2], method='dense')
array([ 1,  2,  3,  2])
>>> rankdata([0, 2, 3, 2], method='ordinal')
array([ 1,  2,  4,  3]) 

はい、一番上の平均 (default is ‘average’) が利用されていました。OKです。

比較する母集団の平均、標準偏差、サンプル数の違いはどのように影響するのでしょうか？簡単な正規分布のプラス側を取り出してチェックしてみました（プロットは正規化済み）。シード値は設けず、何度かテストを繰り返して両側・有意水準5％にかかわる傾向を見てみました。

x1=np.random.normal(1, 10, 1000)
x2=np.random.normal(1, 10, 1000)
x1=x1[x1 > 0]
x2=x2[x2 > 0]   

当然ですが、棄却されない分布です（稀に棄却されるのが気になりますが)。

2群が同じ平均・標準偏差で、一方のサンプルサイズを1000から100000に増やしたケース：検定結果に大きく影響せず（棄却され難い）。

2群が同じ平均・サンプルサイズで、一方の標準偏差を10から20に増やしたケース：検定結果に影響（棄却され易い）。

2群が同じ標準偏差・サンプルサイズで、一方の平均を1から２に増やしたケース：検定結果に大きく影響せず（棄却され難い）。

 

2群が同じ標準偏差・サンプルサイズで、一方の平均を1から3に増やしたケース：検定結果に影響（棄却され易い）。

2群が同じ標準偏差で、一方の平均を1から3、サンプルサイズを1000から100000に増やしたケース：検定結果に影響（棄却され易い）。

ランダムサンプリングをしているため、棄却・棄却されないという結果は固定になりません。それを繰り返し実行しただけですが、比較的容易に傾向が見えてきました。
平均2、３は見た目が微妙ですね。やや大きいかな？というような違いでも「有意差がないとは言えない」という傾向を出せるようです。人によってはやや大きい、あるいは同じように見える、などと意見が分かれるような分布ですが、ある考え方に従い再現性のある答えを出してくれるのはありがたい。

最近では、土研交流研究員の方が景観に配慮した護岸に関する類似性の評価に、この検定を利用されていたように思います。アイデア次第でいろいろな分野に適用できそうです。

好きではなかった統計的手法ですが、その合理性をもっと利用すべきでしょうとも考えるようになってきました。

shpfile in ZIP

先日、防災科研さんの地すべり地形分布図をプロットしようとして、ちょっとした問題に引っかかりました。

1．ZIP ファイルに５つの SHP ファイルが納められている。
2．ZIP ファイルは１５０個。
3．ZIP ファイルから、目的の SHP のみを解凍せずに取り出したい。
4．取り出した SHP を結合し、１つの SHP として保存したい。

QGIS で試みるも、複数の SHP からその一部を機械的に取り出せません。
ArcGIS Pro 2.5 では ZIP を扱えません。技があるのかとサポートに問い合わせましたが、やはり標準では対応しておらず、開発&有償になるとのこと（余談ですが、サポートさんのレスが以前に比べずいぶん遅くなりました）。

それほど大仰なモノでもないでしょう。
ということで、簡素にできたのがコチラ。

import os
import glob
import pandas as pd
import geopandas as gpd

path = './input/*.zip' #folder 'input'
zfiles = glob.glob(path) #zip files in the folder
gdf_300=gpd.GeoDataFrame()

for zfile in zfiles:
    fname=os.path.splitext(os.path.basename(zfile))[0]
    gdf_temp=gpd.read_file('zip://'+str(zfile)+'!ls300_'+fname+'.shp')
    gdf_300=pd.concat([gdf_300, gdf_temp])

gdf_300=gdf_300.reset_index(drop=True)
gdf_300.to_file('./output/300.shp', encoding='utf-8')

GIS のカジュアルユーザーだからでしょうか、Python との併用が合理的に感じます。 

5次メッシュを3次メッシュに集約

AVS30（5次メッシュ）を 3次メッシュ内で平均化したいなあと思案。

全国だと EXCEL では行数が足りないので、正攻法の ArcGISで処理。

時間がかかるので Python でも追っかけ処理。
shpを読んだ後の主な処理はこれだけ。

#AVS30=0(水域)を削除
gdf=gdf[gdf['5th_CODE'] != 0]
#3次メッシュコード作成
gdf['3rd_CODE']=gdf['5th_CODE'].astype(np.int64)//100
#3次メッシュにおける平均を求める
df=gdf.groupby('3rd_CODE').mean()

空間演算が必要ないので、Python の方が早く終わりました。
便利。