Transcribing 'Understanding Statistics through Numerical Simulation' in Python - Chapter 3, Section 3.2 'Expected Value and Variance of Random Variables'

Chapter 3: 'Basics of Random Number Generation Simulation'

Authors: Dr. Koji Kosugi, Dr. Yasunori Kinosada, Dr. Yuji Shimizu

This article covers the 'Understanding Statistics through Numerical Simulation' Chapter 3, 'Basics of Random Number Generation Simulation', Section 3.2, 'Expected Value and Variance of Random Variables', Python transcription activity.

Continuing from last time, we will work on random number generation, which is the core of simulation.
This time, to deepen our familiarity (understanding) with probability distributions, we will verify the expected values and variances of random variables following several probability distributions using simulations and other methods.

R's basic programming syntax is roughly similar to Python's syntax.
Let's trace the details of the code characters and get closer to both R and Python!
Now, let's open the text and set off on a warm-up journey🚀

Introduction

Introduction to the text 'Understanding Statistics through Numerical Simulation'

This series is a Python transcription of the book 'Understanding Statistics through Numerical Simulation: Learning Psychological Statistics by Trying with R' (Gijutsu-Hyohron Co., Ltd., referred to as 'the text').

The text was published in September 2023, and as the subtitle 'Learning Psychological Statistics by Trying with R' suggests, it is a wonderful book that practices numerical simulations useful for understanding psychological statistics using concrete code in R, which has a reputation for statistical processing.

Citation Notation

This article cites text and code published in the book listed in the source, and the published text and code have been modified as appropriate.
[Source]
'Understanding Statistics through Numerical Simulation: Learning Psychological Statistics by Trying with R', First Edition, First Printing, Authors: Koji Kosugi, Yasunori Kinosada, Yuji Shimizu, Gijutsu-Hyohron Co., Ltd.

The illustrations in this article are borrowed from 'Cute Free Material Collection Irasutoya'.
Thank you!

3.2 Expected Value and Variance of Random Variables

As per the text, we will work on the expected value and variance of random variables, and the reproducibility of the normal distribution.

This article writes Python code in Jupyter Notebook format (extension .ipynb).
In general, I use scipy.stats to calculate the characteristic values of probability distributions, and numpy.random.generator to generate random numbers.

I will import the libraries mainly used.

### インポート

# 数値計算
import numpy as np

# 確率・統計
import scipy.stats as stats

# 数値積分
from scipy import integrate

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

3.2.1 Expected Value

■ Expected value of a loaded die
We will calculate the expected value of the loaded die that appeared in section 3.1.1.
The probability of rolling 1 to 4 is 0.1 each, and the probability of rolling 5 and 6 is 0.3 each.

### 83ページ 3.1.1項のイカサマサイコロの期待値
(0.1 * 1) + (0.1 * 2) + (0.1 * 3) + (0.1 * 4) + (0.3 * 5) + (0.3 * 6)

[Execution Result]
It is 4.3.

■ Expected value of a fair die
The expected value of a die where each outcome has a 1/6 probability is the expected value of a random variable following a discrete uniform distribution in the range of 1 to 6.

### 84ページ 公正なサイコロの期待値 a=1, b=6の離散一様分布に従う

# 設定
low = 1    # 下限
high = 6   # 上限
x = np.arange(1, 7) # 出目

# 期待値の計算
sum(stats.randint.pmf(k=x, low=low, high=high + 1) * x)

[Execution Result]
It is 3.5.

■ Expected value of a random variable following a continuous uniform distribution
The expected value of a random variable following a continuous probability distribution seems to require integral calculation.
This is the code for calculating the expected value using numerical integration for a continuous uniform distribution with a lower bound of α=1 and an upper bound of β=6.

### 84ページ 連続一様分布に従う確率変数の期待値 数値積分版
def d_unif_exp(x, alpha=1, beta=6):
    return stats.uniform.pdf(x=x, loc=alpha, scale=beta - alpha) * x

integrate.quad(func=d_unif_exp, a=1, b=6)  # 数値積分の実行(積分値, 誤差推定値)

[Execution Result]
The left side is the expected value. The expected value is 3.5.

Calculate the expected value of a random variable following a continuous uniform distribution using scipy.stats functionality.
Execute mean on the continuous uniform distribution 'uniform'.

### scipy.statsで期待値を算出
# パラメータ設定
alpha = 1
beta = 6
# 期待値算出
stats.uniform.mean(loc=alpha, scale=beta - alpha)

[Execution Result]

■ Expected value of a random variable following a standard normal distribution
Calculate the expected value by following the implementation in the text.

### 84ページ 標準正規分布に従う確率変数の期待値=平均パラメータ
def std_norm_exp(x, mu=0, sigma=1):
    return stats.norm.pdf(x=x, loc=mu, scale=sigma) * x
integrate.quad(func=std_norm_exp, a=-np.inf, b=np.inf)

[Execution Result]
The expected value is 0.
Since the standard normal distribution has a mean parameter of 0, the expected value is 0.

Simulation time!
Simulate the expected value by generating a large number of standard normal random variables and taking their average.

### 85ページ 乱数生成シミュレーションで近似的に期待値を算出
rng = np.random.default_rng(seed=123)
rng.normal(size=100000, loc=0, scale=1).mean()

[Execution Result]
It is almost 0.

Calculate the expected value of a random variable following a standard normal distribution using scipy.stats functionality.
Execute mean on the normal distribution 'norm'.

### scipy.statsで期待値を算出
# パラメータ設定
mu = 0
sigma = 1
# 期待値算出
stats.norm.mean(loc=mu, scale=sigma)

[Execution Result]
It is 0 as expected.

3.2.2 Variance

Next, calculating the variance.

■ Variance of a random variable following a discrete uniform distribution
This is the variance of a discrete uniform distribution with α=1 and β=6.
I will write it following the code in the text.

### 85ページ 離散一様分布に従う確率変数の分散
def d_unif_var(x, alpha=1, beta=6):
    # α = 1, β = 6の離散一様分布に従う確率変数の期待値
    expected_val = np.arange(alpha, beta+1).mean()
    var = stats.randint.pmf(k=x, low=alpha, high=beta+1) * (x - expected_val)**2
    return var

d_unif_var(range(1, 7)).sum()

[Execution Result]

Next, as per the text, I will calculate the variance using the formula for the variance of a random variable following a discrete uniform distribution.

### 86ページ 離散一様分布に従う確率変数の分散（公式利用）
alpha = 1
beta = 6
((beta - alpha + 1)**2 - 1) / 12

[Execution Result]

Of course, you can also find the variance using scipy.stats.

### scipy.statsで分散を算出
# パラメータ設定
alpha = 1
beta = 6
# 分散算出
stats.randint.var(low=alpha, high=beta + 1)

[Execution Result]

■ Variance of a random variable following a standard normal distribution
I will draw it following the code in the text.

### 86ページ 標準正規分布に従う確率変数の分散
def d_norm_var(x, mu=0, sigma=1):
    expected_val = mu   # 標準正規分布に従う確率変数の期待値
    return stats.norm.pdf(x=x, loc=mu, scale=sigma) * (x - expected_val)**2

integrate.quad(func=d_norm_var, a=-np.inf, b=np.inf)

【Execution Result】
The left side is the variance. The variance is 1.
Since the variance parameter of the standard normal distribution is 1, the variance is 1.

Simulation start!
I will perform a simulation to calculate the variance by generating a large number of standard normal random variables.

### 86ページ 乱数生成シミュレーションで近似的に期待値を算出

# 2章で定義した自作関数
def var_p(x):  # xはnumpy配列
    n = len(x)
    mean_x = x.mean()
    var_x = sum((x - mean_x)**2) / n
    return var_x

rng = np.random.default_rng(seed=0)
var_p(rng.normal(size=100000, loc=0, scale=1))

【Execution Result】
It is almost 1.

Of course, you can also find the variance using scipy.stats.

### scipy.statsで分散を算出
# パラメータ設定
mu = 0
sigma = 1
# 期待値算出
stats.norm.var(loc=mu, scale=sigma)

【Execution Result】
It is 1 as expected.

Column: Reproducibility of the Normal Distribution

Following the text, I will perform a simulation of the reproducibility of the normal distribution.

I will create $${Z}$$ by adding normal distribution random variables $${X}$$ and $${Y}$$, and draw a histogram of $${Z}$$.
This histogram is based on the realized values from the simulation.
I will overlay the 'probability density function of a normal distribution with mean $${5}$$ and variance $${125}$$' as the theoretical value on this histogram to check if they match.

### 87ページ 正規分布の再生性のシミュレーション 88ページ図3.15

## パラメータの設定
mu_x = 0
sigma_x = 10
mu_y = 5
sigma_y = 5
n = 20000  # 生成する乱数の個数

## 確率変数x,y,zの算出
rng = np.random.default_rng(seed=0)
# 確率変数x μ = 0, σ = 10の正規分布に従う乱数
x = rng.normal(size=n, loc=mu_x, scale=sigma_x)
# 確率変数y μ = 5, σ = 5の正規分布に従う乱数
y = rng.normal(size=n, loc=mu_y, scale=sigma_y)
# 確率変数z
z = x + y

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# zのヒストグラムの描画
ax.hist(z, bins=30, density=True, color='tab:blue', ec='white', alpha=0.7)
# zが従う正規分布の確率密度関数（理論値）の描画（赤い線）
line_x = np.linspace(-70, 70, 1001)
ax.plot(line_x, stats.norm.pdf(line_x, loc=mu_x+mu_y,
                               scale=np.sqrt(sigma_x**2 + sigma_y**2)),
         color='tab:red', lw=2)
# 修飾
ax.set(xlim=(-70, 70), xlabel='z', ylabel='Density', title='Histogram of z')
ax.grid(lw=0.5);

【Execution Result】
The histogram of variable $${Z}$$ and the probability density function (red line) of the normal distribution with mean $${5}$$ and variance $${125}$$ matched almost perfectly!

I will calculate the mean and variance of variable $${Z}$$.

### 88ページ zの平均（期待値）
z.mean()

【Execution Result】
It approximates the theoretical mean value of 5.

### 88ページ zの分散
var_p(z)  # 分散（2章で定義した自作関数を使用）

【Execution Result】
It approximates the theoretical variance value.

By the way, you can also calculate the variance with numpy.

### numpyで標本分散を算出
z.var()  # ddofのデフォルト値は0

【Execution Result】

■ Normal distribution followed by the sample mean of random variables independently following the same normal distribution
This is an important property of the sample mean.
I will borrow the theorem from the text.

I will verify this theorem with a simulation.
I will create 10,000 sample means $${\bar{X}}$$ of random variables $${X_1, X_2, X_3, X_4}$$ that independently follow a normal distribution with mean $${\mu=50}$$ and variance $${\sigma^2=10^2}$$ by generating random numbers, and draw a histogram of $${\bar{X}}$$.
I will overlay the probability density function of a normal distribution with mean $${50}$$ and variance $${10^2/4=25}$$ (standard deviation $${5}$$) as the theoretical value to confirm its match with the histogram.

### 89ページ 確率変数の平均のヒストグラムの描画 90ページ図3.16

n = 4
mu = 50
sigma = 10
iter = 10000

means = np.zeros(iter)
rng = np.random.default_rng(seed=123)
for i in range(iter):
    means[i] = rng.normal(size=n, loc=mu, scale=sigma).mean()

## 平均のヒストグラムの描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# 平均のヒストグラムの描画
ax.hist(means, bins=30, density=True, color='tab:blue', ec='white', alpha=0.7)
# 平均が従う正規分布の確率密度関数（理論値）の描画（赤い線）
line_x = np.linspace(means.min(), means.max(), 1001)
ax.plot(line_x, stats.norm.pdf(line_x, loc=mu, scale=sigma / np.sqrt(n)),
        color='tab:red', lw=2)
# 修飾
ax.set(xlabel='mean', ylabel='Density', title='Histogram of mean')
ax.grid(lw=0.5);

[Execution Result]
They matched almost perfectly!

That concludes this coding session.

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 explores probability and statistics roughly, using the Statistical Test Grade 2 problem collection as a guide.
Feel free to read it as if it were casual conversation. Please do take a look.
It corresponds to the Statistical Test Grade 2 Official Problem Collection CBT version.
Sample code for Python and EXCEL is also available.

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

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

3. Experiment! Iwanami Data Science 1 Bayesian Modeling 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 Coding: Bayesian, Python, etc.

I will blog about the results of my 'book coding activities' for Bayesian, Python, and others.
I am mainly working on translations into Python.
I hope this serves as sample code for fellow coders 🍀

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

I will practice 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 time series analysis themes!
I realized the depth of time series analysis.
I will enjoy learning time series analysis with my favorite language, Python.

6. Writing about Data Science-like things

I write columns on statistics, data analysis, AI, machine learning, and Python on an irregular basis.
There are many articles related to statistics and data science books.
Series on 'Statistics', 'Python', 'Math and Python', and 'R' have been created.

7. Python Machine Learning Programming Practice Notes

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.

Thank you very much for reading until the end.