foliumを使ってハザードマップを重ねる

１．概要

Pythonの地図ライブラリ「folium」を使ってハザードマップを重ね合わせたいと思います。

ハザードマップは国土交通省国土地理院が公開しているのですが、APIも公開しているので、色々自分でカスタマイズすることもできます。ハザードマップAPIについては以下にまとまってます。
ハザードマップポータルサイト

地図のベースはfoliumを使用します。FoliumはLeaflet.jsライブラリをPythonで元々はJavascriptで作られたライブラリです。サイトとかを作るならLeaflet.jsの方が便利だと思います。 foliumを使うことで、簡単に地図上にマーカー、ライン、ポリゴンなどの地理情報を表示し、これらの要素にカスタマイズしたポップアップ情報を追加することができます。また、地図のズーム、位置、スタイルを制御できます。
GitHub - python-visualization/folium: Python Data. Leaflet.js Maps.

２．foliumについて

まずはfoliumから簡単に触ってみます。

import pandas as pd
import folium

latlons = ['35.658581', '139.745433'] # 東京タワーの緯度経度

fmap = folium.Map(
    location=latlons,
    tiles = "OpenStreetMap",
    zoom_start = 15, 
    width = 500, height = 500
)
folium.Marker(latlons, popup="東京タワー").add_to(fmap)
fmap

上記を実行すると以下のような地図が表示されます。

地図情報はOpenStreetMapというものを使用しています。GoogleMapを使うことも可能ですが、OpenStreetMapは無料なので、個人的にはこちらの方が好きです。 東京タワーの緯度経度情報をインプットに地図を表示しています。少し調べると住所から緯度経度情報を取得する方法も色々あるので、組み合わせてみても良いと思います。

３．ハザードマップを重ねてみる

次にハザードマップを重ねてみようと思います。

import pandas as pd
import folium

latlons = ['35.658581', '139.745433'] # 東京タワーの緯度経度

fmap = folium.Map(
    location=latlons,
    tiles = "OpenStreetMap",
    zoom_start = 15, 
    width = 500, height = 500
)
folium.Marker(latlons, popup="東京タワー").add_to(fmap)

fmap = add_tile_layer(
    fmap=fmap,
    tiles='https://disaportaldata.gsi.go.jp/raster/05_dosekiryukeikaikuiki/{z}/{x}/{y}.png',
    name='土砂災害警戒区域（土石流）'
)

fmap = add_tile_layer(
    fmap=fmap,
    tiles='https://disaportaldata.gsi.go.jp/raster/01_flood_l2_shinsuishin_data/{z}/{x}/{y}.png',
    name='洪水浸水想定区域（想定最大規模）'
)

fmap = add_tile_layer(
    fmap=fmap,
    tiles='https://disaportaldata.gsi.go.jp/raster/04_tsunami_newlegend_data/{z}/{x}/{y}.png',
    name='津波浸水想定'
)
folium.LayerControl().add_to(fmap)
fmap

国土地理院の公開するAPIから土石流、洪水、津波に関するハザードマップを呼び出してます。 LayerControlで重ねたハザードマップのレイヤーを消したり表示したりするコントロールを付与してます。 東京タワー周辺、、結構赤くなってますね。。

他にも公開されている情報はたくさんあるので、色々組み合わせてみると面白いかもしれません。

Pythonで株価データや金融データを取得する方法（pandas_datareader）

pandas-datareaderというpythonライブラリを使用して株価データや金融データを取得する方法についてまとめます。pandas-datareaderとはウェブ上の様々なデータソースにアクセスするライブラリです。今回は以下にアクセスします。

アクセス先	概要
Stooq	ポーランドのサイトで日本の株価データも取り扱っている。
St.Louis FED (FRED)	セントルイス連邦準備銀行の公開している経済統計データ集。 株価データは取り扱っていない。

上記以外のデータソースは以下にまとまっています。
Data Readers — pandas-datareader v0.10.0 documentation

１．pandas-datareaderの使用方法

① インストール

pipコマンドを使用して以下の通りインストール。

pip install pandas-datareader

② Stooqからデータ取得

以下コードでStooqからデンソーの株価データを取得します。

# ライブラリのインポート
import pandas_datareader.stooq as web
from datetime import datetime

# データ取得期間の設定
st = datetime(2015, 1, 1)
ed = datetime(2020, 1, 1)

# 株価価格等の取得
lst1=['6902.JP', #デンソー
    ]
stooq = web.StooqDailyReader(lst1,start=st,end=ed).read()
display(stooq.head())

＜出力結果＞

③ FREDからデータ取得

次にFREDからは株価以外の経済統計データを取得してみます。
試しにResidential Property Prices for Japanを取得してみます。

# ライブラリのインポート
import pandas_datareader.fred as web
from datetime import datetime

# データ取得期間の設定
st = datetime(2015, 1, 1)
ed = datetime(2020, 1, 1)

# Residential Property Prices for Japanを取得
fred = web.FredReader('QJPN368BIS',start=st,end=ed).read()
display(fred.head())

＜出力結果＞

Seleniumエラー：Element is not clickable at point (xxx, xxx)について

Selenium使用時に「Element is not clickable at point (xxx, xxx)」というエラーが出て色々調べたので、その内容についてまとめます。

１．結論

最初に結論だけ述べると、、

• エラー原因：画面外のelementを操作しようとしたため
• 解決方法：対象のelementまでSeleniumでスクロールする

２．事象再現

まずは適当にSeleniumを使って事象を再現してみます。

#モジュールのインポート
import glob
import datetime
import time
from selenium import webdriver
from selenium.webdriver.common.by import By

# ブラウザ起動
driver = webdriver.Chrome()
driver.maximize_window()

# アクセスするURL
TARGET_URL = "https://suumo.jp/chintai/tokyo/ensen/"

# 対象サイトへアクセス
driver.get(TARGET_URL)
time.sleep(2)

#チェックボックスをクリック
elements = driver.find_elements(By.CLASS_NAME, 'js-fr-checkSingle')
elements[10].click()
elements[30].click()

試していただければ分かると思いますが、実はelement[10].click()をコメントアウトすると、element[30]が画面外にあるボタンであるにも関わらず、エラーは起きないです。つまり、より詳細な原因としては、「一つ目のボタンを選択した状態で、画面外のボタンを操作しようとするとエラーになる」ようです。

３．解決方法

解決するために、Seleniumを使って画面スクロールします。

#element[30]が画面の真ん中らへんに来るように調整
window_size = driver.get_window_size()
loc_y =  elements[30].location['y']
if loc_y < window_size['height']/2:
    scroll_y = 0
else:
    scroll_y = loc_y - window_size['height']/2
driver.execute_script(f"window.scrollBy(0, {scroll_y});")
elements[30].click()

以下のようにシンプルに対象elementまでスクロールしちゃってもいいのですが、WEBページによってはヘッダが邪魔で対象elementが隠れちゃうことがあります。隠れちゃうと同じエラーが出るので、少しめんどくさいですが、対象elementが真ん中らへんに来るように調整してます。スマートな実装があれば教えてほしいです。

４．ソース一式

#モジュールのインポート
import glob
import datetime
import time
from selenium import webdriver
from selenium.webdriver.common.by import By

# ブラウザ起動
driver = webdriver.Chrome()
driver.maximize_window()

# アクセスするURL
TARGET_URL = "https://suumo.jp/chintai/tokyo/ensen/"

# 対象サイトへアクセス
driver.get(TARGET_URL)
time.sleep(2)

#チェックボックスをクリック
elements = driver.find_elements(By.CLASS_NAME, 'js-fr-checkSingle')
elements[10].click()

#element[30]が画面の真ん中らへんに来るように調整
window_size = driver.get_window_size()
loc_y =  elements[30].location['y']
if loc_y < window_size['height']/2:
    scroll_y = 0
else:
    scroll_y = loc_y - window_size['height']/2
driver.execute_script(f"window.scrollBy(0, {scroll_y});")
elements[30].click()

タスクスケジューラで陥りやすいトラブル集

タスクスケジューラ関連で個人的に躓いたポイントとその解決策についてまとめておきたいと思います。

１．予定時刻になってもタスクが実行されない問題

① 事象

以下のように繰り返し実行のタスクを作成するが、予定時刻になってもタスクが実行されないという事象。

予定時刻になっても、エラー表示も出ないままタスクが実行されず、勝手に「次回の実行時刻」が延期されます。エラーが出ないので、原因特定に時間がかかりました。

② 原因と対策

どうやら、これは「トリガー」設定に問題があったようです。以下のように、周期を[毎日]と設定したうえで、継続時間を[無期限]と設定するとタスクが正常に動作しないとのことです。Microsoft Japan Windows Technology Support Blogに情報があったので、詳しくはそちらを参照ください。

正確には以下のパターンのときに問題が発生するようです。
× : 毎日 / 毎週 / 毎月 + 繰り返し間隔の継続時間「無期限」

以下の場合はOKとのこと。
○ : 毎日 / 毎週 / 毎月 + 繰り返し間隔の継続時間「1 日間」など
○ : 1 回 + 繰り返し間隔の継続時間「無期限」

２．タスクは実行されるけどバッチが実行されない問題

① 事象

タスクスケジューラ上で「前回の実行時刻」が更新される（タスク自体は実行されている）のに、タスクから呼び出すバッチが起動しない事象。

これもエラーも出ず、ただただバッチが実行されないので、原因がなかなかわかりませんでした。

② 原因と対策

「操作の編集」の「開始（オプション）」が指定されていないのが問題らしいです。タスクスケジューラに関する記事を見ても「開始（オプション）」を必須と紹介する記事は少ないので、設定しなくても問題なく動くパターンもあるんだと思います。

競馬データのスクレイピング

netkeibaから競馬データをスクレイピングする方法についてまとめます。既に同様の記事は世にたくさん出回ってますが、少し改良して効率化してみました。

１．既出のスクレイピング方法の問題点

競馬データのスクレイピングで検索すると、以下のようなソースコードが多いと思います。

for i in range(1, 11):
    for j in range(1, 7):
        for k in range(1, 13):
            for l in range(1, 13):
                Base = "https://race.netkeiba.com/race/result.html?race_id="
                race_id = "2019" + str(i).zfill(2) + str(j).zfill(2) + str(k).zfill(2) + str(l).zfill(2)
                url = Base + race_id
                '''------------------'''
                '''スクレイピング処理'''
                '''------------------'''

上記はrace_idを総当たりスクレイピングする方法ですが、残念ながら実際にはrace_idは飛び飛びで、存在しないrace_idがたくさんあります。よって無駄な検索をしていることになります。スクレイピング中にsleep処理を入れると思いますが、無駄な検索をするたびにsleepが入って、全体の処理時間が長引きます。

２．どうやって効率化するか

netkeibaには1日毎の開催レース一覧をまとめている以下のようなページがあります。
レース一覧 | 2023年5月20日 レース情報(JRA) - netkeiba.com

このページからrace_idの一覧を取得することで、実在するrace_idのリストを作成することができます。ここまで来たら後はあとはrace_idをキーにスクレイピングするだけです。

２．１．開催レース一覧ページからrace_idを取得する方法

開催レース一覧のページは動的にrace_idを設定しているようなので、beautifulsoupは使えません。こういうときはseleniumが効果的です。 以下のような実装でrace_idの一覧を取得できます。

import time
from selenium import webdriver
from selenium.webdriver.chrome.options import Options
from selenium.webdriver.common.by import By

#selenium driverでブラウザ起動
def get_selenium_driver(url):
    options = Options()
    # ヘッドレスモードで実行する場合
    options.add_argument("--headless")
    driver = webdriver.Chrome(options=options)

    # 取得先URLにアクセス
    driver.get(url)
    # コンテンツが描画されるまで待機
    time.sleep(5)
    
    return driver

#指定した開催日に開催されるレースのrace_id一覧を取得
def get_race_ids(kaisai_dates):
    race_ids = []
    for kaisai_date in kaisai_dates:
        url = 'https://race.netkeiba.com/top/race_list.html?kaisai_date=' + kaisai_date
        driver = get_selenium_driver(url)
        elements = driver.find_elements(By.CLASS_NAME,"MyRace_List_Item")
        for element in elements:
            race_id = str(element.get_attribute('id'))
            race_id = race_id.replace("myrace_","")
            race_ids.append(race_id)
    
    return race_ids

race_ids = get_race_ids(["20230528"])

上記で2023年5月28日に開催された全レースの全race_idをリストとして取得できます。

２．２．競馬のレース結果のスクレイピングする方法

これは他にサイトがたくさんあるので、今更紹介する必要はないと思いますが、一応まとめておきます。

from bs4 import BeautifulSoup
import requests
import pandas as pd

def fetch_horse_datas(race_id, horse_datas):
    Base="https://race.netkeiba.com/race/result.html?race_id="
    url = Base + race_id
    kaisai_year = int(race_id[:4])

    request = requests.get(url)
    soup = BeautifulSoup(request.content, 'html.parser')
    horses = soup.find_all(class_='HorseList')

    data1 = soup.find(class_='RaceData01')
    data1 = data1.get_text().replace("/","").split()

    data2 = soup.find(class_='RaceData02')
    data2 = data2.select('span')
    data2 = [t.get_text(strip=True) for t in data2]

    data3 = soup.find(class_='Refundlink')
    data3 = str(data3)
    
    for horse in horses:
        horse_data = {}
        horse_data['start_time']    = data1[0]
        horse_data['race_type']     = data1[1]
        horse_data['ground']        = data1[2]
        horse_data['weather']       = data1[3]
        horse_data['gr_condition']  = data1[4]
        horse_data['place']         = data2[1]
        horse_data['rule']          = data2[3]
        horse_data['grade']         = data2[4]
        horse_data['sex']           = data2[5]
        horse_data['rule6']         = data2[6]
        horse_data['num_horse']     = data2[7]
        horse_data['year']          = kaisai_year
        horse_data['date']          = data3[data3.find("kaisai_date=")+12 : data3.find("kaisai_date=")+20]
        horse_data['race_id']       = race_id
        horse_data['rank']          = horse.find(class_='Rank').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['horse_number']  = horse.find(class_='Txt_C').find('div').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['horse_name']    = horse.find(class_='Horse_Name').find('a').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['seirei']        = horse.find(class_='Lgt_Txt Txt_C').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['jockey_weight'] = horse.find(class_='JockeyWeight').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['jockey_name']   = horse.find(class_='Jockey').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['race_time']     = horse.find(class_='RaceTime').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['popularity']    = horse.find(class_='OddsPeople').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['odds']          = horse.find(class_='Txt_R').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['passage']       = horse.find(class_='PassageRate').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['trainer']       = horse.find(class_='Trainer').find('a').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['weight']        = horse.find(class_='Weight').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_datas.append(horse_data)
    return horse_datas

# レースデータをスクレイピング
horse_datas = []
race_ids = get_race_ids(["20230528"])
for race_id in race_ids:
    try:
        # 出馬表を読み取る    
        horse_datas = fetch_horse_datas(race_id, horse_datas)
        time.sleep(1.5)

    except Exception as e:
        print("error in race_id:",race_id,e)
        pass
                    
df = pd.DataFrame(horse_datas)

上記で2023年5月28日の全レース結果をスクレイピングできます。

３．ソースコード一式

以下にソースコード一式をまとめます。 以下では取得したい「年」を指定することで、その年の全レースデータ結果を取得できるようにしてます。

from bs4 import BeautifulSoup
import requests
import pandas as pd
pd.set_option('display.max_columns', 150)
pd.set_option('display.max_rows', 500)
from tqdm import tqdm
import time
import datetime

from selenium import webdriver
from selenium.webdriver.chrome.options import Options
from selenium.webdriver.common.by import By

# 今日の日付を取得
today = datetime.date.today()

# 読み取り開始年、終了年
st_year = 2015
ed_year = today.year

# 土日リストを取得
def get_weekends():
    weekends = []
    for target_year in range(st_year, ed_year+1):
        baseDate = datetime.date(target_year, 1, 1)
        days = (baseDate - datetime.date(target_year - 1, 1, 1)).days
        weekends.extend([(baseDate + datetime.timedelta(i)).strftime("%Y%m%d") for i in range(0, days) if (baseDate + datetime.timedelta(i)).weekday() >= 5])
    return weekends

def get_kaisai_dates():
    weekends = get_weekends()
    return [weekend for weekend in weekends if (weekend<today.strftime("%Y%m%d"))]

#selenium driverでブラウザ起動
def get_selenium_driver(url):
    options = Options()
    # ヘッドレスモードで実行する場合
    options.add_argument("--headless")
    driver = webdriver.Chrome(options=options)

    # 取得先URLにアクセス
    driver.get(url)
    # コンテンツが描画されるまで待機
    time.sleep(5)
    
    return driver

#指定した開催日に開催されるレースのrace_id一覧を取得
def get_race_ids(kaisai_dates):
    race_ids = []
    for kaisai_date in tqdm(kaisai_dates):
        url = 'https://race.netkeiba.com/top/race_list.html?kaisai_date=' + kaisai_date
        driver = get_selenium_driver(url)
        elements = driver.find_elements(By.CLASS_NAME,"MyRace_List_Item")
        for element in elements:
            race_id = str(element.get_attribute('id'))
            race_id = race_id.replace("myrace_","")
            race_ids.append(race_id)
    
    return race_ids

def fetch_horse_datas(race_id, horse_datas):
    Base="https://race.netkeiba.com/race/result.html?race_id="
    url = Base + race_id
    kaisai_year = int(race_id[:4])

    request = requests.get(url)
    soup = BeautifulSoup(request.content, 'html.parser')
    horses = soup.find_all(class_='HorseList')

    data1 = soup.find(class_='RaceData01')
    data1 = data1.get_text().replace("/","").split()

    data2 = soup.find(class_='RaceData02')
    data2 = data2.select('span')
    data2 = [t.get_text(strip=True) for t in data2]

    data3 = soup.find(class_='Refundlink')
    data3 = str(data3)

    for horse in horses:
        horse_data = {}
        horse_data['start_time']    = data1[0]
        horse_data['race_type']     = data1[1]
        horse_data['ground']        = data1[2]
        horse_data['weather']       = data1[3]
        horse_data['gr_condition']  = data1[4]
        horse_data['place']         = data2[1]
        horse_data['rule']          = data2[3]
        horse_data['grade']         = data2[4]
        horse_data['sex']           = data2[5]
        horse_data['rule6']         = data2[6]
        horse_data['num_horse']     = data2[7]
        horse_data['year']          = kaisai_year
        horse_data['date']          = data3[data3.find("kaisai_date=")+12 : data3.find("kaisai_date=")+20]
        horse_data['race_id']       = race_id
        horse_data['rank']          = horse.find(class_='Rank').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['horse_number']  = horse.find(class_='Txt_C').find('div').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['horse_name']    = horse.find(class_='Horse_Name').find('a').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['seirei']        = horse.find(class_='Lgt_Txt Txt_C').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['jockey_weight'] = horse.find(class_='JockeyWeight').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['jockey_name']   = horse.find(class_='Jockey').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['race_time']     = horse.find(class_='RaceTime').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['popularity']    = horse.find(class_='OddsPeople').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['odds']          = horse.find(class_='Txt_R').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['passage']       = horse.find(class_='PassageRate').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['trainer']       = horse.find(class_='Trainer').find('a').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_data['weight']        = horse.find(class_='Weight').get_text().replace('\n','').replace(' ','').replace('　','')
        horse_datas.append(horse_data)
    return horse_datas

# レースデータをスクレイピング
kaisai_dates = get_kaisai_dates()
race_ids = get_race_ids(kaisai_dates)
horse_datas = []
for race_id in tqdm(race_ids):
    try:
        # 出馬表を読み取る    
        horse_datas = fetch_horse_datas(race_id, horse_datas)
        time.sleep(1.5)

    except Exception as e:
        print("error in race_id:",race_id,e)
        pass
                    
df = pd.DataFrame(horse_datas)

SHAPを使用した回帰問題の機械学習モデルの局所解釈方法

shapライブラリを使用して、回帰問題を解いた機械学習モデルの大局的解釈を行う。

１．SHAPとは

SHAP（SHapley Additive exPlanations）は、機械学習モデルの局所的な解釈可能性を提供するためのフレームワークです。SHAPは、個々の特徴量が予測にどのように寄与しているかを計算することにより、モデルの解釈性を高めます。

２．データセット

今回はsklearnのdiabetes（糖尿病）データセットを例に実装してみる。データセットの概要は以下の通り。

• 目的変数：1年後の糖尿病の進行に関する測定値
• 説明変数：11個の特徴量

特徴量	説明
age	患者の年齢
sex	患者の性別
bmi	患者のBMI
bp	患者の平均血圧
S1	T-Cells (白血球の1種)の活性化に関与する化学物質
S2	インスリンのレセプターの活性に関与する化学物質
S3	グルコースの代謝に関与する化学物質
S4	トリグリセリド (脂質の一種)の代謝に関与する化学物質
S5	血清アルブミン (タンパク質の一種)のレベル
S6	銅の輸送に関与する化学物質

※S1~S6: 血清の6つの化学指標。

３．実装

３．１．Lightgbmを用いた簡易的なモデルを作成

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
import lightgbm as lgb
import pandas as pd
import shap

# データのロード
diabetes = load_diabetes()
x_train, x_test, y_train, y_test = train_test_split(diabetes.data, diabetes.target, test_size=0.3)

# LightGBMモデルのトレーニング
model = lgb.LGBMRegressor().fit(x_train, y_train)

３．２．SHAPの実装

# SHAP値の計算
df_x_test = pd.DataFrame(x_test, columns=diabetes.feature_names)
explainer = shap.Explainer(model.predict, df_x_test)
shap_values = explainer(df_x_test)

# 特定のサンプルのSHAP値を表示（shap_values[x]：x番目のデータのSHAP値を表示、max_display=n：n個の特徴量を描画）
shap.plots.waterfall(shap_values[0], max_display=10)

上記の結果、以下のような図が表示されます。

E[f(x)]は、モデルの期待出力値（期待値）を表しています。つまり、全体的な平均予測値です。 SHAP値は、個々の特徴量がこの期待値からどれだけ上下に影響を与えたかを表します。

PDP／ICEを使用した回帰問題の機械学習モデルの大局的解釈方法

pdpboxライブラリを使用して、回帰問題を解いた機械学習モデルの大局的解釈を行う。

１．データセット

今回はsklearnのボストン住宅価格データセットを例に実装してみる。データセットの概要は以下の通り。

• 目的変数：ボストンの地域別住宅価格
• 説明変数：以下13個の特徴量¹

特徴の名前	説明
CRIM	人口 1 人当たりの犯罪発生数
ZN	25,000 平方フィート以上の住居区画の占める割合＝「広い家の割合」
INDUS	小売業以外の商業が占める面積の割合
CHAS	チャールズ川に関わるダミー変数 (1: 川の周辺, 0: それ以外)
NOX	一酸化窒素の濃度
RM	住居の平均部屋数
AGE	1940 年より前に建てられた物件の割合 ＝「古い家の割合」
DIS	ボストン市の5 つの雇用施設からの距離 (重み付け済) ＝「主要施設への距離」
RAD	高速道路へのアクセスのしやすさ
TAX	$10,000 ドルあたりの固定資産税率
PTRATIO	町毎の生徒と教師の比率
B	町毎の黒人 (Bk) の比率
LSTAT	低所得者人口の割合

２．実装

２．１．Lightgbmを用いた簡易的なモデルを作成

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
import pandas as pd
import lightgbm as lgb

boston = load_boston()
x_train, x_test, y_train, y_test = train_test_split(boston.data, boston.target, test_size=0.3)

#LGB用のデータに変形
lgb_train = lgb.Dataset(x_train, y_train)
lgb_eval = lgb.Dataset(x_test, y_test)

params = {
    'objective': 'regression',
    'boosting_type': 'gbdt',
    'metric': 'mae'
    }

model = lgb.train(
    params = params, 
    train_set = lgb_train,
    valid_sets= [lgb_eval,lgb_train],
)

２．２．PDP／ICEの実装

# pdpboxはpandas.DataFrameしか受け付けない
test_df = pd.DataFrame(x_test, columns=boston.feature_names)

# 特徴量毎にpdp/iceを描画
for col in train_df.columns:
    pdp_boston = pdp.pdp_isolate(
        model=model, dataset=test_df, model_features=test_df.columns, feature=col
    )
    fig, axes = pdp.pdp_plot(pdp_boston, col, plot_lines=True, frac_to_plot=1.0, plot_pts_dist=True)

特徴量毎に以下のような図が描画される。


上図より、一酸化窒素の濃度が0.7までは住宅価格にあまり影響がなく、0.7以上あたりから住宅価格を下げる要因と機械学習モデルが判断していることが分かる。

２．３．PDP（2変数）の実装

ついでに2変数の相互作用を等高線として図示する方法もまとめておく。

import itertools
columns = test_df.columns
combinations = itertools.combinations(columns, 2)

for combination in combinations:
    interaction = pdp.pdp_interact(model=model, dataset=test_df, model_features=test_df.columns, features=combination)
    fig, axes = pdp.pdp_interact_plot(pdp_interact_out=interaction, feature_names=combination)

特徴量の組み合わせ毎に以下のような図が描画される。


上図より、固定資産税率が350を超えたあたりから住宅価格にあまり影響がなく、犯罪発生率の方が影響していることが分かる。

---

1. ボストンデータセットの特徴量の説明：【初心者】ネコでも分かる「scikit-learnのサンプルデータ」まとめ【Python】を引用↩