TRIGRS

なにげなく目に留まった文献の Highlights。

Influence of digital elevation models on the simulation of rainfall-induced landslides in the hillslopes of Guwahati, India
https://www.sciencedirect.com/science/article/abs/pii/S0013795219309603

•Influence of digital elevation model (DEM) on prediction of rainfall-induced landslides.
•Application of TRIGRS, a physically based GIS model.
•Prediction of rainfall-induced landslides for various rainfall events.
•Use of receiver operating characteristics (ROC) to judge suitability of DEMs.
•Higher resolution DEM may not lead to accurate prediction of landslide scenarios.
•Success of ALOS-3D DEM in predicting rainfall-induced landslides in Guwahati city.

ALOS の DEM が選択肢に入るんだなあという点と、 TRIGRS っていうのは何？GIS ベースって簡単そう、といった点に引っ掛かりました。

TRIGRS、調べてみると簡単でした。
document: https://pubs.usgs.gov/of/2008/1159/
code: https://code.usgs.gov/usgs/landslides-trigrs
document: https://pubs.usgs.gov/of/2008/1159/
parallel:https://www.sciencedirect.com/science/article/pii/S1364815216300822?via%3Dihub

準備物は以下の grid のみでOK。GIS で用意するので手軽です。
8の雨量を grid で用意するとなると、東北地方全域のような広域だとやや面倒。5次メッシュか6次メッシュから3次メッシュコードを追加後に grib2 とマージすれば良いでしょうか。そうなると、データ数からして GIS でなく GeoPandas の出番ですね。いえ、そこまで広域に対応させるようなモデルでもないか。

1．DEM grid： dem.asc
2．slope direction index grid:directions.asc (Spatial Analyst, 流向ラスターの作成 (Flow Direction))
3．slope angle grid: slope.asc
4．property zone grid: zones.asc
5．depth grid: zmax.asc
6．initial depth of water table grid: depthwt.asc
7．initial infiltration rate grid: rizero.asc
8．List of file name(s) of rainfall intensity for each period: ri1.asc, ri2.asc

残念ながら、地表流の扱いは以下の通り。モデル化の際に留意すべき点です。

Overland flow between adjacent cells is assumed to occur instantaneously; we have made no attempts to model the rate of overland flow. Consequently, individual storm periods should be long enough to allow surface water to flow to adjacent cells. The routing method enforces mass balance for each time step through the storm; the total precipitation at all the cells must equal the water that infiltrates at all the cells and the water that flows to edges of the problem domain (or any closed depressions) without infiltrating. Thus, excess precipitation that cannot infiltrate at its cell of origin has the opportunity to infiltrate at some point down slope. However, TRIGRS does not carry runoff over from one time step to the next, nor does it track water that enters storm drains. In other words, water that runs off a cell during any given time step will either infiltrate at another cell or reach the edge of the model (or some other sink) within that time step.

Git のコミットを見ると MPI で Parallel 化済み。バイナリは Cygwin 経由だったので、古い2.0のバイナリを起動。テストデータを見ると動きました。簡単。


降雨浸透により地下水が上昇し、斜面崩壊にかかわる安全率が下がる。
昔から言われていますが、地下水の上昇と安定計算の連成、しかも広域の表流水と地下水の収支をチェックしている例をあまり見ません。広域だと表層土砂の層厚分布が不均一なのはもちろんですし、力学・透水パラメータが多様で計算結果を合わせづらくしている点、それらをある一定の精度で把握する術を持っていませんので粗い計算になならざるを得ない点などが理由でしょう。（地すべりなど、単一ブロックを対象にした小さなモデルは容易に作れますので、そちらの方向に行きがちになるのでしょう。）

ただ、粗くても広域を把握する手法があれば、使い道はあります。上記のような計算ではないですが、土壌雨量指数や CL などがその方向性の代表ですし、リスクの顕在化を把握する一つの指標として有用です。
タイムステップの取り方等を工夫し、表流水と疑似三次元、安定計算を組み合わせ Web GIS に集約している地域もあるようです。

残念ながら、私はこの手の計算を実施したことがありませんでした。そもそも、表流水と地下水を連成させるのが GSFLOW まででしたし、それも十分ではありません。
TRIGRS はいくつか問題点もありますが、それを補うツールになり得ます（国内で積極的に使えるとは考えませんが）。

さらなる連成、広域と単一ブロックの中間程度のモデルに適したツールの作成が残る課題となりそうです。
ま、大きな前進ではありませんので、やるべきことはこれまでと同じでしょうか。

Folium

Folium 0.10.1
https://python-visualization.github.io/folium/

GeoDataFrame を Folium で地理院地図ベースに図化する例です。
※見やすさのため不要なインデントをつけています。

import folium

##ポリゴンからコロプレスマップ
  #地物の数が多い（数万）とブラウザがクラッシュ。
  m = folium.Map([36, 140], 
    zoom_start=7, 
    tiles='https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png',
    attr='淡色地図')
  m.choropleth(
    geo_data=gdf,
    name='choropleth',
    data=df,#gdfとdfは結合していない
    columns=['key', 'A'], # dfのkeyと表示データ
    key_on='feature.properties.A', # gdfのキー 
    fill_color='YlGn',# 色パレットを指定
    threshold_scale=[0,2,4,6,8,10,12,14], # 段階
    fill_opacity=0.7,
    line_opacity=0.2, 
    legend_name='aaa')
  m # ポリゴンの多い場合は、コメントアウトして表示させない。
  m.save('./AAA.html')

##ポイント
  #ポイント表示は数千まで?
　#多い場合は表示させずにhtml保存。
  #GeoDataFrameを種類で分割
  gdfA=gdf[(gdf['A']==aaa)]
  gdfB=gdf[(gdf['A']==bbb)]
  #ベースマップ
  m = folium.Map([36, 140], 
    zoom_start=9, 
    tiles='https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png',
    attr='淡色地図')
  #1種類目：赤色でプロット
  for index, row in gdfA.iterrows():
    　　p_txt = '種類:' + row['A']
    　　p_txt = p_txt + '年月日:' +
  　　　　　　　    　　str(row['年'])+'/'+
  　　　　　　　　　　　str(row['月'])+'/'+
 　　　　　　　　　 　　str(row['日'])
    　　popup = folium.Popup(p_txt, max_width=2650)
    　　folium.Circle(#units of meters 
    　　#folium.CircleMarker(#units of pixcels
            　　radius=100,
            　　location=[gdfA.Lat[index], gdfA.Lng[index]],
            　　popup=popup,
            　　color='red',
            　　fill_opacity=0.5, line_opacity=1,
            　　fill=True,
            　　fill_color='red'
        　　　　).add_to(m)
  #2種類目：青色でプロット
  for index, row in gdfB.iterrows():
    　　p_txt = '種類:' + row['A']
    　　p_txt = p_txt + '年月日:' +
  　　　　　　　    　　str(row['年'])+'/'+
  　　　　　　　　　　　str(row['月'])+'/'+
 　　　　　　　　　 　　str(row['日'])
    　　popup = folium.Popup(p_txt, max_width=2650)
    　　folium.Circle(#units of meters 
    　　#folium.CircleMarker(#units of pixcels
            　　radius=100,
            　　location=[gdfB.Lat[index], gdfB.Lng[index]],
            　　popup=popup,
            　　color='blue',
            　　fill_opacity=0.5, line_opacity=1,
            　　fill=True,
            　　fill_color='blue'
        　　　　).add_to(m)
  m # ポイントの多い場合はコメントアウトして表示させない。
  m.save('./AAA.html')

GeoPandas GeoDataFrame 操作

GeoPandas 0.6.0
基本操作です。
※見やすさのため不要なインデントをつけています。

空間演算は 以下を参照。
 spatial overlays
https://github.com/geopandas/geopandas/blob/master/examples/overlays.ipynb
Spatial Joins
https://github.com/geopandas/geopandas/blob/master/examples/spatial_joins.ipynb

import geopandas as gpd
#Ver. 0.6.3

#読み込み
  #SHP、GeoJSON、GeoPackage に対応 
  gdf = gpd.read_file('./A.shp')
  #フォルダ内のSHPファイルを全て読み込み結合
  import glob
  import pandas as pd
  shpfiles = [
   shpfile for shpfile in glob.glob('./shp/*.shp')]
  gdf = pd.concat(
   [gpd.read_file(shpfile)for shpfile in shpfiles])
  #df(csv)からgdfを作成
  #https://gis.stackovernet.com/ja/q/47396
  import pandas as pd
  from shapely.geometry import Point
  geometry = [Point(xy) for xy in zip(df.lon, df.lat)]
  crs = {'init': 'epsg:6668'}
  gdf = gpd.GeoDataFrame(df, crs=crs, geometry=geometry)

#書き込みSHP、GeoJSON、GeoPackage 
  #shpは仕様上、属性255種までの制限あり。型の制限もあり。
  #geojsonはQGISで読みこみ可。Arcはtoolboxから変換。
  #geojsonはdaytime型不可。→事前にstr型変換
  gdf.to_file('./A.shp')
  gdf.to_file('./A.geojson', driver='GeoJSON')

#型変換(pandasと同操作)
  gdf['時刻']=gdf['時刻'].astype(str)     
     
#A列（属性）削除 (pandasと同操作)
  gdf.drop('A',axis=1)

#列名（属性名）変更 (pandasと同操作)
  gdf.columns = ['A','B','geometry'] 

#'都道府県'列（属性）から3県を含む行抽出 (pandasと同操作)
  gdf3=gdf[(gdf['都道府県']=='茨城')|
           (gdf['都道府県']=='埼玉')|
           (gdf['都道府県']=='群馬')|]

#型変換(pandasと同操作)
  gdf['時刻']=gdf['時刻'].astype(str)     
     
#属性の追加、演算(pandasと同操作)
  gdf['C']=gdf['A']-gdf['B']

#県単位でディゾルブ（地物をまとめる）
  gdf=gdf3.dissolve(by='都道府県')
  
#属性結合
  #gdfにA列をkeyとしてDataFrameを属性結合
  #df.merge だと df になる
  gdf=gdf.merge(df, on='A')

#Coordinate Reference Systems
  #CRS確認
  gdf.crs
  #gdf1のCRSをgdf2のCRSに変更（一致）
  #EPSG:4612 : JGD2000
  #EPSG:6668 : JGD2011
  gdf1 = gdf1.to_crs(gdf2.crs)
  
#地図プロット
  #matplotlibベース
  #1レイヤー
  gdf.plot(figsize=(10, 10),
           alpha=0.5,
           edgecolor='k')
     
  #3レイヤー：gdfを重ねて表示  
  　#Choropleth Map(ポリゴン1枚目)
  　base=gdf1.plot(figsize=(10, 10),
      　  column='A',
          cmap='OrRd',
          alpha=0.8,
          legend=True,
          legend_kwds={'label':'aaa',
          orientation:'horizontal'}) 
    #ポリゴンを重ねる（2枚目）
    ax=gdf2.plot(ax=base,
                 facecolor='none',
                 edgecolor='gray',
         　　    alpha=0.8)                
  　#ポイントを重ねる（3枚目）
  　ax=gdf3.plot(ax=base,
          　　   marker='o',
          　　   color='red',
          　　   markersize=5)
　　#仕上げ
    ax.set_xlabel('longitude')
    ax.set_ylabel('latitude')
    ax.grid()
    ax.set_axisbelow(True) #grid moves back
    #保存する場合は import matplotlib.pyplot as plt からの
    plt.savefig("./aaa.png")
  

20200220 Foliumでの図化は移動しました。
https://phreeqc.blogspot.com/2020/02/folium.html


番外編： rasterstats を用いたラスターからの統計量取得。
大量処理は速度面で不可。

#空間演算
  from rasterstats import zonal_stats
  zonal_stats(gdf, './aaa.tif', stats=['max'])
  #フォルダ内のgeotiffを全て読み込み、統計量をgdfに結合
  import glob
  import os
  path = './input/*.tif'
  files= glob.glob(path)
  for i, name in enumerate(files) :
    print(i)
    file_name = os.path.splitext(
                            os.path.basename(files[i]))[0]
    print(file_name)
    stats=pd.DataFrame(zonal_stats(gdf,
                                   files[i],
                                   stats=['max', 'count']))
    gdf_max[file_name]=stats['max']
    gdf_count[file_name]=stats['count']

GeoPandas

GeoPandas 0.6.0
https://geopandas.org/

GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types. Geometric operations are performed by shapely. Geopandas further depends on fiona for file access and descartes and matplotlib for plotting.

pandas の使い勝手が良いので GIS データのハンドリングにこれを選びました。現時点でいくつか問題はあるようですが、pandasに慣れていたら新たに覚えることが少ない、という利点があります。ありがたい。

使い始めたきっかけは速度。
数十万の属性データを ArcGIS Pro 等で操作しようとすると、かなり待たされます。いえ、地物の表示は問題ない速度ですが、属性が遅い。
で、この操作部分を Python で扱おうと。もともと、属性を Python で作っていましたので、その延長で GIS データまで作ってしまえば Arc では表示のみで済むのではないかと考えたわけです。
結果、正解でした。

属性内のデータ検索、抽出を python で実施。CSV データや HDF5 等の部分読み込みからあらたな属性作成・演算・結合も Python で実施。もちろん、CSV データを GIS データに変換可。これ、効率的です。
地物表示もできますが、グリグリ動かすにはブラウザベースだとパワー不足。データ分布の確認と印刷物の作成まででしょうか？

GeoPandas に限りませんが、将来的には Python 等で データ管理と GIS データ作成 を担わせ、閲覧だけを GIS で、といった流れになるのかもしれません。

**********************************
20200216追記
既に PostGIS、PostgreSQL といったものが利用されているようです。
https://postgis.net/
https://www.postgresql.org/

20200217追記
基本操作をまとめました。
https://phreeqc.blogspot.com/2020/02/geopandas-geodataframe.html

pandas DataFrame 操作

先日、pandas が 1.0 に　Ver.UP。
痒いところに手の届く機能が追加されました。ありがたい。

過去、いくつか備忘録としてまとめていますが、コードとテキストが区別されてておらず見難かったため修正しました。このあたり↓
pandas データI/O：https://phreeqc.blogspot.com/2019/02/python3-io.html
df 内容確認：https://phreeqc.blogspot.com/2019/02/python3_62.html
データ型：https://phreeqc.blogspot.com/2019/02/python3.html
ファイル操作：https://phreeqc.blogspot.com/2019/02/python3_54.html

DataFrame 関連でよく使う基本操作を、以下にまとめ直しました。
なお、先頭のインデントは見やすさだけで付けたものです。使用時は外す必要があります。

import pandas as pd
#結合
 # ABでkey列の共通したデータだけ抽出し結合
 df=pd.merge(dfA, dfB, on='key')
 #left指定で最初（左側）のデータフレームは全て残す
 df=pd.merge(dfA, dfB, on='key', how='left')
 #縦に連結（列の重複なし、index重複あり）
 df=pd.concat([dfA, dfB, dfC])
 #横に連結（列の重複あり、index重複なし）
 df=pd.concat([dfA, dfB], axis=1)

#要素数カウント
  df.size
  #A列のユニーク要素数カウント
  df['A'].value_counts()
 
#要素指定
  df.iloc[行番号][列番号]
  df.iloc[:,3:]
  df.loc[行名, 列名]
  df[列名][行名]

#条件検索、操作(行削除、行抽出、要素置換)
  #値を検索して該当行削除・抽出
    #A列が-1の行を削除（-1でない行を抽出）
    df=df[df['A'] != -1]
    #A列が1000より大きな行を削除(1000以下を抽出)
    df=df[df['A'] <= 1000]
    #'くだもの'列から3種指定し行抽出（その他削除）
    df=df[(df['くだもの'] == 'りんご')|
          (df['くだもの'] == 'いちご')|
          (df['くだもの'] == 'めろん')]
  #値を検索して置換
    #A列が0の時は0、それ以外は1で行全体を置換
    df=df.where(df['A'] == 0, 1)
    #A列が0の時は0、それ以外は1でA列のみ置換
    df['A']=df['A'].where(df['A'] == 0, 1)

#Nan
  #A列の要素がNanならTrue判定。B列を追加
  df['B']=df['A'].isnull()
  #A列の要素がNanなら'aaa'で穴埋め
  df['A']=df['A'].fillna('aaa')
 
#要素抽出
  #dfAとdfBの0列を比較し、BにないAの要素を抽出
  dfA[~dfA[0].isin(dfB[0])]
 
#最大最小抽出
  df['A'].idxmax() #A列最大値の行名抽出
  df['A'].idxmin() #A列最小値の行名抽出
  df['A'].max() #A列最大値抽出
  df['A'].min() #A列最小値抽出

#累積
  #行方向累積和-B列追加
  df['B']=df['A'].cumsum()
  #列方向累積和
  df.cumsum(axis=1)
  #行方向最大値維持
  df.cummax()
  #累積積
  df.cumprod()

#行抽出
 df.iloc[行番号] 
 df.loc[行名]
 
#行追加
  df=df.append(aaa, ignore_index=True) #aaa=[・・・]

#行名ふり直し
  #0から連番で振り直し
  df=df.reset_index(drop=True)
  #指定
  #0,3,5→idx=[0,1,2,3,4,5]
  df=df.reindex(idx, fill_value=0) 
#列名取得
  df.columns
  df.columns.tolist()
 
#列名変更
  df.columns = ['A', 'B', 'C']
#列名・行名一部変更
　　df.rename(columns={'A': 'a'},
　　　　　　　　　　　 index={'B': 'b'},
 　　　　　　　　　　　inplace=True)
#列並べ替え
  df=df[['C', 'A', 'B']]

#行削除
  df=df.drop(['A'])
 
#列削除
  df=df.drop(['A'],axis=1)

#集計
  df.values.sum() #全合計

#集約
  pd.crosstab(df['A'], df['B'])
  df.groupby(['A','B']).size().unstack('B', fill_value=0)
  pd.pivot_table(df, index=['A'], columns=['B'],
                 aggfunc=len, fill_value=0)

#時間
  #時間型に変換
  df['date']=pd.to_datetime(df['date'])
  #列を結合し、date列作成
  df['date'] = pd.to_datetime(df['西暦'].astype(str)+'/'+
                              df['月'].astype(str)+'/'+
                              df['日'].astype(str)+'/'+
                              df['時刻'].astype(str),
                              format='%Y/%m/%d/%H:%M:%S')       
  #時刻が書かれていなければ0時を記入
  df['時刻']=df['時刻'].fillna('00:00:00')
  #丸め
    #時間単位でround up
    df['date']=df['date'].dt.ceil("H")
    #30分単位でround
    df['date']=df['date'].dt.round('30min')
    #時間単位でround down
    df['date']=df['date'].dt.floor('30min')
  #日付を四半期表示に変換
  dfQ=pd.PeriodIndex(df['date'], freq='Q')
  #四半期データ作成（期間指定）
  df['date']=pd.period_range(
　　　start=pd.Period('2008Q1', freq='Q'), 
      end=pd.Period('2019Q4', freq='Q'), 
      freq='Q')
#LOOP処理
  #行を順に処理
  for index, row in df.iterrows()
　#iteritems()→カラム名、カラムのデータ（Series形状）を返す
  #iterrows() →ラベル名、ラベルのデータ（Series形状）を返す
 
#DataFrame作成
  #列名を決めた空のdf作成
  df= pd.DataFrame(index=[], columns=['A', 'B', 'C'])
  #データから作成
  import numpy as np
  x, y, z, v = np.random.random((4, 1000))*100
  df=pd.DataFrame({'x':x,'y':y,'z':z,'v':v})
  #クリップボードから作成
  df = pd.read_clipboard()

#不要データフレーム削除
  del(df1,df2)

20200216追加・修正
20200311追加・修正
20200610追加
20210206追加　引用：集約関数

文献整理

たまった文献の整理。
気になったものだけメモ( ..)φ

*****************************************
林ほか「土砂災害に関する大雨の豪雨度と崩壊個数の関係」
砂防学会誌, 72(2019)4,345
広島県呉市
大災害後4年、小災害（1か月複数回）の影響（災害発生数の低下）を除けば、年間崩壊個所数は豪雨度^4.3に比例

菊地ほか「平成23年台風12号における崩壊・非崩壊地に関する地形的特徴の検証」
日本地すべり学会誌,56 巻 (2019) 4 号

深層崩壊あり、なし観点から地形を区分。区分表では代表例を示す。

壱岐「航空レーザ測探機による沿岸の測量」
地質と調査, 2019,2,154
水深10m

国総研資料第686号
「山地河道における流砂水文観測の手引き(案)」
濁度計
ハイドロフォン
国総研作成の流砂量変換プログラム
砂面計
p14【参 考】濁度計の出力値と浮遊砂の容積濃度の関係式

濁度計の出力値（V）から浮遊砂の容積土砂濃度（C）に変換する式は、一般的に次式が用いられる。
C＝AV2+BV
ここで、A、B はキャリブレーションで求める係数。
なお、濁度計に付属するキャリブレーション式は、NTU という単位に変換する式が示されている場合が多く、この式は容積土砂濃度への変換には使用できない。
（参考文献）平成16 年新潟県中越地震後における浮遊砂量の観測、小山内信智・水野秀明・清水武志・沖中健起・原槇利幸、国土技術政策総合研究所資料、第278 号 2006

国総研資料1013号
「ソーシャルメディア分析によるリアルタイム災害発生情報検知手法の確立に関する研究」
twitter
→小林ほか「土砂災害情報収集システムにおける発災推定用キーワードの検証」

日本語は難しい

プログラムを作る際には要求事項、仕様を確認します。

口頭で依頼されたり書面で渡されたりするのですが、わかりにくい場合があります。多くの場合はフローを用いたり数式を用いて、依頼された方の思考過程を追っていくことで解決します。が、それもできず口頭などで説明を繰り返されると理解不能に陥ります。

先日、難しい日本語に出合いました。何を仰りたいのか概要は分かるのですが、具体的にどの数値を、どの順で、どう操作すればよいかまでは変換できない文章です。

この方の場合、話を聞いていくうちに、P(A|B) の条件付確率であれば良いことが分かりました。さらに、求めようとしていた答えよりも良い指標があり、一つ手前に検討すべき内容が抜けていたこともわかりました。シンプルに6文字で表現できていたら、悩まずミスしなかったでしょう。ま、自分の考えた過程を数式等に変換できようであれば、御自分でプログラムを組まれたかもしれません。

思い出したのが中学生の頃の話。言葉は通じなくても H2O と書けば水が欲しいことを伝えられると先生から聞き、なるほどと思ったことがあります。英語の論文でも、化学式や数式だけは読み違えません。容易に通じ合えます。

日本語は難しいのです。日本語を日本語に翻訳するより、数式等で表現した方が間違いなく相手に伝わります。お互いに「誤解を与える表現であった」「そういう意味ではなかった」という逃げも打てません。

自分がどのように考えているのか、その過程を自分で理解し伝えるのは難しいと思います。そこに一つの壁があることは間違いありません。この点では、小学生へのプログラミング教育は有効でしょう。

日本語を鍛えつつ、論理的かつ効率的に考える能力も身に付けたい。相乗効果で伝えやすい文章を作る能力が備わるのだと思われます。

統計学の今後

統計学ではデータセットを training（＋ varidation）, test に区分することはなかったでしょう（況や、CV をや）。

このあたり、機械学習では標準。結果を出すために必要な手順です。
重回帰でもできるのですが、どの統計の図書を見ても理論ばかりで、この流れには触れられていませんでした。予測に対して理論よりも結果（性能）重視の立場が機械学習、理論武装（仕様）重視が統計というような感じを受けます。
単にプログラミングしやすくなったというのも後発に影響しているでしょう。ハードルが低くなり、統計でも CV 等が当たり前になってくるのでしょうか。今、統計のプロはどうされているのでしょう。

機械学習における特徴量選択手法の呼び方 wrapper method は明らかにプログラミングの立場からですね。もちろん、手順に則した呼び方もあります。 Sequential Feature Selection (Sequential Forward Selection , Sequential Backward Elimination, Sequential Forward Floating Selection (SFFS), Sequential Backward Floating Selection (SBFS)), Exhaustive Feature Selection 等。
手元の統計学の図書では SFFS を stepwise と記載されています。手順に関わる呼び方はある程度共通なのかもしれません。が、統計で wrapper method と呼ぶ発想はないでしょう。ある程度の住み分けができているのでしょうね。

どちらの立場からでも良いと思います。両方知っているのが理想でしょう。
PCA、ラフ集合、回帰など、統計の世界にも入り始めましたので、今後も機会があれば学んでいきましょう。

重回帰分析

機械学習側から純然たる統計学に近づいてきました。

プログラミングの側面から見ると、データを作って手法に投げ込む形は同じ。しかもより簡素。これなら理解しようという気になります。
で、regression。いろいろありますね。sklearn のチートでは、特徴量選択を含めると Lasso, Elastic net が選択される流れです。

・重回帰：単回帰の高次元版
・ロバスト回帰：外れ値の影響低減
・Lasso回帰：L1正則化項の追加（係数ノルムの合計×α）・・・次元圧縮
・Ridge回帰：L2正則化項の追加（係数2乗の合計×α/2）
・Elastic Net：L1,L2正則化項の追加

始めたのが重回帰分析。これは単回帰の特徴量を増やしただけなので、答えを出すだけなら新たに理解すべきことはありません。が、特徴量選択や適合性評価については、培われてきた独自の流儀があるようです。

以下、備忘録。
*****************************************
【重回帰モデルのあてはまりの良さの評価】
Adjuted r2：自由度調整済決定係数

AIC：赤池情報量基準

　r：重相関係数

　r2：決定係数　Sr/St=予測値の平方和/実測値の平方和=1-Se/St

　S：平方和　データと平均との差の2乗和=分散をn

【重回帰モデルが予測に役立つかの評価】
F-Value：F値　Vr/Ve=(Sr/p)/(Se/(n-p-1))
　p：特徴量の数
　n：サンプル数

【特徴量選択】

F値、t値、p値
wrapper method、反復特徴量選択：stepwise

【留意点】
サンプル数が少なくても、最小二乗法（誤差合計の式を偏微分）で係数分の式を作るので、係数自体は求まる。しかし、サンプル数Nに比べ求める係数の数Pが多い場合（自由度N-Pが低いとき）、決定係数r2は大きくなりやすい（N=2の単回帰は誤差0、R2=1.0）。
F値、Adjuted r2 等でのモデル評価が重要。