Tạo mô hình hồi quy logistic ba cấp trong pymc3

Question

Tôi đang cố gắng tạo mô hình hồi quy logistic ba cấp trong pymc3. Có một mức cao nhất, mức giữa và mức cá nhân, trong đó các hệ số ở mức trung bình được ước tính từ các hệ số cấp cao nhất. Tôi đang gặp khó khăn trong việc xác định cấu trúc dữ liệu thích hợp cho mức trung bình, tuy nhiên.

Dưới đây là mã của tôi:

with pm.Model() as model: 
    # Hyperpriors 
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.) 
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)  

    # Priors 
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top) 
    mid_level = [pm.Normal('mid_level_{}'.format(j), 
          mu=top_level[mid_to_top_idx[j]], 
          tau=mid_level_tau) 
       for j in range(k_mid)] 

    intercept = pm.Normal('intercept', mu=0., sd=100.) 

    # Model prediction 
    yhat = pm.invlogit(mid_level[mid_to_bot_idx] + intercept) 

    # Likelihood 
    yact = pm.Bernoulli('yact', p=yhat, observed=y) 

Tôi nhận được lỗi "only integer arrays with one element can be converted to an index" (trên đường 16), mà tôi nghĩ có liên quan đến thực tế là biến mid_level là một danh sách, không phải là một pymc container thích hợp. (Tôi cũng không thấy lớp Container trong mã nguồn pymc3).

Bất kỳ trợ giúp nào cũng sẽ được đánh giá cao.

Chỉnh sửa: Thêm một số dữ liệu giả

y = np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0]) 
mid_to_bot_idx = np.array([0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 3, 2, 3, 3, 3, 3, 2, 2, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 2, 2, 3, 2, 2, 3, 3, 3, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 2]) 
mid_to_top_idx = np.array([0, 0, 1, 1]) 
k_top = 2 
k_mid = 4 

Chỉnh sửa # 2:

Dường như có một vài cách khác nhau để giải quyết vấn đề này, mặc dù không hoàn toàn thỏa đáng:

1) Người ta có thể điều chỉnh lại mô hình là:

with pm.Model() as model: 
    # Hyperpriors 
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.) 
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)  

    # Priors 
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top) 
    mid_level = pm.Normal('mid_level', mu=0., tau=mid_level_tau, shape=k_top) 
    intercept = pm.Normal('intercept', mu=0., sd=100.) 

    # Model prediction 
    yhat = pm.invlogit(top_level[top_to_bot_idx] + mid_level[mid_to_bot_idx] + intercept) 

    # Likelihood 
    yact = pm.Bernoulli('yact', p=yhat, observed=y) 

Điều này dường như hoạt động, mặc dù Tôi không thể tìm ra cách mở rộng nó đến trường hợp phương sai mức trung bình không phải là hằng số cho tất cả các nhóm cấp trung bình.

2) Người ta có thể quấn các hệ số giữa mức độ thành một tensor Theano sử dụng theano.tensor.stack: ví dụ,

import theano.tensor as tt 
mid_level = tt.stack([pm.Normal('mid_level_{}'.format(j), 
          mu=top_level[mid_to_top_idx[j]], 
          tau=mid_level_tau) 
       for j in range(k_mid)]) 

Nhưng điều này có vẻ chạy rất chậm trên bộ dữ liệu thực tế của tôi (30k quan sát) , và nó làm cho âm mưu bất tiện (mỗi hệ số mid_level có dấu vết riêng của nó bằng cách sử dụng pm.traceplot).

Dù sao, một số lời khuyên/đầu vào từ các nhà phát triển sẽ được đánh giá cao.

Answer 1

Hóa ra giải pháp rất đơn giản: nó xuất hiện rằng bất kỳ phân phối (như pm.Normal) có thể chấp nhận một vector của các phương tiện như một cuộc tranh cãi, vì vậy thay thế dòng này

mid_level = [pm.Normal('mid_level_{}'.format(j), 
         mu=top_level[mid_to_top_idx[j]], 
         tau=mid_level_tau) 
      for j in range(k_mid)] 

với điều này

mid_level = pm.Normal('mid_level', 
         mu=top_level[mid_to_top_idx], 
         tau=mid_level_tau, 
         shape=k_mid) 

công trình. Phương pháp tương tự cũng có thể được sử dụng để xác định độ lệch chuẩn riêng lẻ cho từng nhóm giữa cấp.

CHỈNH SỬA: Sửa lỗi typo

Answer 2

ít thay đổi (lưu ý tôi đã thay đổi yhat để theta):

theta = pm.Deterministic('theta', pm.invlogit(sum(mid_level[i] for i in mid_to_bot_idx)+intercept)) 
yact = pm.Bernoulli('yact', p=theta, observed=y)