2018-03-14

Pythonで実装しながら緑本を学ぶ (第10章階層ベイズモデル)

Python 統計学読書メモ

データ解析のための統計モデリング入門(通称、緑本)を読み進めています。
述べられている理論を整理しつつ、Rでの実装をPythonに置き換えた際のポイントなども深掘りしていきます。

今回は第10章です。PyMC3を使って階層ベイズモデルを表現します。実装は以下で公開しています。

introduction_to_machine_learning_with_python/chapter10.ipynb at master · ohke/introduction_to_machine_learning_with_python · GitHub

10 GLMのベイズモデル化と事後分布の推定 -GLMMのベイズモデル化-

階層事前分布を使って一般化線形混合モデル(GLMM)を階層ベイズモデルとして扱うことが、本章のテーマです。

10.1 例題：個体差と生存種子数(個体差あり)

1個体あたりの種子数8個を持つ架空植物100個体について、生存種子数を計測した。個体 $i$ における生存種子数を $y_i$ とする。

8個の内 $y_i$ 個というデータですので、二項分布が使えそうです。
しかし個体に由来するランダム効果によって過分散が生じています。平均4.03に対して分散9.93となっており、二項分布が期待する分散2.00(=8×(4.03/8)×(1-(4.03/8)))よりも、かなり大きいです。

f:id:ohke:20180310095005p:plain

10.2 GLMMの階層ベイズモデル化

第7章(こちらの投稿)では、GLMにこうした原因不明の個体差を組み込んだGLMMを導入し、パラメータを最尤推定で求めた。
パラメータが増えると最尤推定では困難になるため、MCMCサンプリングによる事後分布推定を行いたい。そこで、GLMMのベイズモデル化を行う。

まずリンク関数と線形予測子を決めます。二項分布ですので、リンク関数にはロジットを使います(第6章も参照)。 $r_i$ が個体差です。

$$ logit(q_i)=\beta+r_i $$

8章・9章で追ってきたとおり、事後分布は $p({\bf Y} \mid \beta, \{r_i\}) \times \mbox{事前分布}$ に比例します。

データ ${\bf Y }$ が得られる条件付き確率(尤度)は全個体の二項分布の積です。

$$ p({\bf Y} \mid \beta, {r_i})=\prod_i \binom{8}{y_i} q_i^{y_i}(1-q_i)^{8-y_i} $$

$\beta$ の事前分布は、特に制約がないので平均0・標準偏差100の正規分布に従うこととする(無情報事前分布)。

$$ p(\beta)=\frac{1}{\sqrt{2\pi \times 100^2}} \exp(\frac{-\beta^2}{2 \times 100^2}) $$

個体差 $r_i$ は、個体差のばらつきを表す変数 $s$ を導入し、平均0・標準偏差 $s$ の正規分布に従うことにする。

$$ p(r_i \mid s)=\frac{1}{\sqrt{2\pi s^2}} \exp(\frac{-r_i^2}{2s^2}) $$

$s$ は正であれば良いので、0から10000の連続一様分布に従う無情報事前分布とする。

$$ p(s)=U(0, 10^4) $$

このように、別の事前分布 $p(s)$ によって定まる事前分布 $p(r_i \mid s)$ を階層事前分布と呼ばれます。これに対して、 $s$ を超パラメータ、 $p(s)$ を超事前分布)と言います。

階層事前分布を使うベイズモデルが、階層ベイズです。

10.3 階層ベイズモデルの推定・予測

PyMC3を使い、10.2のモデルを構築します。超事前分布 $p(s)$ 、階層事前分布 $p(r \mid s)$ の順に作っていることを確認してください。

Normalで正規分布、Uniformで連続一様分布、Binomialで二項分布を設定できます
- rは個体数分必要ですので、shape=len(data.y)としています
- invlogitでロジット関数の逆関数(つまりロジスティック関数)を二項分布のpに設定しています
NUTSで1600点(最初の100点は捨てる)をサンプリングしたサンプリング列を3つ作る
- 合計4500点をサンプリング

import pymc3 as pm

# モデルの作成
with pm.Model() as model:
    # βの事前分布をN(0, 100)の正規分布で設定(無情報事前分布)
    beta = pm.Normal('beta', mu=0, sd=100)
    
    # 超パラメータsの(超)事前分布をU(0, 10000)の連続一様分布で設定(無情報事前分布)
    s = pm.Uniform('s', lower=0, upper=10000)
    
    # パラメータrの事前分布をN(0, s)の正規分布で設定(階層事前分布)
    r = pm.Normal('r', mu=0, sd=s, shape=len(data.y))
    
    # ロジットリンク関数を設定し、二項分布で推定する
    y = pm.Binomial('y', n=8, p=pm.invlogit(beta + r), observed=data.y.values)

# サンプリング
with model:
    # NUTSで101個目からサンプルを取得するチェインを3つ作る
    trace = pm.sample(1500, step=pm.NUTS(), tune=100, njobs=3, random_seed=0)
    
_ = pm.traceplot(trace)
pm.summary(trace).loc[['beta', 's']]

各パラメータの事後分布(95%信用区間)を見ると $\beta$ は-0.600〜0.723、 $s$ は0.010〜3.802となっており、いずれも $\hat{R}$ は十分1に近い(収束している)。

f:id:ohke:20180310111301p:plain

$r_i$ の事後分布(下図中段左)を見ると、 $r_i$ 毎に異なる分布となっており、個体差が表現されていることがわかる。

f:id:ohke:20180310111413p:plain

パラメータの事後分布から、生存確率を予測してみます。

3つのサンプル列の平均値を求める( $\beta$ と $s$ のそれぞれ1500点ずつ)
$\beta$ と $s$ を1つずつ取り出し、 $s$ から100個体分の $\{r_i \}$ を正規分布で生成する
$\beta$ と $\{r_i \}$ から二項分布で種子数の生存確率を求める
1500サンプル・100個体分の生存確率の平均値を計算する

こうして得られた生存確率から、100個体を調べた時の生存種子数の頻度をプロットしたのが以下となります。過分散を予測できていることがわかります。

f:id:ohke:20180310115057p:plain

10.4 ベイズモデルで使うさまざまな事前分布

事前分布の決め方は3パターン。

ある分布を信じて、主観的に決める (基本的には使わない)
- 他のデータセットや統計モデルで得られた分布を使う場合など
無情報事前分布
階層事前分布

無情報事前分布は、個体ごとに自由に決めてしまうと、フルモデルと同等になるデメリットがある。そのため10.2のように、個体ごとに異なるパラメータなどは階層事前分布とし、階層事前分布の(超)パラメータのみを無情報事前分布とするのが基本戦略となる。
したがって、パラメータの影響が局所的な場合は階層事前分布、大域的な場合は無情報事前分布と使い分ける。

10.5 個体差＋場所差の階層ベイズモデル

個体差に加えて、場所差を組み込んだモデルを見ていきます。

架空植物100個体の種子数を計測し、施肥処理と種子数の影響について調べたい、というケースを想定します。
架空植物は10個の植木鉢A〜Jで10個体ずつ植えられてます。植木鉢A〜Eの50個体は施肥処理無し('C')、F〜Jの50個体は施肥処理有り('T')で育てられています。

個体ごとの種子数をプロットしてみますと、施肥処理無し(左側50個体)の方が若干種子数が多いように見えます。

f:id:ohke:20180310121717p:plain

しかし、植木鉢ごとの種子数を箱ひげ図でプロットしてみますと、種子数の分布が異なっており、植木鉢による差が少なからずあるものと思われます。

f:id:ohke:20180310121731p:plain

個体差と植木鉢差を考慮したモデルを定義していきます。

種子数は離散値・上限なしのため、ポアソン分布を使います。
個体 $i$ におけるパラメータ $\lambda_i$ は、線形予測子と対数リンク関数でモデルを定義します。

$$ \log \lambda_i = \beta_1 + \beta_2 f_i + r_i + r_{j(i)} $$

は、平均0・標準偏差100のの正規分布
- $f_i$ は個体 $i$ の施肥処理の有無(有りの場合に1となる)
$r_i$ は個体差で、平均0・標準偏差 $s$ の正規分布 (階層事前分布)
は植木鉢差で、平均0・標準偏差の正規分布 (階層事前分布)
- パラメータ数は10個( $j$ は0〜9)
$s, s_p$ は0〜10000の連続一様分布

PyMC3での実装は以下のとおりです。

with pm.Model() as model:
    # β1とβ2の事前分布をN(0, 100)の正規分布で設定(無情報事前分布)
    beta1 = pm.Normal('beta1', mu=0, sd=100)
    beta2 = pm.Normal('beta2', mu=0, sd=100)
    
    # sとspの(超)事前分布をU(0, 10000)の連続一様分布で設定(無情報事前分布)
    s = pm.Uniform('s', lower=0, upper=10000)
    sp = pm.Uniform('sp', lower=0, upper=10000)
    
    # 個体差パラメータrの事前分布をN(0, s)の正規分布で設定(階層事前分布)
    r = pm.Normal('r', mu=0, sd=s, shape=len(data.y))
    
    # 場所差パラメータrpの事前分布をN(0, sp)の正規分布で設定(階層事前分布)
    rp = pm.Normal('rp', mu=0, sd=sp, shape=len(data.POT.unique()))
    
    # ログ関数を設定し、ポアソン分布で推定する
    y = pm.Poisson('y', mu=np.exp(beta1 + beta2 * data.F + r + rp[data.POT]), observed=data.y.values)
    
# HMCでサンプリング
with model:
    # 101個目からサンプルを取得するチェインを3つ作る
    trace = pm.sample(1500, step=pm.HamiltonianMC(), tune=100, njobs=3, random_seed=0)

pm.summary(trace).loc[['beta1', 'beta2', 's', 'sp']]

$\beta_2$ に着目すると、95%信用区間で-2.243〜0.714となっており、0を跨いでいるので施肥処理による効果は見られないようです。

f:id:ohke:20180310124626p:plain

2018-03-08

Pythonで実装しながら緑本を学ぶ (第9章 GLMのベイズモデル化と事後分布の推定)

Python 統計学読書メモ

今回は第9章です。PyMC3を使って、GLMをベイズモデルで表現します。実装は以下で公開しています。

introduction_to_machine_learning_with_python/chapter9.ipynb at master · ohke/introduction_to_machine_learning_with_python · GitHub

9 GLMのベイズモデル化と事後分布の推定

9.1 例題：種子数のポアソン回帰(個体差なし)

架空植物20個体において体サイズと種子数を計測した。個体iにおける体サイズを $x_i$ 、種子数を $y_i$ とします。得られたデータから体サイズ $x_i$ と種子数 $y_i$ の関係について調べる、というのが今回の例題です。

計測結果は以下にプロットしています。
なお著者サイトで提供されているサンプルデータはRDataですので、Pythonで扱えるようにするため、PypeRでロード後にDataFrame化しています。PypeRについては以前の投稿も参考にしてください。

f:id:ohke:20180303113433p:plain

種子数は上限が無い離散値ですので、ポアソン分布でばらつきを表現できそうです。
statsmodelsで最尤推定すると、ポアソン分布の平均 $\lambda_i$ は下式で得られます(なお、このデータは作成するときには $\lambda_i=\exp(1.5+0.1x_i)$ から生成されたとのことです)。

$$ \lambda_i = \exp(1.5661 + 0.0833x_i) $$

import statsmodels.formula.api as smf
result = smf.poisson('y ~ x', data=data).fit()
result.summary()

f:id:ohke:20180303120036p:plain

9.2 GLMのベイズモデル化

ベイズモデル化したとしても、中核はポアソン回帰のGLM。

平均 $\lambda_i$ のポアソン分布 $p(y_i \mid \lambda_i)$ に従う
線形予測子と対数リンクを使い、 $\lambda_i =\exp(\beta_1 + \beta_2 x_i)$ で指定する

このモデルでの尤度関数Lは以下となります。 $x_i$ を定数としている(Lのパラメータになっていない)ことに注意してください。

$$ L(\beta_1, \beta_2)=\prod_i p(y_i \mid \lambda_i)=\prod_i p(y_i \mid \beta_1, \beta_2, x_i) $$

ある $\beta_1, \beta_2$ において ${\bf Y}$ ( ${\bf Y}= \{ y_i \}$ )が得られる確率は、尤度関数と一致します(この関係は8.4で導入・説明しています)。

$$ p({\bf Y} \mid \beta_1, \beta_2)=L(\beta_1, \beta_2) $$

ベイズモデルの事後分布は、尤度×事前分布に比例するので、以下の関係が成り立ちます(こちらも8.4を参考にしてください)。

$$ p(\beta_1, \beta_2 \mid {\bf Y}) \propto p({\bf Y} \mid \beta_1, \beta_2)p(\beta_1)p(\beta_2) $$

9.3 無情報事前分布

まずは事前分布 $p(\beta_1)$ と $p(\beta_2)$ (まとめて $p(\beta_{*})$ )を設定する。

データ ${\bf Y}$ が得られていない状態で決める事前分布なので、無情報事前分布と呼ばれます。つまり、どんな値( $[-\infty, +\infty$ ] )でも良いのです。
こうした分布はの生成方法は2つあります。

広い範囲 (例えば[-100000, +100000])の一様分布
平均0で標準偏差が大きい(平べったい)正規分布

分散100の正規分布は以下のようにプロットできます。今回は標準偏差100の正規分布を、事前分として使います。

f:id:ohke:20180303131934p:plain

9.4 ベイズ統計モデルの事後分布の推定

事前分布が定まったので、事後分布 $p(\beta_1, \beta_2 \mid {\bf Y})$ をMCMCサンプリングで推定します。

書籍では、WinBUGS+R2WinBUGSで、RからMCMCサンプリンによる推定が行われています。
ベイズモデルでパラメータ推定できるPythonパッケージとして、PyStanやPyMCが有名です。今回はPyMC3を使います。いつも通り、pip install pymc3でインストールしておきます。

http://docs.pymc.io/

PyMC3によるGLMのパラメータ推定は、3ステップで行います。

モデルを定義する
サンプリングして推定する
得られた結果を確認する

import pymc3 as pm

# モデルを定義する
with pymc3.Model() as model:
    # 事前分布をN(0, 100)の正規分布で設定
    beta1 = pymc3.Normal('beta1', mu=0, sd=100)
    beta2 = pymc3.Normal('beta2', mu=0, sd=100)
    
    # 線形予測子θをβ1+β2xで設定
    theta = beta1 + beta2*data['x'].values
    
    # ログリンク関数(log(μ)=θ⇔μ=exp(θ))を設定し、ポアソン分布で推定する
    y = pymc3.Poisson('y', mu=np.exp(theta), observed=data['y'].values)

# サンプリングして推定する
with model:
    # 101個目から3個置きでサンプルを取得するチェインを3つ作る
    # NUTSではburnとthinが効いていない？
    trace = pymc3.sample(1600, burn=100, thin=100, njobs=3, random_seed=0)

# 得られた結果を確認する
pymc3.traceplot(trace) # サンプリング過程を表示する
pymc3.summary(trace) # 推定結果を表示する

モデルを定義する

まずwith句でモデルクラス(pymc3.Model)を作成し、事前分布、線形予測子、リンク関数、確率分布を指定します。

$\beta_1, \beta_2$ の事前分布には、それぞれ平均0・標準偏差100の正規分布(pymc3.Normal)を使う
線形予測子は $\theta=\beta_1 + \beta_2 x$
ポアソン分布(pymc3.Poisson)で尤度を計算
- ログリンク関数 $\log(\mu)=\theta \Leftrightarrow \mu=\exp(\theta)$ で平均 $\mu$ を設定

import pymc3

# モデルを定義する
with pymc3.Model() as model:
    # 事前分布をN(0, 100)の正規分布で設定
    beta1 = pymc3.Normal('beta1', mu=0, sd=100)
    beta2 = pymc3.Normal('beta2', mu=0, sd=100)
    
    # 線形予測子θをβ1+β2xで設定
    theta = beta1 + beta2*data['x'].values
    
    # ログリンク関数(log(μ)=θ⇔μ=exp(θ))を設定し、ポアソン分布で推定する
    y = pymc3.Poisson('y', mu=np.exp(theta), observed=data['y'].values)

サンプリングしてパラメータ推定する

次に、pymc3.sampleメソッドでサンプラーの定義とサンプリングを行い、パラメータを推定する。

最初の引数(draws)で、サンプル数を指定
stepで、サンプリングアルゴリズムを指定
- Metropolis、HamiltonianMC、NUTSなどが選択できます(デフォルトはNUTS)
tuneで、先頭から捨てるサンプル数を指定する(WinBUGSではn.burninで指定する値)
- サンプルの最初の方は、ランダムに選ばれた初期値の影響を大きく受けるため、捨てた方が良い
njobsで、チェイン数(サンプル列数)を指定
- 3を指定すると、3つの異なる初期値からそれぞれサンプリングが行われる
observedに、観測された従属変数(ここではy)を渡します
2個飛ばしでサンプリング(つまり合計500個)するために、スライス表記[::3]でサンプリング過程を取得

# ハミルトニアンモンテカルロ法
with model:
    # 101個目から3個置きでサンプルを取得するチェインを3つ作る
    trace = pymc3.sample(1500, step=pymc3.HamiltonianMC(), tune=100, njobs=3, random_seed=0)[::3]

得られたサンプリング過程(trace)は、添字でアクセスできます。

print('Trace type:', type(trace)) # Trace type: <class 'pymc3.backends.base.MultiTrace'>
print('Trace length:', len(trace)) # Trace length: 500
print('trace[0]:', trace[0]) # trace[0]: {'beta1': 2.0772965015391716, 'beta2': -0.02971672503615687}

得られた結果を確認する

サンプリング後には、初期値やサンプリング数などの設定値が適切であったかどうかを確認する必要があります。

pymc3.traceplotメソッドで、サンプリング過程がグラフ化できます。

pymc3.traceplot(trace)

メトロポリス法(上図)とハミルトニアンモンテカルロ法(下図、以降はHMCと表記)で並べて比較してみます。

メトロポリス法については、以前の投稿を見てください
HMCについては、こちらの方の記事と、リンクされている伊庭さんのPDFが役に立つかと思います (まだちゃんと調べられてないです)
- ハミルトニアン・モンテカルロ(HMC)法のざっくりとした解説とPythonによる実装 - ともにゃん的データ分析ブログ

右図のパラメータの推移を見ると、HMCではサンプル列同士が近づきつつありますがまだ不安定な状態にある(つまりサンプリング数が不足している)ことがわかります。対して、HMCではいずれのサンプル列同士も十分近づき、類似した波形となっています。
いずれも、3つのサンプル列がプロットされていることや、サンプル数が500(1500から2つ飛ばしでサンプリングしているため)となっていることに注意してください。

メトロポリス法
HMC

また、pymc3.summaryメソッドで、各パラメータの推定値とそれぞれの統計値を確認できます。

pymc3.summary(trace)

こちらもメトロポリス法(上表)とHMC(下表)で並べてみます。
各項目の詳細は次節で見ていきますが、このうちRhatで表記されている $\hat{R}$ 指数はサンプル列間のばらつきを表す値で、パラメータ毎に求められます。この値が1に近いほどサンプル列間のばらつきよりも列内のばらつきが大きくなるので、収束していると言えます。経験的には $\hat{R} \le 1.1$ が1つの目安ですが、 $\hat{R}$ も十分なサンプル数でない場合は安定した結果が得られないので注意が必要です。
メトロポリス法では1.007、HMCでは0.999なので、後者の方が僅かですがより収束していると言えそうです。

メトロポリス法
HMC

9.5 MCMCサンプルから事後分布を推定

得られた推定結果から事後分布を確認します。

HMCで得られたサンプリング過程を再掲します。
左図は、 $\beta_1$ と $\beta_2$ の周辺事後分布で、カーネル密度推定で近似された確率密度関数で表現されています。周辺事後分布は、あるパラメータ1つに関する事後分布で、ここでは $p(\beta_1 \mid {\bf Y})$ と $p(\beta_2 \mid {\bf Y})$ となります。

f:id:ohke:20180304125331p:plain

周辺事後分布に対して、 $p(\beta_1, \beta_2 \mid {\bf Y})$ は同時事後分布と呼ばれます。
パラメータの組み合わせを散布図にプロットします。 $\beta_1$ と $\beta_2$ の相関がかなり強いサンプリングが行われていることがわかります。本書の傾向と大きく異なっており、サンプリングアルゴリズム(本書ではWinBUGSのギブスサンプリング実装)の違いによるものと思われます。

f:id:ohke:20180304141715p:plain

なお、サンプリングした値はtrace(MultiTraceクラス)からget_valuesメソッド(返り値の型はnumpy.ndarray)で取得できます。

1番目の引数(varname)にパラメータ名を指定
chainsオプションに、サンプル列のインデックスを渡す

# サンプル列数分だけ繰り返す
for i in trace.chains:
    # 各サンプル列のパラメータの平均値を計算
    beta1_averages += trace.get_values('beta1', chains=i) / trace.nchains
    beta2_averages += trace.get_values('beta2', chains=i) / trace.nchains

次に、HMCの統計量から事後分布の推定値を確認します。

パラメータは平均1.5599で、95%信用区間は0.8630〜2.2622
- $\beta_1$ は95%の事後確率で0.8630〜2.2622に収まる、と解釈できる
パラメータ $\beta_2$ は平均0.0826で、95%信用区間は0.0017〜0.2145
n_effは有効なサンプルサイズで、サンプル間の相関が高いと、この値が小さくなります

f:id:ohke:20180303232202p:plain

9.6 複数パラメーターのMCMCサンプリング

8章では、1パラメータについてメトロポリス法でMCMCサンプリングする方法を説明しました。

メトロポリス法には、更新前と更新後の値の相関が強く、なかなか収束しないという問題がありました。
この問題を解決するサンプリング方法の1つにギブスサンプリングがあります。ギブスサンプリングは、「新しい値の確率分布を作ってその確率分布からランダムに選択する」という方法で値を更新することで、更新前後の値の相関を弱くします。

さらに今回の例題のように、複数のパラメータ( $\beta_1, \beta_2$ )のMCMCサンプリングを考える必要があります。
こうした場合、全てのパラメータを同時に更新するよりも、 $\beta_1$ と $\beta_2$ を交互に少しずつ更新していく方が簡単です。

2パラメータのギブスサンプリングをまとめると以下のアルゴリズムとなります。

(1) の適当な初期値を設定する
- 例えば ${\beta_1, \beta_2}={1.5, 0}$ と置く
(2) に従う乱数を発生させ、得られた値を新しいとする
- 他の変量を全て定数とする一変量確率分布で、全条件付き分布(FCD)と呼ばれる
- FCDに従う乱数1つを発生させる( が得られたとする)
  - FCDからのサンプリングする方法は理解できませんでした

$$ p(\beta_1 \mid {\bf Y}, \beta_2=0.0) \propto \prod_i \frac{\lambda_i^{y_i} \exp(-\lambda_i)}{y_i !}p(\beta_1) $$

(3) $p(\beta_2 \mid {\bf Y}, \beta_1)$ に従う乱数を発生させ、得られた値を新しい $\beta_2$ とする (2の逆)

$$ p(\beta_1 \mid {\bf Y}, \beta_1=2.052) \propto \prod_i \frac{\lambda_i^{y_i} \exp(-\lambda_i)}{y_i !}p(\beta_2) $$

(4) この新しい ${\beta_1, \beta_2}$ を記録する
(5) 十分なサンプル数が得られるまで(2)〜(4)を繰り返す

2018-03-02

Python: foliumでJupyter Notebookに地図を描画する

Python

Jupyter Notebook上で、緯度経度の情報を地図へ簡単にプロットできる方法を探していたところ、foliumの使い勝手が良かったので紹介します。

folium

PythonからLeaflet.jsで地図をプロットするパッケージです。

https://github.com/python-visualization/folium

使用前にpip install foliumまたはconda install foliumでインストールします。

地図を表示する

まずは、千葉県千葉市の地図を表示します。

Mapで生成できます。

locationオプションに中心の緯度・経度を渡しています
Mapインスタンスを出力すると、leaflet.jsが読み込まれて表示されます
デフォルトではOpenStreetMapの地図が使われます(tilesオプションでMapboxやStamenにも変更できます)
zoom_startで表示する時のズームを変えられます(デフォルト10で、大きくするとより拡大されます)

import folium

chiba_map = folium.Map(location=[35.607451, 140.106340])

chiba_map

f:id:ohke:20180228092918p:plain

マーカをプロットする

次に千葉県の4つの市にマークします(緯度経度は市役所所在地、人口は千葉県のサイトから2018/1/1時点のデータを参照)。

市	緯度	経度	人口
千葉市	35.607451	140.106340	975,535
館山市	34.996596	139.869906	46,349
銚子市	35.734795	140.826926	61,674
浦安市	35.653146	139.902058	168,169

Markerで地図上にマーカをプロットできます。

popupオプションで、マーカをクリックした時表示されるテキストを設定できます
add_toメソッドでマップ上にプロットします

import pandas as pd

chiba_cities = pd.DataFrame({
    'city': ['千葉市', '館山市', '銚子市', '浦安市'],
    'latitude': [35.607451, 34.996596, 35.734795, 35.653146],
    'longtude': [140.106340, 139.869906, 140.826926, 139.902058],
    'population': [975535, 46349, 61674, 168169]
})

chiba_map = folium.Map(location=[35.607451, 140.106340], zoom_start=9)

for i, r in chiba_cities.iterrows():
    folium.Marker(location=[r['latitude'], r['longtude']], popup=r['city']).add_to(chiba_map)
    
chiba_map