[Revised] Bayesian Modeling for Chapter 4 "How Many Others Engage in Sensitive Behaviors" using PyMC Ver. 5

This article is a revised version of the "PyMC article for Chapter 4 'How Many Others Engage in Sensitive Behaviors'" that I posted previously.

You can find the previous article here.

At the time I wrote the previous article, I had a bitter experience because I could not get the posterior distribution to converge.
However, the situation has changed completely.
Thanks to the evolution of generative AI, I have found a solution!

So, take two, let's begin!

Restart

This article is a statistical documentary depicting the process of "experimentally" implementing the Bayesian model from Chapter 4, "How Many Others Engage in Sensitive Behaviors" of the text "Fun Bayesian Modeling 2" using PyMC Ver. 5.

The text covers Bayesian modeling using various "indirect questioning methods" to analyze survey data on "experiences of illegal acts, etc." from Chapter 2 to Chapter 4.
In this chapter, we will infer experiences of sensitive behaviors using the "Aggregated Response Method" (AR method).

I first encountered the AR method in Chapter 1 of the previous book, "Fun Bayesian Models."
At that time, I had a bitter experience because the output of my own model differed from the results in the text.
Now, with a fresh mind, let's enjoy Bayesian modeling with PyMC.

Please refer to the linked article for the introduction of the text, citation format, series preface, and version information for PyMC, etc.

The data used in the text, along with sample scripts for R, Stan, etc., can be obtained by downloading them from the publisher's website.

Summary

Overview of the text

Authors: Dr. Takashi Akiyama, Dr. Hideki Toyoda, Dr. Saori Kubo, Dr. Ritsuko Oka
Model Difficulty: ★★★.. (Average)

Self-Evaluation

Rating

$$
\begin{array}{c:c:c}
Implementation Accuracy & ★★★★★ & GoooD! \\
Result Reproducibility & ★★★★★ & Best! \\
Fun Factor & ★★★★★ & Fun! \\
\end{array}
$$

Evaluation Points

• In my previous article, I couldn't resolve the error occurring with the custom distribution function in the Poisson distribution model, and I ended the performance feeling unsatisfied.

• However, when I simply threw the error message and the model code at a generative AI, it produced code that worked with ease.
  We live in an era where we can easily consult with reliable partners.

Ingenuity, Joy, and Reflection

• I was able to successfully close the remaining homework!
  I did it! ✨️

Model Overview

Overview of Text Research and Experiments

■ Indirect Questioning Method
The indirect questioning method is a technique for obtaining answers to sensitive questions while protecting the respondent's privacy and acquiring the data necessary for analysis.
I am borrowing the definition from Chapter 2 of the text.

The research in this chapter deals with the number of experiences and number of people involved in sensitive behaviors of university students (count data).
Applying this to the definition above, it can be rephrased as "a method to obtain estimated experience levels for an entire group of university students while being considerate so that the researcher cannot know the individual response values of the students even when viewing the results."

■ AR Method
This is one of the techniques of the indirect questioning method.
Please look at the figure below.
Respondents answer the question using "True Value ± Mask Value".

The numerical value added or subtracted from the true value, such as the last digit of a mobile phone number, is called a "mask".
The mask acts as a "concealment" and functions as privacy protection, making it easier for respondents to answer.

■ Research Overview
The survey items are "number of people experienced in sexual activity" (hereinafter abbreviated as S) and "number of times experienced cheating" (hereinafter abbreviated as C).
They created a question list using Google Forms, collected responses from survey participants on the web, and proceeded with the analysis using 130 data points for S and 125 for C.

Modeling in the Text

There is no model included in this chapter, and it is written to refer to Chapter 1 of "Fun Bayesian Modeling".
So, please also refer to the blog for Chapter 1 of "Fun Bayesian Modeling"!
The section on "Modeling in the Text" describes in detail how the mask relates to modeling.

The following modeling content is written with reference to Chapter 1 of "Fun Bayesian Modeling".

■ Objective variable and parameters of interest
The objective variable is the masked response value $${y}$$.
The parameters of interest are the true number of people/times $${{\lambda}}$$ excluding the mask.

■ Bayesian model of the text
I will cite the model from Chapter 1 of the text "Enjoyable Bayesian Modeling".
In this chapter, the assumption that the response value follows a Poisson distribution with an average parameter $${{\lambda}}$$ is reflected in the following mathematical formula.

■ Mask $${m_k}$$, set $${A}$$, and function $${N(k \subset A)}$$
The set $${{\boldsymbol{A}}}$$ of possible values for the mask $${{\boldsymbol{m_k}}}$$ varies depending on the value of the response $${{\boldsymbol{y}}}$$.For example, if the response value is $${y=2}$$, the last digit of the mobile phone number (i.e., the mask) to be added or subtracted ranges from $${-9}$$ (coin tails & last digit 9) to $${+2}$$ (coin heads & last digit 2).

Now, what happens in the case of coin heads & last digit 3?
Since the minimum value for the number of experiences is $${0}$$, the response must be $${3}$$ or greater, so it cannot fit into the response value $${y=2}$$.
Therefore, the mask value cannot be $${+3}$$ or greater.
In summary, for a response value of $${y=2}$$, the set of mask values is $${A=\{ -9, -8, \cdots, -1, -0, +0, 1, 2\}}$$, and the number of elements $${N(k \subset A)}$$ is $${13}$$.

■ Being a mixture distribution
$${f(y-m_k \mid \cdot \ )}$$ indicates that there are as many probability distributions as there are elements in the set $${A}$$ (i.e., the number of masks), which varies depending on the value of $${y}$$$$.
Depending on the value of the response $${{\boldsymbol{y}}}$$, the "number of probability distributions" being mixed changes!
For example, in the case of $${y=2}$$, the mixture distribution is an equal mixture of the probability mass functions of the following $${13}$$ Poisson distributions, each with a weight of $${1/13}$$.

* Six functions in the middle are omitted

■ Analysis and analysis results
I believe the descriptions in the text regarding the analysis methods and numerical values are accurate, so I recommend reading the text.
Please see the "PyMC Implementation" chapter for analysis using inferred values from my own PyMC model.

PyMC Implementation

Let's enjoy PyMC & Python !

Preparation

1. Import

# インポート

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import pymc as pm
import arviz as az

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

S data preparation

Check the response values regarding the number of experiences of sexual activity S.

1. Loading data
Set the response value data described in the R script of the text.

# データの読み込み ※Rスクリプト掲載データを引用（昇順ソートしています）

data = [-9, -9, -9, -8, -8, -8, -8, -8, -8, -8, -7, -7, -7, -7, -6, -6, -6,
        -5, -5, -5, -5, -5, -4, -4, -4, -3, -3, -3, -3, -3, -2, -2, -2, -2, -2,
        -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
        4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5,
        6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
        8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 12, 15, 15]

print('標本サイズ: ', len(data))

[Execution Result]
The sample size is 130.

2. Checking the appearance of the data
First, let's check the summary statistics.

# 要約統計量の表示
pd.DataFrame({'回答値（マスク済み）': data}).describe().round(2)

[Execution Result]
These are response values masked by the last digit of the mobile phone number (0 to 9).
The minimum value of $${-9}$$ is likely the maximum mask value of 9 subtracted, so the number of experienced people is probably 0.
Regarding the maximum value of $${+15}$$, if the mask is +9, the number of experienced people is 6, and if the mask is -9, the number of experienced people is 24.

Let's visualize it with a box plot.
This corresponds to Figure 4.1 in the text.

# 箱ひげ図の描画 ※図4.1に相当

# 描画領域の設定
plt.figure(figsize=(3, 4))
# 箱ひげ図の描画
sns.boxplot(data, color='tab:blue', fill=False)
# スウォームプロットの描画
sns.swarmplot(data, color='blue', alpha=0.3)
plt.ylabel('回答値：経験人数＋マスク');

[Execution Result]
I overlaid a swarm plot to make the distribution of data points easier to see.
It gives the impression of being evenly distributed from about -9 to 9.
Response value = 7 is the most frequent response.

Let's calculate the mean and standard error of the response values $${y}$$ shown in Table 4.1 of the text.

# 平均値と標準誤差の計算 「不偏分散 ÷ 標本サイズ」の平方根 ※表4.1のyに相当
mean1 = np.mean(data)
se1 = np.sqrt(np.var(data, ddof=1) / len(data))
print(f'平均:{mean1:.3f}, 標準誤差:{se1:.3f}')

[Execution Result]

Construction of the Poisson distribution model for S

Mathematical expression of the model
This is a "pseudo-mathematical" notation that mixes in the feel of the PyMC model I want to aim for.

$$
\begin{align*}
\lambda &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
components &= [ \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)-9, \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)-8, \\
&\quad \vdots \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)-0, \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)+0, \\
&\quad \vdots \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)+8, \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)+9, \\
&] \\
likelihood &\sim \text{Mixture}\ (\text{w}=w,\ \text{comp\_dist}=components)\\
\end{align*}
$$

1. Initial value setting
Calculate the mixture ratio $${w}$$ of the mixture distribution.
Find the number of masks $${N}$$ according to the value of the response value $${y}$$.
Since there are 20 possible values for the mask from $${-9}$$ to $${9}$$, the array of weights $${w}$$ consists of 20 elements per row.
In the code, I use a for loop to set $${1/N}$$ for the first $${N}$$ elements of the array.
The remaining elements are set to $${0}$$ so that the mixture distribution with weight 0 is not used.

# 初期値設定

# データのインデックス、カテゴリ変数のインデックスの取得
idx_data = list(range(len(data)))

# 混合分布の混合比率w：観測データ（全130行）ごとにマスク個数20個分の均等割合を算出
# wの初期化
w = np.zeros((len(data), 20))
# 観測値yごとに可能なマスク個数を算出
num_masks = [y+10 if y<=-1 else y+11 if y<=9 else 20 for y in data]
# w(numpy配列)の作成
for i, num_mask in enumerate(num_masks):
    # 配列先頭からマスクの個数までの要素を1/Nにする
    w[i, :num_mask] = 1 / num_mask
# wの表示
print('混合比率w: shape=', w.shape)
print(w)

[Execution Result]
A portion of the mixture ratio $${w}$$ is displayed.

2. Model definition
I will describe the model faithfully to the mathematical expression.
The Poisson distribution + mask m is implemented using the custom distribution class pm.CustomDist.dist.

# モデルの定義

# logp（対数尤度）関数の定義 ★変更点
def custom_logp1(value, lam, m):
    '''
    value: 観測値 (PyTensor変数が渡される)
    lam: パラメータ
    m: マスク
    '''
    # 「観測値 value から m を引いた値がポアソン分布(λ)に従う」という設定
    return pm.logp(pm.Poisson.dist(mu=lam), value - m)

# coordsの定義
coords = {'data': idx_data}

# モデルの定義
with pm.Model(coords=coords) as model1:

    # Dataの定義：目的変数
    y = pm.Data('y', value=data, dims='data')
    
    # 事前分布の定義
    lam = pm.Uniform('lam', lower=0, upper=100)
    
    # 混合分布の構成要素の定義
    # マスクのfull集合 [-9, -8, …, 0, 0, …, 9]
    masks = np.hstack([np.arange(-9, 1), np.arange(0, 10)])
    # リスト内包表記で20個(マスク個数分)のカスタムポアソン分布を生成 ★変更点
    components = [
        pm.CustomDist.dist(
            lam, m, 
            logp=custom_logp1,            # logp関数を指定
            dist=lambda lam, m, size: (   # 標本サンプリング関数を指定
                pm.Poisson.dist(mu=lam, size=size) + m
            )
        ) for m in masks
    ]
    
    # 尤度：混合分布
    likelihood = pm.Mixture('likelihood', w=w, comp_dists=components,
                            observed=y, dims='data'))

[Model Annotations]

• Definition of coord
  You can name coordinates and set the values those coordinates can take.
  This time, I set the following one.

  • Coordinate of the data row: name "data", value "row index"

• Definition of data
  I set the response value $${y}$$.

• Prior distribution of parameters

  • The prior distribution for the number of experienced people $${lambda}$$ regarding the content of the key item is a uniform distribution in the interval $${[0,100]}$$.} , {

• Likelihood

  • It is a mixture distribution of 20 Poisson distributions.

  • I implemented the mixture distribution as follows.

    • Specify the mixture distribution "pm.Mixture" in the likelihood.
      Refer to "components" for the composition of the mixture distribution.
      Specify the weight (mixing ratio) of each distribution with w.

    • Using the following settings in the components list, I generated 20 custom probability distributions "pm.Custom.dist()" using list comprehensions.

      • In the log probability logp, specify the Poisson distribution and parameter $${λ}$$ for the "difference between the response value $${y}$$ and the mask m".

      • In the distribution dist, use a Python anonymous function lambda to specify the sample sampling function for the response value $${y}$$ based on "Poisson distribution + mask m".

By using logp, I was able to avoid the error that occurred in the previous article!
And by obtaining logp, the expressive power of the custom distribution has been powered up!

Note that the estimated value of the parameter $${λ}$$ is the "difference between the response value $${y}$$ and the mask m", that is, the average value of "the number of people S has experienced".

3. Checking the model appearance

# モデルの表示
model1

[Execution result]
The right side of the mixture distribution in the second line is cut off.
It is composed of 20 custom distributions (Poisson distribution + mask).

# モデルの可視化
pm.model_to_graphviz(model1)

[Execution result]
The figure has a very simple structure.

4. Sampling from the posterior distribution
Execute MCMC.
I used "nutpie" for the NUTS sampler.
nutpie official site▶ https://pymc-devs.github.io/nutpie/
The processing time was about 50 seconds.

# 事後分布からのサンプリング ※NUTSサンプラーにnutpieを使用
# テキスト：iter=21000, warmup=1000, chains=5
with model1:
    idata1 = pm.sample(draws=20000, tune=1000, chains=5,
                       nuts_sampler='nutpie', random_seed=1234)

[Execution result]
In the previous article, an error occurred during this MCMC sampling process...
This time, there were no errors, no divergences, and the process completed successfully!

5. Checking sampling data
Check $${R̂}$$, summary statistics of the posterior distribution, and trace plots.
Following the text, I set $${R̂≤1.1}$$ for checking the convergence of the posterior distribution.

# r_hat>1.01の確認

# 設定
idata_in = idata1        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

[Execution result]
For now, I checked with a threshold of $${R̂>1.01}$$.
There are 0 parameters with $${R̂>1.01}$$.
I was able to confirm that all parameters are $${R̂ ≤ 1.01}$$.

I will roughly check the summary statistics and trace plots of the posterior distribution samples.

# 推論データの要約統計情報の表示
var_names = ['lam']
pm.summary(idata1, hdi_prob=0.95, var_names=var_names, round_to=3)

[Execution result]
The posterior mean of the average parameter $${λ}$$ for the number of people S has experienced is $${1.547}$$ people.

# トレースプロットの表示
pm.plot_trace(idata1, var_names=var_names, compact=False, figsize=(8, 2))
plt.tight_layout();

[Execution result]
The graph on the right is drawn evenly, and it feels like it has converged.
The five colors in the graph on the left = the five Markov chains overlap, and they are producing almost identical inference values.

Let's visualize the posterior distribution of the average parameter $${λ}$$.
Draw a histogram with the posterior distribution sampling data.
Create a graph with a KDE curve using seaborn's histplot.

# パラメータの事後分布の可視化 ※図4.3に相当

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))

# λのヒストグラムの描画
sns.histplot(idata1.posterior.lam.data.flatten(), bins=30, kde=True,
             stat='density', color='steelblue', edgecolor='white', alpha=0.5,
             ax=ax)
ax.set(title='λの事後分布', xlabel='λ', ylabel='density')
plt.show()

[Execution result]
It is an almost symmetrical distribution!

Analysis of S's Poisson distribution model

Regarding the number of people S has experienced, I will proceed with the analysis following the text.

1. Posterior distribution
Display the posterior statistics of the average parameter $${λ}$$.
This corresponds to Table 4.1 in the text.

# 推論データの事後統計量の算出 ※表4.1に相当

# 要約統計量計算関数の定義
def calc_stats(x):
    ci = np.quantile(x, q=[0.025, 0.975])
    hdi = az.hdi(x, hdi_prob=0.95)
    return np.mean(x), np.std(x), ci[0], ci[1], hdi[0], hdi[1]

# 要約統計量データフレームの作成
stats_df1 = pd.DataFrame(
    {'λ': calc_stats(idata1.posterior.lam.data.flatten())},
    index = ['EAP', 'post.sd', 'CI2.5%', 'CI97.5%', 'HDI2.5%', 'HDI97.5%']
    ).T

# データフレームの表示
display(stats_df1.round(3))

[Execution result]
The result was almost the same as the text!

[Analysis]
The mean value (EAP) of the parameter $${λ}$$ is $${1.547}$$, the posterior standard deviation (post.sd) is $${0.297}$$, and the 95% credible interval is $${[1.008,\ 2.171]}$$.
Based on the posterior distribution, we can say with 95% confidence that the number of experienced people is between $${1.008}$$ and $${2.171}$$.

2. Posterior Predictive Distribution
Next, let's look at the posterior predictive distribution.
In this case, we use MCMC samples to generate posterior predictive samples of the number of experienced people $${x^*}$$ using Poisson random numbers from numpy.
$${x^*}$$ includes the uncertainty (variance) of the Poisson distribution.

We will plot the posterior predictive distribution corresponding to Figure 4.4.

# 事後予測分布の可視化 ※図4.4に相当

# 推定された lam のサンプリング結果（全サンプル）を取得
post_lam = idata1.posterior['lam'].values.flatten()

# 推定された lam を使って、マスクなしの純粋なポアソン乱数を生成
rng = np.random.default_rng(seed=1234)  # 乱数生成器の初期化
pure_y_samples = rng.poisson(post_lam)

# 事後予測分布のヒストグラムの描画
plt.hist(pure_y_samples, bins=range(0, max(pure_y_samples) + 3),
         rwidth=0.3, align='left', color='tab:blue')
plt.title('Sの経験人数の事後予測分布')
plt.xlabel('経験人数の事後予測値', fontsize=12)
plt.ylabel('頻度', fontsize=12)
plt.show()

[Execution Result]
We can see that the number of experienced people is most often between $${0}$$ and $${2}$$!

Let's calculate the posterior predictive distribution and cumulative distribution from Table 4.2.

# 事後予測分布と累積分布 ※表4.2に相当
unique, counts = np.unique(pure_y_samples, return_counts=True)
ratios = counts / counts.sum()
pd.DataFrame(
    {'x*': ratios, '累積分布': np.cumsum(ratios)}, index=unique
).T.round(2)

[Execution Result]

We will check the mean and unbiased variance of $${x^*}$$.

# x*の平均と不偏分散の算出
print(f'x*の平均　　: {np.mean(pure_y_samples):.3f}')
print(f'x*の不偏分散: {np.var(pure_y_samples, ddof=1):.3f}')

[Execution Result]
The unbiased variance feels slightly larger compared to the text.
(The unbiased variance in the text is $${1.278}$$)

3. PHC Curve: Probability that the research hypothesis is correct
Finally, we will draw a PHC curve using the MCMC samples of the posterior distribution.
The PHC curve places the threshold $${c}$$ on the horizontal axis and plots the ratio of MCMC samples that are greater than $${c}$$.
It shows the probability that the average number of experienced people $${λ}$$ is greater than $${c}$$ people.

We will define a function to draw the PHC curve.

# 推論データidataを受け取ってPHC曲線を描画する関数

def phc_plot(data, xrange=10):
    
    ## データ加工
    # x軸目盛値となるcの値を生成
    x = np.linspace(0, xrange, 1001)
    # λ, μ > c の割合を計算
    y = [(data > c).sum() / len(data) for c in x]
    
    ## 描画処理
    plt.plot(x, y)
    plt.xlabel('c')
    plt.ylabel('確率:PHC')

We will draw the PHC curve.
This corresponds to Figure 4.5.

# PHC曲線の描画 ※図4.5に相当
plt.figure(figsize=(6, 4))
phc_plot(idata1.posterior.lam.data.flatten(), xrange=5)
plt.title(f'「平均値 λ > c 」の確率：PHC曲線')
plt.legend(['経験人数のPHC'])
plt.grid(lw=0.3);

[Execution Result]
The probability drops sharply when the number of experienced people is between 1 and 2.

Following the text, we will calculate the probability of more than $${c=1}$$: 1 person, and the probability of more than $${c=2}$$: 2 people.

# Sの経験人数の平均が1より大きい（p>1）というPHCの算出
p = idata1.posterior.lam.data.flatten()
c = 1
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

[Execution Result]
The probability that the average number of experienced people is greater than 1 is $${97.7\%}$$.

# Sの経験人数の平均が2より大きい（p>2）というPHCの算出
p = idata1.posterior.lam.data.flatten()
c = 2
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

[Execution Result]
The probability that the average number of experienced people is greater than 2 is $${7.1\%}$$.

Visually, we can promote an intuitive understanding that "there is likely at least one experienced person"!
(The posterior mean EAP of the average number of experienced people was $${1.547}$$ people.)

That concludes the analysis of S.

Data preparation for C

Next, we will proceed to modeling the number of times cheating C was experienced.
Following the text, we will work on the "Poisson distribution model" and the "Normal distribution model".

Note that for the AR method of the number of times C was experienced, the mask is set as the "birth month".
The elements of the mask are $${-12, -11, \cdots, -1, 1, \cdots, 11, 12}$$, and the number of elements is 24.

1. Loading the data
Set the response value data described in the text's R script.

# データの読み込み ※Rスクリプト掲載データを引用（昇順ソートしています）

data2 = [-11, -8, -7, -7, -7, -7, -7, -7, -7, -6, -5, -5, -5, -5, -5,
         -4, -4, -4, -4, -4, -4, -3, -3, -3, -3, -3, -3, -3,
         -2, -2, -2, -2, -2, -1, -1, -1, -1, -1, -1,
         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,
         2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
         4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6,
         7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10,
         11, 12, 16, 18, 19, 19, 20, 24]

print('標本サイズ: ', len(data2))

[Execution Result]
The sample size is 125.

2. Checking the appearance of the data
First, let's check the summary statistics.

# 要約統計量の表示
pd.DataFrame({'回答値（マスク済み）': data2}).describe().round(2)

[Execution Result]
These are response values masked by birth month (1 to 12).
The minimum value of $${-11}$$ is likely 11 or 12 subtracted from the mask, so the number of experiences is probably 0 or 1.
Regarding the maximum value of $${+24}$$, if the mask is +12, the number of experiences is 12, and if the mask is -12, the number of experiences is 36.

Let's visualize it with a box plot.
This corresponds to Figure 4.6 in the text.

# 箱ひげ図の描画 ※図4.6に相当

# 描画領域の設定
plt.figure(figsize=(3, 4))
# 箱ひげ図の描画
sns.boxplot(data2, color='tab:blue', fill=False)
# スウォームプロットの描画
sns.swarmplot(data2, color='blue', alpha=0.3)
plt.ylabel('回答値：経験回数＋マスク');

[Execution Result]
I overlaid a swarm plot to make the distribution of data points visible.
It is generally contained between about -8 and 12, and some outlier-like values can be seen outside the whiskers.
A response value of 0 is the most frequent answer.
Cases where the response value is 0 correspond to instances where the number of experiences is 1 to 12, since the mask is $${-1}$$ to $${-12}$$.

Let's calculate the mean and standard error of the response values $${y}$$ shown in Table 4.3 of the text.

# 平均値と標準誤差の計算 「不偏分散 ÷ 標本サイズ」の平方根 ※表4.3のyに相当
mean2 = np.mean(data2)
se2 = np.sqrt(np.var(data2, ddof=1) / len(data2))
print(f'平均:{mean2:.3f}, 標準誤差:{se2:.3f}')

[Execution Result]
The mean of $${2.272}$$ is the estimated number of experiences using the conventional method.

Model consideration for C

Next, I will repurpose the S model and arrange it for C.

Mathematical expression of the model
This is a "pseudo-mathematical" notation that mixes in the feel of the PyMC model I want to aim for.
I am using a normal distribution instead of a Poisson distribution.

$$
\begin{align*}
\mu &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
\sigma &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
components &= [ \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)-12, \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)-11, \\
&\quad \vdots \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)-1, \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)+1, \\
&\quad \vdots \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)+11, \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)+12, \\
&] \\
likelihood &\sim \text{Mixture}\ (\text{w}=w,\ \text{comp\_dist}=components)\\
\end{align*}
$$

Initial value setting
Calculate the mixing ratio $${w}$$ of the mixture distribution.
Find the number of masks $${N}$$ according to the value of the response $${y}$$.
Since there are 24 possible values for the mask from $${-12}$$ to $${12}$$, the array of weights $${w}$$ consists of 24 elements per row.
In the code, I set $${1/N}$$ for the first $${N}$$ elements of the array using a for loop.
The remaining elements are set to $${0}$$ so that mixture distributions with a weight of 0 are not used.

# 初期値設定

# データのインデックス、カテゴリ変数のインデックスの取得
idx_data = list(range(len(data2)))

# 混合分布の混合比率w：観測データ（全125行）ごとにマスク個数24個分の均等割合を算出
# wの初期化
w = np.zeros((len(data2), 24))
# 観測値yごとに可能なマスク個数を算出
num_masks = [y+13 if y<=-1 else y+12 if y<=12 else 24 for y in data2]
# w(numpy配列)の作成
for i, num_mask in enumerate(num_masks):
    # 配列先頭からマスクの個数までの要素を1/Nにする
    w[i, :num_mask] = 1 / num_mask
# wの表示
print('混合比率w: shape=', w.shape)
print(w)

[Execution Result]
A portion of the mixing ratio $${w}$$ is displayed.

Construction of the Poisson distribution model for C

First, I will perform the modeling for the Poisson distribution model.
I will create it by repurposing the code for the Poisson distribution model of S from earlier.

1. Definition of the model

# モデルの定義 ポアソン分布モデル

# logp（対数尤度）関数の定義
def custom_logp2(value, lam, m):
    '''
    value: 観測値 (PyTensor変数が渡される)
    lam: パラメータ
    m: マスク
    '''
    # 「観測値 value から m を引いた値がポアソン分布(λ)に従う」という設定
    return pm.logp(pm.Poisson.dist(mu=lam), value - m)

# coordsの定義
coords = {'data': idx_data}

# モデルの定義
with pm.Model(coords=coords) as model2:

    # Dataの定義：目的変数
    y = pm.Data('y', value=data2, dims='data')
    
    # 事前分布の定義
    lam = pm.Uniform('lam', lower=0, upper=100)
    
    # 混合分布の構成要素の定義
    # マスクのfull集合 [-12, -1, …, -1, 1, …, 12]
    masks = np.delete(np.arange(-12, 13), 12)
    # リスト内包表記で24個(マスク個数分)のカスタムポアソン分布を生成
    components = [
        pm.CustomDist.dist(
            lam, m, 
            logp=custom_logp2,            # logp関数を指定
            dist=lambda lam, m, size: (  # 標本サンプリング関数を指定
                pm.Poisson.dist(mu=lam, size=size) + m
            )
        ) for m in masks
    ]
    
    # 尤度：混合分布
    likelihood = pm.Mixture('likelihood', w=w, comp_dists=components,
                            observed=y, dims='data')

2. Checking the appearance of the model

# モデルの表示
model2

[Execution Result]
The mixture distribution on the second line is cut off on the right.
It is composed of 24 custom distributions (Poisson distribution + mask).

# モデルの可視化
pm.model_to_graphviz(model2)

[Execution Result]
It has the same appearance as the Poisson distribution of S from earlier.

3. Sampling from the posterior distribution
The processing time was approximately 1 minute.

# 事後分布からのサンプリング ※NUTSサンプラーにnutpieを使用
# テキスト：iter=21000, warmup=1000, chains=5
with model2:
    idata2 = pm.sample(draws=20000, tune=1000, chains=5,
                       nuts_sampler='nutpie', random_seed=1234)

[Execution Result]
There are 0 divergences.

4. Checking sampling data
We will check the R-hat, summary statistics of the posterior distribution, and trace plots.
Following the text, we will use R-hat <= 1.1 to confirm convergence of the posterior distribution.

# r_hat>1.01の確認

# 設定
idata_in = idata2        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

[Execution Result]
For now, I checked with a threshold of R-hat > 1.01.
There are 0 parameters with R-hat > 1.01.
I was able to confirm that all parameters satisfy R-hat <= 1.01.

I will roughly check the summary statistics and trace plots of the posterior distribution sampling data.

# 推論データの要約統計情報の表示
var_names = ['lam']
pm.summary(idata2, hdi_prob=0.95, var_names=var_names, round_to=3)

[Execution Result]
The posterior mean (EAP) of the average parameter lambda for the number of experiences C is 5.008.

# トレースプロットの表示
pm.plot_trace(idata2, var_names=var_names, compact=False, figsize=(8, 2))
plt.tight_layout();

[Execution Result]
The graph on the right is drawn evenly, and it feels like it has converged.
The five colors on the left graph, representing the five Markov chains, overlap and produce almost identical inferred values.

Analysis of the Poisson distribution model for C

1. Posterior distribution
First, let's look at the summary statistics of the posterior distribution, which corresponds to Table 4.3 in the text.

# 推論データの事後統計量の算出 ※表4.3に相当

# 要約統計量データフレームの作成
stats_df2 = pd.DataFrame(
    {'λ': calc_stats(idata2.posterior.lam.data.flatten())},
    index = ['EAP', 'post.sd', 'CI2.5%', 'CI97.5%', 'HDI2.5%', 'HDI97.5%']
    ).T

# データフレームの表示
display(stats_df2.round(3))

[Execution Result]

[Analysis]
The mean value (EAP) of parameter lambda is 5.008, the posterior standard deviation (post.sd) is 0.537, and the 95% credible interval is [3.972, 6.079].
Based on the posterior distribution, we can say with 95% confidence that the number of experiences is between 3.972 and 6.079.

Let's visualize the posterior distribution.

# パラメータの事後分布の可視化

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))

# λのヒストグラムの描画
sns.histplot(idata2.posterior.lam.data.flatten(), bins=30, kde=True,
             stat='density', color='steelblue', edgecolor='white', alpha=0.5,
             ax=ax)
ax.set(title='λの事後分布', xlabel='λ', ylabel='density')
plt.show()

[Execution Result]
It's a beautiful bell shape!

2. Posterior predictive distribution
Using the MCMC samples, let's generate and plot samples for the posterior predictive distribution of the number of experiences x*.

# 事後予測分布の可視化

# 推定された lam のサンプリング結果（全サンプル）を取得
post_lam2 = idata2.posterior['lam'].values.flatten()

# 推定された lam を使って、マスクなしの純粋なポアソン乱数を生成
rng = np.random.default_rng(seed=1234)  # 乱数生成器の初期化
pure_y_samples2 = rng.poisson(post_lam2)

# 事後予測分布のヒストグラムの描画
plt.hist(pure_y_samples2, bins=range(0, max(pure_y_samples2) + 2),
         rwidth=0.3, align='left', color='tab:blue')
plt.title('Cの経験人数の事後予測分布')
plt.xlabel('経験人数の事後予測値', fontsize=12)
plt.ylabel('頻度', fontsize=12)
plt.show()

[Execution Result]
There is a peak between 4 and 5 times.

Let's check the mean and unbiased variance of the posterior predictive distribution of x*.

# 事後予測分布と累積分布
unique, counts = np.unique(pure_y_samples2, return_counts=True)
ratios = counts / counts.sum()
pd.DataFrame(
    {'x*': ratios, '累積分布': np.cumsum(ratios)}, index=unique
).T.round(2)

[Execution Result]

3. PHC curve
Although not in the text, let's look at the PHC curve for the Poisson distribution model.

# PHC曲線の描画
plt.figure(figsize=(6, 4))
phc_plot(idata2.posterior.lam.data.flatten(), xrange=10)
plt.title(f'「平均値 λ > c 」の確率：PHC曲線')
plt.legend(['経験回数のPHC'])
plt.grid(lw=0.3);

[Execution Result]
The probability drops sharply when the number of experiences is between 4 and 5.

Calculate the probability of more than c=4 experiences and the probability of more than c=6 experiences.

# Cの経験回数の平均が4より大きい（p>4）というPHCの算出
p = idata2.posterior.lam.data.flatten()
c = 4
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

[Execution Result]
The probability that the average number of experiences is greater than 4 is $${97.2\%}$$.

# Cの経験回数の平均が6より大きい（p>6）というPHCの算出
p = idata2.posterior.lam.data.flatten()
c = 6
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

[Execution Result]
The probability that the average number of experiences is greater than 6 is $${3.4\%}$$.

Visually, we can infer that the number of experiences is likely at least 4!
(The posterior mean EAP of the average number of experiences was $${5.008}$$ times.)

C. Construction of the Normal Distribution Model

Next, we will perform modeling using a normal distribution model.
We will reuse the code pattern from the Poisson distribution model and rewrite it for the normal distribution.

1. Model Definition

# モデルの定義 正規分布モデル

# logp（対数尤度）関数の定義
def custom_logp3(value, mu, sigma, m):
    '''
    value: 観測値 (PyTensor変数が渡される)
    mu: 平均
    sigma: 標準偏差
    m: マスク
    '''
    # 「観測値 value から m を引いた値が正規分布(μ, σ)に従う」という設定
    return pm.logp(pm.Normal.dist(mu=mu, sigma=sigma), value - m)

# coordsの定義
coords = {'data': idx_data}

# モデルの定義
with pm.Model(coords=coords) as model3:
    
    # Dataの定義
    y = pm.Data('y', value=data2, dims='data')
  
    # 事前分布の定義
    mu = pm.Uniform('mu', lower=0, upper=100)
    sigma = pm.Uniform('sigma', lower=0, upper=100)
    
    # 混合分布の構成要素の定義
    # マスクのfull集合 [-12, -1, …, -1, 1, …, 12]
    masks = np.delete(np.arange(-12, 13), 12)
    # リスト内包表記で24個(マスク個数分)のカスタム正規分布を生成
    components = [
        pm.CustomDist.dist(
            mu, sigma, m,
            logp=custom_logp3,                   # logp関数を指定
            dist=lambda mu, sigma, m, size: (    # 標本サンプリング関数を指定
                pm.Normal.dist(mu=mu, sigma=sigma, size=size) + m
            )
        ) for m in masks
    ]
    
    # 尤度：混合分布
    likelihood = pm.Mixture('likelihood', w=w, comp_dists=components,
                             observed=y, dims='data')

[Model Annotation]

• Definition of coords
  You can assign names to coordinates and set the values those coordinates can take.
  This time, I set the following one.

  • Data row coordinate: name 'data', value 'row index'

• Definition of data
  I set the response value $${y}$$.

• Prior distribution of parameters

  • The prior distributions for the mean $${\mu}$$ and standard deviation $${\sigma}$$ of the number of experiences related to the content of the key item are uniform distributions in the interval $${[0,100]}$$.

• Likelihood

  • This is a mixture distribution of 24 normal distributions.

  • The mixture distribution was implemented as follows.

    • Specify the mixture distribution 'pm.Mixture' for the likelihood.
      Refer to 'components' for the composition of the mixture distribution.
      The weights (mixing ratios) for each distribution are specified by w.

    • Using the settings below in the components list, generate 24 custom probability distributions 'pm.Custom.dist()' using list comprehension.

      • In the log probability logp, specify the normal distribution and parameters $${\mu,\ \sigma}$$ for the 'difference between the response value $${y}$$ and the mask m'.

      • In the distribution dist, use the Python anonymous function lambda to specify the sample sampling function for the response value $${y}$$ based on 'normal distribution + mask m'.

2. Checking the Model Appearance

# モデルの表示
model3

[Execution Result]
The mixture distribution on the second line is cut off on the right.
It is composed of 24 custom distributions (normal distribution + mask).

# モデルの可視化
pm.model_to_graphviz(model3)

[Execution Result]
Is there a Hidden Mickey!?

3. Sampling from the Posterior Distribution
The processing time was approximately 1 minute.

# 事後分布からのサンプリング ※NUTSサンプラーにnutpieを使用
# テキスト：iter=21000, warmup=1000, chains=5
with model3:
    idata3 = pm.sample(draws=20000, tune=1000, chains=5,
                       nuts_sampler='nutpie', random_seed=1234)

[Execution Result]
There are 0 divergences.

4. Checking the Sampled Data
Check $${\hat{R}}$$, the summary statistics of the posterior distribution, and the trace plot.
Following the text, we will use $${\hat{R} \leq 1.1}$$ to confirm the convergence of the posterior distribution.

# r_hat>1.01の確認
# 設定
idata_in = idata3        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

[Execution Result]
For now, I checked with a threshold of $${\hat{R} > 1.01}$$.
There are 0 parameters with $${\hat{R} > 1.01}$$.
I was able to confirm that all parameters satisfy $${\hat{R} \leq 1.01}$$.

I will roughly check the summary statistics and trace plot of the posterior distribution sampled data.

# 推論データの要約統計情報の表示
var_names = ['mu', 'sigma']
pm.summary(idata3, hdi_prob=0.95, var_names=var_names, round_to=3)

[Execution Result]
The mean of the mean $${\mu}$$ of the number of experiences of C (EAP in the text) is 5.340.

# トレースプロットの表示
pm.plot_trace(idata3, var_names=var_names, compact=False, figsize=(8, 4))
plt.tight_layout();

[Execution Result]
The graph on the right is drawn evenly, and it feels like it has converged.
The five colors in the graph on the left = the five Markov chains overlap, producing almost identical inferred values.

Analysis of C's Normal Distribution Model

Regarding the number of experiences of C, I will proceed with the analysis based on the posterior distribution of the normal distribution model.

1. Posterior Distribution
Display the posterior statistics for the mean $${\mu}$$ and standard deviation $${\sigma}$$.
This corresponds to Table 4.4 in the text.

# 推論データの事後統計量の算出 ※表4.4に相当

# 要約統計量データフレームの作成
stats_df2 = pd.DataFrame(
    {'μ': calc_stats(idata3.posterior.mu.data.flatten()),
     'σ': calc_stats(idata3.posterior.sigma.data.flatten())},
    index = ['EAP', 'post.sd', 'CI2.5%', 'CI97.5%', 'HDI2.5%', 'HDI97.5%']
    ).T

# データフレームの表示
display(stats_df2.round(3))

[Execution Result]
The results are almost the same as the text!

[Analysis]
The mean value (EAP) of the mean parameter $${\mu}$$ is $${5.340}$$, and the posterior standard deviation (post.sd) is $${0.516}$$.

Compare point estimation and interval estimation.

$$
\begin{array}{c|ccc}
Model & Mean & Std. Dev./Error & 95\% \text{CI Lower} & 95\% \text{CI Upper} \\
\hline
Normal Distribution & 5.340 & 0.516 & 4.353 & 6.376 \\
Poisson Distribution & 5.008 & 0.537 & 3.972 & 6.079 \\
Conventional Method & 2.272 & 0.551 & - & - \\
\end{array}
$$

• Regarding the mean value

  • The normal distribution model and the Poisson distribution model are similar.

  • On the other hand, the estimated values of the two models are far from those of the conventional method.

  • The mean value of the conventional method, $${2.272}$$, is not included in the 95% CI interval of the normal distribution model and the Poisson distribution model.

• Regarding standard deviation, error, and 95% credible interval

  • The standard deviations of the normal distribution model and the Poisson distribution model are smaller than the standard error of the conventional method, $${0.551}$$.

  • The standard deviation and the width of the 95% credible interval of the normal distribution model are smaller and narrower than those of the Poisson distribution model.

Visualize the posterior distribution of the normal distribution model.

# パラメータの事後分布の可視化

# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharey=True)
# μのヒストグラムの描画
sns.histplot(idata3.posterior.mu.data.flatten(), bins=30, kde=True,
             stat='density', color='steelblue', edgecolor='white', alpha=0.5,
             ax=ax1)
ax1.set(title='μの事後分布', xlabel='μ', ylabel='density')
# σのヒストグラムの描画
sns.histplot(idata3.posterior.sigma.data.flatten(), bins=30, kde=True,
             stat='density', color='tomato', edgecolor='white', alpha=0.5,
             ax=ax2)
ax2.set(title='σの事後分布', xlabel='σ')
# 最終化
plt.tight_layout()
plt.show()

[Execution Result]
It is a beautiful bell shape!

2. Posterior Predictive Distribution
Let's use the MCMC samples to generate and plot samples of the posterior predictive distribution for the number of experiences $${x^*}$$.

# 事後予測分布の可視化

# 推定された mu, sigma のサンプリング結果（全サンプル）を取得
post_mu = idata3.posterior['mu'].values.flatten()
post_sigma = idata3.posterior['sigma'].values.flatten()

# 推定された mu, sigma を使って、マスクなしの純粋な正規乱数を生成
rng = np.random.default_rng(seed=1234)  # 乱数生成器の初期化
pure_y_samples3 = rng.normal(loc=post_mu, scale=post_sigma)

# 事後予測分布のヒストグラムの描画
plt.hist(pure_y_samples3, bins=30, align='left', color='tab:blue',
         edgecolor='white')
plt.title('Cの経験人数の事後予測分布')
plt.xlabel('経験人数の事後予測値', fontsize=12)
plt.ylabel('頻度', fontsize=12)
plt.show()

[Execution Result]
There is a peak between 4 and 6.
Negative values are occurring, which is a bit disappointing...

Check the mean and unbiased variance of the posterior predictive of $${x^*}$$.

# x*の平均と不偏分散の算出
print(f'x*の平均　　: {np.mean(pure_y_samples3):.3f}')
print(f'x*の不偏分散: {np.var(pure_y_samples3, ddof=1):.3f}')

[Execution Result]

3. PHC Curve
Next, we will plot the PHC curve.
This corresponds to Figure 4.7 in the text, and it is a curve representing the probability that the research hypothesis that "the mean number of experiences of C is greater than $${c}$$" is correct, assuming a normal distribution model.

# PHC曲線の描画 ※図4.7に相当
plt.figure(figsize=(6, 4))
phc_plot(idata3.posterior.mu.data.flatten())
plt.title(f'「平均値 μ > c 」の確率：PHC曲線')
plt.legend(['経験回数のPHC'])
plt.grid(lw=0.3);

[Execution Result]
We obtained almost the same results as the text.
It can be seen that the probability fluctuates significantly (drops sharply) between 4 and 6 times.

Let's calculate the PHC that the mean number of experiences is greater than 4 ($${p>4}$$) and the PHC that it is greater than 6 ($${p>6}$$).

# Cの経験回数の平均が4より大きい（p>4）というPHCの算出
p = idata3.posterior.mu.data.flatten()
c = 4
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

[Execution Result]
The probability of exceeding 4 times is almost 100%.
The population mean of C's experiences is likely at least 4 times.

# Cの経験回数の平均が6より大きい（p>6）というPHCの算出
p = idata3.posterior.mu.data.flatten()
c = 6
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

[Execution Result]
The probability of exceeding 6 times has decreased, but it seems to be $${10\%}$$.

Saving the idata file is omitted.

This concludes Chapter 4.

Conclusion

The origins of the three indirect questioning method articles

We worked on three indirect questioning methods in Chapters 2 through 4.
Personally, the randomized response technique in Chapter 3 wassimpleand easy to understand.
However, it seems that the randomized response technique can only be applied when the response values are binary.
Wondering if there were any othersimple methods(other than the item count technique and the AR method) that could handle count data or continuous values, I searched the web and...
I discovered the following site!

https://www.waseda.jp/sem-toyoda-lab/backbook/gform.html

The faces of the authors of Chapters 2 through 4 are also there.
It seems they were not using Bayesian methods at the time.
For your reference.

🍀🍀🍀

Library Version Upgrade

This revised article was written on March 9, 2026.
Actually, at the beginning of March, ArviZ reached Ver. 1.0.0, which was a major update (or rather, the release of the official version).
What concerns me is that there have been significant changes.

• Split into three libraries

  • arviz-base

  • arviz-stats

  • arviz-plots

• The InferenceData format has been phased out and replaced with the xarray.DataTree format

With this change, it seems PyMC will also be undergoing updates to accommodate it...

I wonder what I should do about the PyMC & ArviZ code I've written (and uploaded to my blog) so far...
It's a year-end full of worries.

Series Articles

Next Article

Previous Article

Table of Contents

Blog Introduction

I am writing seven series of articles on note.
Please come and take a look!

1. Relaxed Statistics

This is a blog that roughly explores probability and statistics using the Statistical Test Grade 2 problem collection as a guide.
Feel free to read it like casual conversation. Please come and take a look.
It corresponds to the Statistical Test Grade 2 Official Problem Collection CBT version.
There is also a distribution of sample code for Python and EXCEL.

2. Experiment! Fun Bayesian Modeling 1 & 2 with PyMC Ver. 5

I will draw and analyze the Bayesian models used in the psychology research of the books 'Fun Bayesian Modeling' and 'Fun Bayesian Modeling 2' using PyMC Ver. 5.
Like this book, many Bayesian models are written in R + Stan.
I will strive to explore the possibilities of PyMC and make it easy to practice Bayesian modeling.
Since these are familiar and easy-to-visualize themes, please try running them with PyMC and let's enjoy it together!

3. Experiment! Bayesian Modeling from Iwanami Data Science 1 with PyMC Ver. 5

I will draw and analyze the Bayesian models by four Bayesians from the book 'Experiment! Iwanami Data Science Vol. 1' using PyMC Ver. 5.
This book is a great resource for learning the basics of Bayesian programming.
I feel like I've become friends with Bayesian methods by happily running PyMC models.
Everyone, please try running them with PyMC and let's play and learn together!

4. Fun Copying: Bayesian, Python, etc.

I am blogging about the results of my "book transcription activities" involving Bayesian methods, Python, and other topics.
I am primarily working on translating them into Python.
I hope this serves as sample code for fellow learners who are also transcribing books.🍀

5. Introduction to Time Series Analysis for Psychology with R and Stan, implemented in Python and PyMC Ver. 5

I am practicing the time series analysis from the book "Introduction to Time Series Analysis for Psychology with R and Stan" using Python and PyMC Ver. 5.
This book is packed with themes on time series analysis!
I have truly realized the depth of time series analysis.
I will enjoy learning time series analysis using my favorite language, Python.

6. Writing about data science-related topics

I write columns on statistics, data analysis, AI, machine learning, and Python on an irregular basis.
Many of my articles are related to statistics and data science books.
I have created series on "Statistics," "Python," "Mathematics and Python," and "R."

7. Notes on practicing Python machine learning programming

I have written articles about my various thoughts while studying the book "Python Machine Learning Programming: PyTorch & scikit-learn Edition."
This book is a textbook for scikit-learn and PyTorch.
Please feel free to give it a try if you are interested.

Thank you very much for reading until the end.