1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@chicken_data_analyst(ひよっこ アナライズ)

鳥取にスタバを増やしたい：One-Class SVMで「隠れたカフェ立地候補」を見つける

1
Last updated at Posted at 2026-07-14

はじめに

「スタバはないが、日本一の砂場（砂丘）はある」

いつかの鳥取県知事のこの発言が記憶に残っている。

現在スタバは複数店舗あるようだ。とはいえ、やはり他府県と比べると少ないらしく、「「Q. スタバがない県ってある？」「A. ない」でも“ほとんどない” 県は存在した…！」によると、2024年11月時点では鳥取・島根が最も店舗数の少ない県らしい。

「うーーん、もっと出店を増やすためにはどこがカフェ適地なのかがわかればいいのでは？」と考えた。（余計なお世話）

そこで、既存カフェの立地パターンをデータから学習し、"カフェが出店しそうな場所"を機械学習で予測してみることにした。

使うのは異常検知手法として知られるOne-Class SVM(OCSVM)。今回は「異常」を探すのではなく「正常（=カフェらしい場所）」をキャッチするために使う。

全体の処理をまとめると、人流・土地利用・カフェ立地データをメッシュ単位で結合し、既存カフェが存在するメッシュ（Presence）の特徴量からOne-Class SVMで「カフェらしい場所」のパターンを学習する。カフェが存在しないメッシュ（Background）とのROC-AUCでハイパーパラメータを最適化し、ROC曲線から決めた閾値で二値化したモデルを鳥取県の土地利用データに適用することで、実在しないが立地条件的に「カフェらしい」エリアを推定・可視化する。
OCSVM_pipeline.png

コードは一通り生成AIにまとめてもらいコメントも書いてもらった。間違ったコメントとかあったらサーセン。

ちなみに今回の結果は以下。赤い点が濃いほどスタバの立地適正が高い。(赤い点に隠れているけど)青いX印は既存のスタバ。
image.png

参考

import

必要ライブラリ

import glob

import numpy as np
import pandas as pd
import geopandas as gpd
import osmnx as ox
from tqdm import tqdm

import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import mpl_toolkits.axes_grid1
import seaborn as sns
import cartopy.crs as ccrs
import contextily as cx
from matplotlib_scalebar.scalebar import ScaleBar

from sklearn.svm import OneClassSVM
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import KFold
from sklearn.metrics import pairwise_distances
from scipy.stats import wasserstein_distance
from scipy.stats import wasserstein_distance_nd
from sklearn.model_selection import KFold
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import roc_auc_score

import jismesh.utils as ju

import japanize_matplotlib
# 日本語フォント読み込み
japanize_matplotlib.japanize()
jpn_font = japanize_matplotlib.get_font_ttf_path()
prop = fm.FontProperties(fname=jpn_font)
sns.set()
plt.rcParams['font.family'] = prop.get_name() #全体のフォントを設定

データ

人流データの取得

全国の人流オープンデータ（1kmメッシュ、市区町村単位発地別）から人流データを../dataにダウンロードする。（今回は中国地方のみ）

../data直下にある.zipファイルをすべて同じ../dataフォルダに展開。人流データの配布形式では、まず年度や地域の大分類ごとにzipが分かれているので、この段階では「1段目の圧縮」を解いているだけ。

# ============================================================
# 0-1. ダウンロードした人流データ（zip）を解凍
#      ../data 直下にある zip を、直下に展開する
# ============================================================
import zipfile

for dr in glob.glob('../data/*.zip'):
    with zipfile.ZipFile(dr) as zf:
        zf.extractall('../data')

人流データは「フォルダ1（例：都道府県）／フォルダ2（例：年）／フォルダ3（例：月）」という3階層となっている。

  • glob.glob(flow_dir + '/[0-9]*[!zip.geojson.csv]')で、数字始まりかつ拡張子がzip/geojson/csvではないフォルダ（＝1段目で展開されたディレクトリ）だけを抽出
  • その中をdr1→dr2→dr3と3段階掘り進み、各dr3フォルダ内にある.zipを見つけて、同じ階層（dr3）に展開
%%time
# ============================================================
# 0-2. 年フォルダ配下にネストされた zip（都道府県・メッシュ単位など）を再帰的に解凍
#      構造想定: flow_dir/<フォルダ1>/<フォルダ2>/<フォルダ3>/*.zip
# ============================================================
flow_dir = '../data'
# 数字始まりのフォルダのみを対象（"zip", "geojson", "csv" で終わるものは除外）
flow_folders = sorted(glob.glob(flow_dir + '/[0-9]*[!zip.geojson.csv]'))

for dr1 in tqdm(flow_folders):
    print('#######', dr1, '#######')
    for dr2 in glob.glob(dr1 + '/*'):
        for dr3 in tqdm(glob.glob(dr2 + '/*')):
            # print(dr3)  # デバッグ用（通常は不要）
            for dr4 in glob.glob(dr3 + '/*.zip'):
                # print(dr4)  # デバッグ用（通常は不要）
                with zipfile.ZipFile(dr4) as zf:
                    zf.extractall(dr3)  # 同じ階層に展開

解凍済みの生CSV群を読み込み、分析用に使いやすい1つの集計テーブルへまとめる。

  • 同じフォルダ構造をもう一度走査し、2019年・2020年・2021年フォルダ内の全CSVをpd.concatで縦結合
  • dayflag == 2（全日）かつtimezone == 0（昼間）の条件でレコードを絞り込み、カフェを利用する機会が多い「全日昼間」の人流のみを対象とする
  • groupby(['mesh1kmid', 'prefcode', 'citycode', 'year'])で、1kmメッシュ×都道府県×市区町村×年の単位ごとにpopulationのmin・median・maxを計算し、各地点の人口変動の代表値を取得
  • 全フォルダ分の集計結果をpersonflow_dに積み上げ、最終的に../data/人流データ.csvとして保存
# ============================================================
#    解凍済みCSV（2019〜2021年の人流データ）を1つのDataFrameに集約
#    メッシュ×都道府県×市区町村×年 単位で population の min/median/max を計算
# ============================================================
flow_dir = '../data'
flow_folders = sorted(glob.glob(flow_dir + '/[0-9]*[!zip.geojson.csv]'))

personflow_d = pd.DataFrame()
for dr1 in flow_folders:
    print('#######', dr1, '#######')
    tmp1 = pd.DataFrame()

    # 2019, 2020, 2021 年フォルダ内のCSVをすべて読み込んで縦結合
    for dr2 in glob.glob(dr1 + '/2019/*') + glob.glob(dr1 + '/2020/*') + glob.glob(dr1 + '/2021/*'):
        tmp2 = pd.DataFrame()
        for dr3 in glob.glob(dr2 + '/*.csv'):
            tmp3 = pd.read_csv(dr3)
            tmp2 = pd.concat([tmp2, tmp3])
        tmp1 = pd.concat([tmp1, tmp2])

    # 全日昼間（dayflag==2, timezone==0）に絞り、メッシュ×属性×年別に統計量を集計
    tmp1_mean = (
        tmp1[(tmp1['dayflag'] == 2) & (tmp1['timezone'] == 0)]
        .groupby(['mesh1kmid', 'prefcode', 'citycode', 'year'], as_index=False)['population']
        .agg(['min', 'median', 'max'])
    )
    personflow_d = pd.concat([personflow_d, tmp1_mean])

personflow_d = personflow_d.reset_index(drop=True)
personflow_d.to_csv("../data/人流データ.csv", index=False)  # 集計結果をキャッシュとして保存
display(personflow_d)

image.png

縦持ちの人流データをpivot_tableで年ごとに横展開し、min/median/maxの3テーブルをmesh1kmidで結合することにより「1メッシュ1行」の横持ちデータに変換。2019年のmedian欠損メッシュを除外した上で、以降の分析では2019年のmin/median/maxを候補列とし、代表値として2019_medianを採用。（コロナ禍を除外）

# ============================================================
#    年ごとに縦持ちだったデータを横持ち（mesh1kmid × 年_統計量）に変換
#    2019/2020/2021 × min/median/max の組み合わせを1行1メッシュにまとめる
# ============================================================
personflow_d_min = personflow_d.pivot_table(index='mesh1kmid', columns='year', values='min').reset_index()
personflow_d_median = personflow_d.pivot_table(index='mesh1kmid', columns='year', values='median').reset_index()
personflow_d_max = personflow_d.pivot_table(index='mesh1kmid', columns='year', values='max').reset_index()

# 列名末尾に _min/_median/_max を付けながら mesh1kmid をキーに結合
personflow_d_agg = (
    personflow_d_min
    .merge(personflow_d_median, on=['mesh1kmid'], suffixes=[None, '_median'])
    .merge(personflow_d_max, on=['mesh1kmid'], suffixes=['_min', '_max'])
    .dropna(subset=['2019_median'])  # 分析の基準年である2019年median欠損のみ除外
)
display(personflow_d_agg)

# 分析で使う人流の統計量列（2019年の min/median/max）と、代表値として使う列を指定
flow_cols = ['2019_min', '2019_median', '2019_max']
selected_flow = '2019_median'

image.png

カフェ立地データの取得

鳥取県のスタバの立地候補を見つけたいので、中国地方の立地状況を学習させたい。
OpenStreetMap（OSM）から中国地方5県のスターバックス店舗座標を取得し、GeoJSONとして保存する。

  • osmnx（ox）のfeatures_from_place関数を使い、鳥取・島根・岡山・広島・山口の5県を対象にOSM上の地物データを取得
    • tagsでbrandまたはnameタグが「Starbucks」「スターバックス」「スターバックス コーヒー」のいずれかに一致する地物を絞り込み条件として指定
    • OSMでは店舗タグの表記が英語・日本語・カタカナ表記などで揺れることがあるため、複数の表記パターンを列挙
  • 取得した地物の座標形式を整え、必要な列だけを抽出して保存
    • OSM上の店舗は建物のポリゴン（多角形）として登録されている場合があるため、representative_point()でポリゴンの内部にある代表点を1つ選び、扱いやすい単一のPoint（点）ジオメトリに変換
    • id・amenity・brand・name・name:en・name:ja・geometryの列に絞り込み、drop_duplicates()で重複データを除去した上で、../data/スタバ_中国地方.geojsonとして保存
# ============================================================
# 1. データ取得: 中国地方5県のスターバックスcafe座標
# ============================================================
places = [
    "Tottori Prefecture, Japan",
    "Shimane Prefecture, Japan",
    "Okayama Prefecture, Japan",
    "Hiroshima Prefecture, Japan",
    "Yamaguchi Prefecture, Japan"
]

tags = {
    "brand": ["Starbucks", "スターバックス", "スターバックス コーヒー"],
    "name": ["Starbucks", "スターバックス", "スターバックス コーヒー"]
}

print("中国地方のデータを取得中...")
starbucks_gdf = ox.features_from_place(places, tags).reset_index()
starbucks_gdf["geometry"] = starbucks_gdf.geometry.representative_point()  # Polygon -> Point に変換

starbucks_gdf[["id", "amenity", "brand", "name", "name:en", "name:ja", "geometry"]].drop_duplicates().to_file("../data/スタバ_中国地方.geojson")
display(starbucks_gdf)

image.png

土地利用データの取得

国土数値情報ダウンロードサイト TOP > 国土数値情報 > 土地利用3次メッシュデータから3次メッシュの土地利用データを../data/land_use/にダウンロードする。
../data/land_use/配下の各geojsonを読み込んでCRSを"EPSG:4326"に統一し、pd.concatで結合することで、都道府県別などに分かれていた土地利用データを扱いやすい単一テーブルにまとめる。

# ============================================================
# 2. 土地利用データの読み込み・結合（geojsonを1つに統合）
# ============================================================
# 国土数値情報の1kmメッシュ土地利用
luse_files = glob.glob("../data/land_use/*.geojson")
print(luse_files)

lu_list = []
for f in tqdm(luse_files, total=len(luse_files)):
    lu = gpd.read_file(f)
    lu = lu.to_crs("EPSG:4326")  # CRS を統一
    lu_list.append(lu)

landuse_merged = gpd.GeoDataFrame(pd.concat(lu_list, ignore_index=True), crs="EPSG:4326")

# 単位をヘクタールに変換
lu_cols = ["田", "その他の農用地", "森林", "荒地", "建物用地", "道路", "鉄道",
           "その他の用地", "河川地及び湖沼", "海浜", "海水域", "ゴルフ場", "解析範囲外"]
landuse_merged[lu_cols] = landuse_merged[lu_cols] / 10000  # ヘクタール
display(landuse_merged)

image.png

各データを結合

カフェ店舗座標に土地利用属性と人流データを結合し、OCSVMの学習に使う特徴量テーブルを作る。まずsjoinでカフェ座標に土地利用ポリゴンの属性を空間結合し、続いて各座標を1kmメッシュコードに変換して人流データ（personflow_d_agg）とメッシュ単位で結合。
さらに「田」と「その他の農用地」を「農地」として1つの列に統合し、森林・農地・荒地・その他の用地・建物用地・道路・鉄道という7つの土地利用面積列（areas_cols）を分析用の特徴量として定義。最後に、これらの面積列と人流の代表値（selected_flow）にlog1p変換を適用する。

# ============================================================
# 3. cafe座標に土地利用属性を空間結合（1回だけ実施）
# ============================================================
all_coords_geos = starbucks_gdf.copy()
all_coords_with_luse = all_coords_geos.sjoin(landuse_merged, how="inner")

# 人流データと結合
all_coords_with_luse['mesh1kmid'] = all_coords_with_luse.geometry.map(lambda geo: ju.to_meshcode(geo.y, geo.x, 3))
all_coords_with_luse = all_coords_with_luse.merge(personflow_d_agg, on=['mesh1kmid'], how='inner')

# 田 + その他の農用地 → 農地 として統合
landuse_merged["農地"] = landuse_merged["田"] + landuse_merged["その他の農用地"]
all_coords_with_luse["農地"] = all_coords_with_luse["田"] + all_coords_with_luse["その他の農用地"]

# 分析に使う特徴量列（対数変換版も用意）
areas_cols = ["森林", "農地", "荒地", "その他の用地", "建物用地", "道路", "鉄道"]
all_coords_with_luse_log = all_coords_with_luse.copy()
all_coords_with_luse_log[areas_cols+[selected_flow]] = np.log1p(all_coords_with_luse_log[areas_cols+[selected_flow]])

display(all_coords_with_luse_log)

image.png

データ確認

土地利用データと人流データの傾向を確認する。バイオリンプロットで可視化。

# ============================================================
# 4. 既存Cafeの土地利用カテゴリ別 面積分布を可視化（violinplot）
# ============================================================
violin_df = pd.concat([
    pd.DataFrame({"cat": np.full(len(all_coords_with_luse), c), "area": all_coords_with_luse[c]})
    for c in areas_cols
]).reset_index(drop=True)

plt.rcParams["font.family"] = prop.get_name()
plt.figure(figsize=(12, 5))

sns.violinplot(
    data=violin_df, x="cat", y="area", inner="quart", hue="cat",
    palette="BrBG_r", hue_order=areas_cols, native_scale=True, alpha=0.8,
)
sns.stripplot(
    data=violin_df, x="cat", y="area", hue="cat",
    palette="BrBG_r", hue_order=areas_cols, linewidth=1,
)

plt.xlabel("土地利用カテゴリ")
plt.ylabel("面積 [ha]")
plt.title("土地利用カテゴリ別面積の分布（既存Cafe）")
plt.tight_layout()
plt.show()


# ============================================================
# 4. 既存Cafeの人流分布を可視化（violinplot）
# ============================================================
violin_df = pd.concat([
    pd.DataFrame({"cat": np.full(len(all_coords_with_luse), c), "area": all_coords_with_luse[c]})
    for c in flow_cols
]).reset_index(drop=True)

plt.rcParams["font.family"] = prop.get_name()
plt.figure(figsize=(12, 5))

sns.violinplot(
    data=violin_df, x="cat", y="area", inner="quart", hue="cat",
    palette="Greys", hue_order=flow_cols, native_scale=True, alpha=0.8,
)
sns.stripplot(
    data=violin_df, x="cat", y="area", hue="cat",
    palette="Greys", hue_order=flow_cols, linewidth=1,
)

plt.xlabel("統計量")
plt.ylabel("人流 [人]")
plt.title("人流の分布（既存Cafe）")
plt.tight_layout()
plt.show()

image.png
image.png

当たり前だが市街地にあることがほとんどなので、既存カフェの1kmメッシュ内は建物用地面積が大きい。意外と森林(というか山間部？)が近くにある店舗もあるようだ。人流はばらつきがかなり大きそう。中心地と郊外の店舗でもかなり変わりそうなのでそりゃそうか。今回は2019_medianのみ使用する。

OCSVMによる立地候補抽出モデル作成

共通関数定義

正常/異常ラベルの比率を円グラフで表示する共通関数。

def plot_label_ratio_pie(pred_labels, title="OCSVM 判定結果の比率"):
    """正常/異常ラベルの比率を円グラフで表示する共通関数（重複していた描画処理をまとめた）"""
    unique, counts = np.unique(pred_labels, return_counts=True)
    label_map = {1: "正常 (1)", -1: "異常 (-1)"}
    labels = [label_map[u] for u in unique]
    colors = ["#EF553B" if u == -1 else "#636EFA" for u in unique]

    fig, ax = plt.subplots(figsize=(6, 5))
    ax.pie(
        counts, labels=labels, colors=colors, wedgeprops=dict(width=0.6), startangle=90,
        autopct=lambda p: f"{p:.1f}%\n({int(round(p/100*sum(counts)))}件)",
    )
    ax.set_title(title)
    plt.tight_layout()
    plt.show()

バックグラウンドデータの取得

指定メッシュコードの周囲nマス分（n=5なら一辺11マス、計121マス）のメッシュコードを算出する関数。メッシュの南西・北東端の座標から1マス分の緯度経度幅を求め、その幅を使って周囲のメッシュ中心座標を計算し直している。この関数を使い、既存スタバ店舗が存在する全メッシュについて周囲5マス分をinclude_meshesとして集め、店舗が実在するメッシュ自体はinclude_meshesから取り除いている。この処理の目的は、後述の背景データ（Background）を作成することである。

def get_surrounding_meshcodes(meshcode, level=3, n=1):
    """指定した次数の3次メッシュコードについて、自身と周囲 n マス分
    （一辺 2n+1 マスの正方形、計 (2n+1)^2 マス）のメッシュコードを返す

    n=1 なら隣接8マス+自身の計9マス（従来と同じ挙動）
    n=2 なら周囲2マス分（一辺5マス、計25マス）
    """
    # メッシュの南西端・北東端座標から、1マス分の緯度経度幅を計算
    lat_sw, lon_sw = ju.to_meshpoint(meshcode, 0, 0)
    lat_ne, lon_ne = ju.to_meshpoint(meshcode, 1, 1)
    lat_size = lat_ne - lat_sw
    lon_size = lon_ne - lon_sw

    # 自身の中心点座標
    lat_c, lon_c = ju.to_meshpoint(meshcode, 0.5, 0.5)

    neighbor_codes = []
    for dlat in range(-n, n + 1):
        for dlon in range(-n, n + 1):
            lat_n = lat_c + dlat * lat_size
            lon_n = lon_c + dlon * lon_size
            neighbor_codes.append(ju.to_meshcode(lat_n, lon_n, level))

    return list(set(neighbor_codes))  # 重複除去（境界付近での丸め誤差対策）

include_meshes = set()
for m in all_coords_with_luse_log['mesh1kmid'].unique():
    include_meshes.update(get_surrounding_meshcodes(m, level=3, n=5))

include_meshes = [i for i in include_meshes if i not in all_coords_with_luse_log['mesh1kmid'].unique()]
print(f"バックグラウンド対象メッシュ数: {len(include_meshes)}")
# >> バックグラウンド対象メッシュ数: 2705

ハイパーパラメータの探索

カフェの「在データ（Presence）」と「背景データ（Background）」を比較して、ハイパーパラメータを探索したい。

まず、背景データ（Background）のサンプリング設計だが、「店舗の近隣だが店舗そのものではない」メッシュ群を背景データの候補としている。この設計により、背景データが店舗周辺の似た環境から選ばれるため、都市/郊外のような粗い違いだけでAUCが不当に高くなることを防ぎ、より微妙な立地差をモデルに学習させる狙いがある。（要は「山や耕作地じゃなくて市街地が候補だよ！」というような単純なモデルになってしまうことを避けたい。）

evaluate_ocsvm_auc関数は、層化K-Fold（5分割）で在データ・背景データを分割し、各foldの学習時には在データ（y=1）のみでOneClassSVMの境界を学習し、評価時には在データ・背景データが混在するテストセットに対するdecision_functionのAUCを計算する。

gamma_grid（scale・autoおよび0.01から10.0までの13パターン）とnu_grid（0.01から0.5までの5パターン）の全組み合わせについてAUCの平均・標準偏差を計算し、results_dfをAUC平均の降順に並べて上位60件を表示することで、最も汎化性能の高いgamma、nuの組み合わせを確認。

# ============================================================
# 5. OneClassSVM: gamma / nu のグリッドサーチ（AUCで評価）
#   ※ recall や nu の内部閾値(predict)ではなく、
#     decision_function の「順序付け能力」自体をAUCで評価する。
# ============================================================

# 使用する特徴量（土地利用面積 + 人流の代表値）
features = areas_cols + [selected_flow]

# --- 在データ（Presence）: 既存スタバ店舗の特徴量 ---
X_Presence = (
    all_coords_with_luse_log[features]
    .replace([np.inf, -np.inf], np.nan).dropna()  # inf/nanを除外
)

# --- 背景データ（Background）: 在データの周辺メッシュ（include_meshes）からサンプリング ---
#   在データそのものではなく「近隣の似た環境」を背景にすることで、
#   「都市 vs 郊外」のような粗い分離ではなく、微妙な立地差を学習させる
#   （全域からランダム抽出すると、都市/郊外の違いだけでAUCが不当に高くなるため）
X_Background = (
    landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))][features]
    .replace([np.inf, -np.inf], np.nan).dropna()
    .sample(n=min(10000, len(landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))])), random_state=0)
    # 背景データが多すぎる場合に備え、最大10000件にランダムサンプリングして件数を制限
)

# 在データ・背景データを結合し、1=在データ, 0=背景データのラベルを付与
x_all = pd.concat([X_Presence, X_Background], ignore_index=True)
y_all = np.concatenate([np.ones(len(X_Presence)), np.zeros(len(X_Background))])

# OneClassSVM(rbfカーネル)は距離ベースの手法のため、
# 特徴量間のスケール差（面積[ha] vs 人流[人]など）を揃えるために標準化する
scaler = StandardScaler().fit(x_all)
x_all = pd.DataFrame(scaler.transform(x_all), columns=x_all.columns)


def evaluate_ocsvm_auc(x, y, gamma, nu, n_splits=5, random_state=0):
    """
    在データ・背景データを層化K-Foldで分割し、
    「学習に使っていないfold」に対するAUCを計算して汎化性能を評価する関数。

    - 学習(fit)は各foldの「在データ(正常, y=1)のみ」を使う
      （OneClassSVMは教師なし学習のため、背景データは学習には使わない）
    - 評価(predict)は在データ・背景データ両方を含むtestデータに対して行う
    """
    skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=random_state)
    aucs = []
    for train_idx, test_idx in skf.split(x, y):
        # 学習用インデックスの中から「在データ(y=1)のみ」を抽出
        train_Presence_idx = train_idx[y[train_idx] == 1]

        model = OneClassSVM(kernel="rbf", gamma=gamma, nu=nu, max_iter=-1)  # tolはデフォルト固定
        model.fit(x.values[train_Presence_idx])  # 正常データのみで境界を学習

        # decision_function: 境界からの距離（値が大きいほど「在データらしい」）
        y_score = model.decision_function(x.values[test_idx])

        # test側は「在データ」「背景データ」が混在しているため、
        # このスコアで両者をどれだけ正しく順位付けできるかをAUCで測る
        aucs.append(roc_auc_score(y[test_idx], y_score))
    return np.mean(aucs), np.std(aucs)


# ハイパーパラメータの探索範囲
gamma_grid = ["scale", "auto", 0.01, 0.1, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 5.0, 8.0, 10.0]
nu_grid = [0.01, 0.05, 0.1, 0.2, 0.5]

# 全組み合わせについてAUCを計算し、比較用テーブルを作成
results = []
for gamma in gamma_grid:
    for nu in nu_grid:
        auc_mean, auc_std = evaluate_ocsvm_auc(x_all, y_all, gamma, nu)
        results.append({
            "gamma": gamma,
            "nu": nu,
            "auc_mean": auc_mean,
            "auc_std": auc_std,
        })

# AUC平均が高い順に並べ、最も汎化性能の高いgammaを確認
results_df = pd.DataFrame(results).sort_values("auc_mean", ascending=False)
display(results_df.head(60))

正直各組合せはそこまで大きな差が無い。なのでThe one standard error (1SE) rule of thumbに則って、シンプルなモデルとなる組み合わせを選択する。gamma = 0.01、nu = 0.5
image.png
~~
image.png

判定の閾値を決定

OneClassSVM.predict()を使わず、背景データとクロスバリデーションを使って閾値を求める。

  • KFold(n_splits=5)で在データを5分割し、各foldでは「そのfoldのテスト対象を除いた在データ」のみを使ってOneClassSVMを学習。学習に使われていないテスト対象の店舗に対するdecision_functionの値をoof_scoresに記録する
  • 同時に、同じfoldで学習したモデルを使って背景データ（X_Background_s_fixed）のスコア(decision_functionの値)も計算し、Background_scores_per_foldに記録
  • 5つのfoldモデルによる背景スコアを平均したBackground_scores_fixedを、在データのoof_scoresと結合し、ROC曲線（roc_curve）とAUC（roc_auc_score）を計算。最適な閾値は「感度＋特異度が最大になる点」（tpr - fprが最大になるインデックス）として選ばれ、このoptimal_thresholdを基準に、在データ全体のスコア（oof_scores）が閾値以上であれば1（正常＝カフェらしい）、そうでなければ-1（カフェ立地候補ではない）というラベル（pred_oof）を最終的に付与
# ============================================================
# 6. OOF（Out-of-Fold）方式でPresence・Backgroundを同一基準で評価し、
#    閾値をMTSSと同じ考え方（在データ vs 背景データの分離最大化）で決定する
#
#   注意: OneClassSVM.predict()が返すラベルは、nuに基づく内部閾値(rho)を
#   使っており、これは学習fold内でのみ最適化された値。
#   fold間で閾値の厳しさが揃っていない
# ============================================================

# グリッドサーチ結果から選択したハイパーパラメータ
gamma = 0.01
nu = 0.5

# 在データ全体（既存スタバ店舗）を取得。特徴量にinf/nanが無い行のみを対象にする
df = all_coords_with_luse_log.copy()
mask = np.isfinite(df[areas_cols+[selected_flow]].values).all(axis=1)
X_labeled_normal = df.loc[mask, areas_cols+[selected_flow]].values

# スケーラーは在データ全体（54件）で作成し、以降すべての変換に共通で使う
# （foldごとに別のスケーラーを作ると、背景データとのスケール比較ができなくなるため統一）
scaler = StandardScaler().fit(X_labeled_normal)

# 背景データ（在データ周辺メッシュからサンプリング）を取得
# ※ このデータは学習には使わず、「モデルが背景をどう評価するか」の基準点としてのみ使う
X_Background = (
    landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))][features]
    .replace([np.inf, -np.inf], np.nan).dropna()
    .sample(n=min(10000, len(landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))])), random_state=0)
)
X_Background_s_fixed = scaler.transform(X_Background.values)

# OOF予測・スコアを格納する配列を用意
oof_scores = np.zeros(len(X_labeled_normal))  # decision_function値（採用する）
# 背景データに対するスコアも、foldごとのモデルで計算し後で平均する
# （Presence側と同じ「学習に使われていないモデル」で評価しスケールを揃えるため）
Background_scores_per_fold = np.zeros((5, len(X_Background_s_fixed)))

kf = KFold(n_splits=5, shuffle=True, random_state=0)

for fold_i, (train_idx, test_idx) in enumerate(kf.split(X_labeled_normal)):
    # 学習データ（そのfoldのtest対象を除いた在データ）を標準化
    X_train_s = scaler.transform(X_labeled_normal[train_idx])
    X_test_s = scaler.transform(X_labeled_normal[test_idx])

    # 「テスト対象の店舗を含まない」データだけでモデルを学習（リーク防止）
    m = OneClassSVM(kernel="rbf", gamma=gamma, nu=nu, max_iter=-1)
    m.fit(X_train_s)

    # テスト対象の店舗（未学習データ扱い）に対する予測・スコアを記録
    oof_scores[test_idx] = m.decision_function(X_test_s)  # 値が大きいほど正常

    # 同じfoldモデル（学習に使われていない視点）で背景データも評価
    # → Presenceと同じ「厳しさ」のモデルで背景を測ることでスケールを揃える
    Background_scores_per_fold[fold_i] = m.decision_function(X_Background_s_fixed)

# 5つのfoldモデルによる背景スコアを平均し、Presence(oof_scores)と比較可能な基準値にする
Background_scores_fixed = Background_scores_per_fold.mean(axis=0)

# Presence(在データ)・Background(背景データ)のスコアとラベルを結合
scores_all = np.concatenate([oof_scores, Background_scores_fixed])
labels_all = np.concatenate([np.ones(len(oof_scores)), np.zeros(len(Background_scores_fixed))])

# --- AUC・最適閾値の計算 ---
print('ROC-AUC', f'{roc_auc_score(labels_all, scores_all):.4f}')
fpr, tpr, thresholds = roc_curve(labels_all, scores_all)
optimal_idx = np.argmax(tpr - fpr)  # 感度+特異度が最大になる点
optimal_threshold = thresholds[optimal_idx]
print(f"最適な異常判定閾値(OOF基準・背景データ対比): {optimal_threshold:.4f}")


# --- 最終的な閾値決定: optimal_thresholdベースでラベル付け(リークなし・基準統一) ---
pred_oof = np.where(oof_scores >= optimal_threshold, 1, -1)
all_coords_with_luse_log.loc[mask, "ocsvm_label"] = pred_oof
all_coords_with_luse_log.loc[mask, "ocsvm_score"] = oof_scores
plot_label_ratio_pie(pred_oof)

image.png

精度確認

バックグラウンドデータのラベル1と-1でどのくらい特徴量に差があるのか確認する。

データの用意。

# ============================================================
# 検証: 背景データ（既存カフェ周辺メッシュ）に最終モデルを適用
#      → 「都市 vs 郊外」のような粗い違いではなく、
#        既存カフェとその近傍という似た環境内でモデルが
#        微妙な立地差をどれだけ判別できているかを確認する
# ============================================================
backg_df = landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))][areas_cols+[selected_flow]].copy()
backg_df_val = backg_df.values

# inf/nanを含む行を除外（マスクで有効行のみ後段の処理対象にする）
mask = np.isfinite(backg_df_val).all(axis=1)
backg_df_val = backg_df_val[mask]

# ステップ6で在データ全体でfitしたscalerを再利用（学習時と同じスケール基準で変換）
backg_df_scaled = scaler.transform(backg_df_val)

# 最終モデル(ocsvm)とOOFで決めたoptimal_thresholdで、背景データを判定
scores_bg = ocsvm.decision_function(backg_df_scaled)
pred_labels_bg = np.where(scores_bg >= optimal_threshold, 1, -1)

backg_df.loc[mask, "ocsvm_label"] = pred_labels_bg
backg_df.loc[mask, "ocsvm_score"] = scores_bg

plot_label_ratio_pie(pred_labels_bg)
display(backg_df)

439件の候補地がある。
image.png

既存カフェ（all_coords_with_luse_log）、OCSVMで正常:立地候補（ocsvm_label==1）と判定されたメッシュ、異常:候補外（ocsvm_label!=1）と判定されたメッシュの3種類のデータの特徴量の分布をバイオリンプロットにより確認。

# ============================================================
# 8. 既存cafe vs 立地候補(正常) vs 候補外(異常) の面積分布比較（violinplot）
# ============================================================

violin_df = pd.concat([
    pd.DataFrame({"cat": np.full(len(all_coords_with_luse_log), c), "area": all_coords_with_luse_log[c]})
    for c in areas_cols
]).reset_index(drop=True)

violin_df_test = pd.concat([
    pd.DataFrame({
        "cat": np.full(len(backg_df[backg_df["ocsvm_label"] == 1]), c),
        "area": backg_df[backg_df["ocsvm_label"] == 1][c],
    }) for c in areas_cols
]).reset_index(drop=True)

violin_df_test_except = pd.concat([
    pd.DataFrame({
        "cat": np.full(len(backg_df[backg_df["ocsvm_label"] != 1]), c),
        "area": backg_df[backg_df["ocsvm_label"] != 1][c],
    }) for c in areas_cols
]).reset_index(drop=True)

plt.rcParams["font.family"] = prop.get_name()
fig = plt.figure(figsize=(15, 10))

# 上段: 既存cafe(青) vs 立地候補(緑)
ax = fig.add_subplot(2, 1, 1)
violin_df_all = pd.concat([violin_df.assign(group="立地"), violin_df_test.assign(group="立地候補")])
sns.violinplot(
    data=violin_df_all, x="cat", y="area", hue="group", split=True, inner="quart",
    palette={"立地": "steelblue", "立地候補": "seagreen", "候補外": "gold"},
    density_norm="width", ax=ax,
)
ax.set_xlabel("土地利用カテゴリ"); ax.set_ylabel("面積")
ax.set_title("土地利用カテゴリ別面積分布（青:立地 / 緑:立地候補）")
ax.legend(title="グループ", loc="upper right")

# 下段: 既存cafe(青) vs 候補外(金)
ax = fig.add_subplot(2, 1, 2)
violin_df_all = pd.concat([violin_df.assign(group="立地"), violin_df_test_except.assign(group="候補外")])
sns.violinplot(
    data=violin_df_all, x="cat", y="area", hue="group", split=True, inner="quart",
    palette={"立地": "steelblue", "立地候補": "seagreen", "候補外": "gold"},
    density_norm="width", ax=ax,
)
ax.set_xlabel("土地利用カテゴリ"); ax.set_ylabel("面積")
ax.set_title("土地利用カテゴリ別面積分布（青:立地 / 金:候補外）")
ax.legend(title="グループ", loc="upper right")

plt.tight_layout()
plt.show()

上段では既存カフェ（青）と立地候補（緑）を、下段では既存カフェ（青）と候補外（金）を、土地利用カテゴリごとにsplit=Trueのviolinplotで左右に分割表示し、分布の形状を直接比較できるようにしている。青の形状に近い＝既存カフェの立地の特徴に近いと言える。ラベル1=立地候補のメッシュの方が既存カフェの特徴に近いことが確認できる。
image.png


同様に、人流の傾向。 
violin_df = pd.concat([
    pd.DataFrame({"cat": np.full(len(all_coords_with_luse_log), c), "area": all_coords_with_luse_log[c]})
    for c in [selected_flow]
]).reset_index(drop=True)

violin_df_test = pd.concat([
    pd.DataFrame({
        "cat": np.full(len(backg_df[backg_df["ocsvm_label"] == 1]), c),
        "area": backg_df[backg_df["ocsvm_label"] == 1][c],
    }) for c in [selected_flow]
]).reset_index(drop=True)

violin_df_test_except = pd.concat([
    pd.DataFrame({
        "cat": np.full(len(backg_df[backg_df["ocsvm_label"] != 1]), c),
        "area": backg_df[backg_df["ocsvm_label"] != 1][c],
    }) for c in [selected_flow]
]).reset_index(drop=True)

plt.rcParams["font.family"] = prop.get_name()
fig = plt.figure(figsize=(15, 10))

# 上段: 既存cafe(青) vs 立地候補(緑)
ax = fig.add_subplot(1, 2, 1)
violin_df_all = pd.concat([violin_df.assign(group="立地"), violin_df_test.assign(group="立地候補")])
sns.violinplot(
    data=violin_df_all, x="cat", y="area", hue="group", split=True, inner="quart",
    palette={"立地": "steelblue", "立地候補": "seagreen", "候補外": "gold"},
    density_norm="width", ax=ax,
)
ax.set_xlabel("年(統計量)"); ax.set_ylabel("人流[人]")
ax.set_title("人流分布（青:立地 / 緑:立地候補）")
ax.legend(title="グループ", loc="upper right")

# 下段: 既存cafe(青) vs 候補外(金)
ax = fig.add_subplot(1, 2, 2)
violin_df_all = pd.concat([violin_df.assign(group="立地"), violin_df_test_except.assign(group="候補外")])
sns.violinplot(
    data=violin_df_all, x="cat", y="area", hue="group", split=True, inner="quart",
    palette={"立地": "steelblue", "立地候補": "seagreen", "候補外": "gold"},
    density_norm="width", ax=ax,
)
ax.set_xlabel("年(統計量)"); ax.set_ylabel("人流[人]")
ax.set_title("人流分布（青:立地 / 金:候補外）")
ax.legend(title="グループ", loc="upper right")

plt.tight_layout()
plt.show()

ラベル1=立地候補のメッシュの方が既存カフェの特徴に近いことが確認できる。（とはいってもやはり既存カフェの人流とは差がある。）
image.png


土地利用カテゴリごとに既存カフェ分布との近さをWasserstein距離という分布間距離で定量化し比較する。
土地利用面積列と人流の分布を対象に既存カフェと、立地候補・候補外のWasserstein距離を計算。この距離は値が小さいほど両者の分布が近いことを意味し、緑（既存カフェvs立地候補）と金（既存カフェvs候補外）に分けて算出している。

# ============================================================
# 9. カテゴリ別 Wasserstein距離: 青 vs 緑（立地候補）/ 青 vs 金（候補外）
# ============================================================

violin_df = pd.concat([
    pd.DataFrame({"cat": np.full(len(all_coords_with_luse_log), c), "area": all_coords_with_luse_log[c]})
    for c in areas_cols+[selected_flow]
]).reset_index(drop=True)

violin_df_test = pd.concat([
    pd.DataFrame({
        "cat": np.full(len(backg_df[backg_df["ocsvm_label"] == 1]), c),
        "area": backg_df[backg_df["ocsvm_label"] == 1][c],
    }) for c in areas_cols+[selected_flow]
]).reset_index(drop=True)

violin_df_test_except = pd.concat([
    pd.DataFrame({
        "cat": np.full(len(backg_df[backg_df["ocsvm_label"] != 1]), c),
        "area": backg_df[backg_df["ocsvm_label"] != 1][c],
    }) for c in areas_cols+[selected_flow]
]).reset_index(drop=True)

cats = violin_df["cat"].unique()
dist_per_cat_green, dist_per_cat_gold = {}, {}

for c in cats:
    x_blue = violin_df.loc[violin_df["cat"] == c, "area"].dropna().values
    x_green = violin_df_test.loc[violin_df_test["cat"] == c, "area"].dropna().values
    x_gold = violin_df_test_except.loc[violin_df_test_except["cat"] == c, "area"].dropna().values

    dist_per_cat_green[c] = wasserstein_distance(x_blue, x_green)
    dist_per_cat_gold[c] = wasserstein_distance(x_blue, x_gold)

df_dist = pd.DataFrame({
    "cat": list(dist_per_cat_green.keys()),
    "W_bl_green": list(dist_per_cat_green.values()),
    "W_bl_gold": list(dist_per_cat_gold.values()),
})

plt.rcParams["font.family"] = prop.get_name()
fig, ax = plt.subplots(figsize=(10, 4))

df_long = df_dist.melt(id_vars="cat", value_vars=["W_bl_green", "W_bl_gold"],
                        var_name="group", value_name="distance")
sns.barplot(
    data=df_long, x="cat", y="distance", hue="group", ax=ax,
    palette={"W_bl_green": "green", "W_bl_gold": "gold"},
)
ax.set_xlabel("土地利用カテゴリ"); ax.set_ylabel("Wasserstein距離")
ax.set_title("青分布とのWasserstein距離（緑:立地候補 / 金:候補外）")

handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, ["青 vs 緑", "青 vs 金"], title="比較ペア", loc="upper right")
plt.tight_layout()
plt.show()

立地候補のメッシュの方が既存カフェの特徴に近いことが確認できる。
image.png


Wasserstein距離を多次元に拡張し、既存カフェ・立地候補・候補外の3グループ間の分布の近さを、複数特徴量を同時に考慮した形で定量化。
土地利用面積と人流の特徴量をStandardScalerにより標準化した上で、既存カフェと立地候補、既存カフェと候補外のペアについて多次元Wasserstein距離を500回のランダムサンプリングで繰り返し計算し、偶然の偏りを抑えた平均・標準偏差を求めている。
最後に500回分の距離をKDEで可視化し、立地候補が既存カフェの多次元的な特徴パターンに候補外よりも近いかを、特徴量全体の組み合わせとして確認している。

# ============================================================
# 9.5. 多次元Wasserstein距離（サンプリング500回×KDEで比較）
# ============================================================

features = areas_cols + [selected_flow]

# 前処理: StandardScaler で標準化してから距離計算（スケールの偏りを防ぐため）
X_blue = all_coords_with_luse_log[features].values
X_blue = X_blue[np.isfinite(X_blue).all(axis=1)]
scaler_wnd = StandardScaler().fit(X_blue)
X_blue_scaled = scaler_wnd.transform(X_blue)

X_green = backg_df[backg_df['ocsvm_label'] == 1][features].values
X_green = X_green[np.isfinite(X_green).all(axis=1)]
X_green_scaled = scaler_wnd.transform(X_green)

X_gold = backg_df[backg_df['ocsvm_label'] != 1][features].values
X_gold = X_gold[np.isfinite(X_gold).all(axis=1)]
X_gold_scaled = scaler_wnd.transform(X_gold)

n_sample = min(len(X_blue_scaled), len(X_green_scaled), len(X_gold_scaled), 1000)

wnd_green_list = []
wnd_gold_list = []

iter = 500
for ix in tqdm(range(iter)):
    idx_b = np.random.choice(len(X_blue_scaled), n_sample, replace=False)
    idx_g = np.random.choice(len(X_green_scaled), n_sample, replace=False)
    idx_gold = np.random.choice(len(X_gold_scaled), n_sample, replace=False)

    wnd_green_list.append(
        wasserstein_distance_nd(X_blue_scaled[idx_b], X_green_scaled[idx_g])
    )
    wnd_gold_list.append(
        wasserstein_distance_nd(X_blue_scaled[idx_b], X_gold_scaled[idx_gold])
    )

wnd_green_arr = np.array(wnd_green_list)
wnd_gold_arr = np.array(wnd_gold_list)

print("青-緑 多次元Wasserstein距離: mean=%.4f, std=%.4f" % (wnd_green_arr.mean(), wnd_green_arr.std()))
print("青-金 多次元Wasserstein距離: mean=%.4f, std=%.4f" % (wnd_gold_arr.mean(), wnd_gold_arr.std()))

df_wnd = pd.DataFrame({
    'iteration': range(iter),
    'W_nd_green': wnd_green_arr,
    'W_nd_gold': wnd_gold_arr,
})

display(df_wnd)

fig, ax = plt.subplots(figsize=(8, 6))
sns.kdeplot(wnd_green_arr, ax=ax, color='green', label='青 vs 緑', fill=True, alpha=0.3)
sns.kdeplot(wnd_gold_arr, ax=ax, color='goldenrod', label='青 vs 金', fill=True, alpha=0.3)
ax.set_xlabel('多次元Wasserstein距離')
ax.set_ylabel('density（確率密度）')
ax.set_title('サンプリングによる多次元Wasserstein距離の分布')
ax.legend()
plt.tight_layout()
plt.show()

# >> 青-緑 多次元Wasserstein距離: mean=2.0840, std=0.0852
# >> 青-金 多次元Wasserstein距離: mean=6.4242, std=0.3796

立地候補のメッシュの方が既存カフェの特徴に近いことが確認できる。
image.png


特徴量空間でのユークリッド距離でも既存カフェとの類似性を確認する。

既存カフェ・立地候補・候補外の3グループから同じ件数をランダムに取り出し、既存カフェの各点と立地候補の各点、既存カフェの各点と候補外の各点、というように全ての点同士の組み合わせでユークリッド距離を計算する処理を200回繰り返す。
前段の多次元Wasserstein距離では「既存カフェの分布」と「立地候補の分布」を比べて、その差を1つの数値に要約していた。一方この処理では、そうした要約はせず、点と点の距離を1つずつすべて求め、その距離の値がどのように分布しているか（近い距離が多いのか、遠い距離も多く混じっているのか）をKDEで描いている。

# ============================================================
# 10. 多次元空間でのユークリッド距離分布（サンプリング200回×KDEで比較）
# ============================================================
green_dist, gold_dist = {}, {}

X_blue = all_coords_with_luse_log[areas_cols+[selected_flow]].values
X_blue = X_blue[np.isfinite(X_blue).all(axis=1)]
scaler_dist = StandardScaler().fit(X_blue)
X_blue_scaled = scaler_dist.transform(X_blue)

X_green = backg_df[backg_df["ocsvm_label"] == 1][areas_cols+[selected_flow]].values
X_green = X_green[np.isfinite(X_green).all(axis=1)]
X_green_scaled = scaler_dist.transform(X_green)

X_gold = backg_df[backg_df["ocsvm_label"] != 1][areas_cols+[selected_flow]].values
X_gold = X_gold[np.isfinite(X_gold).all(axis=1)]
X_gold_scaled = scaler_dist.transform(X_gold)

n_sample = min(len(X_blue_scaled), len(X_green_scaled), len(X_gold_scaled), 1000)

for ix in range(200):
    idx_b = np.random.choice(len(X_blue_scaled), n_sample, replace=False)
    idx_g = np.random.choice(len(X_green_scaled), n_sample, replace=False)
    idx_gold = np.random.choice(len(X_gold_scaled), n_sample, replace=False)

    green_dist[ix] = pairwise_distances(X_blue_scaled[idx_b], X_green_scaled[idx_g], metric="euclidean").ravel()
    gold_dist[ix] = pairwise_distances(X_blue_scaled[idx_b], X_gold_scaled[idx_gold], metric="euclidean").ravel()

fig, ax = plt.subplots(figsize=(8, 6))

# 緑（立地候補）の距離分布を緑系グラデーションでKDE描画
cmap_g = plt.cm.Greens
colors_g = [cmap_g(i / max(len(green_dist) - 1, 1)) for i in range(len(green_dist))]
for i, (k, v) in enumerate(green_dist.items()):
    sns.kdeplot(v, ax=ax, color=colors_g[i], fill=False, alpha=0.2, linewidth=2)

# 金（候補外）の距離分布を暖色系グラデーションでKDE描画（同じaxに重ねる）
cmap_gold = plt.cm.autumn
colors_gold = [cmap_gold(i / max(len(gold_dist) - 1, 1)) for i in range(len(gold_dist))]
for i, (k, v) in enumerate(gold_dist.items()):
    sns.kdeplot(v, ax=ax, color=colors_gold[i], fill=False, alpha=0.2, linewidth=2)

ax.set_xlabel("euclidean距離")
ax.set_ylabel("density（確率密度）")
ax.set_title("青-緑・青-金間のユークリッド距離分布（200回サンプリング）")
plt.tight_layout()
plt.show()

立地候補のメッシュの方が既存カフェの特徴に近いことが確認できる。
image.png

OCSVMモデルを全体に適用

学習済みOCSVMを適用する対象範囲として、鳥取県を含む中国地方5県の土地利用ポリゴン全体に人流データを結合し、対数変換まで済ませた予測用データセットを準備する。
まず都道府県境界データ（prefectures.geojson）と土地利用データをsjoinで空間結合し、鳥取・島根・岡山・広島・山口の5県に絞り込んだうえで、各ポリゴンの代表点座標から1kmメッシュコードを算出し、人流データ（personflow_d_agg）とメッシュ単位で結合。
最後に、土地利用面積列と人流の代表値にlog1p変換を適用し、学習時に在データへ施した処理と同じ形式に整えることで、次のステップで学習済みモデルにそのまま入力できる予測用データ（landuse_merged_log）を作成。

# ============================================================
# 7. 中国地方の土地利用ポリゴン全体に、学習済みモデルを適用（立地候補の判定）
# ============================================================
pref = gpd.read_file("../../basic_data/prefectures.geojson")

landuse_merged_pref = landuse_merged.sjoin(pref[["name", "geometry"]], how="inner").drop(columns=["index_right"])
landuse_merged_pref = landuse_merged_pref[
    landuse_merged_pref["name"].isin(["鳥取県", "島根県", "岡山県", "広島県", "山口県"])
].copy()

# 人流データと結合
landuse_merged_pref['mesh1kmid'] = landuse_merged_pref.geometry.representative_point().map(lambda geo: ju.to_meshcode(geo.y, geo.x, 3))
landuse_merged_pref = landuse_merged_pref.merge(personflow_d_agg, on=['mesh1kmid'], how='inner')

landuse_merged_log = landuse_merged_pref.copy()
landuse_merged_log[areas_cols+[selected_flow]] = np.log1p(landuse_merged_log[areas_cols+[selected_flow]])

display(landuse_merged_pref)
display(landuse_merged_log)

image.png

これまでの検証で決めたgamma・nuと最適閾値を使い、最終的なOCSVMモデルを在データ全件で学習し、中国地方全域のメッシュに対してカフェ立地の適性を予測する。
予測対象となる中国地方全域の特徴量は、学習時に使ったスケーラーで変換することによりスケール基準を統一。
判定にはモデル内部のnuベース閾値ではなく、計算した前述のoptimal_thresholdを使ってdecision_functionのスコアを二値化している。
予測ラベルとスコアをlanduse_merged_logに書き込み、円グラフで正常(＝カフェらしい)/異常(=カフェ立地候補ではない)の比率を可視化することで、中国地方全域における「実在しないが立地条件的にカフェらしい」候補メッシュを最終的に抽出している。

# 最終モデルは全在データで学習（もう検証目的ではないため、全件を使ってよい）
X_train_s = scaler.transform(X_labeled_normal)
ocsvm = OneClassSVM(kernel="rbf", gamma=gamma, nu=nu, max_iter=-1)
ocsvm.fit(X_train_s)

# 中国地方全域メッシュ（候補地点）に対して適性を予測
X_test = landuse_merged_log[areas_cols+[selected_flow]].values
mask = np.isfinite(X_test).all(axis=1)
X_test = X_test[mask]
X_test_scaled = scaler.transform(X_test)  # 既存cafeでfitしたscalerを再利用（スケール基準を統一）

# pred_labels = ocsvm.predict(X_test_scaled)
scores = ocsvm.decision_function(X_test_scaled)
pred_labels = np.where(scores >= optimal_threshold, 1, -1)

landuse_merged_pref.loc[mask, "ocsvm_label"] = pred_labels
landuse_merged_pref.loc[mask, "ocsvm_score"] = scores
landuse_merged_log.loc[mask, "ocsvm_label"] = pred_labels
landuse_merged_log.loc[mask, "ocsvm_score"] = scores

plot_label_ratio_pie(pred_labels)
display(landuse_merged_log)

最終的に688メッシュの立地候補が抽出されたことになる。
image.png

立地候補地図可視化

精度の確認で、既存カフェの特徴に近いメッシュが抽出できていることがわかった。
では中国地方全体と、鳥取県全体で立地候補の場所を確認してみる。
すでにモデルを全体に適用して、688メッシュの立地候補は出せているので、それを可視化する。

# ============================================================
# 11. 地図描画: 立地候補スコアと既存cafeを衛星写真ベースマップ上にプロット
# ============================================================
def vec_cbar(ax):
    """cartopy軸にカラーバー用の通常軸(cax)を追加する関数"""
    divider = mpl_toolkits.axes_grid1.make_axes_locatable(ax)
    cax = divider.append_axes("right", "4%", pad=0.1, axes_class=plt.Axes)  # axes_class を明示してKeyError回避
    return cax


except_mesh = all_coords_with_luse_log[['geometry']].sjoin(landuse_merged_log[['メッシュ', 'geometry']].drop_duplicates(), 
                                                           how='inner'
                                                          ).drop(columns=['index_right'])['メッシュ'].unique()
landuse_merged_log_point = landuse_merged_log.copy()
landuse_merged_log_point = landuse_merged_log_point[~landuse_merged_log_point['メッシュ'].isin(except_mesh)]
landuse_merged_log_point.geometry = landuse_merged_log_point.geometry.representative_point()  # Polygon -> Point

lon_min_broad, lat_min_broad, lon_max_broad, lat_max_broad = landuse_merged_pref.total_bounds

plt.rcParams["font.family"] = prop.get_name()
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.grid(False)
ax.set_extent(
    [lon_min_broad - 0.05, lon_max_broad + 0.05, lat_min_broad - 0.05, lat_max_broad + 0.05],
    crs=ccrs.PlateCarree(),
)

cx.add_basemap(ax, crs=ccrs.PlateCarree(), source=cx.providers.Esri.WorldImagery, zoom=10, zorder=1, alpha=0.8)

cax = vec_cbar(ax)

# 立地候補（正常判定）のスコアを色の濃淡で表示
landuse_merged_log_point[landuse_merged_log_point["ocsvm_label"] == 1].plot(
    column="ocsvm_score", ax=ax, legend=True, cmap="Reds",
    marker=".", markersize=5 ** 2, linewidth=0, cax=cax, zorder=4, alpha=0.9,
)

# 既存cafeをシアン×黒フチのXマークで強調表示
all_coords_with_luse_log.plot(
    color="c", ax=ax, legend=True, markersize=10 ** 2,
    marker="X", linewidth=1.2, edgecolor="k", alpha=1.0, zorder=3,
)

# 県境を薄いグレーの細線で表示（ノイズを抑える）
pref.plot(ax=ax, facecolor="None", linewidth=0.6, edgecolor="lightgray", alpha=0.7, zorder=2)

lat_center = (lat_min_broad + lat_max_broad) / 2
dx = 111320 * np.cos(np.radians(lat_center))  # 緯度に応じた1度あたりの距離(m)を近似計算

scalebar = ScaleBar(
    dx=dx, units="m", length_fraction=0.2, location="lower right",
    box_alpha=0.7, rotation="horizontal-only", font_properties={"size": 10},
)
ax.add_artist(scalebar)

ax.set_aspect("equal")
plt.show()

青×が既存カフェ、赤い点が立地候補。赤が濃いほど既存カフェと類似していることを示している。
やはり既存カフェが立地している市街地に候補地が多く出る。これは仕方ないが、既存カフェとの距離とかを特徴量に入れたら改善するかもしれない。
ただ、既存カフェが立地していない場所にも赤い点があり、こういったところは今から出店するチャンスかもしれない。
image.png

鳥取県に絞ってみる。

# ============================================================
# 12. 鳥取地図描画: 立地候補スコアと既存cafeを衛星写真ベースマップ上にプロット
# ============================================================

landuse_merged_log_point = landuse_merged_log.copy()
landuse_merged_log_point = landuse_merged_log_point[~landuse_merged_log_point['メッシュ'].isin(except_mesh)]
landuse_merged_log_point.geometry = landuse_merged_log_point.geometry.representative_point()  # Polygon -> Point

lon_min_broad, lat_min_broad, lon_max_broad, lat_max_broad = landuse_merged_pref[landuse_merged_pref['name']=='鳥取県'].total_bounds

plt.rcParams["font.family"] = prop.get_name()
fig = plt.figure(figsize=(14, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.grid(False)
ax.set_extent(
    [lon_min_broad - 0.05, lon_max_broad + 0.05, lat_min_broad - 0.05, lat_max_broad + 0.05],
    crs=ccrs.PlateCarree(),
)

cx.add_basemap(ax, crs=ccrs.PlateCarree(), source=cx.providers.Esri.WorldImagery, zoom=11, zorder=1, alpha=0.8)

cax = vec_cbar(ax)

# 立地候補（正常判定）のスコアを色の濃淡で表示
landuse_merged_log_point[landuse_merged_log_point["ocsvm_label"] == 1].plot(
    column="ocsvm_score", ax=ax, legend=True, cmap="Reds",
    marker="o", markersize=6 ** 2, linewidth=0, cax=cax, zorder=4, alpha=0.9,
)

# 既存cafeをシアン×黒フチのXマークで強調表示
all_coords_with_luse_log.plot(
    color="c", ax=ax, legend=True, markersize=10 ** 2,
    marker="X", linewidth=1.2, edgecolor="k", alpha=1.0, zorder=3,
)

# 県境を薄いグレーの細線で表示（ノイズを抑える）
pref.plot(ax=ax, facecolor="None", linewidth=0.6, edgecolor="lightgray", alpha=0.7, zorder=2)


lat_center = (lat_min_broad + lat_max_broad) / 2
dx = 111320 * np.cos(np.radians(lat_center))  # 緯度に応じた1度あたりの距離(m)を近似計算

scalebar = ScaleBar(
    dx=dx, units="m", length_fraction=0.2, location="lower right",
    box_alpha=0.7, rotation="horizontal-only", font_properties={"size": 10},
)
ax.add_artist(scalebar)

ax.set_aspect("equal")
plt.show()

鳥取県の中央部と、最北西端あたりは既存カフェもなく、立地候補なのでは！？実際にGoogle Mapで見てもこのあたりにスタバは無かった。
image.png


拡大。

landuse_merged_log_point = landuse_merged_log.copy()
landuse_merged_log_point = landuse_merged_log_point[~landuse_merged_log_point['メッシュ'].isin(except_mesh)]
landuse_merged_log_point.geometry = landuse_merged_log_point.geometry.representative_point()  # Polygon -> Point

row_geometry_y, row_geometry_x = 35.454166665,133.84375

lon_min_broad, lat_min_broad, lon_max_broad, lat_max_broad = row_geometry_x-0.001, row_geometry_y-0.001, row_geometry_x+0.001, row_geometry_y+0.001

plt.rcParams["font.family"] = prop.get_name()
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.grid(False)
ax.set_extent(
    [lon_min_broad - 0.05, lon_max_broad + 0.05, lat_min_broad - 0.05, lat_max_broad + 0.05],
    crs=ccrs.PlateCarree(),
)

cx.add_basemap(ax, crs=ccrs.PlateCarree(), source=cx.providers.Esri.WorldImagery, zoom=14, zorder=1, alpha=0.8)

cax = vec_cbar(ax)

# 立地候補（正常判定）のスコアを色の濃淡で表示
landuse_merged_log_point[landuse_merged_log_point["ocsvm_label"] == 1].plot(
    column="ocsvm_score", ax=ax, legend=True, cmap="Reds",
    marker="o", markersize=14 ** 2, linewidth=1, edgecolor='k', cax=cax, zorder=4, alpha=0.9,
)

# 既存cafeをシアン×黒フチのXマークで強調表示
all_coords_with_luse_log.plot(
    color="c", ax=ax, legend=True, markersize=12 ** 2,
    marker="X", linewidth=1.2, edgecolor="k", alpha=1.0, zorder=3,
)

# 県境を薄いグレーの細線で表示（ノイズを抑える）
pref.plot(ax=ax, facecolor="None", linewidth=0.6, edgecolor="lightgray", alpha=0.7, zorder=2)


lat_center = (lat_min_broad + lat_max_broad) / 2
dx = 111320 * np.cos(np.radians(lat_center))  # 緯度に応じた1度あたりの距離(m)を近似計算

scalebar = ScaleBar(
    dx=dx, units="m", length_fraction=0.2, location="lower right",
    box_alpha=0.7, rotation="horizontal-only", font_properties={"size": 10},
)
ax.add_artist(scalebar)

ax.set_aspect("equal")
plt.show()


landuse_merged_log_point = landuse_merged_log.copy()
landuse_merged_log_point = landuse_merged_log_point[~landuse_merged_log_point['メッシュ'].isin(except_mesh)]
landuse_merged_log_point.geometry = landuse_merged_log_point.geometry.representative_point()  # Polygon -> Point

row_geometry_y, row_geometry_x = 35.545833335,133.24375

lon_min_broad, lat_min_broad, lon_max_broad, lat_max_broad = row_geometry_x-0.001, row_geometry_y-0.001, row_geometry_x+0.001, row_geometry_y+0.001

plt.rcParams["font.family"] = prop.get_name()
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.grid(False)
ax.set_extent(
    [lon_min_broad - 0.05, lon_max_broad + 0.05, lat_min_broad - 0.05, lat_max_broad + 0.05],
    crs=ccrs.PlateCarree(),
)

cx.add_basemap(ax, crs=ccrs.PlateCarree(), source=cx.providers.Esri.WorldImagery, zoom=14, zorder=1, alpha=0.8)

cax = vec_cbar(ax)

# 立地候補（正常判定）のスコアを色の濃淡で表示
landuse_merged_log_point[landuse_merged_log_point["ocsvm_label"] == 1].plot(
    column="ocsvm_score", ax=ax, legend=True, cmap="Reds",
    marker="o", markersize=14 ** 2, linewidth=1, edgecolor='k', cax=cax, zorder=4, alpha=0.9,
)

# 既存cafeをシアン×黒フチのXマークで強調表示
all_coords_with_luse_log.plot(
    color="c", ax=ax, legend=True, markersize=12 ** 2,
    marker="X", linewidth=1.2, edgecolor="k", alpha=1.0, zorder=3,
)

# 県境を薄いグレーの細線で表示（ノイズを抑える）
pref.plot(ax=ax, facecolor="None", linewidth=0.6, edgecolor="lightgray", alpha=0.7, zorder=2)


lat_center = (lat_min_broad + lat_max_broad) / 2
dx = 111320 * np.cos(np.radians(lat_center))  # 緯度に応じた1度あたりの距離(m)を近似計算

scalebar = ScaleBar(
    dx=dx, units="m", length_fraction=0.2, location="lower right",
    box_alpha=0.7, rotation="horizontal-only", font_properties={"size": 10},
)
ax.add_artist(scalebar)

ax.set_aspect("equal")
plt.show()

鳥取県中央部。鳥取県倉吉市市街地の天神川沿い。
image.png

鳥取県最北西端。鳥取県境港市。
image.png

おわりに

というわけで、鳥取県のどこに"隠れスタバ適地"があるかはモデル上では判明した。あとはスターバックスさんの管轄である。
ちなみに私は鳥取県に行ったことはない。（←え？）
機会があれば今回の立地適地に行ってみて街並みを確認するのも面白いかもしれない。

以上！

おまけ

Isolation Forest

Isolation Forestでも同様にやってみた。以下のようにn_estimators、max_samplesのグリッドサーチを実施。（contaminationもやるべきだったか…）

# ============================================================
# (おまけ) IsolationForest: n_estimators / max_samples のグリッドサーチ（AUCで評価）
#   OCSVMと同じ在データ・背景データ・評価方針を使い、
#   別の異常検知手法でも同様の性能が出るかを比較検証する
# ============================================================

# 在データ・背景データを結合してラベル付け
features = areas_cols + [selected_flow]

# --- 在データ（Presence）: 既存スタバ店舗の特徴量 ---
X_Presence = (
    all_coords_with_luse_log[features]
    .replace([np.inf, -np.inf], np.nan).dropna()  # inf/nanを除外
)
# --- 背景データ（Background）: 在データ周辺メッシュ（include_meshes）からサンプリング ---
#   OCSVMのときと同じ考え方で、都市/郊外の粗い違いではなく
#   近傍内の微妙な立地差を評価対象にする
X_Background = (
    landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))][features]
    .replace([np.inf, -np.inf], np.nan).dropna()
    .sample(n=min(10000, len(landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))])), random_state=0)
    # 背景データが多すぎる場合に備え、最大10000件にランダムサンプリングして件数を制限
)

# 在データ・背景データを結合し、1=在データ, 0=背景データのラベルを付与
x_all = pd.concat([X_Presence, X_Background], ignore_index=True)
y_all = np.concatenate([np.ones(len(X_Presence)), np.zeros(len(X_Background))])


def evaluate_iforest_auc(x, y, n_estimators, max_samples, n_splits=5, random_state=0):
    """
    OCSVMのevaluate_ocsvm_aucと同じ考え方で、
    在データ・背景データを層化K-Foldで分割し、
    「学習に使っていないfold」に対するAUCで汎化性能を評価する関数。
    IsolationForestも教師なし学習の異常検知手法のため、学習(fit)は
    各foldの「在データ(正常, y=1)のみ」を使う。
    """
    skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=random_state)
    aucs = []
    for train_idx, test_idx in skf.split(x, y):
        # 学習は「在データ（正常）のみ」で行う（IsolationForestは教師なしのため）
        train_Presence_idx = train_idx[y[train_idx] == 1]

        model = IsolationForest(
            n_estimators=n_estimators,
            max_samples=max_samples,
            random_state=random_state,
        )
        model.fit(x.values[train_Presence_idx])

        # decision_function: 値が大きいほど「正常（在データらしい）」
        # OCSVMと同じくスコアの「順序付け能力」自体をAUCで評価する
        y_score = model.decision_function(x.values[test_idx])
        aucs.append(roc_auc_score(y[test_idx], y_score))
    return np.mean(aucs), np.std(aucs)


# ハイパーパラメータの探索範囲
# n_estimators: 木の数、max_samples: 各木の学習に使うサンプル数/割合
results = []
for n_estimators in [100, 200, 300]:
    for max_samples in ["auto", 0.5, 0.8]:
        auc_mean, auc_std = evaluate_iforest_auc(x_all, y_all, n_estimators, max_samples)
        results.append({
            "n_estimators": n_estimators,
            "max_samples": max_samples,
            "auc_mean": auc_mean,
            "auc_std": auc_std,
        })

# AUC平均が高い順に並べ、OCSVMのAUCと比較できるようにする
results_df = pd.DataFrame(results).sort_values("auc_mean", ascending=False)
display(results_df)

# グリッドサーチ結果から選択したハイパーパラメータ
# n_estimators, max_samples = 300, 0.5

候補地の結果は以下。まあだいたいOCSVMと同じような感じ。
image.png
image.png
image.png
image.png

MaxEnt

今回実施した「バックグラウンドデータを活用する」というアイデアは、MaxEnt（最大エントロピー法によるSDM: Species Distribution Model）で使われる手法を参考にしたものである。

MaxEntは、真の「不在（Absence）」データが得られない生物種の分布データに対応できる手法である。生物の調査では、「その種を観測できなかった」ことが必ずしも「その種がそこにいない」ことを意味しない（探索不足や、見つけにくい種である可能性があるため）。そのため、基本的には在データ（Presence）しか得られないことが多い。

真の不在データがない代わりに、MaxEntでは調査対象領域全体からランダムに点をサンプリングし、それを「Background（背景）」として扱う。そして、Presence地点とBackground地点の環境条件（特徴量）の違いから、その種にとっての生息適地を統計的に推定する。

今回のOCSVMによるアプローチとMaxEntによるアプローチは、Presence-Background設計という点では共通するが、モデルの学習方法は異なる。

MaxEntは、PresenceとBackgroundを比較し、「どの環境でその種が背景より相対的に多く現れているか」をスコア化する。このスコアが高い環境ほど、その種にとって生息適地であるとみなせる。モデルは、最大エントロピー原理に従ってできるだけ単純な形を保ちながら、"Presenceの環境分布"と"Presence/Backgroundの偏り（分布比）を説明するための指数型モデルの分布"のずれ（KLダイバージェンス）を小さくするように学習される。言い換えれば、「PresenceがBackgroundと比べてどこにどれだけ偏っているか」を、データに忠実な形で再現できるようにパラメータが選ばれる。その結果得られたスコアは0〜1の存在確率っぽい値に変換されて出力される。

一方OCSVMはPresenceデータだけで「その分布の外側」を決める境界（サポートベクターによる決定境界）を学習する。今回Backgroundは主に閾値決定やモデル評価の基準としてだけ利用した。(Isolation Forestも)

MaxEntの考え方は、今回のスタバの立地適地を推定するという問題にもそのまま当てはめられると思ったので、MaxEntを用いた立地適地の推定もやってみた。
MaxEntには、Javaで書かれた本家のオープンソースGUIツール（MaxEnt）や、それをR経由で使えるライブラリ（SDMtune）などがある。Python経由の場合は、純Pythonで書かれたelapidを使うしかない。純Pythonであるため、大きなデータではやや扱いづらいが、今回程度の規模であれば問題ない。

在データと背景データ（PresenceとBackground）を定義。

# ============================================================
# (おまけ) MaxEnt (elapid) による立地適性モデリング
#   OCSVMと同じPresence-Background設計を、MaxEntの枠組みでも試す
#   在データ (Presence) = 既存スタバ店舗の特徴量
#   背景データ (Background) = カフェ周辺メッシュ（include_meshes）から
#     ランダムサンプリング（真の「不在」データではなく、MaxEnt流の
#     疑似不在点として扱う）
# ============================================================
import elapid as ela

# 使用する特徴量列（OCSVMと同じ土地利用面積 + 人流の代表値を使う）
features = areas_cols + [selected_flow]

# --- 在データ（Presence）: 既存スタバ店舗 ---
X_Presence = (
    all_coords_with_luse_log[features]
    .replace([np.inf, -np.inf], np.nan)  # inf/nanを除外
    .dropna()
    .reset_index(drop=True)
)

# --- 背景データ（Background）: カフェ周辺メッシュ（include_meshes）からサンプリング ---
#   OCSVMのときと同じ考え方で、都市/郊外のような粗い違いではなく
#   近傍内の微妙な立地差をMaxEntにも学習させる
X_Background = landuse_merged_log[landuse_merged_log['mesh1kmid'].isin(list(include_meshes))].copy()
n_Background = min(10000, len(X_Background))  # 件数が多すぎる場合は上限を設定
X_Background = (
    X_Background[features]
    .replace([np.inf, -np.inf], np.nan)
    .dropna()
    .sample(n=n_Background, random_state=0)  # 上限件数までランダムサンプリング
    .reset_index(drop=True)
)

print(f"在データ件数: {len(X_Presence)}, 背景データ件数: {len(X_Background)}")

# --- x（特徴量）と y（1=在, 0=背景）を作成 ---
x = pd.concat([X_Presence, X_Background], ignore_index=True)
y = pd.Series(
    np.concatenate([np.ones(len(X_Presence)), np.zeros(len(X_Background))]),
    name="class"
)

elapidによるモデリング。AUCでハイパーパラメータ探索。

# ============================================================
# (おまけ) MaxEnt: feature_types / beta_multiplier のグリッドサーチ（AUCで評価）
#   OCSVMのgamma/nuグリッドサーチと同じ考え方で、MaxEntの
#   モデル複雑さ(feature_types)と正則化の強さ(beta_multiplier)を
#   総当たりで探索し、CV AUCが最も高い組み合わせを確認する
# ============================================================

# feature_types: モデルに使う応答関数の種類（項が増えるほど非線形な表現が可能になる）
#   linear: 線形項のみ（最も単純）
#   linear+quadratic: 線形+二次項（応答曲線に山/谷を許容）
#   linear+quadratic+hinge: さらにヒンジ項を加え、より柔軟な非線形応答を許容
feature_type_grid = [
    ["linear"],
    ["linear", "quadratic"],
    ["linear", "quadratic", "hinge"],
]

# beta_multiplier: 正則化の強さ（値が大きいほど正則化が強く、モデルが滑らかで
#   過学習しにくくなる。値が小さいほどPresenceデータへの適合度が上がる）
beta_grid = [0.5, 1.0, 1.5, 2.0, 3.0]

results = []
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)

# feature_types × beta_multiplier の全組み合わせについてCV AUCを計算
for ft in tqdm(feature_type_grid):
    for beta in beta_grid:
        aucs = []
        for train_idx, test_idx in skf.split(x, y):
            # Presence(y=1)とBackground(y=0)を使ってMaxEntを学習
            m = ela.MaxentModel(feature_types=ft, transform="cloglog", beta_multiplier=beta)
            m.fit(x.iloc[train_idx], y.iloc[train_idx])

            # predict: 値が大きいほど「Presenceらしい（立地適性が高い）」
            y_sc = m.predict(x.iloc[test_idx])
            aucs.append(roc_auc_score(y.iloc[test_idx], y_sc))
        results.append({"feature_types": "+".join(ft), "beta": beta, "auc_mean": np.mean(aucs), "auc_std": np.std(aucs)})

# AUC平均が高い順に並べ、OCSVM・IsolationForestのAUCと比較できるようにする
results_df = pd.DataFrame(results).sort_values("auc_mean", ascending=False)
display(results_df)

image.png

クロスバリデーションで評価。

# ============================================================
# (おまけ) MaxEnt: 層化K-FoldによるAUC評価（線形項のみ・正則化を強めた設定）
#   OCSVMのevaluate_ocsvm_aucと同じ評価方針で、MaxEntの
#   汎化性能（未学習データへの順位付け能力）をAUCで確認する
# ============================================================
def evaluate_maxent_auc(x, y, n_splits=5, random_state=0):
    skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=random_state)
    aucs = []
    for train_idx, test_idx in skf.split(x, y):
        # MaxEntはPresence(y=1)とBackground(y=0)の両方を使って学習する
        # （OCSVMと異なり、教師なしではなくPresence-Background判別として学習する）
        x_tr, x_te = x.iloc[train_idx], x.iloc[test_idx]
        y_tr, y_te = y.iloc[train_idx], y.iloc[test_idx]

        # feature_types=['linear']: 環境応答を線形項のみで表現し、モデルを単純化
        # transform="cloglog": 出力を生息確率に近いスケールに変換
        # beta_multiplier=0.5: 正則化を弱めに設定し、Presenceデータへの適合度を上げる
        #   （値が小さいほど正則化が弱まり、過学習方向に振れやすくなる）
        m = ela.MaxentModel(feature_types=['linear'], transform="cloglog", beta_multiplier=0.5)
        m.fit(x_tr, y_tr)

        # predict: 値が大きいほど「Presenceらしい（立地適性が高い）」
        y_sc = m.predict(x_te)
        aucs.append(roc_auc_score(y_te, y_sc))
    return np.mean(aucs), np.std(aucs)


# 5-Fold CVでAUCの平均・標準偏差を算出し、feature_types/beta_multiplierを変えた
# 他の設定（例: linear+quadratic+hinge, beta_multiplier=1.5等）と比較する
auc_mean, auc_std = evaluate_maxent_auc(x, y)
print(f"CV AUC: {auc_mean:.3f} ± {auc_std:.3f}")

# >> CV AUC: 0.931 ± 0.023

全域に適用。

# 全候補地点への適用
model_final = ela.MaxentModel(feature_types=['linear'], transform="cloglog", beta_multiplier=0.5)
model_final.fit(x, y)

X_full = landuse_merged_log[features].replace([np.inf, -np.inf], np.nan)
mask_full = X_full.notna().all(axis=1)

landuse_merged_log.loc[mask_full, "maxent_score"] = model_final.predict(X_full.loc[mask_full])

display(landuse_merged_log[["maxent_score"] + features].sort_values("maxent_score", ascending=False))

応答曲線で、「他の変数を平均的な値に固定した状態で、1つの変数だけを変化させたときに、モデルの予測スコアがどう変わるか」を可視化。

plt.rcParams["font.family"] = prop.get_name()

fig, ax = model_final.partial_dependence_plot(x, labels=features, dpi=100, figsize=(10, 6))
plt.tight_layout()
plt.show()

image.png

在不在推定の閾値の計算。（sensitivity(tpr) + specificity(1-fpr) が最大になる点を探す）

# 学習データ全体（在データ + 背景データ）でスコアを予測
y_score_all = model_final.predict(x)

# ROC曲線を計算（fpr, tpr, thresholds はすべて同じ長さの配列）
fpr, tpr, thresholds = roc_curve(y, y_score_all)

# sensitivity(tpr) + specificity(1-fpr) が最大になる点を探す
sss = tpr + (1 - fpr)  # sensitivity + specificity
optimal_idx = np.argmax(sss)

mtss_threshold = thresholds[optimal_idx]
print(f"Maximum training sensitivity plus specificity threshold: {mtss_threshold:.4f}")
print(f"  Sensitivity: {tpr[optimal_idx]:.3f}, Specificity: {1 - fpr[optimal_idx]:.3f}")
# >> Maximum training sensitivity plus specificity threshold: 0.1035
# >>   Sensitivity: 0.944, Specificity: 0.803

654件の立地候補。

landuse_merged_log['maxent_label_mtss'] = (
    landuse_merged_log['maxent_score'] >= mtss_threshold
).astype(int)

print(landuse_merged_log['maxent_label_mtss'].value_counts())

# maxent_label_mtss
# 0    16489
# 1      654
# Name: count, dtype: int64

結果はOCSVMやIsolation Forestとおおよそ似ている。
image.png
image.png
image.png
image.png

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?