Copying "Introduction to Statistical Analysis" in Python Vol. 7 - Chapter 3 "First Statistical Estimation" Part 3: Sample Size, Maximum Likelihood Estimation

Chapter 3 "First Statistical Estimation"

Book author: Dr. Sadao Ishimura

This article covers the "Introduction to Statistical Analysis" Chapter 3 "First Statistical Estimation" Python copying activity.

This is a copying series that calmly converts the book's figures, tables, and calculations into Python.
In this article, among the statistical estimation themes in Chapter 3, we will tackle sample size calculation and parameter estimation by maximum likelihood method.
We will continue to utilize ChatGPT!

Now, let's open the book and set off on a journey of statistical analysis 🚀

Introduction

This blog series introduces the "joy of statistical analysis" gained through copying the Python code from the book "Introduction to Statistical Analysis" (Tokyo Tosho, referred to as the "Text").

The introduction of the book and citation notation are posted in the linked article.

Chapter 3 First Statistical Estimation

This article covers the following sections of Chapter 3.

The data used in the article cites the data published in the text itself.
For items with a small number of data points, the data is registered in the code, and for items with a large number of data points, the data is read by creating a CSV file.

Import the libraries used in Chapter 3.

### インポート

# 数値計算
import math                      # python標準ライブラリ
import numpy as np
import pandas as pd
from scipy.special import gamma  # ガンマ関数

# 統計
import scipy.stats as stats

# ユーティリティ
import collections               # カウンター

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

Introduction

ChatGPT gave a rough overview of this theme.

✏️ Sample size calculation

"How many people should I ask to get 'reliable results' for a survey?"
"With how many monitors should I conduct a product test?"

Actually, when deciding on such a "number of people to investigate," sample size calculation plays an active role.
To reach a statistically reliable conclusion, it is very important to think about "how much data is needed?"

✏️ Parameter estimation by maximum likelihood method

'What is the most likely reason this data was created?'

Maximum Likelihood Estimation is a method to answer such a question.
It finds the 'probability mechanism' that best fits the observed data, and searches for which numbers (parameters) are most plausible within that framework
.

'Maximum Likelihood Estimation' is a difficult concept, but let's do our best to tackle it!
Let's get started!

Section 3.6 Deciding the Sample Size?!

Sample size refers to the 'number of data points N'.
The motivation for wanting to decide the sample size is to calculate the 'required number of data points' in advance before conducting surveys such as questionnaires.
We will consider a sample size that 'keeps data collection costs down' while also 'enabling statistically reliable analysis'.

The text begins with a formula for determining sample size by focusing on the 'error of the confidence interval'.

◆ ◆ ◆

■ Formula for sample size N in interval estimation of population proportion p.128
Think of surveys such as usage rates, viewership ratings, approval ratings, or the percentage of 'yes' answers in a questionnaire.
'Interval estimation of population proportion' seems useful here!

I will borrow the formula for interval estimation of population proportion from the text.
The sample size N required to keep the error of the 100(1-α)% confidence interval for the population proportion within E can be calculated as follows.
For a 95% confidence interval, α=0.05.

p is the assumed value of the population proportion, such as usage rate or approval rate.
z(α/2) is the 100·(α/2) percentile of the standard normal distribution.
For a 95% confidence interval, α=0.05, so z(α/2) is the 2.5% point.

If the population proportion p cannot be predicted, use the following formula.

To be safe, we assume p=0.5, where p(1-p) is maximized, using 0.5 * (1-0.5) = 1/2 * 1/2 = 1/4.

The calculated sample size N is determined so that the width of the confidence interval is

.

By the way, this formula assumes that "all data is collected."
In the case of questionnaire surveys, there may be cases where people do not respond to the distributed questionnaires.
It is also possible to adjust the sample size by the questionnaire response rate, etc., to determine the "number of questionnaires to distribute."

◆

We will define a function to calculate the sample size according to the formula.
We will use scipy.stats to calculate the % point of the standard normal distribution.

### 母比率の区間推定のサンプルサイズNの算出関数 p.128～130
# E: 誤差（比率）、alpha: 1-信頼係数
# p: 想定母比率（想定できないときはサンプルサイズ最大になるp=0.5を使用）

def calc_sample_size_pop_ratio(E, alpha, p=0.5):
    # 標準正規分布の上側100*(α/2)%点の取得 scipy.stats利用
    z = stats.norm.ppf(q=1 - alpha/2, loc=0, scale=1)
    return (z / E)**2 * p * (1 - p)

Let's calculate the sample size for the population proportion using test data with an error of $${E=0.02}$$ (error ratio 0.02) and a confidence coefficient of $${100(1-\alpha)\%}$$ where $${\alpha=0.05}$$.

# テスト
calc_sample_size_pop_ratio(E=0.02, alpha=0.05)

[Execution Result]
Rounding up the fractional part of the calculation result, the minimum sample size is $${2401}$$.

📢 By the way, for your information 📢
There is a "Calculation form for the required sample size in interval estimation of population proportion" in an article on TokeiWEB.
Let's give it a try!

I will quote the screen showing the result calculated using the test data as input.

◆ ◆ ◆

■ Example of sample size for population proportion p.130
We will calculate the sample size for the population proportion based on the assumed population proportion $${p=0.34}$$, confidence interval 95% ($${\alpha=0.05}$$), and error $${E=0.02}$$ in the example on p.130.

First, let's quickly calculate the sample size using the sample size calculation function from earlier.

### 関数利用 p.130
calc_sample_size_pop_ratio(p=0.34, E=0.02, alpha=0.05)

[Execution Result]
A sample size of $${2156}$$ or more is required.

Next, let's calculate the sample size step by step following the calculation procedure on p.130 of the text.

(1) Calculation of $${\alpha/2}$$

### テキストの計算を辿る p.130

# 誤差Eと想定母比率p
E, p = 0.02, 0.34

# α
alpha = 0.05
alpha / 2

[Execution Result]

(2) Calculation of the upper $${100 \cdot \alpha/2 \%}$$ point of the standard normal distribution

# 標準正規分布の95%点
z = stats.norm.ppf(q=1 - alpha/2, loc=0, scale=1)
z

[Execution Result]

(3) Calculation of sample size

# サンプルサイズの算出
(z / E)**2 * p * (1 - p)

[Execution Result]
A sample size of $${2156}$$ or more is required.

■ Formula for sample size $${N}$$ in interval estimation of population mean p.130
Next, let's move on to the topic of "population mean"!
I will borrow the formula for interval estimation of the population mean from the text.

The sample size $${N}$$ required to keep the error of the $${100(1-\alpha)\%}$$ confidence interval for the population mean within $${E}$$ can be calculated as follows.
Information regarding the assumed population variance $${\sigma^2}$$ is required.

$${z (\tfrac{\alpha}{2})}$$ is the $${100\cdot\tfrac {\alpha}{2}}$$ percentile of the standard normal distribution.

The calculated sample size $${N}$$ is determined so that the width of the confidence interval becomes

.

We will define a function to calculate the sample size according to the formula.
We use scipy.stats to calculate the percentiles of the standard normal distribution.

### 母平均の区間推定のサンプルサイズNの算出関数 p.130
# E: 誤差（個、cmなどのサンプルの単位）、alpha: 1-信頼係数, sigma2: 想定母分散 

def calc_sample_size_pop_mean(E, alpha, sigma2):
    # 標準正規分布の上側100*(α/2)%点の取得 scipy.stats利用
    z = stats.norm.isf(q=alpha/2, loc=0, scale=1)
    return (z / E * np.sqrt(sigma2))**2

Let's calculate the sample size for the population mean using test data with an error of $${E=1}$$ (e.g., 1mm error), a confidence coefficient of $${100(1-\alpha)\%}$$ where $${\alpha=0.05}$$, and a population variance of $${\sigma^2=3}$$.

# テスト
calc_sample_size_pop_mean(E=1, alpha=0.05, sigma2=3)

[Execution Result]
Rounding up the fractional part of the calculation result, the minimum sample size is $${12}$$$$.

In determining the sample size, we focused on $${N}$$.
In the next topic, "Maximum Likelihood Estimation," the focus will shift to a different subject.

Section 3.7 What is Maximum Likelihood Estimation?

■ Parameter Estimation and Maximum Likelihood Estimation
When assuming that data follows a certain probability distribution, "trying to estimate the characteristics (i.e., parameters) of that distribution from the data" is what we call parameter estimation.
Parameters are, for example...

• In the case of a binomial distribution $${\text{Bin}(n, p)}$$

  • The parameter to be estimated is the population proportion $${p}$$

• In the case of a normal distribution $${\mathcal{N}(\mu, \sigma^2)}$$

  • The parameters to be estimated are the population mean $${\mu}$$ and the population variance $${\sigma^2}$$

Parameters are important values that determine the "center" or "dispersion" of a probability distribution.

Maximum likelihood estimation is one of the representative methods for estimating parameters.
By the way, "point estimation" and "interval estimation," which we covered up to the previous article, are also methods for estimating parameters, which are characteristics of the population.

◆ ◆ ◆

■ Definition of Maximum Likelihood Estimation p.137
I will borrow the definition from the text.

A likelihood function is a function that calculates the "likelihood" (=likelihood) of each case by tentatively applying various parameters, assuming that "this data might have been observed because the parameters were this way."the plausibility (=likelihood) of eachis calculated.

Specifically, you tentatively decide the value of the parameter $${\theta}$$, and multiply (take the product of) theprobability (or probability density)of each piece of data $${\{x_1, x_2, \ldots, x_N\}}$$ occurring:

The fact that the result of this multiplication—the likelihood—is large means that "it is highly probable that this data would appear under that parameter." The goal of the maximum likelihood method is to apply various parameter values, calculate the likelihood, and
choose the parameter value that yields the highest likelihood.

◆ ◆ ◆

【Aside】
The "introductory part of the maximum likelihood method" (excluding the parts quoted from the text) was woven together with ChatGPT.
I want to use ChatGPT's ability to generate "intuitive" text to help improve the quality of my article drafts!
Now, back to the main topic.

◆ ◆ ◆

■ When the distribution type of the population is known to be a binomial distribution p.133
I will borrow the formula for the likelihood function of a binomial distribution from the text.

I will define the likelihood function according to the formula using [Python standard library only].

### 二項分布の尤度関数の定義

# 二項分布の尤度関数
def binom_likelihood_func(x_list, n, p):
    # 二項係数部分の掛け算
    binom_coef = math.prod(
        [math.factorial(n) / (math.factorial(x) * math.factorial(n - x))
         for x in x_list])
    # pの指数expo1と1-pの指数expo2のｎ算出
    expo1 = sum(x_list)
    expo2 = sum([n - x for x in x_list])
    # 尤度の計算
    llf = binom_coef * p**expo1 * (1 - p)**expo2
    # 戻り値：尤度
    return llf

I will check the likelihood function for the case where the parameters are $${n=8, p=0.32}$$ using test data.

# テスト
x_list = [4, 2, 3]
binom_likelihood_func(x_list, n=8, p=0.32)

【Execution Result】
The likelihood function is $${0.012}$$.

I will verify the calculation using scipy.stats.
Using stats.binom.pmf(), I calculate the binomial distribution probability of the test data for the specified parameters, and multiply the probabilities using np.prod().

# scipy.statsで答え合わせ
np.prod(stats.binom.pmf(x_list, n=8, p=0.32))

【Execution Result】
It was correct!

■ Example: When the population is a binomial distribution $${Bin(8,p)}$$: Maximum likelihood method for population proportion $${p}$$ p.134
Using the sample $${{4, 2, 3}}$$ and the trial count parameter $${n=8}$$ from the text, I will estimate the population proportion parameter $${p}$$ using the maximum likelihood method.
I cannot calculate the likelihood function without setting a tentative value for $${p}$$.
Here, I have set the tentative values for $${p}$$—the search range—from $${0.32}$$ to $${0.42}$$ in increments of $${0.01}$$.

### 関数利用 p.134

## 設定
x_list = [4, 2, 3]  # 標本
n = 8               # 試行回数

## 初期値
# 探索するパラメータpの値
p_valus = np.arange(0.32, 0.43, 0.01)
# 結果を格納するデータフレーム
result_df = pd.DataFrame(columns=['パラメータp', '尤度関数L'])

## パラメータpの尤度関数の算出
for i, p in enumerate(p_valus):
    L = binom_likelihood_func(x_list, n=n, p=p)
    result_df.loc[i, :] = p, L

## 結果の表示:尤度最大の行をハイライト
(result_df
 .style
 .set_properties(**{'background-color': 'yellow'}, 
                 subset=pd.IndexSlice[result_df['尤度関数L'].argmax(), :])
 .format({'パラメータp': '{:.2f}', '尤度関数L': '{:.8f}'}))

【Execution Result】
These are the search values for parameter $${p}$$ and the corresponding values of the likelihood function $${L}$$.
The value $${p=0.38}$$, where the likelihood function value is the largest, is the estimate of the maximum likelihood method (maximum likelihood estimator).

◆ ◆ ◆

■ When the distribution type of the population is known to be a normal distribution p.135
I will borrow the formula for the likelihood function of a normal distribution from the text.

I will define the likelihood function according to the formula using [Python standard library only].

### 正規分布の尤度関数の定義

# 正規分布の尤度関数
def norm_likelihood_func(list, mu, sigma):
    # 標本サイズの算出
    N = len(list)
    # eの指数部の計算
    expo = -1/2 * sum([((x - mu)/sigma)**2 for x in list])
    # 尤度の計算
    llf = (1 / (math.sqrt(2 * math.pi) * sigma))**N * math.e**expo
    # 戻り値：尤度
    return llf

I will check the likelihood function for test data with parameters $${\mu=4, \sigma=3}$$.
$${\sigma}$$ is the standard deviation.

# テスト
x_list = [5, 3, 4]
norm_likelihood_func(x_list, mu=4, sigma=3)

[Execution Result]
The likelihood function is $${0.0021}$$.

I will verify the calculation using scipy.stats.
I am using stats.norm.pdf() to calculate the probability density of the normal distribution for the test data with the specified parameters, and multiplying the probabilities using np.prod().

# scipy.statsで答え合わせ
np.prod(stats.norm.pdf(x_list, loc=4, scale=3))

[Execution Result]
It matched!

◆ ◆ ◆

■ Example: Maximum likelihood method for the population mean $${\mu}$$ when the population is a normal distribution $${N(\mu,3^2)}$$ p.136
Using the sample $${{5, 3, 4}}$$ and variance parameter $${\sigma^2=3^2}$$ from the text, I will estimate the population mean parameter $${\mu}$$ using the maximum likelihood method.
I cannot calculate the likelihood function without setting a tentative value for $${\mu}$$.
Here, I set the tentative values for $${\mu}$$ (the search range) from $${3.40}$$ to $${4.60}$$ in increments of $${0.02}$$.

### 関数利用 p.136

## 設定
x_list = [5, 3, 4]  # 標本
sigma = 3           # 標準偏差

## 初期値
# 探索するパラメータμの値
mu_valus = np.arange(3.4, 4.7, 0.2)
# 結果を格納するデータフレーム
result_df2 = pd.DataFrame(columns=['パラメータμ', '尤度関数L'])

## パラメータμの尤度関数の算出
for i, mu in enumerate(mu_valus):
    L = norm_likelihood_func(x_list, mu=mu, sigma=sigma)
    result_df2.loc[i, :] = mu, L

## 結果の表示:尤度最大の行をハイライト
(result_df2
 .style
 .set_properties(**{'background-color': 'yellow'}, 
                 subset=pd.IndexSlice[result_df2['尤度関数L'].argmax(), :])
 .format({'パラメータμ': '{:.3f}', '尤度関数L': '{:.8f}'}))

[Execution Result]
These are the search values for parameter $${\mu}$$ and the corresponding likelihood function $${L}$$ values.
The value $${\mu=4.000}$$, where the likelihood function value is the largest, is the estimate (maximum likelihood estimator) from the maximum likelihood method.

■ Calculating likelihood maximization using scipy's optimization function
In the maximum likelihood method used so far, I set a limited "parameter search range" to calculate the likelihood.
Setting a search range is a bit troublesome...

So, let's use scipy's optimization function, minimize().
It will estimate the parameters that maximize the likelihood without needing to define a search range.
In practice, it performs estimation by minimizing the negative likelihood.

I referred to this web article.
Thank you!

I assume a population that follows a normal distribution.
I will import minimize() from scipy.

# 追加インポート
from scipy.optimize import minimize

◆ ◆ ◆

First, I will estimate the population mean $${\mu}$$ using the maximum likelihood method with a known population variance.
Sample $${{5, 3, 4}}$$ and variance parameter $${\sigma^2=3^2}$$ from the text

#### 尤度最大化をscipyで最適化（最小化）してみる 【平均μの推定】
# 参考サイト: https://note.com/united_code/n/nd72d381f7617

# 尤度関数の定義
def likelihood_function_mu(params, list, sigma):
    # パラメータの初期値の取り出し
    mu = params
    # 標本サイズの算出
    N = len(list)
    # eの指数部の計算
    expo = -1/2 * sum([((x - mu)/sigma)**2 for x in list])
    # 尤度の計算
    llf = (1 / (math.sqrt(2 * math.pi) * sigma))**N * math.e**expo
    # 戻り値：尤度のマイナス（最小化するので）
    return -llf

# データの設定
x_list = [5, 3, 4]
sigma = 3           # 標準偏差

# パラメータの初期値の設定
x_bar = sum(x_list) / len(x_list)  # 標本平均
init_params = [x_bar]

# 尤度最大となるパラメータμ,σをscipyで探索
result = minimize(likelihood_function_mu, init_params, args=(x_list, sigma, ))

# パラメータ推定値の表示
print(f'平均μの推定値: {result.x[0]}')
print(f'尤度　　　　 : {-result.fun:.8f}')

[Execution Result]
It matched the estimate of the population mean $${\mu}$$ and the likelihood calculated using the formula in the text.

(Reference) The estimate of the population mean $${\mu}$$ and the likelihood calculated using the formula in the text

◆ ◆ ◆

The second case assumes the population variance is unknown, and we estimate the population mean $$\mu$$ and population variance $$\sigma^2$$ (actually the standard deviation $$\sigma$$ of the population variance) using the maximum likelihood method.

### 尤度最大化をscipyで最適化（最小化）してみる 【平均μと標準偏差σの推定】
# 参考サイト: https://note.com/united_code/n/nd72d381f7617

# 尤度関数の定義
def likelihood_function_mu_sigma(params, list):
    # パラメータの初期値の取り出し
    mu, sigma = params
    # 標本サイズの算出
    N = len(list)
    # eの指数部の計算
    expo = -1/2 * sum([((x - mu)/sigma)**2 for x in list])
    # 尤度の計算
    llf = (1 / (math.sqrt(2 * math.pi) * sigma))**N * math.e**expo
    # 戻り値：尤度のマイナス（最小化するので）
    return -llf

# データの設定
x_list = [5, 3, 4]

# パラメータの初期値の設定
N = len(x_list)                                                   # 標本サイズ
x_bar = sum(x_list) / N                                           # 標本平均
x_std = math.sqrt(sum([(x-x_bar)**2 for x in x_list]) / (N - 1))  # 標本標準偏差
init_params = [x_bar, x_std]

# 尤度最大となるパラメータμ,σをscipyで探索
result = minimize(likelihood_function_mu_sigma, init_params, args=(x_list,))

# パラメータ推定値の表示
mu_hat, sigma_hat = result.x
print(f'平均μの推定値　　: {mu_hat:.4f}')
print(f'標準偏差σの推定値: {sigma_hat:.4f}')
print(f'尤度　　　　　　 : {-result.fun:.8f}')

[Execution Results]
Compared to the $${3}$$ in the case where the population variance is known, the estimated value for the standard deviation $$\sigma$$ is smaller here.
The sample variance and standard deviation of the three data points are also $${1}$$.
The likelihood has increased.

Maximum likelihood estimation was interesting, wasn't it!

Let's have ChatGPT wrap up the end of the article!

📘 A word from ChatGPT:

That is all for this session of copying code.

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 treat it like casual conversation. Please come and take a look.
It corresponds to the Statistical Test Grade 2 Official Problem Collection CBT Edition.
Python and EXCEL sample code are also available for distribution.

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.
Starting with these books, many Bayesian models are written in R language + Stan.
I strive to explore the possibilities of PyMC and make it easy to practice Bayesian modeling.
These are familiar and easy-to-visualize themes, so 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" for Bayes, Python, and others.
I am mainly working on translating them into Python.
I hope this will serve as sample code for fellow learners who are also transcribing.

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

I will practice the time series analysis from the book "Introduction to Time Series Analysis for Psychology using R and Stan" using Python and PyMC Ver. 5.
This book is packed with themes on time series analysis!
I have 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," "Mathematics and Python," and "R" have been created.

7. Python Machine Learning Programming Practice Log

I wrote articles about my various thoughts when I studied 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 like.

Thank you very much for reading until the end.