使用確定性分佈的PyMC3代碼中的AssertionError

我一直試圖在過去的日子裏習慣於PyMC最終從我有直接代碼（我最終對參數估計感興趣）的某些模型中進行一些MCMC分佈採樣。使用確定性分佈的PyMC3代碼中的AssertionError

據我所知，沒有那麼多的例子顯示他們的代碼（如外部C或FORTRAN代碼），他們成功地使用PyMC3工作。到目前爲止，我發現了here或here賄賂。因此，從簡單的問題開始，我嘗試用PyMC3複製現有Python代碼中的一些結果，這些結果來自於使用PyMC的「複雜」（閱讀：比doc示例多）發行版的MCMC，即運行於此的簡單result我的筆記本電腦上10000個樣本需要2秒鐘。

要定義PyMC3中的確定性變量，應該使用@theano.compile.ops.as_op裝飾器（它在PyMC中是@deterministic，現在在PyMC3中已棄用），這就是我所做的。

這是我的代碼（I使用的Spyder等IPython的），由PyMC documentation，其中我遇到AssertionError第二小區之後激發（在估計過程，所以採樣之前）。我一直在試圖解決這個過去兩天，但我真的不明白錯誤可能是什麼。我相信它應該是一種PyMC或Theano的伎倆，我還沒有趕上，因爲我相信我真的很接近。

#%% Define model 
import numpy,math 
import matplotlib.pyplot as plt 
import random as random 

import theano.tensor as t 
import theano 

random.seed(1) # set random seed 

# copy-pasted function from the specific model used in the source and adapted with as_op 
@theano.compile.ops.as_op(itypes=[t.iscalar, t.dscalar, t.fscalar, t.fscalar],otypes=[t.dvector]) 
def sampleFromSalpeter(N, alpha, M_min, M_max):   
    log_M_Min = math.log(M_min) 
    log_M_Max = math.log(M_max) 
    maxlik = math.pow(M_min, 1.0 - alpha) 
    Masses = [] 
    while (len(Masses) < N):    
     logM = random.uniform(log_M_Min,log_M_Max) 
     M = math.exp(logM)    
     likelihood = math.pow(M, 1.0 - alpha)    
     u = random.uniform(0.0,maxlik) 
     if (u < likelihood): 
      Masses.append(M) 
    return Masses 

# SAME function as above, used to make test data (so no Theano here) 
def sampleFromSalpeterDATA(N, alpha, M_min, M_max):   
    log_M_Min = math.log(M_min) 
    log_M_Max = math.log(M_max)   
    maxlik = math.pow(M_min, 1.0 - alpha)   
    Masses = []   
    while (len(Masses) < N):    
     logM = random.uniform(log_M_Min,log_M_Max) 
     M = math.exp(logM)    
     likelihood = math.pow(M, 1.0 - alpha)    
     u = random.uniform(0.0,maxlik) 
     if (u < likelihood): 
      Masses.append(M) 
    return Masses 

# Generate toy data. 
N  = 1000000 # Draw 1 Million stellar masses. 
alpha = 2.35 
M_min = 1.0 
M_max = 100.0 
Masses = sampleFromSalpeterDATA(N, alpha, M_min, M_max) 

#%% Estimation process 
import pymc3 as pm 
basic_model = pm.Model() 
with basic_model: 

    # Priors for unknown model parameters 
    alpha2 = pm.Normal('alpha2', mu=3, sd=10) 

    N2=t.constant(1000000) 
    M_min2 = t.constant(1.0) 
    M_max2 = t.constant(100.0) 

    # Expected value of outcome 
    m = sampleFromSalpeter(N2, alpha2, M_min2, M_max2) 

    # Likelihood (sampling distribution) of observations 
    Y_obs = pm.Normal('Y_obs', mu=m, sd=10, observed=Masses) 

#%% Sample 
with basic_model: 
    step = pm.Metropolis() 
    trace = pm.sample(10000, step=step) 

#%% Plot posteriors 
_ = pm.traceplot(trace) 
pm.summary(trace)

來源

2017-09-13 viiv

我剛剛編輯了一個錯字：估計_process_而不是估計_problem_。問題仍然有效 – viiv

我在Discourse of PyMC得到了答案。我確實非常接近，因爲我不得不改變函數定義的最後一行來返回[numpy.array（Masses）]。

現在，這段代碼似乎需要大約300小時才能運行（而不是在沒有PyMC的情況下在原始實現中需要2秒），並且我想知道是否有人知道爲什麼它會如此慢？我實際上已經嘗試過很簡單的發行版，並且使用as_op，並且注意到一個相當大的放緩，但甚至沒有那麼多。我想原始代碼直接在logp上工作，它是什麼使整體差異？

來源

2017-09-15 07:43:20 viiv

使用確定性分佈的PyMC3代碼中的AssertionError

回答

相關問題