평균, 중앙값, 최빈값¶

import matplotlib
matplotlib.use("Agg")
import pandas as pd
import matplotlib.pyplot as plt

# 그림에 한글을 쓰므로 한글 글꼴을 지정한다. 맥이면 'Apple SD Gothic Neo',
# 윈도우면 'Malgun Gothic', 리눅스면 'NanumGothic' 정도가 무난하다.
# 글꼴을 바꾸면 마이너스 기호가 깨지므로 unicode_minus 도 함께 꺼 준다.
plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
loans_data = pd.read_csv(url)

mean_income = loans_data['x'].mean()

fig, ax = plt.subplots(figsize=(12, 3))
# density=True 라서 y축이 밀도다. 소득 자료의 밀도는 1e-5 규모이므로
# 아래 세로선의 높이 1.6e-5 도 그 눈금에 맞춘 값이다.
ax.hist(loans_data['x'], bins=20, density=True, alpha=0.35,
        color="#1565C0", edgecolor="white")
# 평균 위치에 세로선을 긋는다. (x, x)와 (0, 높이)를 이어 그린 것이다.
ax.plot([mean_income, mean_income], [0, 1.6e-5], "--",
        color="#E65100", lw=2, label="평균")
ax.legend()
ax.set_title("소득 자료의 히스토그램과 평균")
ax.set_xlabel("소득 (달러)")
ax.set_ylabel("밀도")
ax.spines[["top", "right"]].set_visible(False)
fig.savefig("mean_median_mode_44.png", dpi=170, facecolor="white",
            bbox_inches="tight")

import numpy as np
import pandas as pd

url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists'
       '/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv')
income = pd.read_csv(url)['x'].values.astype(float)
m = income.mean()
dev = income - m

print(f"n = {len(income):,},  평균 {m:,.2f}")
print(f"편차의 합 = {dev.sum():.3e}   (0 이어야 한다)")
print(f"  자료 규모로 나눈 상대오차 = {abs(dev.sum()) / np.abs(income).sum():.1e}")
print(f"왼쪽 토크  sum(m - x),  x < m : {(m - income[income < m]).sum():>15,.0f}")
print(f"오른쪽 토크 sum(x - m),  x > m : {(income[income > m] - m).sum():>15,.0f}")

below, above = np.sum(income < m), np.sum(income > m)
print(f"\n평균 미만 {below:,}명 ({below / len(income):.2%})")
print(f"평균 초과 {above:,}명 ({above / len(income):.2%})")
print(f"중앙값 {np.median(income):,.0f} 은 평균보다 {m - np.median(income):,.0f} 작다")

n = 50,000,  평균 68,760.52
편차의 합 = -6.217e-08   (0 이어야 한다)
  자료 규모로 나눈 상대오차 = 1.8e-17
왼쪽 토크  sum(m - x),  x < m :     640,946,246
오른쪽 토크 sum(x - m),  x > m :     640,946,246

평균 미만 28,811명 (57.62%)
평균 초과 21,189명 (42.38%)
중앙값 62,000 은 평균보다 6,761 작다

import matplotlib
matplotlib.use("Agg")
import pandas as pd
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
loans_data = pd.read_csv(url)

# 같은 자료에 평균과 중앙값을 함께 표시한다.
# 오른쪽으로 치우친 분포에서는 평균이 중앙값보다 오른쪽에 놓인다.
# 긴 오른쪽 꼬리가 평균만 끌어당기기 때문이다.
mean_income = loans_data['x'].mean()
median_income = loans_data['x'].median()

fig, ax = plt.subplots(figsize=(12, 3))
ax.hist(loans_data['x'], bins=20, density=True, alpha=0.35,
        color="#1565C0", edgecolor="white")
ax.plot([mean_income, mean_income], [0, 1.6e-5], "--",
        color="#E65100", lw=2, label="평균")
ax.plot([median_income, median_income], [0, 1.6e-5], "--",
        color="#33691E", lw=2, label="중앙값")
ax.legend()
ax.set_title("소득 자료의 히스토그램과 평균, 중앙값")
ax.set_xlabel("소득 (달러)")
ax.set_ylabel("밀도")
ax.spines[["top", "right"]].set_visible(False)
fig.savefig("mean_median_mode_87.png", dpi=170, facecolor="white",
            bbox_inches="tight")

print(f"평균   {mean_income:>9,.0f}")
print(f"중앙값 {median_income:>9,.0f}")
print(f"차이   {mean_income - median_income:>9,.0f}  (양수 = 오른쪽 치우침)")

평균      68,761
중앙값    62,000
차이       6,761  (양수 = 오른쪽 치우침)

import numpy as np
import pandas as pd
from scipy import stats

url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists'
       '/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv')
income = pd.read_csv(url)['x'].values.astype(float)
mu, me, sd = income.mean(), np.median(income), income.std(ddof=1)
print(f"소득 자료 : 평균 {mu:,.2f}  중앙값 {me:,.2f}  표준편차 {sd:,.2f}")
print(f"            |평균 - 중앙값| / 표준편차 = {abs(mu - me) / sd:.6f}   (1 이하여야 한다)")

dists = {
    "정규(0,1)": stats.norm(),
    "지수(1)": stats.expon(),
    "로그정규(1)": stats.lognorm(s=1),
    "카이제곱(1)": stats.chi2(1),
    "파레토(3)": stats.pareto(3),
    "감마(0.2)": stats.gamma(0.2),
    "베타(0.1, 5)": stats.beta(0.1, 5),
}
print(f"\n{'분포':>12}{'평균':>12}{'중앙값':>12}{'표준편차':>12}{'|mu-m|/sigma':>14}")
for name, d in dists.items():
    mu_, me_, sd_ = d.mean(), d.median(), np.sqrt(d.var())
    print(f"{name:>12}{mu_:>12.6f}{me_:>12.6f}{sd_:>12.6f}{abs(mu_ - me_) / sd_:>14.6f}")

소득 자료 : 평균 68,760.52  중앙값 62,000.00  표준편차 32,872.04
            |평균 - 중앙값| / 표준편차 = 0.205662   (1 이하여야 한다)

          분포          평균         중앙값        표준편차  |mu-m|/sigma
     정규(0,1)    0.000000    0.000000    1.000000      0.000000
       지수(1)    1.000000    0.693147    1.000000      0.306853
     로그정규(1)    1.648721    1.000000    2.161197      0.300168
     카이제곱(1)    1.000000    0.454936    1.414214      0.385418
      파레토(3)    1.500000    1.259921    0.866025      0.277219
     감마(0.2)    0.200000    0.020746    0.447214      0.400823
  베타(0.1, 5)    0.019608    0.000130    0.056137      0.346967

import matplotlib
matplotlib.use("Agg")
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
loans_data = pd.read_csv(url)
income_data = loans_data['x'].values

# --- 1부: 원자료 ---
mean_income = income_data.mean()
median_income = np.median(income_data)

fig, (hist_ax, box_ax) = plt.subplots(1, 2, figsize=(12, 3))
fig.suptitle("원자료", fontsize=16)

n, bin_edges, _ = hist_ax.hist(income_data, bins=20, density=True, alpha=0.5,
                                color="#DCEBFB", edgecolor="#1565C0")
hist_ax.plot([mean_income, mean_income], [0, n.max()], "--",
             color="#E65100", lw=2, label="평균")
hist_ax.plot([median_income, median_income], [0, n.max()], "--",
             color="#33691E", lw=2, label="중앙값")
hist_ax.legend()
hist_ax.set_title("히스토그램")
hist_ax.set_xlabel("소득 (달러)")
hist_ax.set_ylabel("밀도")

box_ax.boxplot(income_data, vert=False, patch_artist=True)
box_ax.set_title("상자그림")
box_ax.set_xlabel("소득 (달러)")

for ax in (hist_ax, box_ax):
    ax.spines[["top", "right"]].set_visible(False)

fig.tight_layout()
fig.savefig("mean_median_mode_117_0.png", dpi=170, facecolor="white",
            bbox_inches="tight")

# --- 2부: 이상치를 넣는다 ---
# 5만 개 자료에 2천만 달러짜리 20개를 더한다. 전체의 0.04%에 불과하다.
# 그런데도 평균은 크게 밀리고 중앙값은 사실상 그대로다.
outliers = np.array([20_000_000] * 20)
data_with_outliers = np.concatenate((income_data, outliers))

mean_outliers = data_with_outliers.mean()
median_outliers = np.median(data_with_outliers)

fig, (hist_ax, box_ax) = plt.subplots(1, 2, figsize=(12, 3))
fig.suptitle("이상치를 넣은 뒤", fontsize=16)

# 구간 경계를 원자료와 똑같이 맞춘다. 그래야 두 히스토그램을 나란히 견줄 수 있다.
# 이상치는 이 구간 밖이라 막대로는 보이지 않고, 평균선의 이동으로만 드러난다.
n, bin_edges, _ = hist_ax.hist(data_with_outliers, bins=bin_edges,
                                density=True, alpha=0.5,
                                color="#DCEBFB", edgecolor="#1565C0")
hist_ax.plot([mean_outliers, mean_outliers], [0, n.max()], "--",
             color="#E65100", lw=2, label="평균")
hist_ax.plot([median_outliers, median_outliers], [0, n.max()], "--",
             color="#33691E", lw=2, label="중앙값")
hist_ax.legend()
hist_ax.set_title("히스토그램")
hist_ax.set_xlabel("소득 (달러)")
hist_ax.set_ylabel("밀도")

box_ax.boxplot(data_with_outliers, vert=False, patch_artist=True)
box_ax.set_title("상자그림")
box_ax.set_xlabel("소득 (달러)")

for ax in (hist_ax, box_ax):
    ax.spines[["top", "right"]].set_visible(False)

fig.tight_layout()
fig.savefig("mean_median_mode_117_1.png", dpi=170, facecolor="white",
            bbox_inches="tight")

# 이상치 20개(전체의 0.04%)가 두 측도를 각각 얼마나 움직였는가
print(f"{'':10}{'원자료':>14}{'이상치 추가 후':>18}{'변화율':>10}")
print(f"{'평균':10}{mean_income:>14,.0f}{mean_outliers:>18,.0f}"
      f"{(mean_outliers/mean_income - 1):>9.1%}")
print(f"{'중앙값':10}{median_income:>14,.0f}{median_outliers:>18,.0f}"
      f"{(median_outliers/median_income - 1):>9.1%}")

                     원자료          이상치 추가 후       변화율
평균                68,761            76,730    11.6%
중앙값               62,000            62,000     0.0%

import numpy as np
import pandas as pd

url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists'
       '/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv')
income = pd.read_csv(url)['x'].values.astype(float)
n, k, M = len(income), 20, 20_000_000.0
xbar = income.mean()
shifted = np.concatenate([income, np.full(k, M)])

pred = k * (M - xbar) / (n + k)
print(f"예측 이동량 k(M - xbar)/(n+k) = {pred:,.4f}")
print(f"실제 이동량                  = {shifted.mean() - xbar:,.4f}")
print(f"새 평균 {shifted.mean():,.4f}  ({shifted.mean() / xbar - 1:.2%} 증가)")

# 이상치 크기를 키우면 이동량도 비례해서 커진다.
print(f"\n{'M':>14}{'예측 이동량':>16}{'실제 이동량':>16}{'중앙값':>12}")
for Mi in (2e7, 2e8, 2e9, 1e30):
    q = np.concatenate([income, np.full(k, Mi)])
    print(f"{Mi:>14.1e}{k * (Mi - xbar) / (n + k):>16.4e}{q.mean() - xbar:>16.4e}"
          f"{np.median(q):>12,.0f}")

# 중앙값은 색인만 옮겨 간다.
srt, srt2 = np.sort(income), np.sort(shifted)
print(f"\n원자료 n = {n}:   x_({n // 2}) = {srt[n // 2 - 1]:,.0f},"
      f"  x_({n // 2 + 1}) = {srt[n // 2]:,.0f}   -> 중앙값 {np.median(income):,.0f}")
m2 = len(shifted)
print(f"이후   n = {m2}: x_({m2 // 2}) = {srt2[m2 // 2 - 1]:,.0f},"
      f"  x_({m2 // 2 + 1}) = {srt2[m2 // 2]:,.0f}   -> 중앙값 {np.median(shifted):,.0f}")
print(f"62,000 달러인 관측이 {np.sum(income == 62000):,}명이나 된다")

예측 이동량 k(M - xbar)/(n+k) = 7,969.3081
실제 이동량                  = 7,969.3081
새 평균 76,729.8265  (11.59% 증가)

             M          예측 이동량          실제 이동량         중앙값
       2.0e+07      7.9693e+03      7.9693e+03      62,000
       2.0e+08      7.9941e+04      7.9941e+04      62,000
       2.0e+09      7.9965e+05      7.9965e+05      62,000
       1.0e+30      3.9984e+26      3.9984e+26      62,000

원자료 n = 50000:   x_(25000) = 62,000,  x_(25001) = 62,000   -> 중앙값 62,000
이후   n = 50020: x_(25010) = 62,000,  x_(25011) = 62,000   -> 중앙값 62,000
62,000 달러인 관측이 429명이나 된다

import pandas as pd
from scipy.stats import trim_mean

# 미국 50개 주의 인구와 살인율. 오른쪽으로 크게 치우친 전형적인 자료다.
url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
state = pd.read_csv(url)

# 보통의 평균. 캘리포니아 같은 큰 주 하나에 끌려 올라간다.
mean_pop = state['Population'].mean()
print(f"평균       : {mean_pop:,.0f}")

# 절사평균. scipy 의 0.1 은 "양끝에서 각각 10%씩" 이다 (아래 (2) 참조).
trimmed_mean_pop = trim_mean(state['Population'], 0.1)
print(f"10% 절사평균: {trimmed_mean_pop:,.0f}")

# 중앙값은 순서만 보므로 꼬리가 아무리 길어도 흔들리지 않는다.
median_pop = state['Population'].median()
print(f"중앙값     : {median_pop:,.0f}")

평균       : 6,162,876
10% 절사평균: 4,783,697
중앙값     : 4,436,370

import numpy as np
import pandas as pd
from scipy.stats import trim_mean

url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists'
       '/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
pop = np.sort(pd.read_csv(url)['Population'].values.astype(float))
n = len(pop)
print(f"n = {n}")

# scipy 는 양끝에서 각각 proportiontocut 만큼 버린다.
print(f"\n{'인수':>6}{'양끝에서 버리는 개수':>22}{'남는 개수':>10}{'값':>14}")
for q in (0.05, 0.06, 0.10, 0.25):
    cut = int(n * q)
    print(f"{q:>6.2f}{cut:>22}{n - 2 * cut:>10}{trim_mean(pop, q):>14,.0f}")
    assert np.isclose(trim_mean(pop, q), pop[cut:n - cut].mean())

# 책의 정의 1 은 각 꼬리에서 ceil(p*n/2) 개를 버린다 -- 전체 p 다.
print(f"\n{'p':>6}{'ceil(pn/2)':>12}{'남는 개수':>10}{'정의 1 의 값':>16}{'scipy trim_mean(x, p)':>24}")
for p in (0.05, 0.10, 0.20, 0.40):
    c = int(np.ceil(p * n / 2))
    print(f"{p:>6.2f}{c:>12}{n - 2 * c:>10}{pop[c:n - c].mean():>16,.0f}"
          f"{trim_mean(pop, p):>24,.0f}")

print(f"\n평균 {pop.mean():,.0f},  중앙값 {np.median(pop):,.0f}")

n = 50

    인수           양끝에서 버리는 개수     남는 개수             값
  0.05                     2        46     5,316,412
  0.06                     3        44     5,102,369
  0.10                     5        40     4,783,697
  0.25                    12        26     4,334,488

     p  ceil(pn/2)     남는 개수        정의 1 의 값   scipy trim_mean(x, p)
  0.05           2        46       5,316,412               5,316,412
  0.10           3        44       5,102,369               4,783,697
  0.20           5        40       4,783,697               4,413,916
  0.40          10        30       4,413,916               4,281,384

평균 6,162,876,  중앙값 4,436,370

import pandas as pd
import numpy as np

# 미국 50개 주의 인구와 살인율. 오른쪽으로 크게 치우친 전형적인 자료다.
url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
state = pd.read_csv(url)

# 가중하지 않은 평균은 주 50개를 똑같이 한 표씩 센다.
# 인구 60만의 와이오밍과 인구 3900만의 캘리포니아가 같은 무게를 갖는다.
unweighted_mean = state['Murder.Rate'].mean()
print(f"단순평균 살인율: {unweighted_mean:.3f}")

# 인구로 가중하면 사람 한 명이 한 표가 된다. "미국 사람이 겪는 평균"에 가깝다.
weighted_mean = np.average(state['Murder.Rate'], weights=state['Population'])
print(f"가중평균 살인율: {weighted_mean:.3f}")

단순평균 살인율: 4.066
가중평균 살인율: 4.446

import numpy as np
import pandas as pd

url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists'
       '/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
state = pd.read_csv(url)
pop = state['Population'].values.astype(float)
rate = state['Murder.Rate'].values.astype(float)

uw = rate.mean()
w = np.average(rate, weights=pop)
cov = np.mean((pop - pop.mean()) * (rate - rate.mean()))
print(f"단순평균 {uw:.6f}")
print(f"가중평균 {w:.6f}")
print(f"차       {w - uw:.6f}")
print(f"Cov(w, x) / wbar = {cov / pop.mean():.6f}   <- 같아야 한다")
print(f"상관계수 corr(인구, 살인율) = {np.corrcoef(pop, rate)[0, 1]:.4f}")

print(f"\n인구 상위 5 개 주")
top = np.argsort(pop)[::-1][:5]
for j in top:
    print(f"  {state['State'][j]:>12}  인구 {pop[j]:>12,.0f}  살인율 {rate[j]:>5.1f}"
          f"  가중치 {pop[j] / pop.sum():.4f}")
print(f"  상위 5 개 주가 가중치의 {pop[top].sum() / pop.sum():.1%} 를 갖는다")

print(f"\n총 살인 건수 / 총 인구 (10만 명당) = {np.sum(rate * pop) / np.sum(pop):.6f}")

단순평균 4.066000
가중평균 4.445834
차       0.379834
Cov(w, x) / wbar = 0.379834   <- 같아야 한다
상관계수 corr(인구, 살인율) = 0.1821

인구 상위 5 개 주
    California  인구   37,253,956  살인율   4.4  가중치 0.1209
         Texas  인구   25,145,561  살인율   4.4  가중치 0.0816
      New York  인구   19,378,102  살인율   3.1  가중치 0.0629
       Florida  인구   18,801,310  살인율   5.8  가중치 0.0610
      Illinois  인구   12,830,632  살인율   5.3  가중치 0.0416
  상위 5 개 주가 가중치의 36.8% 를 갖는다

총 살인 건수 / 총 인구 (10만 명당) = 4.445834

import pandas as pd

# 미국 50개 주의 인구와 살인율. 오른쪽으로 크게 치우친 전형적인 자료다.
url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
state = pd.read_csv(url)

# 가중하지 않은 중앙값: 주를 크기와 상관없이 한 표씩 센다.
# 즉 캘리포니아(3900만 명)와 와이오밍(56만 명)이 같은 무게를 갖는다.
unweighted_median = state['Murder.Rate'].median()
print(f"Unweighted Median: {unweighted_median:.1f}")


def weighted_median(values, weights):
    """누적 가중치가 전체의 절반에 도달하는 값을 찾는다.

    가중중앙값의 정의 그대로다. 외부 패키지 없이 세 줄로 구현된다.
    """
    d = pd.DataFrame({"v": values, "w": weights}).sort_values("v")
    cum = d["w"].cumsum() / d["w"].sum()      # 누적 가중치 비율
    return d.loc[cum >= 0.5, "v"].iloc[0]     # 0.5를 처음 넘는 값

# 인구로 가중한 중앙값: 사람 한 명씩을 세는 것과 같다.
# "미국의 중앙값 시민이 사는 주의 살인율"이라고 읽으면 된다.
wm = weighted_median(state['Murder.Rate'], state['Population'])
print(f"Weighted Median: {wm:.1f}")

Unweighted Median: 4.0
Weighted Median: 4.4

import numpy as np
import pandas as pd

url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists'
       '/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
state = pd.read_csv(url)
pop = state['Population'].values.astype(float)
rate = state['Murder.Rate'].values.astype(float)

grid = np.linspace(rate.min(), rate.max(), 20001)
wl1 = np.array([np.sum(pop * np.abs(rate - c)) for c in grid])
wl2 = np.array([np.sum(pop * (rate - c) ** 2) for c in grid])
l1 = np.array([np.sum(np.abs(rate - c)) for c in grid])

print(f"{'목적함수':>22}{'격자 최소점':>14}{'닫힌 꼴의 답':>16}")
print(f"{'sum w|x - c|':>22}{grid[wl1.argmin()]:>14.4f}{4.4:>16.4f}")
print(f"{'sum w(x - c)^2':>22}{grid[wl2.argmin()]:>14.4f}"
      f"{np.average(rate, weights=pop):>16.4f}")
print(f"{'sum |x - c| (가중 없음)':>22}{grid[l1.argmin()]:>14.4f}{np.median(rate):>16.4f}")
print(f"격자 간격 {grid[1] - grid[0]:.6f}")

# 누적 가중치가 0.5 를 넘는 자리
d = pd.DataFrame({"v": rate, "w": pop}).sort_values("v")
cum = (d["w"].cumsum() / d["w"].sum()).values
v = d["v"].values
j = int(np.argmax(cum >= 0.5))
print(f"\n누적 가중치가 0.5 를 처음 넘는 자리: 값 {v[j]}, 누적 {cum[j]:.6f}")
print(f"  바로 앞 값 {v[j - 1]}, 누적 {cum[j - 1]:.6f}")

                  목적함수        격자 최소점         닫힌 꼴의 답
          sum w|x - c|        4.4001          4.4000
        sum w(x - c)^2        4.4457          4.4458
   sum |x - c| (가중 없음)        4.0001          4.0000
격자 간격 0.000470

누적 가중치가 0.5 를 처음 넘는 자리: 값 4.4, 누적 0.563949
  바로 앞 값 4.4, 누적 0.443051

import statistics

# 최빈값은 가장 자주 나온 값이다. 여기서는 2 만 두 번 나온다.
data = [4, 1, 2, 2, 3, 5]
mode = statistics.mode(data)
print(f"{mode = }")  # mode = 2

mode = 2

import statistics

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = rng.normal(0, 1, 1000)          # 참 최빈값은 0
print(f"n = {len(x)},  서로 다른 값의 개수 = {len(set(x))}   (모두 다르다)")

r = stats.mode(x, keepdims=False)
print(f"\nscipy.stats.mode  -> 값 {r.mode:.6f},  도수 {r.count}")
print(f"  자료의 최솟값     =  {x.min():.6f}   <- 같은 값이다")
print(f"statistics.mode   -> 값 {statistics.mode(x):.6f}")
print(f"  자료의 첫 원소    =  {x[0]:.6f}   <- 같은 값이다")
print(f"statistics.multimode 가 돌려주는 값의 개수 = {len(statistics.multimode(list(x)))}")

# 구간을 잡아 주면? 어떻게 잡느냐에 따라 답이 달라진다.
print(f"\n{'반올림 자리':>12}{'최빈값':>10}{'그 도수':>9}")
for d in (0, 1, 2):
    xr = np.round(x, d)
    vals, cnts = np.unique(xr, return_counts=True)
    print(f"{d:>12}{vals[cnts.argmax()]:>10.2f}{cnts.max():>9}")

n = 1000,  서로 다른 값의 개수 = 1000   (모두 다르다)

scipy.stats.mode  -> 값 -3.899422,  도수 1
  자료의 최솟값     =  -3.899422   <- 같은 값이다
statistics.mode   -> 값 0.125730
  자료의 첫 원소    =  0.125730   <- 같은 값이다
statistics.multimode 가 돌려주는 값의 개수 = 1000

      반올림 자리       최빈값     그 도수
           0      0.00      389
           1     -0.30       50
           2      0.36       11

import statistics

# 이번에는 2 와 3 이 나란히 두 번씩 나온다. 최빈값이 둘인 자료다.
data = [4, 1, 2, 2, 3, 3, 5]

# mode 는 그중 먼저 나온 것 하나만 돌려준다. 나머지 하나가 조용히 감춰진다.
mode = statistics.mode(data)
print(f"{mode = }")

# 그래서 최빈값이 여럿일 수 있는 자료에는 multimode 를 쓴다.
modes = statistics.multimode(data)
print(f"{modes = }")  # [2, 3]

mode = 2
modes = [2, 3]

import statistics

import numpy as np
from scipy import stats

# 평균 > 중앙값 인데 왼쪽으로 치우친 자료
a = np.array([0, 7, 7, 7, 7, 9, 9, 9, 9], dtype=float)
print(f"자료 {a.astype(int)}")
print(f"평균 {a.mean():.4f}   중앙값 {np.median(a):.1f}   최빈값 {statistics.multimode(list(a))}")
print(f"평균 - 중앙값 = {a.mean() - np.median(a):+.4f}   -> 어림은 '오른쪽 치우침' 이라 한다")
print(f"표본왜도 g1   = {stats.skew(a):+.4f}            -> 실제로는 왼쪽 치우침")
print(f"  (불편보정 G1 = {stats.skew(a, bias=False):+.4f})")

print(f"\n{'c':>6}{'sum|x - c|':>12}{'sum(x - c)^2':>14}")
for c in (6.0, 7.0, 7.111111, 8.0, 9.0):
    print(f"{c:>6.2f}{np.abs(a - c).sum():>12.4f}{((a - c) ** 2).sum():>14.4f}")

# 거울상을 보면 부호가 모두 뒤집힌다
b = -a
print(f"\n거울상 -x : 평균 {b.mean():.4f}  중앙값 {np.median(b):.1f}  왜도 {stats.skew(b):+.4f}")

자료 [0 7 7 7 7 9 9 9 9]
평균 7.1111   중앙값 7.0   최빈값 [7.0, 9.0]
평균 - 중앙값 = +0.1111   -> 어림은 '오른쪽 치우침' 이라 한다
표본왜도 g1   = -1.9092            -> 실제로는 왼쪽 치우침
  (불편보정 G1 = -2.3143)

     c  sum|x - c|  sum(x - c)^2
  6.00     22.0000       76.0000
  7.00     15.0000       65.0000
  7.11     15.1111       64.8889
  8.00     16.0000       72.0000
  9.00     17.0000       97.0000

거울상 -x : 평균 -7.1111  중앙값 -7.0  왜도 +1.9092

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
base = rng.normal(50, 10, 100)

print(f"{'오염 비율':>10}{'평균':>12}{'10% 절단':>12}{'25% 절단':>12}{'중앙값':>11}")
for frac in (0.0, 0.05, 0.10, 0.20, 0.30, 0.49):
    x = base.copy()
    k = int(frac * 100)
    if k:
        x[:k] = 1e6                      # 극단값으로 오염시킨다
    print(f"{frac:>10.2f}{x.mean():>12.1f}{stats.trim_mean(x, 0.10):>12.2f}"
          f"{stats.trim_mean(x, 0.25):>12.2f}{np.median(x):>11.2f}")

     오염 비율          평균      10% 절단      25% 절단        중앙값
      0.00        50.8       50.73       50.38      50.73
      0.05     50048.3       51.69       51.21      51.75
      0.10    100045.7       52.77       52.08      52.17
      0.20    200041.2   125046.93       55.36      54.12
      0.30    300036.2   250040.68   100051.97      57.39
      0.49    490025.8   487527.68   480030.67      69.12

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
B, n = 200_000, 25

print(f"{'분포':>7}{'Var(평균)':>14}{'Var(중앙값)':>13}{'Var(25% 절단)':>15}{'중앙값 효율':>13}")
for label, draw in [("정규", lambda: rng.normal(0, 1, (B, n))),
                    ("t(3)", lambda: rng.standard_t(3, (B, n))),
                    ("코시", lambda: rng.standard_cauchy((B, n)))]:
    x = draw()
    v_mean = x.mean(1).var()
    v_med = np.median(x, axis=1).var()
    v_trim = stats.trim_mean(x, 0.25, axis=1).var()
    print(f"{label:>7}{v_mean:>14.4f}{v_med:>13.4f}{v_trim:>15.4f}"
          f"{v_mean / v_med:>13.4f}")

print(f"\n정규분포에서 중앙값의 점근 상대효율 이론값 2/pi = {2 / np.pi:.4f}")

     분포       Var(평균)     Var(중앙값)    Var(25% 절단)       중앙값 효율
     정규        0.0400       0.0618         0.0471       0.6483
   t(3)        0.1203       0.0751         0.0626       1.6004
     코시   145970.0658       0.1113         0.1242 1311827.9195

정규분포에서 중앙값의 점근 상대효율 이론값 2/pi = 0.6366

import numpy as np
from scipy import stats

# (1) 곱해지는 양: 수익률
r = np.array([0.50, -0.30, 0.40, -0.20, 0.25])
growth = np.prod(1 + r)
n = len(r)
print("연간 수익률", r)
print(f"  산술평균 {r.mean():.4%}  → 5년 뒤 {(1 + r.mean()) ** n:.4f}배  (틀림)")
print(f"  기하평균 {growth ** (1 / n) - 1:.4%}  → 5년 뒤 {growth:.4f}배")
print(f"  실제 누적 {growth:.4f}배")

# (2) 나누어지는 양: 속력
d, v = 100.0, np.array([60.0, 40.0])
print(f"\n같은 거리 {d}km 를 각각 {v[0]}, {v[1]} km/h 로 달리면")
print(f"  산술평균 {v.mean():.2f} km/h  (틀림)")
print(f"  조화평균 {len(v) / np.sum(1 / v):.2f} km/h")
print(f"  실제: {2 * d}km 를 {d / v[0] + d / v[1]:.4f}시간 → "
      f"{2 * d / (d / v[0] + d / v[1]):.2f} km/h")

x = np.array([1.0, 2.0, 4.0, 8.0])
print(f"\nAM ≥ GM ≥ HM:  {x.mean():.4f} ≥ {stats.gmean(x):.4f} ≥ {stats.hmean(x):.4f}")

연간 수익률 [ 0.5  -0.3   0.4  -0.2   0.25]
  산술평균 13.0000%  → 5년 뒤 1.8424배  (틀림)
  기하평균 8.0099%  → 5년 뒤 1.4700배
  실제 누적 1.4700배

같은 거리 100.0km 를 각각 60.0, 40.0 km/h 로 달리면
  산술평균 50.00 km/h  (틀림)
  조화평균 48.00 km/h
  실제: 200.0km 를 4.1667시간 → 48.00 km/h

AM ≥ GM ≥ HM:  3.7500 ≥ 2.8284 ≥ 2.1333

import numpy as np
from scipy.stats import gaussian_kde

rng = np.random.default_rng(1)
n, B = 300, 2000
grid = np.linspace(-5, 5, 1000)

mean = np.empty(B); med = np.empty(B)
mode_hist = np.empty(B); mode_kde = np.empty(B)
for b in range(B):
    d = rng.normal(0, 1, n)                  # 참 최빈값은 0
    mean[b] = d.mean()
    med[b] = np.median(d)
    h, e = np.histogram(d, bins=20)
    mode_hist[b] = ((e[:-1] + e[1:]) / 2)[h.argmax()]
    mode_kde[b] = grid[gaussian_kde(d)(grid).argmax()]

print(f"N(0,1) 에서 n={n}, 반복 {B}회 (참값 0)")
for label, v in [("표본평균", mean), ("표본중앙값", med),
                 ("최빈값 (히스토그램 20구간)", mode_hist), ("최빈값 (KDE)", mode_kde)]:
    print(f"  {label:<26} 평균 {v.mean():+.4f}   표준편차 {v.std():.4f}")

print(f"\n최빈값(히스토그램)은 표본평균보다 {mode_hist.std() / mean.std():.1f}배 흔들린다")
print(f"최빈값(KDE)은          {mode_kde.std() / mean.std():.1f}배")

N(0,1) 에서 n=300, 반복 2000회 (참값 0)
  표본평균                       평균 -0.0014   표준편차 0.0590
  표본중앙값                      평균 -0.0010   표준편차 0.0714
  최빈값 (히스토그램 20구간)           평균 -0.0260   표준편차 0.3227
  최빈값 (KDE)                  평균 +0.0059   표준편차 0.1981

최빈값(히스토그램)은 표본평균보다 5.5배 흔들린다
최빈값(KDE)은          3.4배