RのサンプルデータをPythonで使う～標準ベイズ統計学第8.3節を題材として～

はじめに

ベイズ統計学の基礎をあらためて学ぶために，ピーター・D・ホフ著，入江薫・菅澤翔之助・橋本真太郎訳「標準ベイズ統計学」を読むことにした。

本記事では，第8.3節を題材とした，RのサンプルデータをPythonで使う方法をまとめた。

本記事のハイライト

「標準ベイズ統計学」のサンプルデータのありかを説明した。
RのサンプルコードをPythonで移植する際に利用する，.RDataや.rdaファイルを読み込むためのパッケージを紹介した。

はじめに
- 本記事のハイライト
- 目次
第8章　グループ比較と階層モデリング
- データを探せ！～「標準ベイズ統計学」のサンプルデータ～
  - 学校データを探す
RのデータをPythonで使う
- rdataの利用例
- pyreadrの利用例
まとめと感想
- 参考サイト

第8章　グループ比較と階層モデリング

データを探せ！～「標準ベイズ統計学」のサンプルデータ～

「標準ベイズ統計学」を読み進めていると，サンプルコードが出てくる。
サンプルコードがあると，数式をプログラムに実装するための方法もわかるし，細かな挙動が分かるので理解が深まる。

「標準ベイズ統計学」の紹介ページを見ていると，「関連情報」に『標準ベイズ統計学』原著者サポートサイトのリンクがあり，これをクリックすると原著である"A First Course in Bayesian Statistical Methods"へのリンクが貼られている。

原著のサポートサイトを見ると，

という2つのリンクがあり，ここにサンプルコードやデータが保管されていた。

本書中に書かれているサンプルコードは，Data and code to replicate figures and numerical results の方に書かれている。
しかしコードを細かく見ていると，Chapter8.Rの115～118行目を見ていると，

#### Put data into list form. This is different from in the book. 
Y<-list()
YM<-NULL
J<-max(Y.school.mathscore[,1])

のように，急にY.school.mathscoreというデータが出てきた。

どのデータなのかよくわからず，いろいろ探していたところ，もう1つのリンクであるData and code for inline examples において，Y.school.mathscoreというデータがあった。
そのため，本書のサンプルコードを読んでいて，よくわからないデータが出てきたら，サポートサイトのリンクを細かく探してみることが重要だと気付いた。

学校データを探す

8.4節では，上記のY.school.mathscore を使うことになる。この内容を確認すると，以下のようなテキストデータになっていた。

テキストでは書かれているものの，Pythonで扱えるようにするためには，数値データや列名を読み解き，pandas.DataFrameやnumpy.arrayで扱えるように整形する必要がある。

上記のデータを整形するためのコードを，生成AIと相談しながら作成した。なおY.school.mathscoreは，ダウンロードした後に拡張しを追加して，Y.school.mathscore.txt にファイル名を変更した。

データ整形用サンプルコードはクリックで展開

import re
import numpy as np
import pandas as pd


def read_r_dput_matrix(filename):
    """
    R の dput()/dget() 形式で保存された
    structure(c(...), .Dim=..., .Dimnames=...)
    を pandas.DataFrame として読み込む
    """

    with open(filename, "r", encoding="utf-8") as f:
        txt = f.read()

    # .Dim を取得
    dim_match = re.search(
        r"\.Dim\s*=\s*c\(\s*(\d+)L?\s*,\s*(\d+)L?\s*\)",
        txt,
        re.DOTALL,
    )

    if dim_match is None:
        raise ValueError(".Dim が見つかりません")

    nrow = int(dim_match.group(1))
    ncol = int(dim_match.group(2))

    # c(...) の中身を取得
    data_match = re.search(
        r"structure\s*\(\s*c\((.*?)\)\s*,\s*\.Dim",
        txt,
        re.DOTALL,
    )

    if data_match is None:
        raise ValueError("データ部分が見つかりません")

    data_str = data_match.group(1)

    values = np.fromstring(data_str, sep=",")

    # Rは列優先格納
    arr = values.reshape((nrow, ncol), order="F")

    # 列名取得
    col_match = re.search(
        r'c\("([^"]+)"\s*,\s*"([^"]+)"\)',
        txt
    )

    if col_match:
        columns = [col_match.group(1), col_match.group(2)]
    else:
        columns = [f"V{i+1}" for i in range(ncol)]

    return pd.DataFrame(arr, columns=columns)

df = read_r_dput_matrix("../r_codes/Y.school.mathscore.txt")
df["school"] = df["school"].astype(int)
df.to_csv("../r_codes/Y_school_mathscore.csv", index=False)

上記のサンプルコードでは，

元データのうち，scool列は整数型なので，整数に変換した。
後で使いやすくするよう，元のデータを"Y_school_mathscore.csv" という名称のCSVファイルにした。

のような処理を行なっている。

RのデータをPythonで使う

Rで書かれたサンプルコードを扱っていると，入力データの拡張子が.RDataや.rdaのような，R用のデータの拡張子になっていることがある。
.csvや.txtであればPythonでもすぐに扱えるが，これらR用のデータの拡張子ではひと工夫がいる。

このような場合には，

rdataパッケージ
pyreadrパッケージ

を用いればよい。

これらのパッケージの利用方法は，以下の過去記事にサンプルコードを示した。

rdataの利用例

参考：「データ解析のための統計モデリング入門」を読む～第11章空間構造のある階層ベイズモデル～ - jiku log

import rdata
parsed = rdata.parser.parse_file('/content/drive/My Drive/Colab Notebooks/StatModeling/Y.RData')
converted = rdata.conversion.convert(parsed)
print(converted)

import pandas as pd
d = pd.DataFrame(converted)

pyreadrの利用例

参考：馬場ベイズ本をNumPyroで実装～第5部・第8章 DGLM(二項分布)～ - jiku log

# ============================================
# データの読み込み
# ============================================

# KFASパッケージのboatデータ
# https://github.com/helske/KFAS/blob/main/data/boat.rda をダウンロードする

import pyreadr

result = pyreadr.read_r("boat.rda")

print(result.keys())

boat = result["boat"]

print(boat.head())