0
을 pymc하기 은하계의 종류 인 경우, LF는 감마 함수 일뿐입니다. 은하 조사가 대개 목표의 상당 부분, 특히 낮은 광도의 자료를 놓치기 때문에 데이터 절단의 잠재적 인 원인이됩니다. 이 모델에서 나는 아래 모든 것을 놓친다 lmin추가 측정 오차 내가 pymc2에서 다음 모델 한 모델
이 방법의 세부 사항은 this paper by Kelly et al에서 찾을 수있다.
이 모델
작동 : 나는 모델에MAP 및
MCMC을 실행하고
lmin가 성장함에 따라 내가 증가 불확실성, 매개 변수 내 시뮬레이션 데이터
sample에서
alpha 및
scale을 복구 할 수 있습니다.
이제 가우스 측정 오류를 삽입하고 싶습니다. 단순화를 위해 모든 데이터는 동일한 정밀도를가집니다. 나는 또한 오류를 포함시킬 가능성을 수정하지는 않습니다.
alpha = pymc.Uniform('alpha', 0.01, 2.0)
scale = pymc.Uniform('scale',1.0, 4.0)
sig = 0.1
tau = math.pow(sig, -2.0)
@pymc.deterministic(plot=False)
def beta(scale=scale):
return 1.0/scale
@pymc.potential
def p_factor(alpha=alpha, scale=scale, lmin=lmin, n=len(sample)):
dist = gamma(alpha, loc=0., scale=scale)
fp = 1.0 - dist.cdf(lmin)
return -(n+1) * np.log(fp)
dist = pymc.Gamma("dist", alpha=alpha, beta=beta)
obs = pymc.Normal("obs", mu=dist, tau=tau, value=sample, observed=True)
하지만이 모델이 작동하지 않기 때문에 분명히 여기서 뭔가 잘못하고 있습니다. 나는이 모델에 pymc.MAP을 실행하면 내가 pymc.MCMC, alpha 및 beta이 더 전혀 추적하는 실행되지 때 나는 초기 alpha의 가치와 scale
vals = {'alpha': alpha, 'scale': scale, 'beta': beta,
'p_factor': p_factor, 'obs': obs, 'dist': dist}
M2 = pymc.MAP(vals)
M2.fit()
print M2.alpha.value, M2.scale.value
>>> (array(0.010000000006018368), array(1.000000000833973))
을 복구 할 수 있습니다.
M = pymc.MCMC(vals)
M.sample(10000, burn=5000)
...
M.stats()['alpha']
>>> {'95% HPD interval': array([ 0.01000001, 0.01000502]),
'mc error': 2.1442678276712383e-07,
'mean': 0.010001588137798096,
'n': 5000,
'quantiles': {2.5: 0.0100000088679046,
25: 0.010000382359859467,
50: 0.010001100377476166,
75: 0.010001668672799679,
97.5: 0.0100050194240779},
'standard deviation': 2.189828287191421e-06}
다시 초기 값. 실제로 alpha을 0.02로 시작하면 alpha의 복구 값은 0.02입니다.
이것은 a notebook with the working model plus simulated data입니다.
이것은 a notebook with the error model plus simulated data입니다.
이 작업에 대한 지침은 정말 감사하겠습니다.