Saturation curve — 광고비를 더 부어도 안 늘어나는 지점 찾기

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def hill(x, L, k):
    """Hill saturation 함수"""
    return x**k / (x**k + L**k)

# 채널별 주간 지출과 관측 기여도 (MMM 회귀 후 채널 기여도 시계열)
spend_meta   = np.array([20e6, 40e6, 60e6, 80e6, 100e6, 120e6, 150e6])
contrib_meta = np.array([0.18, 0.32, 0.43, 0.51, 0.57,  0.61,  0.65])

# 초기값과 bounds 설정
p0     = [80e6, 1.5]                          # L 초기값, k 초기값
bounds = ([10e6, 0.5], [500e6, 5.0])          # (L_min, k_min), (L_max, k_max)

popt, pcov = curve_fit(hill, spend_meta, contrib_meta,
                        p0=p0, bounds=bounds, maxfev=5000)
L_fit, k_fit = popt
L_std, k_std = np.sqrt(np.diag(pcov))

print(f"절반 포화 지점 L = {L_fit/1e6:.1f}M ± {L_std/1e6:.1f}M")
print(f"곡선 가파름     k = {k_fit:.2f} ± {k_std:.2f}")

x_plot = np.linspace(5e6, 200e6, 300)
plt.plot(x_plot / 1e6, hill(x_plot, L_fit, k_fit), label="Hill fit")
plt.scatter(spend_meta / 1e6, contrib_meta, color='red', zorder=5, label="관측값")
plt.axvline(L_fit / 1e6, color='gray', ls='--', label=f"L = {L_fit/1e6:.0f}M")
plt.xlabel("주간 지출 (백만 원)"); plt.ylabel("기여도 (정규화)")
plt.title("Meta 채널 Hill Saturation Fit"); plt.legend(); plt.show()

def hill_deriv(x, L, k):
    """Hill 함수의 도함수 = 한계 saturation 기울기"""
    return k * L**k * x**(k - 1) / (x**k + L**k)**2

# 채널별 파라미터 (posterior 평균 또는 curve_fit 결과)
channels = ['Meta', 'Naver', 'TV']
beta  = np.array([2.5, 1.8, 3.2])       # 회귀 계수
L_arr = np.array([80e6, 40e6, 150e6])   # 절반 포화 지점
k_arr = np.array([1.5, 1.2, 2.0])       # 곡선 가파름

# 현재 분기 실제 지출
current_spend = np.array([95e6, 35e6, 140e6])

mroas = beta * hill_deriv(current_spend, L_arr, k_arr)
avg_roas_proxy = beta * hill(current_spend, L_arr, k_arr) / (current_spend / 1e8)

print(f"{'채널':<8} {'현재지출':>10} {'평균ROAS':>10} {'한계ROAS':>10} {'판단':>12}")
for i, ch in enumerate(channels):
    pos = "L 우측(줄여야)" if current_spend[i] > L_arr[i] else "L 좌측(더 가능)"
    print(f"{ch:<8} {current_spend[i]/1e6:>8.0f}M {avg_roas_proxy[i]:>10.2f}x "
          f"{mroas[i]:>10.2f}x {pos:>14}")

채널     현재지출    평균ROAS    한계ROAS          판단
Meta       95M       3.45x       1.28x    L 우측(줄여야)
Naver      35M       2.91x       3.47x    L 좌측(더 가능)
TV        140M       2.18x       0.91x    L 우측(줄여야)

import pymc as pm
import arviz as az

# spend_meta, contrib_meta: 위와 동일

with pm.Model() as saturation_model:
    # L prior: 현재 지출 평균의 1~5배 사이
    L = pm.HalfNormal("L", sigma=150e6)
    # k prior: 곡선 가파름 (0.5~4 사이)
    k = pm.Gamma("k", alpha=2, beta=1)   # 평균=2, 분산=2

    # Hill saturation
    mu = (spend_meta**k) / (spend_meta**k + L**k)

    # 관측 오차
    sigma_obs = pm.HalfNormal("sigma_obs", sigma=0.05)
    obs = pm.Normal("obs", mu=mu, sigma=sigma_obs,
                    observed=contrib_meta)

    idata_sat = pm.sample(1000, tune=500, target_accept=0.9,
                          random_seed=42, return_inferencedata=True)

# 수렴 진단
summary = az.summary(idata_sat, var_names=["L", "k"])
print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]])

# posterior 곡선 시각화
x_plot = np.linspace(5e6, 200e6, 300)
posterior_L = idata_sat.posterior["L"].values.flatten()
posterior_k = idata_sat.posterior["k"].values.flatten()

curves = np.array([
    hill(x_plot, l, kk) for l, kk in zip(posterior_L[:500], posterior_k[:500])
])

plt.figure(figsize=(9, 5))
plt.fill_between(x_plot / 1e6,
                 np.percentile(curves, 5, axis=0),
                 np.percentile(curves, 95, axis=0),
                 alpha=0.25, color='steelblue', label='90% 신용구간')
plt.plot(x_plot / 1e6, np.median(curves, axis=0),
         color='steelblue', lw=2, label='posterior 중앙값')
plt.scatter(spend_meta / 1e6, contrib_meta, color='red', zorder=5, label='관측값')
plt.xlabel("주간 지출 (백만 원)"); plt.ylabel("기여도 (정규화)")
plt.title("Hill Saturation — PyMC posterior 90% 신용구간")
plt.legend(); plt.tight_layout(); plt.show()
# → 신용구간이 넓은 지출 구간 = 데이터 부족, 의사결정 신중히