PyMC3에서 두 개의 시프트 된 감마 분포의 혼합 모델을 맞 춥니 다.

Question

저는 파이썬과 MCMC 기법에 익숙하지 않으며 PyMC3에 대해 작업하고 있습니다. PyMC3에 익숙해지기 위해 생성 된 데이터에 두 개의 시프트 된 감마 분포의 혼합 모델을 맞추고 싶습니다. 다음 단계로 나는 이것을 임의의 수의 이동 된 감마 (gamma)로 확장하는 것을 좋아할 것이다. 그러나 한 번에 한 걸음 씩 나아 간다. 전체 코드 추가 링크

포인터는 다음 PyMC3 GitHub의 문제에 대한 내 댓글에서 찾을 수 있습니다 여기에있다

with pm.Model() as model: 
    p = pm.Beta("p", 1., 1., testval=0.5) #this is the fraction that come from mean1 vs mean2 

    alpha = pm.Uniform('alpha', 0, 10) # == k 
    beta = pm.Uniform('beta' , 0, 1) # == theta 
    shifts = pm.Normal('shifts', mu=20, sd=50, shape=2) 

    gamma = ShiftedGamma.dist(alpha, beta, c=shifts) 
    process = pm.Mixture('obs', tt.stack([p, 1-p]), gamma, observed=data) 

with model: 
    step = pm.Metropolis() 
    trace = pm.sample(10000, step=step) 

을 그리고 다음은 https://github.com/pymc-devs/pymc3/issues/864#issuecomment-285864556

내가 구현하려고 모델입니다 내가 해낸 ShiftedGamma의 구현 :

class ShiftedLog(pm.distributions.continuous.transforms.ElemwiseTransform): 
    name = "log" 
    def __init__(self, c=0.0): 
     self.c = c 

    def backward(self, x): 
     return tt.exp(x + self.c) 

    def forward(self, x): 
     return tt.log(x - self.c) 

class ShiftedGamma(pm.distributions.continuous.Continuous): 
    def __init__(self, alpha=None, beta=None, mu=None, sd=None, c=0.0, *args, **kwargs): 
     super(ShiftedGamma, self).__init__(transform=ShiftedLog(c), *args, **kwargs) 
     alpha, beta = self.get_alpha_beta(alpha, beta, mu, sd) 
     self.alpha = alpha 
     self.beta = beta 
     self.mean = alpha/beta + c 
     self.mode = tt.maximum((alpha - 1)/beta, 0) + c 
     self.variance = alpha/beta**2 
     self.c = c 
     # self.testval = kwargs.pop('testval', None) 
     # if not self.testval: 
     #  self.testval = c + 10.0 

     pm.distributions.continuous.assert_negative_support(alpha, 'alpha', 'Gamma') 
     pm.distributions.continuous.assert_negative_support(beta, 'beta', 'Gamma') 

    def get_alpha_beta(self, alpha=None, beta=None, mu=None, sd=None): 
     if (alpha is not None) and (beta is not None): 
      pass 
     elif (mu is not None) and (sd is not None): 
      alpha = mu**2/sd**2 
      beta = mu/sd**2 
     else: 
      raise ValueError('Incompatible parameterization. Either use ' 
          'alpha and beta, or mu and sd to specify ' 
          'distribution.') 

     return alpha, beta 

    def random(self, point=None, size=None, repeat=None): 
     alpha, beta, c = draw_values([self.alpha, self.beta, self.c], point=point) 
     return generate_samples(stats.gamma.rvs, alpha, scale=1./beta, dist_shape=self.shape, size=size, loc=c) 

    def logp(self, value): 
     alpha = self.alpha 
     beta = self.beta 
     c = self.c 
     x = value - c 
     return bound(-gammaln(alpha) + logpow(beta, alpha) - beta * x + logpow(x, alpha - 1), x >= 0, alpha > 0, beta > 0) 

는 기본적으로 내가이 질문이 :

1) 구현이 "좋음"/ "정확합니까?" 나는 특히 ShiftedLog 변환을 구현 한 변환 매개 변수에 대해 잘 모릅니다. 그리고 임의의 장소에서 random() 메서드를 사용합니까? 여기서 중단 점을 설정하면 절대로 호출되지 않습니다.

2)이 모델은 기본적으로 내가 사전 또는 샘플링 방법에 변경 한 사항에 비해 매우 취약합니다. PyMC3의 결과가 나올 때까지 요지에서 많은 의견을보고 실험을 반복 해보십시오.

포인트 2)는 포인트 1에서 실수를 저질렀거나 StanCon에서 Michael Betancourt가 말한 마술 사유 중 하나 일 수 있습니다. "Markov Chain Monte에 대해 알아야 할 모든 것 카를로 "?

내 질문 1) 데이터에 시프트 된 gammas가있는 모델을 피팅하는 것과 동일한 효과를 얻을 수있는 더 좋은 방법이나 쉬운 방법이 있는지 알고 싶습니다. ShiftedGamma를 구현하는 것이 실제로 필요한지 아니면 pm.Deterministic() 또는 이와 유사한 방법으로 동일한 효과를 얻을 수 있습니까?

도움과 통찰력을 가져 주셔서 감사합니다.

내가 원래 GitHub의 문제에 게시

Answer 1

기독교

,하지만 난이 더 좋은 장소 같아요.

변환에 오타가있는 것 같습니다. 역순으로 tt.exp (x) + self.c가되어서는 안됩니까? 이 경우 pm.distributions.transforms.lowerbound와 동일합니다. 그렇지 않으면 접근 방법이 잘된 것 같습니다. 혼합물이 분배를 필요로하기 때문에 나는 deteministics로 이것을하는 어떤 방법을 모른다. 거의 모든 다른 상황에서 ShiftedGamma는 필요하지 않습니다.

난 당신이 뭔가에 다른 코드를 단순화 수 있다고 생각 : 초기화 및 샘플링 어려움에 대해

class ShiftedGamma(pm.Gamma): 
    def __init__(self, alpha, beta, shift, *args, **kwargs): 
     transform = pm.distributions.transforms.lowerbound(shift) 
     super().__init__(alpha=alpha, beta=beta, *args, **kwargs, 
         transform=transform) 
     self.shift = shift 
     self.mean += shift 
     self.mode += shift 

    def random(self): 
     return super().random() + self.shift 

    def logp(self, x): 
     return super().logp(x - self.shift) 


with pm.Model() as model: 
    p = pm.Beta("p", 1., 1., testval=0.5) 

    alpha = pm.Uniform('alpha', 0, 10) # == k 
    beta = pm.Uniform('beta' , 0, 1) # == theta 
    shifts = pm.Normal('shifts', mu=20, sd=50, shape=2, testval=floatX([0., 0.])) 

    gamma = ShiftedGamma.dist(alpha=alpha, beta=beta, shift=shifts) 

    pm.Mixture('y', tt.stack([p, 1-p]), gamma, observed=data)) 

(... 내가 이것을 제대로 테스트하지 않았다, 다시 한 번 확인하시기 바랍니다) : 당신은 대부분 MAPI를 사용하면 초기화에 advi를 사용하는 것이 일반적으로 더 좋습니다. 당신이 그런 어려움을 겪는 이유는 많은 shift 값들이 -inf의 logp가되는 것, 즉 교대 최소값이 데이터 셋의 최소값보다 큰 경우라고 생각합니다. 어떻게 든 다시 매개 변수화를 시도하여 이러한 일이 발생하지 않도록 할 수 있습니다.아마도 다음과 같은 것일 수 있습니다 :

shift1_raw = pm.HalfCauchy('shift1_raw', beta=1) 
shift2_raw = pm.HalfCauchy('shift2_raw', beta=1) 
shift1 = data.min() - shift1_raw 
shift2 = shift1 + shift2_raw 
shift = tt.stack([shift1, shift2]) 
shift = pm.Deterministic('shift', shift) 

그러나 이것은 분명히 더 많은 생각을 필요로합니다. 이 버전은 이전 버전을 변경하고이 새로운 이전 버전은 데이터 세트에 따라 달라지며 shift1_raw 및 shift2_raw는 아마도 상관 관계가 높아 샘플러를 불만족스럽게 만듭니다.

또한 혼합 모델은 일반적으로 멀티 모달로 끝나며 NUTS는 일반적으로 잘 작동하지 않습니다.

편집 pymc3에 대한 전반적인 이해를 원하고 Mixture 모델에 특별한 관심이 없다면 그것들로 시작하지 않을 수도 있습니다. 그들은 문제가되는 경향이 있습니다.