pymc와statsmodels로 정적 분포 파라미터를 추정하는 노트

3295 단어 statsmodelspymcPython
이 글은 정적 분포의 매개 변수를 추산한다statsmodels라는 두 가지 방법을 제시했다.
학생은 필기를 대신해서 쓴 것이니 무슨 잘못이 있으면 지적해 주십시오.

데이터 세트 만들기


어쨌든 이번에 준비한 것은 간단하고 소음이 없는 데이터 집합(정규 분포된 소음을 넣으면 전체적으로 특정한 정규 분포에 따라 분포되기 때문)이다.
샘플 수량이 증가했기 때문에 샘플 수량은 10010000000이다.
import numpy as np

MU = 10.0
SIGMA = 1.0

y0 = np.random.normal(MU, SIGMA, size=100)  
y1 = np.random.normal(MU, SIGMA, size=1000) 
y2 = np.random.normal(MU, SIGMA, size=10000) 
y = [y1, y2, y3]
얻은 데이터를 다음과 같이 그려라.

Y2가 되면 정규 분포에 가까워진다.

statsmodels에 의한 추정


그럼 statsmodels로 추측해 봅시다.
import statsmodels.api as sm

for i in range(3):
    size = 10**(i + 2)
    x = [1 for j in range(size)]
    glm = sm.GLM(np.array(y[i]), np.array(x), family=sm.families.Gaussian())
    res = glm.fit()
    print(res.summary())
결과는 다음과 같다.
샘플 번호
균등하다
편차
95% 신뢰 구간(하한)
95% 신뢰 구간(상한)
0
9.936
0.874
9.752
10.119
1
9.981
1.010
9.919
10.043
2
9.995
1.016
9.976
10.015
샘플 물량이 늘면서 신뢰 구간도 95% 축소됐다.

Pymc의 추측에 의하면


pymc3


드디어 Pymc로 추정하고 싶습니다.
코드가 아까와 다르니 샘플에 대해 실행하니 주의하세요.
또한pymc는tau를 입력했기 때문에 분산된 값이 아니라는 것을 주의하십시오.
import statsmodels.api as sm
from scipy import optimize
import pymc3 as pm3

with pm3.Model() as model:
    # mu, tauの事前分布は一様分布
    mu = pm3.Uniform("mu", upper= 10**2, lower= -(10**2))
    tau = pm3.Uniform("tau", upper= 10**2, lower= 0)
    Y_obs = pm3.Normal("Y_obs", mu=mu, sd=tau, observed=np.array(y3)) 

with model:
    # MAP推定量の計算
    start = pm3.find_MAP(fmin=optimize.fmin_powell)
    # メトロポリス法のサンプリング
    # NUTS, Metropolis, Slice, HamiltonianMCなどがある
    step = pm3.Metropolis()
    # 5000回サンプリング
    trace = pm3.sample(5000, start=start, step=step)
    # 2000個を捨てる
    pm3.traceplot(trace[2000:])
    plt.show()
    pm3.summary(trace[2000:])
샘플 번호
mu
tau
95% 신용 구간(하한선,mu,tau순)
95% 신용 구간(상한선,mu,tau순)
0
9.931
0.948
9.742, 0.811
10.091, 1.080
1
9.979
1.007
9.915, 0.966
10.040, 1.056
2
9.995
1.008
9.975, 0.993
10.016, 1.022
또한 y2에 대한 샘플링 결과를 게재합니다.

pymc2


pymc2 코드도 넣고.
# モデルの作成
mu = pm.Uniform("mu", upper= 10**2, lower= -(10**2))
tau = pm.Uniform("tau", upper= 10**2, lower= 0)
obs = pm.Normal("obs", mu=mu, tau=tau, value=y3, observed=True)
model = pm.Model([obs, mu, tau])
# MAP推定量の計算
map_ = pm.MAP(model)
map_.fit()
# MCMCの計算
mcmc = pm.MCMC(model)
mcmc.sample(5000, 2000)

# 可視化
pm.Matplot.plot(mcmc.trace("mu"), common_scale=False)
pm.Matplot.plot(mcmc.trace("tau"), common_scale=False)
plt.show()

끝말


주파수론적 추정과 베이시즘적 추정 기법을 시도했다.
statsmodels와pymc 둘 다 처리하는 페이지가 많지 않은 것 같아서 필기를 겸했습니다. MCMC의 각 샘플링 기법에 대한 내용을 잘 모르기 때문에 앞으로 배울 것이 있습니다.

좋은 웹페이지 즐겨찾기