히스토그램과 밀도 그림¶

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def plot_histogram_with_density(data, figsize=(12, 3)):
    """히스토그램에 적합한 정규 확률밀도함수를 겹쳐 그린다.

    자료의 평균과 표준편차로 만든 곡선이므로 위치와 척도는 이미 맞아 있다.
    남는 차이는 오직 모양이다.

    매개변수
    --------
    data : 그릴 자료
    figsize : 그림 크기 (가로, 세로)
    """
    fig, ax = plt.subplots(figsize=figsize)

    # density=True 로 넓이의 합을 1 로 맞춰야 확률밀도함수와 겹칠 수 있다.
    _, bins, _ = ax.hist(data, bins=20, density=True, alpha=0.5, label="Data Histogram")

    mu = data.mean()
    sigma = data.std()
    pdf = stats.norm(loc=mu, scale=sigma).pdf(bins)

    ax.plot(bins, pdf, "--r", label="Normal PDF")

    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    ax.set_title('Histogram with Density Plot')
    ax.set_xlabel('Value')
    ax.set_ylabel('Density')
    ax.legend()

    plt.show()

if __name__ == "__main__":
    # 정규자료: 곡선이 히스토그램에 잘 얹힌다.
    np.random.seed(0)
    sample_data = np.random.normal(loc=0, scale=1, size=1000)
    plot_histogram_with_density(sample_data)

import numpy as np
from scipy import stats

np.random.seed(0)
x = np.random.normal(loc=0, scale=1, size=1000)

# 쪽의 함수가 쓰는 것과 같은 설정: bins=20, density=True, sigma = std(ddof=0).
mu, sigma = x.mean(), x.std()
dens, bins = np.histogram(x, bins=20, density=True)
cnt, _ = np.histogram(x, bins=20)
k = int(dens.argmax())

print(f"mu = {mu:.4f}   sigma(ddof=0) = {sigma:.4f}   sigma(ddof=1) = {x.std(ddof=1):.4f}")
print(f"구간 폭 = {bins[1] - bins[0]:.4f}   구간 수 = 20")
print(f"최고 막대 밀도 = {dens.max():.4f} @ 구간중심 {(bins[k] + bins[k + 1]) / 2:+.4f}")
print(f"정규 pdf 최대  = {stats.norm(mu, sigma).pdf(mu):.4f} @ mu = {mu:+.4f}")
print(f"양끝 구간 도수: {cnt[:3]} ... {cnt[-3:]}")


def n_peaks(counts):
    """막대 높이에서 국소최대의 개수를 센다."""
    pad = np.r_[-1.0, counts.astype(float), -1.0]
    return int(np.sum((pad[1:-1] > pad[:-2]) & (pad[1:-1] >= pad[2:])))


print("\n구간 수를 바꿔 가며 봉우리를 센다 (자료는 그대로다):")
for m in (5, 10, 20, 40, 80, 160):
    c, b = np.histogram(x, bins=m)
    print(f"  구간 {m:3d}개  폭 {b[1] - b[0]:.4f}  봉우리 {n_peaks(c):2d}개"
          f"  최고 도수 {c.max():3d}  빈 구간 {int((c == 0).sum()):2d}개")

n = len(x)
q1, q3 = np.percentile(x, [25, 75])
h_fd = 2 * (q3 - q1) * n ** (-1 / 3)
h_sc = 3.49 * x.std(ddof=1) * n ** (-1 / 3)
rng_ = x.max() - x.min()
print(f"\nSturges  k = 1 + log2(n) = {1 + np.log2(n):.2f}")
print(f"Scott    h = {h_sc:.4f}  ->  k = {rng_ / h_sc:.2f}")
print(f"FD       h = {h_fd:.4f}  ->  k = {rng_ / h_fd:.2f}")

mu = -0.0453   sigma(ddof=0) = 0.9870   sigma(ddof=1) = 0.9875
구간 폭 = 0.2903   구간 수 = 20
최고 막대 밀도 = 0.4237 @ 구간중심 +0.2920
정규 pdf 최대  = 0.4042 @ mu = -0.0453
양끝 구간 도수: [3 6 9] ... [16 11  4]

구간 수를 바꿔 가며 봉우리를 센다 (자료는 그대로다):
  구간   5개  폭 1.1611  봉우리  1개  최고 도수 456  빈 구간  0개
  구간  10개  폭 0.5805  봉우리  1개  최고 도수 239  빈 구간  0개
  구간  20개  폭 0.2903  봉우리  1개  최고 도수 123  빈 구간  0개
  구간  40개  폭 0.1451  봉우리 10개  최고 도수  64  빈 구간  0개
  구간  80개  폭 0.0726  봉우리 27개  최고 도수  37  빈 구간  5개
  구간 160개  폭 0.0363  봉우리 53개  최고 도수  25  빈 구간 23개

Sturges  k = 1 + log2(n) = 10.97
Scott    h = 0.3446  ->  k = 16.84
FD       h = 0.2611  ->  k = 22.24

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def plot_histogram_with_density(data, figsize=(12, 3)):
    """히스토그램에 적합한 정규 확률밀도함수를 겹쳐 그린다.

    자료의 평균과 표준편차로 만든 곡선이므로, 위치와 척도는 이미 맞아
    있다. 남는 차이는 오직 모양이다.
    """
    fig, ax = plt.subplots(figsize=figsize)
    _, bins, _ = ax.hist(data, bins=20, density=True, alpha=0.5, label="Data Histogram")

    mu = data.mean()
    sigma = data.std()
    pdf = stats.norm(loc=mu, scale=sigma).pdf(bins)
    ax.plot(bins, pdf, "--r", label="Normal PDF")

    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.set_title('Histogram with Density Plot')
    ax.set_xlabel('Value')
    ax.set_ylabel('Density')
    ax.legend()
    plt.show()

if __name__ == "__main__":
    # 지수자료: 곡선이 왼쪽에서 음수 구간까지 뻗어 나가고 봉우리 자리도
    # 어긋난다. 정규로 볼 수 없다는 것이 한눈에 보인다.
    np.random.seed(0)
    sample_data = np.random.exponential(scale=1, size=1000)
    plot_histogram_with_density(sample_data)

import numpy as np
from scipy import stats

np.random.seed(0)
x = np.random.exponential(scale=1, size=1000)

mu, sigma = x.mean(), x.std()
dens, bins = np.histogram(x, bins=20, density=True)
cnt, _ = np.histogram(x, bins=20)
fit = stats.norm(mu, sigma)
k = int(dens.argmax())

print(f"mu = {mu:.4f}   sigma = {sigma:.4f}   (Exp(1) 의 참값은 둘 다 1)")
print(f"구간 폭 = {bins[1] - bins[0]:.4f}")
print(f"최고 막대 밀도 = {dens.max():.4f} @ 구간중심 {(bins[k] + bins[k + 1]) / 2:.4f}")
print(f"정규 pdf 최대  = {fit.pdf(mu):.4f} @ mu = {mu:.4f}   (높이 비 {dens.max() / fit.pdf(mu):.2f})")
print(f"자료 최소 = {x.min():.4f},  Exp 의 참 최빈값 = 0")
print(f"왼쪽 세 구간 도수 = {cnt[:3]},  오른쪽 세 구간 = {cnt[-3:]}")

# 적합한 정규가 음수 구간에 얼마나 질량을 놓는가. 지수자료는 거기 있을 수 없다.
print(f"P(적합 정규 < 0) = {fit.cdf(0):.6f} = {100 * fit.cdf(0):.2f}%")
print(f"실제로 음수인 관측값 = {int((x < 0).sum())}개")
print(f"g1 = {stats.skew(x):.4f}  (Exp 의 이론 왜도 = 2)")
print(f"g2 = {stats.kurtosis(x):.4f}  (Exp 의 이론 초과첨도 = 6)")

mu = 1.0035   sigma = 1.0291   (Exp(1) 의 참값은 둘 다 1)
구간 폭 = 0.4280
최고 막대 밀도 = 0.8387 @ 구간중심 0.2146
정규 pdf 최대  = 0.3877 @ mu = 1.0035   (높이 비 2.16)
자료 최소 = 0.0005,  Exp 의 참 최빈값 = 0
왼쪽 세 구간 도수 = [359 226 143],  오른쪽 세 구간 = [0 0 1]
P(적합 정규 < 0) = 0.164738 = 16.47%
실제로 음수인 관측값 = 0개
g1 = 2.0526  (Exp 의 이론 왜도 = 2)
g2 = 6.4761  (Exp 의 이론 초과첨도 = 6)

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def plot_histogram_with_density(data, figsize=(12, 3)):
    """히스토그램에 적합한 정규 확률밀도함수를 겹쳐 그린다.

    자료의 평균과 표준편차로 만든 곡선이므로, 위치와 척도는 이미 맞아
    있다. 남는 차이는 오직 모양이다.
    """
    fig, ax = plt.subplots(figsize=figsize)
    _, bins, _ = ax.hist(data, bins=20, density=True, alpha=0.5, label="Data Histogram")

    mu = data.mean()
    sigma = data.std()
    pdf = stats.norm(loc=mu, scale=sigma).pdf(bins)
    ax.plot(bins, pdf, "--r", label="Normal PDF")

    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.set_title('Histogram with Density Plot')
    ax.set_xlabel('Value')
    ax.set_ylabel('Density')
    ax.legend()
    plt.show()

if __name__ == "__main__":
    # 카이제곱(자유도 10): 지수보다 훨씬 가깝지만 왼쪽 어깨가 조금 어긋난다.
    np.random.seed(0)
    sample_data = np.random.chisquare(df=10, size=1000)
    plot_histogram_with_density(sample_data)

import numpy as np
from scipy import stats

np.random.seed(0)
x = np.random.chisquare(df=10, size=1000)

mu, sigma = x.mean(), x.std()
dens, bins = np.histogram(x, bins=20, density=True)
fit = stats.norm(mu, sigma)
k = int(dens.argmax())
ctr = (bins[:-1] + bins[1:]) / 2

print(f"mu = {mu:.4f}   sigma = {sigma:.4f}   (chi2_10 의 참값: 10 과 sqrt(20) = {np.sqrt(20):.4f})")
print(f"구간 폭 = {bins[1] - bins[0]:.4f}")
print(f"최고 막대 밀도 = {dens.max():.4f} @ 구간중심 {ctr[k]:.4f}")
print(f"정규 pdf 최대  = {fit.pdf(mu):.4f} @ mu = {mu:.4f}")
print(f"chi2_10 의 참 최빈값 = d - 2 = 8")
print(f"P(적합 정규 < 0) = {fit.cdf(0):.6f} = {100 * fit.cdf(0):.2f}%")

# 막대와 곡선의 차이를 구간마다 재서 가장 크게 어긋난 자리를 찾는다.
diff = dens - fit.pdf(ctr)
j = int(np.abs(diff).argmax())
print(f"가장 크게 어긋난 구간: 중심 {ctr[j]:.3f}, 막대 {dens[j]:.4f} - 곡선 {fit.pdf(ctr[j]):.4f}"
      f" = {diff[j]:+.4f}  (곡선의 {100 * diff[j] / fit.pdf(ctr[j]):+.1f}%)")
print("구간중심 / 막대 / 곡선 / 차:")
for i in range(0, 20, 2):
    print(f"  {ctr[i]:6.2f}  {dens[i]:.4f}  {fit.pdf(ctr[i]):.4f}  {diff[i]:+.4f}")
print(f"g1 = {stats.skew(x):.4f}  (이론 sqrt(8/10) = {np.sqrt(0.8):.4f})")
print(f"g2 = {stats.kurtosis(x):.4f}  (이론 12/10 = 1.2)")

mu = 9.9224   sigma = 4.2793   (chi2_10 의 참값: 10 과 sqrt(20) = 4.4721)
구간 폭 = 1.2419
최고 막대 밀도 = 0.0974 @ 구간중심 7.2509
정규 pdf 최대  = 0.0932 @ mu = 9.9224
chi2_10 의 참 최빈값 = d - 2 = 8
P(적합 정규 < 0) = 0.010205 = 1.02%
가장 크게 어긋난 구간: 중심 6.009, 막대 0.0854 - 곡선 0.0614 = +0.0240  (곡선의 +39.1%)
구간중심 / 막대 / 곡선 / 차:
    2.28  0.0072  0.0189  -0.0117
    4.77  0.0580  0.0451  +0.0129
    7.25  0.0974  0.0767  +0.0207
    9.73  0.0950  0.0931  +0.0019
   12.22  0.0636  0.0807  -0.0171
   14.70  0.0370  0.0500  -0.0129
   17.19  0.0169  0.0221  -0.0052
   19.67  0.0089  0.0070  +0.0019
   22.15  0.0032  0.0016  +0.0017
   24.64  0.0040  0.0003  +0.0038
g1 = 0.8519  (이론 sqrt(8/10) = 0.8944)
g2 = 0.8951  (이론 12/10 = 1.2)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

rng = np.random.default_rng(42)
data = rng.normal(5, 2, 500)

fig, ax = plt.subplots()
ax.hist(data, bins=30, density=True, alpha=0.5, edgecolor="black", label="Histogram")

kde = gaussian_kde(data)
x_grid = np.linspace(data.min() - 1, data.max() + 1, 200)
ax.plot(x_grid, kde(x_grid), "r-", lw=2, label="KDE")

ax.set_xlabel("Value")
ax.set_ylabel("Density")
ax.legend()
plt.show()