統計学メモ

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()