봉우리 수¶

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

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

np.random.seed(0)

# 중심이 0 과 6 으로 떨어진 정규분포 둘에서 각각 1000개씩 뽑는다.
data_normal_1 = stats.norm().rvs(1_000)
data_normal_2 = stats.norm(loc=6).rvs(1_000)

# 둘을 한 자료로 섞는다. 어느 값이 어느 쪽에서 왔는지는 지워진다.
combined_data = np.concatenate((data_normal_1, data_normal_2))

# 그래도 히스토그램에는 봉우리가 둘로 남는다. 섞인 흔적이 모양에 드러난 것이다.
fig, ax = plt.subplots(figsize=(12, 3))
ax.hist(combined_data, bins=30, color='skyblue', edgecolor='black')
ax.set_title("이봉 분포의 히스토그램")
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()

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

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


def mixture(mu, sigma=1.0, m=200_001):
    """중심이 ±mu 인 등가중 정규혼합의 밀도를 촘촘한 격자에서 잰다."""
    g = np.linspace(-(mu + 6 * sigma), mu + 6 * sigma, m)
    d = 0.5 * (stats.norm.pdf(g, -mu, sigma) + stats.norm.pdf(g, mu, sigma))
    return g, d


def peaks_and_dip(mu, sigma=1.0, w=200):
    """봉우리 위치와 골의 깊이 1 - f(0)/f(최대) 를 돌려준다."""
    g, d = mixture(mu, sigma)
    out = []
    for i in range(w, len(g) - w):
        # 창 안에서 최대이고 양쪽 끝보다 높아야 봉우리로 센다.
        # 평탄부를 봉우리 여럿으로 세지 않기 위한 장치다.
        if d[i] == d[i - w:i + w + 1].max() and d[i] > d[i - w] and d[i] > d[i + w]:
            if not out or abs(g[i] - out[-1]) > 0.05:
                out.append(g[i])
    dip = 1 - d[np.argmin(np.abs(g))] / d.max()
    return out, dip


print(f"{'mu/sigma':>9}{'mu^2-sigma^2':>14}{'봉우리 수':>10}{'봉우리 위치':>22}{'골의 깊이':>12}")
for mu in (0.9, 1.0, 1.1, 1.2, 1.5, 2.0, 3.0):
    pk, dip = peaks_and_dip(mu)
    print(f"{mu:>9.2f}{mu * mu - 1:>+14.4f}{len(pk):>10d}"
          f"{str([round(q, 3) for q in pk]):>22}{dip:>12.6f}")

# mu = 3 (위 코드의 분리 6) 에서 봉우리가 성분 평균에 얼마나 가까운가.
a = 3.0
dfdx = lambda x: -((x - a) * stats.norm.pdf(x - a) + (x + a) * stats.norm.pdf(x + a))
root = brentq(dfdx, 1.0, 5.0, xtol=1e-14)
print(f"\nmu = 3 의 오른쪽 봉우리 = {root:.9f}   (어림 mu - 2mu e^(-2mu^2) = "
      f"{a - 2 * a * np.exp(-2 * a * a):.9f})")

fig, (ax, az) = plt.subplots(1, 2, figsize=(12, 3.4))

for mu, c, ls in ((0.9, "#1565C0", "-"), (1.0, "#37474F", "--"), (1.1, "#E65100", "-")):
    g, d = mixture(mu)
    ax.plot(g, d, color=c, ls=ls, lw=1.8, label=f"$\\mu/\\sigma = {mu}$")
ax.set_xlim(-2.6, 2.6)
ax.set_xlabel("$x$")
ax.set_ylabel("밀도")
ax.set_title("임계값 근처의 혼합밀도 — $\\mu/\\sigma = 1$ 에서 갈린다", fontsize=11)
ax.legend(fontsize=9)

rs = np.linspace(0.5, 3.0, 120)
dips = [peaks_and_dip(r)[1] for r in rs]
az.plot(rs, dips, color="#1565C0", lw=1.8)
az.axvline(1.0, color="#D32F2F", ls="--", lw=1.4)
az.axhline(0.0, color="#90A4AE", lw=0.8)
for r in (1.1, 1.2, 1.5, 2.0, 3.0):
    az.plot([r], [peaks_and_dip(r)[1]], "o", ms=5, color="#E65100")
az.set_xlabel("$\\mu/\\sigma$")
az.set_ylabel("골의 깊이 $1 - f(0)/f_{\\max}$")
az.set_title("이봉이 된 뒤에도 골은 한참 얕다", fontsize=11)

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

plt.tight_layout()
plt.show()

 mu/sigma  mu^2-sigma^2     봉우리 수                봉우리 위치       골의 깊이
     0.90       -0.1900         1                 [0.0]    0.000000
     1.00       +0.0000         1                 [0.0]    0.000000
     1.10       +0.2100         2       [-0.737, 0.737]    0.025963
     1.20       +0.4400         2       [-1.001, 1.001]    0.089434
     1.50       +1.2500         2       [-1.463, 1.463]    0.358217
     2.00       +3.0000         2       [-1.999, 1.999]    0.729420
     3.00       +8.0000         2           [-3.0, 3.0]    0.977782

mu = 3 의 오른쪽 봉우리 = 2.999999909   (어림 mu - 2mu e^(-2mu^2) = 2.999999909)

import numpy as np
from scipy import stats

def analyse(s):
    a = s / 2
    f = lambda x: 0.5 * (stats.norm.pdf(x, -a) + stats.norm.pdf(x, a))
    g = np.linspace(-6, 6, 200_001)
    d = f(g)
    w = 50                                   # 평탄부에 속지 않도록 창으로 비교
    peaks = [g[i] for i in range(w, len(g) - w)
             if d[i] == d[i - w:i + w + 1].max() and d[i] > d[i - w] and d[i] > d[i + w]]
    merged = []
    for pk in peaks:
        if not merged or abs(pk - merged[-1]) > 0.05:
            merged.append(pk)
    return a * a - 1, merged

print(f"{'분리 s':>8}{'a^2 - 1':>11}{'봉우리 위치':>26}{'판정':>8}")
for s in (1.0, 1.5, 1.9, 1.99, 2.0, 2.01, 2.1, 2.5, 3.0):
    sign, peaks = analyse(s)
    print(f"{s:>8.2f}{sign:>+11.4f}{str([round(q, 3) for q in peaks]):>26}"
          f"{('이봉' if len(peaks) > 1 else '단봉'):>8}")

    분리 s    a^2 - 1                    봉우리 위치      판정
    1.00    -0.7500                     [0.0]      단봉
    1.50    -0.4375                     [0.0]      단봉
    1.90    -0.0975                     [0.0]      단봉
    1.99    -0.0100                     [0.0]      단봉
    2.00    +0.0000                    [-0.0]      단봉
    2.01    +0.0100           [-0.173, 0.173]      이봉
    2.10    +0.1025           [-0.534, 0.534]      이봉
    2.50    +0.5625               [-1.1, 1.1]      이봉
    3.00    +1.2500           [-1.463, 1.463]      이봉

import numpy as np
from scipy.stats import gaussian_kde

rng = np.random.default_rng(0)

def n_modes(x):
    g = np.linspace(x.min() - 1, x.max() + 1, 1000)
    d = gaussian_kde(x)(g)                      # 스콧 기본 대역폭
    return sum(1 for i in range(1, len(d) - 1) if d[i] > d[i - 1] and d[i] > d[i + 1])

R = 300
print(f"'봉우리 2개 이상'으로 판정할 확률")
print(f"{'n':>7}{'단봉 자료':>12}{'2.5 sigma':>12}{'3.0 sigma':>12}{'4.0 sigma':>12}")
for n in (50, 200, 1000, 5000):
    fp = sum(n_modes(rng.normal(0, 1, n)) >= 2 for _ in range(R)) / R
    row = f"{n:>7}{fp:>12.3f}"
    for s in (2.5, 3.0, 4.0):
        hits = sum(n_modes(np.concatenate([rng.normal(-s / 2, 1, n // 2),
                                           rng.normal(s / 2, 1, n // 2)])) >= 2
                   for _ in range(R))
        row += f"{hits / R:>12.3f}"
    print(row)

'봉우리 2개 이상'으로 판정할 확률
      n       단봉 자료   2.5 sigma   3.0 sigma   4.0 sigma
     50       0.220       0.463       0.733       0.993
    200       0.257       0.670       0.950       1.000
   1000       0.330       0.970       1.000       1.000
   5000       0.533       1.000       1.000       1.000

import numpy as np
from scipy.stats import gaussian_kde

rng = np.random.default_rng(4)

def n_modes(x, bw="scott"):
    g = np.linspace(np.min(x) - 1, np.max(x) + 1, 1500)
    d = gaussian_kde(x, bw_method=bw)(g)
    return sum(1 for i in range(1, len(d) - 1) if d[i] > d[i - 1] and d[i] > d[i + 1])

# (1) 두 시점·두 장비에서 잰 자료를 합쳤다 — 보정 오차가 있다
a = rng.normal(100, 3, 3000)                  # 장비 A
b = rng.normal(100, 3, 3000) + 7              # 장비 B: +7 만큼 치우쳐 보정됨
pooled = np.concatenate([a, b])
print(f"장비별로는 단봉: A {n_modes(a)}개, B {n_modes(b)}개")
print(f"합치면:          {n_modes(pooled)}개   ← 자료가 아니라 보정 오차")

# (2) 절단된 자료를 이어 붙였다 — 특정 구간이 기록되지 않았다
z = rng.normal(0, 1, 20_000)
censored = z[(z < -0.35) | (z > 0.35)]        # 중앙 구간이 누락됨
print(f"\n원자료: {n_modes(z)}개")
print(f"중앙 구간이 누락된 자료: {n_modes(censored)}개   "
      f"(관측 {len(censored)}/{len(z)}개만 남음)")

장비별로는 단봉: A 2개, B 1개
합치면:          2개   ← 자료가 아니라 보정 오차

원자료: 1개
중앙 구간이 누락된 자료: 2개   (관측 14563/20000개만 남음)

import numpy as np
from sklearn.mixture import GaussianMixture

rng = np.random.default_rng(5)

def fit_mixture(sep, n=4000, w=0.5):
    k = int(n * w)
    x = np.concatenate([rng.normal(0, 1, k), rng.normal(sep, 1, n - k)])
    g = GaussianMixture(2, random_state=0, n_init=5).fit(x.reshape(-1, 1))
    order = np.argsort(g.means_.ravel())
    return (g.means_.ravel()[order], np.sqrt(g.covariances_.ravel()[order]),
            g.weights_[order])

print(f"{'참 분리':>8}{'추정 평균':>22}{'추정 표준편차':>20}{'추정 가중치':>20}")
for sep in (4.0, 3.0, 2.0, 1.0):
    m, s, w = fit_mixture(sep)
    print(f"{sep:>8.1f}{str(np.round(m, 2)):>22}{str(np.round(s, 2)):>20}"
          f"{str(np.round(w, 2)):>20}")
print("\n참값: 평균 [0, 분리],  표준편차 [1, 1],  가중치 [0.5, 0.5]")

    참 분리                 추정 평균             추정 표준편차              추정 가중치
     4.0           [0.04 4.05]         [0.99 1.  ]           [0.5 0.5]
     3.0         [-0.01  3.01]         [1.01 1.01]           [0.5 0.5]
     2.0         [-0.06  2.  ]         [0.95 0.98]         [0.49 0.51]
     1.0         [-0.25  1.21]         [0.81 0.84]         [0.49 0.51]

참값: 평균 [0, 분리],  표준편차 [1, 1],  가중치 [0.5, 0.5]