統計学メモ


```
import numpy as np
from scipy import stats

# 表示桁数の指定（オプション）
np.set_printoptions(precision=3)

# データの準備
fish_data = np.array([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])

# 合計の計算
total_sum = np.sum(fish_data)
print(f"合計: {total_sum}")

# サンプルサイズの計算
sample_size = len(fish_data)
print(f"サンプルサイズ: {sample_size}")

# 平均値の計算
mean_value = np.mean(fish_data)
print(f"平均値: {mean_value}")

# 標本分散の計算
sample_variance = np.var(fish_data, ddof=0)
print(f"標本分散: {sample_variance}")

# 不偏分散の計算
unbiased_variance = np.var(fish_data, ddof=1)
print(f"不偏分散: {unbiased_variance}")

# 標準偏差の計算
std_deviation = np.std(fish_data, ddof=1)
print(f"標準偏差: {std_deviation}")

# データの標準化
standardized_data = (fish_data - mean_value) / std_deviation
print(f"標準化されたデータ: {standardized_data}")
print(f"標準化されたデータの平均値: {np.mean(standardized_data)}")
print(f"標準化されたデータの標準偏差: {np.std(standardized_data, ddof=1)}")

# 最大値、最小値、中央値の計算
max_value = np.max(fish_data)
min_value = np.min(fish_data)
median_value = np.median(fish_data)
print(f"最大値: {max_value}")
print(f"最小値: {min_value}")
print(f"中央値: {median_value}")

# 外れ値を含むデータの平均値と中央値の比較
outlier_data = np.append(fish_data, [100])
mean_with_outlier = np.mean(outlier_data)
median_with_outlier = np.median(outlier_data)
print(f"外れ値を含むデータの平均値: {mean_with_outlier}")
print(f"外れ値を含むデータの中央値: {median_with_outlier}")

# 四分位点の計算
quartiles = np.percentile(fish_data, [25, 50, 75])
print(f"四分位点: {quartiles}")


```



```
import pandas as pd
import scipy as sp

# 表示桁数の指定
%precision 3

# --- グループ別の統計量 ---

# データの作成
fish_multi_data = {
    "species": ["A", "A", "A", "B", "B", "B"],
    "length": [2, 3, 4, 6, 8, 10]
}
fish_multi = pd.DataFrame(fish_multi_data)
print("魚のデータ:\n", fish_multi)

# グループごとの平均値を計算
group = fish_multi.groupby("species")
print("\n平均値:\n", group.mean())

# グループごとの標準偏差を計算
print("\n標準偏差:\n", group.std(ddof=1))

# グループごとの詳細な統計量を計算
print("\n詳細な統計量:\n", group.describe())

# --- 実装：クロス集計表 ---

# データの作成
shoes_data = {
    "store": ["tokyo", "tokyo", "osaka", "osaka"],
    "color": ["blue", "red", "blue", "red"],
    "sales": [10, 15, 13, 9]
}
shoes = pd.DataFrame(shoes_data)
print("\n靴のデータ:\n", shoes)

# クロス集計表を作成
cross = pd.pivot_table(
    data=shoes,
    values="sales",
    aggfunc="sum",
    index="store",
    columns="color"
)
print("\nクロス集計表:\n", cross)

# --- 共分散を計算するためのデータの作成 ---

# データの作成
cov_data = {
    "x": [18.5, 18.7, 19.1, 19.7, 21.5, 21.7, 21.8, 22.0, 23.4, 23.8],
    "y": [34, 39, 41, 38, 45, 41, 52, 44, 44, 49]
}
cov_df = pd.DataFrame(cov_data)
print("\n共分散データ:\n", cov_df)

# データの取り出し
x = cov_df["x"]
y = cov_df["y"]

# サンプルサイズ
N = len(cov_df)

# 平均値の計算
mu_x = sp.mean(x)
mu_y = sp.mean(y)

# 標本共分散
cov_sample = sum((x - mu_x) * (y - mu_y)) / N
print("\n標本共分散:", cov_sample)

# 共分散
cov = sum((x - mu_x) * (y - mu_y)) / (N - 1)
print("共分散:", cov)

# --- 実装：分散共分散行列 ---

# 標本共分散行列
cov_matrix_sample = sp.cov(x, y, ddof=0)
print("\n標本共分散行列:\n", cov_matrix_sample)

# 不偏共分散行列
cov_matrix_unbiased = sp.cov(x, y, ddof=1)
print("\n不偏共分散行列:\n", cov_matrix_unbiased)

# --- 実装：ピアソンの積率相関係数 ---

# 分散の計算
sigma_2_x = sp.var(x, ddof=1)
sigma_2_y = sp.var(y, ddof=1)

# 相関係数
rho = cov / sp.sqrt(sigma_2_x * sigma_2_y)
print("\n相関係数:", rho)

# サンプル相関係数
sigma_2_x_sample = sp.var(x, ddof=0)
sigma_2_y_sample = sp.var(y, ddof=0)
rho_sample = cov_sample / sp.sqrt(sigma_2_x_sample * sigma_2_y_sample)
print("サンプル相関係数:", rho_sample)

# 相関行列
corr_matrix = sp.corrcoef(x, y)
print("\n相関行列:\n", corr_matrix)


```




```

# 数値計算に使うライブラリ
import numpy as np
import pandas as pd

# 表示桁数の指定
%precision 3

# グラフを描画するライブラリ
from matplotlib import pyplot as plt

# グラフをjupyter Notebook内に表示させるための指定
%matplotlib inline

# データの準備
x = np.array([0,1,2,3,4,5,6,7,8,9])
y = np.array([2,3,4,3,5,4,6,7,4,8])

# pyplotによる折れ線グラフ
plt.plot(x, y, color='black')
plt.title("lineplot matplotlib")
plt.xlabel("x")
plt.ylabel("y")

# seabornのセットアップ
import seaborn as sns
sns.set()

# seaborn + pyplotによる折れ線グラフ
plt.plot(x, y, color='black')
plt.title("lineplot seaborn")
plt.xlabel("x")
plt.ylabel("y")

# seabornによるヒストグラム
fish_data = np.array([2,3,3,4,4,4,4,5,5,6])
sns.distplot(fish_data, bins=5, color='black', kde=False)

# カーネル密度推定によるヒストグラム平滑化
sns.distplot(fish_data, bins=1, color='black', kde=False)
sns.distplot(fish_data, color='black')

# 2変量データに対するヒストグラム
fish_multi = pd.read_csv("3-3-2-fish_multi_2.csv")
length_a = fish_multi.query('species == "A"')["length"]
length_b = fish_multi.query('species == "B"')["length"]

sns.distplot(length_a, bins=5, color='black', kde=False)
sns.distplot(length_b, bins=5, color='gray', kde=False)

# 箱髭図
sns.boxplot(x="species", y="length", data=fish_multi, color='gray')

# 散布図
cov_data = pd.read_csv("3-2-3-cov.csv")
sns.jointplot(x="x", y="y", data=cov_data, color='black')

# ペア・プロット
iris = sns.load_dataset("iris")
sns.pairplot(iris, hue="species")

```






```
# ライブラリのインポート
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
from matplotlib import pyplot as plt
import seaborn as sns

sns.set()

# 5尾の魚しかいない湖からの標本抽出
fish_5 = np.array([2, 3, 4, 5, 6])
print("魚の標本（5尾の魚）:", fish_5)

# 母集団からのランダムサンプリング
np.random.seed(1)  # シードを設定
sample1 = np.random.choice(fish_5, size=1, replace=False)
print("サンプル1:", sample1)

sample2 = np.random.choice(fish_5, size=3, replace=False)
print("サンプル2:", sample2)

sample3 = np.random.choice(fish_5, size=3, replace=False)
print("サンプル3:", sample3)

# 乱数の種を指定すると、同じ乱数が何度も出る
np.random.seed(1)
sample4 = np.random.choice(fish_5, size=3, replace=False)
print("サンプル4（シード1）:", sample4)

# 標本平均を計算するシミュレーション
np.random.seed(1)
sample_mean = np.mean(np.random.choice(fish_5, size=3, replace=False))
print("標本平均（シード1）:", sample_mean)

# もっとたくさんの魚がいる湖からの標本抽出
try:
    fish_100000 = pd.read_csv("3-4-1-fish_length_100000.csv")["length"]
    print("魚の標本（10万尾の魚）の先頭5行:\n", fish_100000.head())
    print("標本の長さ:", len(fish_100000))
    
    # サンプリングと標本平均の計算
    sampling_result = np.random.choice(fish_100000, size=10, replace=False)
    print("サンプリング結果:", sampling_result)
    
    sample_mean_100000 = np.mean(sampling_result)
    print("標本平均（10万尾の魚）:", sample_mean_100000)

    # 母集団分布
    population_mean = np.mean(fish_100000)
    print("母集団平均:", population_mean)

    population_std = np.std(fish_100000, ddof=0)
    print("母集団標準偏差:", population_std)

    population_var = np.var(fish_100000, ddof=0)
    print("母集団分散:", population_var)

    # データの可視化
    sns.histplot(fish_100000, kde=False, color='black')
    plt.title("魚の長さのヒストグラム（10万尾の魚）")
    plt.xlabel("魚の長さ")
    plt.ylabel("頻度")
    plt.show()

    # 正規分布との比較
    x = np.arange(start=1, stop=7.1, step=0.1)
    pdf = stats.norm.pdf(x=x, loc=4, scale=0.8)
    
    plt.plot(x, pdf, color='black', label='正規分布')
    sns.histplot(fish_100000, kde=False, stat='density', color='black', label='ヒストグラム')
    plt.title("魚の長さのヒストグラムと正規分布の比較")
    plt.xlabel("魚の長さ")
    plt.ylabel("密度")
    plt.legend()
    plt.show()

except FileNotFoundError:
    print("ファイル '3-4-1-fish_length_100000.csv' が見つかりませんでした。")


```







```
# Import Libraries
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns

# Create Population
mu, sigma = 4, 0.8
population = np.random.normal(mu, sigma, 10000)  # Large sample to represent the population

# Calculate Sample Means
sample_size = 30
num_samples = 1000
sample_means = np.zeros(num_samples)

for i in range(num_samples):
    sample = np.random.choice(population, sample_size, replace=False)
    sample_means[i] = np.mean(sample)

# Analyze Sample Means
mean_sample_means = np.mean(sample_means)
std_sample_means = np.std(sample_means)

plt.figure(figsize=(10, 6))
sns.histplot(sample_means, kde=True)
plt.title('Distribution of Sample Means')
plt.xlabel('Sample Mean')
plt.ylabel('Frequency')
plt.show()

# Effect of Sample Size
sample_sizes = [10, 30, 50, 100]
means_by_size = []

for size in sample_sizes:
    temp_means = [np.mean(np.random.choice(population, size, replace=False)) for _ in range(num_samples)]
    means_by_size.append(np.mean(temp_means))

plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, means_by_size, marker='o')
plt.title('Sample Mean vs. Sample Size')
plt.xlabel('Sample Size')
plt.ylabel('Mean of Sample Means')
plt.show()

# Function for Sample Mean
def compute_sample_means(sample_size, num_trials):
    means = [np.mean(np.random.choice(population, sample_size, replace=False)) for _ in range(num_trials)]
    return means

# Distribution of Sample Means with Varying Sample Sizes
all_means = []
for size in sample_sizes:
    means = compute_sample_means(size, num_samples)
    all_means.append(means)

plt.figure(figsize=(10, 6))
sns.violinplot(data=all_means, inner='quartile')
plt.title('Distribution of Sample Means for Different Sample Sizes')
plt.xlabel('Sample Size')
plt.ylabel('Sample Mean')
plt.xticks(ticks=range(len(sample_sizes)), labels=sample_sizes)
plt.show()

# Standard Error
standard_errors = [np.std(np.random.choice(population, size, replace=False)) / np.sqrt(size) for size in sample_sizes]

plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, standard_errors, marker='o', label='Standard Error')
plt.title('Standard Error vs. Sample Size')
plt.xlabel('Sample Size')
plt.ylabel('Standard Error')
plt.show()

# Sample Variance
variances = [np.var(np.random.choice(population, size, replace=False), ddof=0) for size in sample_sizes]

plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, variances, marker='o', label='Sample Variance')
plt.title('Sample Variance vs. Sample Size')
plt.xlabel('Sample Size')
plt.ylabel('Sample Variance')
plt.show()

# Unbiased Sample Variance
unbiased_variances = [np.var(np.random.choice(population, size, replace=False), ddof=1) for size in sample_sizes]

plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, unbiased_variances, marker='o', label='Unbiased Sample Variance')
plt.title('Unbiased Sample Variance vs. Sample Size')
plt.xlabel('Sample Size')
plt.ylabel('Unbiased Sample Variance')
plt.show()


```





```
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

%precision 3
%matplotlib inline

# 円周率
print(sp.pi)

# 指数関数
print(sp.exp(1))

# 「平均4、分散0.64(標準偏差0.8)の正規分布」における、確率変数が3であるときの確率密度
x = 3
mu = 4
sigma = 0.8
density = 1 / (sp.sqrt(2 * sp.pi * sigma**2)) * sp.exp(- ((x - mu)**2) / (2 * sigma**2))
print(density)

# scipyのnorm.pdf関数を使用して確率密度関数を計算
print(stats.norm.pdf(loc = 4, scale = 0.8, x = 3))

# normクラスのインスタンスを用いて確率密度関数を計算
norm_dist = stats.norm(loc = 4, scale = 0.8)
print(norm_dist.pdf(x = 3))

# 正規分布の確率密度関数のグラフを描画
x_plot = np.arange(start = 1, stop = 7.1, step = 0.1)
plt.plot(
    x_plot, 
    stats.norm.pdf(x = x_plot, loc = 4, scale = 0.8),
    color = 'black'
)
plt.title('Normal Distribution PDF')
plt.xlabel('x')
plt.ylabel('Probability Density')
plt.grid(True)
plt.show()

# 標本の生成
np.random.seed(1)
simulated_sample = stats.norm.rvs(
    loc = 4, scale = 0.8, size = 100000
)

# 標本中で値が3以下のデータの数をカウント
count = sp.sum(simulated_sample <= 3)
print(count)

# 標本全体に対する割合を計算
proportion = sp.sum(simulated_sample <= 3) / len(simulated_sample)
print(proportion)

# 累積分布関数 (CDF) の計算
cdf_3 = stats.norm.cdf(loc = 4, scale = 0.8, x = 3)
cdf_4 = stats.norm.cdf(loc = 4, scale = 0.8, x = 4)
print(cdf_3)
print(cdf_4)

# パーセント点 (分位点) の計算
percentile_025 = stats.norm.ppf(loc = 4, scale = 0.8, q = 0.025)
print(percentile_025)

# 確率0.025に対応するパーセント点
sitagawa = stats.norm.cdf(loc = 4, scale = 0.8, x = 3)
percentile_sitagawa = stats.norm.ppf(loc = 4, scale = 0.8, q = sitagawa)
print(percentile_sitagawa)

# 確率0.5（中央値）に対応するパーセント点
percentile_50 = stats.norm.ppf(loc = 4, scale = 0.8, q = 0.5)
print(percentile_50)

# t値の標本分布のシミュレーション
np.random.seed(1)
t_value_array = np.zeros(10000)
norm_dist = stats.norm(loc = 4, scale = 0.8)
for i in range(10000):
    sample = norm_dist.rvs(size = 10)
    sample_mean = sp.mean(sample)
    sample_std = sp.std(sample, ddof = 1)
    sample_se = sample_std / sp.sqrt(len(sample))
    t_value_array[i] = (sample_mean - 4) / sample_se

# t値のヒストグラムと標準正規分布の確率密度関数のプロット
sns.distplot(t_value_array, color = 'black')

x = np.arange(start = -8, stop = 8.1, step = 0.1)
plt.plot(x, stats.norm.pdf(x = x), 
         color = 'black', linestyle = 'dotted')
plt.title('T-value Distribution and Standard Normal Distribution')
plt.xlabel('x')
plt.ylabel('Probability Density')
plt.grid(True)
plt.show()

# t分布と標準正規分布の比較プロット
plt.plot(x, stats.norm.pdf(x = x), 
         color = 'black', linestyle = 'dotted')
plt.plot(x, stats.t.pdf(x = x, df = 9), 
         color = 'black')
plt.title('Comparison of T-distribution and Standard Normal Distribution')
plt.xlabel('x')
plt.ylabel('Probability Density')
plt.grid(True)
plt.show()

# 正規化されたt値のヒストグラムとt分布のプロット
sns.distplot(t_value_array, 
             color = 'black', norm_hist = True)
plt.plot(x, stats.t.pdf(x = x, df = 9), 
         color = 'black', linestyle = 'dotted')
plt.title('Normalized T-value Histogram and T-distribution')
plt.xlabel('x')
plt.ylabel('Probability Density')
plt.grid(True)
plt.show()


```



```
import pandas as pd
import numpy as np
import statsmodels.api as sm
import seaborn as sns
import matplotlib.pyplot as plt

# 仮想データの生成
np.random.seed(42)
n_samples = 100

# 遺伝子型データ（SNPs）の生成
genotype_data = np.random.randint(0, 3, size=(n_samples, 5))
genotype_df = pd.DataFrame(genotype_data, columns=[f'SNP_{i+1}' for i in range(genotype_data.shape[1])])

# 表現型データの生成
phenotype_data = np.random.normal(loc=genotype_data.sum(axis=1), scale=2)
phenotype_df = pd.DataFrame(phenotype_data, columns=['Phenotype'])

# データの結合
data = pd.concat([genotype_df, phenotype_df], axis=1)

# データの可視化
sns.pairplot(data)
plt.show()

# 遺伝子型データと表現型データの関係を調べる
X = sm.add_constant(genotype_df)  # 定数項を追加
y = phenotype_df

model = sm.OLS(y, X).fit()
print(model.summary())

# 特定のSNPの効果を調べる
snp = 'SNP_1'
sns.boxplot(x=data[snp], y=data['Phenotype'])
plt.title(f'Effect of {snp} on Phenotype')
plt.show()

# 遺伝子型-表現型の関係の可視化
for snp in genotype_df.columns:
    sns.boxplot(x=data[snp], y=data['Phenotype'])
    plt.title(f'Effect of {snp} on Phenotype')
    plt.show()


```

```

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 仮想データの生成
np.random.seed(42)

# 全人口における疾患の事前確率
P_D = 0.01  # 疾患の有病率（1%）

# 症状がある場合の疾患の確率
P_S_given_D = 0.8  # 疾患がある場合の症状の確率

# 症状がない場合の疾患の確率
P_S_given_not_D = 0.1  # 疾患がない場合の症状の確率

# 疾患のない場合の確率
P_not_D = 1 - P_D

# 症状がある全体の確率
P_S = P_S_given_D * P_D + P_S_given_not_D * P_not_D

# ベイズの定理を用いた事後確率の計算
P_D_given_S = (P_S_given_D * P_D) / P_S

print(f"症状がある場合の疾患の確率（事後確率）: {P_D_given_S:.4f}")

# 疾患リスクのプロット
labels = ['Disease', 'No Disease']
sizes = [P_D_given_S, 1 - P_D_given_S]
colors = ['lightcoral', 'lightskyblue']
explode = (0.1, 0)  # "explode" the 1st slice (i.e. 'Disease')

fig1, ax1 = plt.subplots()
ax1.pie(sizes, explode=explode, labels=labels, colors=colors,
        autopct='%1.1f%%', shadow=True, startangle=140)
ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.

plt.title('Probability of Disease Given Symptom')
plt.show()

```



```


```