Model comparison issue: "Found several log likelihood arrays var_name cannot be None"

[mservantufc]

Hi community, I have two models that I would like to compare using LOO and WAIC.

Model 1:

with pm.Model() as Model_1:
    α_v = pm.HalfCauchy('α_v',2)
    α_a = pm.HalfCauchy('α_a', 2) 
    α_max_ampl = pm.HalfCauchy('α_max_ampl',2)  
    α_tau = pm.HalfCauchy('α_tau', 2)     
    α_ter = pm.HalfCauchy('α_ter', 2)  

    β = pm.Normal('β', 0.334, 0.1)

    γ_v = pm.Normal('γ_v', params_group[-1,0], .1)
    γ_a = pm.Normal('γ_a', params_group[-1,1], .1) 
    γ_ter = pm.Normal('γ_ter', params_group[-1, 2], .1)     
    γ_max_ampl = pm.Normal('γ_max_ampl', params_group[-1,3], .1)
    γ_tau = pm.Normal('γ_tau', params_group[-1,4], .1)     

    μ_v = -α_v*pm.math.exp(-β*Age_years_group) + γ_v
    μ_a = α_a*pm.math.exp(-β*Age_years_group) + γ_a      
    μ_max_ampl = α_max_ampl*pm.math.exp(-β*Age_years_group) + γ_max_ampl
    μ_tau = α_tau*pm.math.exp(-β*Age_years_group) + γ_tau 
    μ_ter = α_ter*pm.math.exp(-β*Age_years_group) + γ_ter     

    σ_v = pm.HalfNormal('σ_v', .1)
    σ_a = pm.HalfNormal('σ_a', .1)
    σ_max_ampl = pm.HalfNormal('σ_max_ampl', .1)
    σ_tau = pm.HalfNormal('σ_tau', .1)      
    σ_ter = pm.HalfNormal('σ_ter', .1)     

    y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_group[:, 0])    
    y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_group[:, 1])     
    y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_group[:, 2]) 
    y_max_ampl = pm.Normal('y_max_ampl', μ_max_ampl, σ_max_ampl, observed=params_group[:, 3])     
    y_tau = pm.Normal('y_tau', μ_tau, σ_tau, observed=params_group[:, 4])     

Model_2:

with pm.Model() as Model_2:
    α_v = pm.HalfCauchy('α_v',2)
    α_a = pm.HalfCauchy('α_a', 2) 
    α_max_ampl = pm.HalfCauchy('α_max_ampl',2)  
    α_tau = pm.HalfCauchy('α_tau', 2)     
    α_ter = pm.HalfCauchy('α_ter', 2)  

    β_v = pm.Normal('β_v', 0.334, 0.1)#Kail (1991)
    β_a = pm.Normal('β_a', 0.334, 0.1)#Kail (1991)
    β_max_ampl = pm.Normal('β_max_ampl', 0.334, 0.1)#Kail (1991)
    β_tau = pm.Normal('β_tau', 0.334, 0.1)#Kail (1991)
    β_ter = pm.Normal('β_ter', 0.334, 0.1)#Kail (1991)

    #asymptote (get the mean value of adults)
    γ_v = pm.Normal('γ_v', params_group[-1,0], .1)
    γ_a = pm.Normal('γ_a', params_group[-1,1], .1) 
    γ_ter = pm.Normal('γ_ter', params_group[-1, 2], .1)     
    γ_max_ampl = pm.Normal('γ_max_ampl', params_group[-1,3], .1)
    γ_tau = pm.Normal('γ_tau', params_group[-1,4], .1)     

    μ_v = -α_v*pm.math.exp(-β_v*Age_years_group) + γ_v
    μ_a = α_a*pm.math.exp(-β_a*Age_years_group) + γ_a      
    μ_max_ampl = α_max_ampl*pm.math.exp(-β_max_ampl*Age_years_group) + γ_max_ampl
    μ_tau = α_tau*pm.math.exp(-β_tau*Age_years_group) + γ_tau 
    μ_ter = α_ter*pm.math.exp(-β_ter*Age_years_group) + γ_ter     

    σ_v = pm.HalfNormal('σ_v', .1)
    σ_a = pm.HalfNormal('σ_a', .1)
    σ_max_ampl = pm.HalfNormal('σ_max_ampl', .1)
    σ_tau = pm.HalfNormal('σ_tau', .1)      
    σ_ter = pm.HalfNormal('σ_ter', .1)     

    y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_group[:, 0])    
    y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_group[:, 1])     
    y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_group[:, 2]) 
    y_max_ampl = pm.Normal('y_max_ampl', μ_max_ampl, σ_max_ampl, observed=params_group[:, 3])     
    y_tau = pm.Normal('y_tau', μ_tau, σ_tau, observed=params_group[:, 4])           

I thus tried to compute:

compare_dict = {Model_1: trace, Model_2: trace_2}
compare_LOO = az.compare(compare_dict) 

and got the following error:

Traceback (most recent call last):

  File "<ipython-input-35-26e3f30f8929>", line 1, in <module>
    compare_LOO = az.compare(compare_dict)

  File "C:\Users\mservant\Anaconda3\lib\site-packages\arviz\stats\stats.py", line 211, in compare
    ics = ics.append([ic_func(dataset, pointwise=True, scale=scale)])

  File "C:\Users\mservant\Anaconda3\lib\site-packages\arviz\stats\stats.py", line 646, in loo
    log_likelihood = _get_log_likelihood(inference_data, var_name=var_name)

  File "C:\Users\mservant\Anaconda3\lib\site-packages\arviz\stats\stats_utils.py", line 412, in get_log_likelihood
    "Found several log likelihood arrays {}, var_name cannot be None".format(var_names)

TypeError: Found several log likelihood arrays ['y_v', 'y_a', 'y_ter', 'y_max_ampl', 'y_tau'], var_name cannot be None

Could you please help me solve this issue ? I am confused because I could compute LOO and WAIC using a previous version of arviz.

[OriolAbril]

Hi, when there are several likelihoods present there are generally several technically correct ways of calculating loo/waic, and choosing the right one considering the problem and modeling question at hand can't be done automatically.

We have an in progress guide on ArviZ resources (link is to my fork with the open PR) about how to tell ArviZ how to calculate the desired loo/waic.

I hope the linked notebook is clear enough, but don't hesitate to ask any question! I'll try to answer more clearly and improve the notebook at the same time with the feedback received :smile:

[mservantufc]

Hi, Thank you very much for your help. Your notebook is very clear but I am still a little confused about model selection, and the best strategy to solve the question I am trying to answer.. Let me briefly explain what I would like to do. Each of the 5 subplots below represents the development of a particular brain process as a function of age (fits are from Model_1 explained below):

The development of each process has traditionally been modeled using an exponential function of the following form: α*pm.math.exp(-β*Age) + γ My question is very simple: I would like to determine whether the 5 brain processes develop at the same rate (β in the above formula) or not. What I did so far : (1) Fit Model_1 to data (one β parameter for the 5 variables) (2) Fit Model_2 to data (one β parameter for each of the 5 variables) Model_1 and Model_2 are specified in my previous post. I then wanted to compare the 2 models using LOO or WAIC to answer my question. This model comparison appears quite different from the examples in your notebook and I must say I am a bit lost here (my apologies, my Bayesian skills are not very strong). Thanks for your help !

[OriolAbril]

I understand that in this case one observation would be having a value for each of these 5 variables?

[mservantufc]

Yes that is correct.

[OriolAbril]

Then you would have to do something similar to the "Predicting the outcome of a match" section. Your observations are 5d values (like they are 2d values in the match section).

log_lik["y"] = log_lik.y_v + log_lik.y_a + log_lik.y_ter ... 
az.loo(idata, var_name="y")

I am not sure how well will the PSIS approximation work in a 5d case though.

[mservantufc]

Alright, so I have a few issues here. 1) I updated pymc3 and arviz to be able to use 'return_inferencedata=True' and get the log likelihood. I am then trying to save the trace:

pm.save_trace(trace)

but I get the following error: AttributeError: 'InferenceData' object has no attribute '_straces'

I then provided a directory:

folder = 'C:\\Users\\mservant\\Desktop\\Research\\Simon_developpement_behav_EMG\\Modeling\\DMC_ster\\exponential_modeling\\Chains\\group_Model1'
pm.save_trace(trace, folder)

and got the following error: OSError: Cautiously refusing to overwrite the already existing C:\Users\mservant\Desktop\Research\Simon_developpement_behav_EMG\Modeling\DMC_ster\exponential_modeling\Chains\group_Model1\! Please supply a different directory, or setoverwrite=True``

If I add overwrite = True, pymc3 retuèrns: 'this file is used by another process.

[OriolAbril]

You should use to_netcdf to store InferenceData objects, pm.save_trace expects a pm.MultiTrace object. InferenceData follows the netCDF convention with its groups, variables and coords which is quite common and popular in many areas and has interfaces available in most programming languages.

[mservantufc]

Hi, I am still facing a few remaining technical issues regarding this model comparison project, and would appreciate some help/guidance. I have now fit the two exponential models to each individual dataset. PPC from Model 1, with the exponential rate parameter constrained to be similar across the 5 brain processes:

PPC from Model 2, with the exponential rate parameter constrained to be similar across the 5 brain processes:

Graphically, the performance of the 2 models appear very similar, so I would expect a big win from Model 1 (due to the fewer number of parameters). For each model, I then summed the 5 likelihoods:

log_lik_model_1["y"] = log_lik_model_1.y_v + log_lik_model_1.y_a + log_lik_model_1.y_ter + log_lik_model_1.y_max_ampl + log_lik_model_1.y_tau
log_lik_model_2["y"] = log_lik_model_2.y_v + log_lik_model_2.y_a + log_lik_model_2.y_ter + log_lik_model_2.y_max_ampl + log_lik_model_2.y_tau

Issue n°1: I get a memory error, which prevents me from computing a LOO or WAIC: MemoryError: Unable to allocate 128. GiB for an array with shape (4, 2000, 129, 129, 129) and data type float64 Any idea how to solve it ?

Issue n°2: I reduced the number of chains and samples, and did some preliminary tests with LOO and WAIC. However, I get warnings with both computations, e.g. For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. In addition, LOO and WAIC appear to favor Model 2, which doesn't make sense given the plots above (and return a large SE. Is there a more robust way to perform this model comparison? Thank you so much for the precious help.

[OriolAbril]

Issue n°1: I get a memory error, which prevents me from computing a LOO or WAIC: MemoryError: Unable to allocate 128. GiB for an array with shape (4, 2000, 129, 129, 129) and data type float64 Any idea how to solve it ?

Your log likelihood should have shape 4, 2000, 129, the 2 extra 129 are probably due to not having the correct dimensions in the variables of the log likelihood group. xarray broadcasts automatically based on the dimension names, thus, if you have 2 DataArrays with dimensions chain, draw, obs_id and chain, draw, obs_id the result will have shape chain, draw, obs_id but if the 2nd DataArray had dimensions chain, draw, observation instead the result would be chain, draw, obs_id, observation. This multiplies several times the size of the resulting array, ending up in the memory error.

I am not sure it makes sense to look at issue 2 before solving this.

[mservantufc]

Oh ok, that makes sense! I solved the issue by doing:

import xarray as xr
a = np.array(log_lik_model_1.y_v)
b = np.array(log_lik_model_1.y_a)
c = np.array(log_lik_model_1.y_ter)
d = np.array(log_lik_model_1.y_max_ampl)
e = np.array(log_lik_model_1.y_tau)
f = a+b+c+d+e
g = xr.DataArray(f, dims = ("chain","draw","y_dim"))
log_lik_model_1["y"] = g

and did the same for model_2. I then computed loo for each model:

LOO_model_1 = az.loo(trace, var_name="y")
LOO_model_2 = az.loo(trace2, var_name="y")

and got this warning:

UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  "Estimated shape parameter of Pareto distribution is greater than 0.7 for "

Looking at the results, it seems that the two models cannot be dissociated on the basis of these data: Model_1:

         Estimate       SE
elpd_loo   997.86    50.23
p_loo       32.30        -

There has been a warning during the calculation. Please check the results.

Model_2:

         Estimate       SE
elpd_loo   997.66    51.27
p_loo       34.58        -

There has been a warning during the calculation. Please check the results.

Given the posterior predictive checks, I don't understand why the two models have a very similar LOO, and suspect a problem in the LOO estimation (as suggested by the warning).

[OriolAbril]

How many points have bad shape parameter?

You can use LOO_model_1 = az.loo(trace, var_name="y", pointwise=True) to get a more informative repr. If there are only a handful it could make sense to use reloo, but being 5d observations there will probably be many bad shape parameters present. I think this is the part that is still not finished on the doc I shared in one of my first comments. PSIS struggles with these kind of problems, here there is some more guidance and a useful reference to use marginalization to solve this issue (not sure if it can be applied to your problem).

TL;DR: whenever PSIS works we can have reliable estimates of cross validation results without having to actually run cross validations, but there are cases where PSIS LOO simply can't be applied because the approximation does not hold. In such cases one must fall back to actual cross validation (or partial via reloo) or marginalization.

[mservantufc]

Thanks, your feedback is very useful. Please find below the results: Model 1:

         Estimate       SE
elpd_loo   997.86    50.23
p_loo       32.30        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)      126   97.7%
 (0.5, 0.7]   (ok)          2    1.6%
   (0.7, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    1    0.8%

Model 2:

         Estimate       SE
elpd_loo   997.66    51.27
p_loo       34.58        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)      126   97.7%
 (0.5, 0.7]   (ok)          2    1.6%
   (0.7, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    1    0.8%

Based on these results, it seems to me that Pareto k diagnostics are good for both models, which suggests reliable PSIS LOO estimates. However, there might be highly influential observations, which may completely mess up the PSIS LOO estimates. Is it possible to determine whether this is the case based on PSIS LOO diagnostics? (if not, then how can we be confident in PSIS LOO results?) Computing reloo should solve this issue, but I am struggling with the syntax (what is the 'wrapper' in the arguments?). Could you please help me implement it?

[ahartikainen]

-Pareto k diagnostics are good for both models, which indicates reliable PSIS LOO estimates.

You have on influential point, so you should do reloo.

Why would overlapping results be surprising?

[mservantufc]

The ppc look very similar between the 2 models, which suggests a win from Model 1. Of course it is just a graphical representation, but it gives a sense of the relative performance of each model no? I have difficulties to implement the reloo function. I tried: exact_LOO_model_1 = az.reloo(trace) which returns: AttributeError: 'InferenceData' object has no attribute 'check_implemented_methods' Looks like I need to transform trace to be compatible with the wrapper argument of the function, but I don't know how to proceed and can't find examples. Could you please help me implement reloo?

[OriolAbril]

You should take a look at the notebooks in https://github.com/arviz-devs/arviz/pull/1373 to get an idea about how to use reloo.

I have to admit that I expected a much worse outcome, only 1 bad khat is quite amazing :ok_person:

[mservantufc]

Thanks. Please find below my implementation if the wrapper thus far, and the errors returned by arviz

#https://github.com/arviz-devs/arviz/blob/9fa8246378ce16d2881960b4021f3f299e213e99/doc/notebooks/pymc3_refitting_xr_lik.ipynb
def compile_Model_1(Age_years_indiv, params_indiv):
    with pm.Model() as Model_1:
        α_v = pm.HalfCauchy('α_v',2)
        α_a = pm.HalfCauchy('α_a', 2) 
        α_max_ampl = pm.HalfCauchy('α_max_ampl',2)  
        α_tau = pm.HalfCauchy('α_tau', 2)     
        α_ter = pm.HalfCauchy('α_ter', 2)  

        β = pm.Normal('β', 0.334, 0.1)#Kail (1991)

        #asymptote (get the mean value of adults)
        γ_v = pm.Normal('γ_v', params_group[-1,0], .1)
        γ_a = pm.Normal('γ_a', params_group[-1,1], .1) 
        γ_ter = pm.Normal('γ_ter', params_group[-1, 2], .1)     
        γ_max_ampl = pm.Normal('γ_max_ampl', params_group[-1,3], .1)
        γ_tau = pm.Normal('γ_tau', params_group[-1,4], .1)     

        μ_v = -α_v*pm.math.exp(-β*Age_years_indiv) + γ_v
        μ_a = α_a*pm.math.exp(-β*Age_years_indiv) + γ_a      
        μ_max_ampl = α_max_ampl*pm.math.exp(-β*Age_years_indiv) + γ_max_ampl
        μ_tau = α_tau*pm.math.exp(-β*Age_years_indiv) + γ_tau 
        μ_ter = α_ter*pm.math.exp(-β*Age_years_indiv) + γ_ter     

        σ_v = pm.HalfNormal('σ_v', .1)
        σ_a = pm.HalfNormal('σ_a', .1)
        σ_max_ampl = pm.HalfNormal('σ_max_ampl', .1)
        σ_tau = pm.HalfNormal('σ_tau', .1)      
        σ_ter = pm.HalfNormal('σ_ter', .1)     

        y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_indiv[:, 0])    
        y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_indiv[:, 1])     
        y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_indiv[:, 2]) 
        y_max_ampl = pm.Normal('y_max_ampl', μ_max_ampl, σ_max_ampl, observed=params_indiv[:, 3])     
        y_tau = pm.Normal('y_tau', μ_tau, σ_tau, observed=params_indiv[:, 4])           

    return Model_1   
sample_kwargs = {"draws": 2000, "tune": 2000, "chains": 4, "target_accept": .95, init: "adapt_diag"}
with compile_Model_1(Age_years_indiv, params_indiv) as Model_1:
    trace = pm.sample(**sample_kwargs)

which returns this error: NameError: name 'adapt_diag' is not defined I need to use adapt_diag for this specific problem (otherwise the sampling crashes). So I did:

sample_kwargs = {"draws": 2000, "tune": 2000, "chains": 4, "target_accept": .95}
with compile_Model_1(Age_years_indiv, params_indiv) as Model_1:
    trace = pm.sample(**sample_kwargs,init = "adapt_diag")

But I am sure whether adapt_diag will be used when the model will be called later. Then I did:

dims = {"y": ["params_indiv"], "x": ["Age_years_indiv"]}
idata_kwargs = {
    "dims": dims,
    "log_likelihood": False,
}

idata = az.from_pymc3(trace, model=Model_1, **idata_kwargs)

def calculate_log_lik(Age_years_indiv, params_indiv, α_v, α_a, α_max_ampl, α_tau, α_ter, β, γ_v,γ_a,γ_max_ampl,γ_tau,γ_ter,σ_v,σ_a,σ_max_ampl,σ_tau,σ_ter):
    μ_v = -α_v*pm.math.exp(-β*Age_years_indiv) + γ_v
    μ_a = α_a*pm.math.exp(-β*Age_years_indiv) + γ_a      
    μ_max_ampl = α_max_ampl*pm.math.exp(-β*Age_years_indiv) + γ_max_ampl
    μ_tau = α_tau*pm.math.exp(-β*Age_years_indiv) + γ_tau 
    μ_ter = α_ter*pm.math.exp(-β*Age_years_indiv) + γ_ter     
    a = stats.norm(mu, σ_v).logpdf(params_indiv[:, 0])
    b = stats.norm(mu, σ_a).logpdf(params_indiv[:, 1])
    c = stats.norm(mu, σ_max_ampl).logpdf(params_indiv[:, 2])
    d = stats.norm(mu, σ_tau).logpdf(params_indiv[:, 3])
    e = stats.norm(mu, σ_ter).logpdf(params_indiv[:, 4])
    f = a+b+c+d+e
    return f

log_lik = xr.apply_ufunc(
    calculate_log_lik,
    idata.constant_data["Age_years_indiv"],
    idata.observed_data["params_indiv"],
    idata.posterior["α_v"],
    idata.posterior["α_a"],
    idata.posterior["α_max_ampl"],
    idata.posterior["α_tau"],
    idata.posterior["α_ter"],
    idata.posterior["β"],
    idata.posterior["γ_v"],
    idata.posterior["γ_a"],
    idata.posterior["γ_max_ampl"],
    idata.posterior["γ_tau"],
    idata.posterior["γ_ter"],
    idata.posterior["σ_v"],
    idata.posterior["σ_a"],
    idata.posterior["σ_max_ampl"],
    idata.posterior["σ_tau"],
)
idata.add_groups(log_likelihood=log_lik)

which returns this error: AttributeError: 'InferenceData' object has no attribute 'constant_data' I guess you have implemented this recently, but how should I update arviz to get the latest developments? In addition, I have a doubt with regards to the implementation of the calculate_log_lik function. Am I right to sum the 5 likelihoods?

[OriolAbril]

sample_kwargs = {"draws": 2000, "tune": 2000, "chains": 4, "target_accept": .95, init: "adapt_diag"} NameError: name 'adapt_diag' is not defined

The init key is missing the " so this should not work, not sure I completely understand how this relates to the error message you shared though.

Then I did:

def calculate_log_lik(Age_years_indiv, params_indiv, α_v, α_a, α_max_ampl, α_tau, α_ter, β, γ_v,γ_a,γ_max_ampl,γ_tau,γ_ter,σ_v,σ_a,σ_max_ampl,σ_tau,σ_ter):
    μ_v = -α_v*pm.math.exp(-β*Age_years_indiv) + γ_v
    ...  
    a = stats.norm(mu, σ_v).logpdf(params_indiv[:, 0])

Two comments here: 1) you can use np.exp (I think pymc will still work, but not completely sure so wanted to be extra sure), and 2) here you are using mu as the mean of the normal, but in the model it is μ_v, and the same for the other terms.

AttributeError: 'InferenceData' object has no attribute 'constant_data'

You have to use pm.Data to get data added to the constant_data group, otherwise you'd have to add is manually somehow. I think this blogpost will help to see how are InferenceData objects created. I should also help with dims and coordinates:

{"y": ["params_indiv"], "x": ["Age_years_indiv"]}

Won't work, as you'll see, the way that dims are defined, its keys are model variables and its values are dimension names. You'll probably want something like {"γ_v": ["indiv"], "γ_a": ["indiv"], "Age_years_indiv": ["indiv"]} or any other dimension name. I remember you had some issues with the dimension names that resulted in 129, 129, 129 shapes, you'll have to make sure that all the variables involved in the log likelihood computation have the right dimension names.

Note that for reloo you should not use coordinate values (thus not coords argument passed to pm.Model nor to from_pymc3) because each iteration would have different values for them causing the program to error out.

I guess you have implemented this recently, but how should I update arviz to get the latest developments?

The constant_data group integration was added a couple of version back, however the wrapper is yet to be merged, you should install ArviZ from the branch of this PR: pip install git+git://github.com/arviz-devs/arviz.git@refitting should do the trick.

In addition, I have a doubt with regards to the implementation of the calculate_log_lik function. Am I right to sum the 5 likelihoods?

You are right to sum the 5 likelihoods yes, otherwise you'd be calculating reloo on a different question than the original loo calculation

[mservantufc]

Alright, thanks, I think I am very close. I just have a problem regarding the implementation of the final SamplingWrapper. So far I 've done the following:

def compile_Model_1(Age_years_indiv, params_indiv):
    with pm.Model() as Model_1:
        x = pm.Data("x", Age_years_indiv)
        α_v = pm.HalfCauchy('α_v',2)
        α_a = pm.HalfCauchy('α_a', 2) 
        α_max_ampl = pm.HalfCauchy('α_max_ampl',2)  
        α_tau = pm.HalfCauchy('α_tau', 2)     
        α_ter = pm.HalfCauchy('α_ter', 2)  

        β = pm.Normal('β', 0.334, 0.1)#Kail (1991)

        #asymptote (get the mean value of adults)
        γ_v = pm.Normal('γ_v', params_group[-1,0], .1)
        γ_a = pm.Normal('γ_a', params_group[-1,1], .1) 
        γ_ter = pm.Normal('γ_ter', params_group[-1, 2], .1)     
        γ_max_ampl = pm.Normal('γ_max_ampl', params_group[-1,3], .1)
        γ_tau = pm.Normal('γ_tau', params_group[-1,4], .1)     

        μ_v = -α_v*pm.math.exp(-β*x) + γ_v
        μ_a = α_a*pm.math.exp(-β*x) + γ_a      
        μ_max_ampl = α_max_ampl*pm.math.exp(-β*x) + γ_max_ampl
        μ_tau = α_tau*pm.math.exp(-β*x) + γ_tau 
        μ_ter = α_ter*pm.math.exp(-β*x) + γ_ter     

        σ_v = pm.HalfNormal('σ_v', .1)
        σ_a = pm.HalfNormal('σ_a', .1)
        σ_max_ampl = pm.HalfNormal('σ_max_ampl', .1)
        σ_tau = pm.HalfNormal('σ_tau', .1)      
        σ_ter = pm.HalfNormal('σ_ter', .1)     

        y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_indiv[:, 0])    
        y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_indiv[:, 1])     
        y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_indiv[:, 2]) 
        y_max_ampl = pm.Normal('y_max_ampl', μ_max_ampl, σ_max_ampl, observed=params_indiv[:, 3])     
        y_tau = pm.Normal('y_tau', μ_tau, σ_tau, observed=params_indiv[:, 4])           

    return Model_1

Note that I now use pm.Data for my age (x) variable.

sample_kwargs = {"draws": 2000, "tune": 2000, "chains": 4, "target_accept": .95, "init": "adapt_diag"}
with compile_Model_1(Age_years_indiv, params_indiv) as Model_1:
    trace = pm.sample(**sample_kwargs)

dims = {"x":["indiv"], "y_v":["indiv"], "y_a":["indiv"], "y_max_ampl":["indiv"], "y_tau":["indiv"],"y_ter":["indiv"]}

idata_kwargs = {
    "dims": dims,
    "log_likelihood": False,
}

idata = az.from_pymc3(trace, model=Model_1, **idata_kwargs)

def calculate_log_lik(x, y_v,y_a,y_max_ampl,y_tau,y_ter,  α_v, α_a,α_max_ampl, α_tau, α_ter,  β,   γ_v,γ_a,γ_max_ampl,γ_tau,γ_ter,  σ_v,σ_a,σ_max_ampl,σ_tau,σ_ter):
    μ_v = -α_v*np.exp(-β*x) + γ_v
    μ_a = α_a*np.exp(-β*x) + γ_a      
    μ_max_ampl = α_max_ampl*np.exp(-β*x) + γ_max_ampl
    μ_tau = α_tau*np.exp(-β*x) + γ_tau 
    μ_ter = α_ter*np.exp(-β*x) + γ_ter         

    a = stats.norm(μ_v, σ_v).logpdf(y_v)
    b = stats.norm(μ_a, σ_a).logpdf(y_a)
    c = stats.norm(μ_max_ampl, σ_max_ampl).logpdf(y_max_ampl)
    d = stats.norm(μ_tau, σ_tau).logpdf(y_tau)
    e = stats.norm(μ_ter, σ_ter).logpdf(y_ter)
    f = a+b+c+d+e

    return f

Given that I have several dependent variables (observed_data), I then did:

import xarray as xr
log_lik = xr.apply_ufunc(
    calculate_log_lik,
    idata.constant_data["x"],
    idata.observed_data["y_v"],
    idata.observed_data["y_a"],
    idata.observed_data["y_max_ampl"],
    idata.observed_data["y_tau"],
    idata.observed_data["y_ter"],    
    idata.posterior["α_v"],
    idata.posterior["α_a"],
    idata.posterior["α_max_ampl"],
    idata.posterior["α_tau"],
    idata.posterior["α_ter"],
    idata.posterior["β"],
    idata.posterior["γ_v"],
    idata.posterior["γ_a"],
    idata.posterior["γ_max_ampl"],
    idata.posterior["γ_tau"],
    idata.posterior["γ_ter"],
    idata.posterior["σ_v"],
    idata.posterior["σ_a"],
    idata.posterior["σ_max_ampl"],
    idata.posterior["σ_tau"],
    idata.posterior["σ_ter"]
)
idata.add_groups(log_likelihood=log_lik)

Am I correct? The issue now concerns the final phase i.e the definition of the SamplingWrapper. I must say I am completely lost here, and I don't understand how to modify your example below to the present problem, and would greatly appreciate your help to finish this model comparison!

class PyMC3LinRegWrapper(az.SamplingWrapper):        
    def sample(self, modified_observed_data):
        with self.model(*modified_observed_data) as linreg_model:
            idata = pm.sample(
                **self.sample_kwargs, 
                return_inferencedata=True, 
                idata_kwargs=self.idata_kwargs
            )
        return idata

    def get_inference_data(self, idata):
        return idata

    def sel_observations(self, idx):
        xdata = self.idata_orig.constant_data["x"]
        ydata = self.idata_orig.observed_data["y"]
        mask = np.isin(np.arange(len(xdata)), idx)
        data__i = [ary[~mask] for ary in (xdata, ydata)]
        data_ex = [ary[mask] for ary in (xdata, ydata)]
        return data__i, data_ex 

[OriolAbril]

Am I correct?

Up until here everything looks good

The issue now concerns the final phase i.e the definition of the SamplingWrapper. I must say I am completely lost here, and I don't understand how to modify your example below to the present problem, and would greatly appreciate your help to finish this model comparison!

Both sample and get_inference_data should not be modified, only sel_observations which is the model dependent part of the wrapper, depending on both the model and the predictive question at hand it will vary. I think something similar to this will work, but some work on the shapes/stacks may be required:

Note that data__i will be passed as is to wrapper.sample (and thus as *data__i to compile_Model_1) so it should be a tuple/list of 2 elements, age and paramsindiv. We can recreate these stacking the `y*variables (or also by addingparams_indivtoidata_orig`).

data_ex will be passed to log_likelihood__i instead (which is not seen in the notebook, it's part of the base sampling wrapper and is thus inherited), so we want it to be a tuple/list of x, y_v,y_a,y_max_ampl,y_tau,y_ter

    def sel_observations(self, idx):
        x = self.idata_orig.constant_data["x"]
        y_list = [self.idata_orig.observed_data[var_name] for var_name in (y_v,y_a,y_max_ampl,y_tau,y_ter)] # add "
        mask = np.isin(np.arange(len(x)), idx)
        y__i = [ary[~mask] for ary in (y_v,y_a,y_max_ampl,y_tau,y_ter)] # add "
        y_ex = [ary[mask] for ary in (y_v,y_a,y_max_ampl,y_tau,y_ter)] # add "
        data__i = (x[~mask], np.hstack(y__i)) # maybe vstack?
        data_ex = [x[mask], *y_ex]
        return data__i, data_ex 

[mservantufc]

I see, thanks. In your notebook example, you then call:

pymc3_wrapper = PyMC3LinRegWrapper(
    model=compile_linreg_model, 
    log_lik_fun=calculate_log_lik, 
    posterior_vars=("b0", "b1", "sigma_e"),
    idata_orig=idata,
    sample_kwargs=sample_kwargs, 
    idata_kwargs=idata_kwargs,
)

I am confused about log_lik_fun, which is not defined in your class (and thus returns an error in my case: __init__() got an unexpected keyword argument 'log_lik_fun'

[OriolAbril]

The class has no init method defined at all so it inherits the init method of the base class. The updated base class is only available in the PR, it has not been merged yet which is why it won't work unless you install the branch of the refitting PR.

I'd recommend making sure that the example notebook runs to make sure you correctpy installed the branch instead of pypi or github master versions. I shared the install command above

[mservantufc]

I did: pip install git+git://github.com/arviz-devs/arviz.git@refitting but your example notebook returns the very same error:

pymc3_wrapper = PyMC3LinRegWrapper(
    model=compile_linreg_model, 
    log_lik_fun=calculate_log_lik, 
    posterior_vars=("b0", "b1", "sigma_e"),
    idata_orig=idata,
    sample_kwargs=sample_kwargs, 
    idata_kwargs=idata_kwargs,
)
TypeError: __init__() got an unexpected keyword argument 'log_lik_fun'

[OriolAbril]

This is probably because it's not being installed (i.e. the version number is the same as the latest release so pip thinks it the same code). Try a --force-reinstall or downgrading and then installing from the branch

[mservantufc]

Alright thanks, this issue is solved. I think we are very close. Running: loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1) returns the following error:

C:\Users\mservant\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py:99: UserWarning: reloo is an experimental and untested feature
  warnings.warn("reloo is an experimental and untested feature", UserWarning)
arviz.stats.stats_refitting - INFO - Refitting model excluding observation 82
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~\Anaconda3\lib\site-packages\xarray\core\dataset.py in _construct_dataarray(self, name)
   1171         try:
-> 1172             variable = self._variables[name]
   1173         except KeyError:

KeyError: y_v

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
<ipython-input-35-dcd11ffc8f30> in <module>
----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose)
    103             if verbose:
    104                 _log.info("Refitting model excluding observation %d", idx)
--> 105             new_obs, excluded_obs = wrapper.sel_observations(idx)
    106             fit = wrapper.sample(new_obs)
    107             idata_idx = wrapper.get_inference_data(fit)

<ipython-input-29-2c5a4c4ebe9e> in sel_observations(self, idx)
     14     def sel_observations(self, idx):
     15         x = self.idata_orig.constant_data["x"]
---> 16         y_list = [self.idata_orig.observed_data[var_name] for var_name in (y_v,y_a,y_max_ampl,y_tau,y_ter)] # add "
     17         mask = np.isin(np.arange(len(x)), idx)
     18         y__i = [ary[~mask] for ary in (y_v,y_a,y_max_ampl,y_tau,y_ter)] # add "

<ipython-input-29-2c5a4c4ebe9e> in <listcomp>(.0)
     14     def sel_observations(self, idx):
     15         x = self.idata_orig.constant_data["x"]
---> 16         y_list = [self.idata_orig.observed_data[var_name] for var_name in (y_v,y_a,y_max_ampl,y_tau,y_ter)] # add "
     17         mask = np.isin(np.arange(len(x)), idx)
     18         y__i = [ary[~mask] for ary in (y_v,y_a,y_max_ampl,y_tau,y_ter)] # add "

~\Anaconda3\lib\site-packages\xarray\core\dataset.py in __getitem__(self, key)
   1270 
   1271         if hashable(key):
-> 1272             return self._construct_dataarray(key)
   1273         else:
   1274             return self._copy_listed(np.asarray(key))

~\Anaconda3\lib\site-packages\xarray\core\dataset.py in _construct_dataarray(self, name)
   1173         except KeyError:
   1174             _, name, variable = _get_virtual_variable(
-> 1175                 self._variables, name, self._level_coords, self.dims
   1176             )
   1177 

~\Anaconda3\lib\site-packages\xarray\core\dataset.py in _get_virtual_variable(variables, key, level_vars, dim_sizes)
    154 
    155     if not isinstance(key, str):
--> 156         raise KeyError(key)
    157 
    158     split_key = key.split(".", 1)

KeyError: y_v

[OriolAbril]

All the (y_v,y_a,y_max_ampl,y_tau,y_ter) tuples have to be tuples of strings, you have to add the " to each variable name for it to work

[mservantufc]

Good point, but there is still a problem:

class PyMC3Model_1Wrapper(az.SamplingWrapper):        
    def sample(self, observed_data):
        with self.model(*observed_data) as Model_1:
            idata = pm.sample(
                **self.sample_kwargs, 
                return_inferencedata=True, 
                idata_kwargs=self.idata_kwargs
            )
        return idata

    def get_inference_data(self, idata):
        return idata

    def sel_observations(self, idx):
        x = self.idata_orig.constant_data["x"]
        y_list = [self.idata_orig.observed_data[var_name] for var_name in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
        mask = np.isin(np.arange(len(x)), idx)
        y__i = [ary[~mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
        y_ex = [ary[mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
        data__i = (x[~mask], np.hstack(y__i)) # maybe vstack?
        data_ex = [x[mask], *y_ex]
        return data__i, data_ex 

pymc3_wrapper = PyMC3Model_1Wrapper(
    model=compile_Model_1, 
    log_lik_fun=calculate_log_lik, 
    posterior_vars=("α_v", "α_a","α_max_ampl","α_tau","α_ter","α_tau","β","γ_v","γ_a","γ_max_ampl","γ_tau","γ_ter","σ_v","σ_a","σ_max_ampl","σ_tau","σ_ter"),
    idata_orig=idata,
    sample_kwargs=sample_kwargs, 
    idata_kwargs=idata_kwargs,
)

loo_Model_1 = az.loo(idata, pointwise=True)

loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-33-dcd11ffc8f30> in <module>
----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose)
    103             if verbose:
    104                 _log.info("Refitting model excluding observation %d", idx)
--> 105             new_obs, excluded_obs = wrapper.sel_observations(idx)
    106             fit = wrapper.sample(new_obs)
    107             idata_idx = wrapper.get_inference_data(fit)

<ipython-input-29-fd08f94b51d7> in sel_observations(self, idx)
     16         y_list = [self.idata_orig.observed_data[var_name] for var_name in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     17         mask = np.isin(np.arange(len(x)), idx)
---> 18         y__i = [ary[~mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     19         y_ex = [ary[mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     20         data__i = (x[~mask], np.hstack(y__i)) # maybe vstack?

<ipython-input-29-fd08f94b51d7> in <listcomp>(.0)
     16         y_list = [self.idata_orig.observed_data[var_name] for var_name in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     17         mask = np.isin(np.arange(len(x)), idx)
---> 18         y__i = [ary[~mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     19         y_ex = [ary[mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     20         data__i = (x[~mask], np.hstack(y__i)) # maybe vstack?

TypeError: only integer scalar arrays can be converted to a scalar index

[ahartikainen]

Hi, can you print both the inputs on that line and also type(...) for inputs?

I wonder if there is an empty selection which goes to float?

[mservantufc]

Hi, do you mean line 18? I did:

    def sel_observations(self, idx):
        x = self.idata_orig.constant_data["x"]
        y_list = [self.idata_orig.observed_data[var_name] for var_name in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
        mask = np.isin(np.arange(len(x)), idx)
        for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter"):
            print (ary)
            print (type(ary))**
        y__i = [ary[~mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
        y_ex = [ary[mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
        data__i = (x[~mask], np.hstack(y__i)) # maybe vstack?
        data_ex = [x[mask], *y_ex]
        return data__i, data_ex 

which returns:

y_v
<class 'str'>
y_a
<class 'str'>
y_max_ampl
<class 'str'>
y_tau
<class 'str'>
y_ter
<class 'str'>
y_v
<class 'str'>
y_a
<class 'str'>
y_max_ampl
<class 'str'>
y_tau
<class 'str'>
y_ter
<class 'str'>
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-33-dcd11ffc8f30> in <module>
----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose)
    103             if verbose:
    104                 _log.info("Refitting model excluding observation %d", idx)
--> 105             new_obs, excluded_obs = wrapper.sel_observations(idx)
    106             fit = wrapper.sample(new_obs)
    107             idata_idx = wrapper.get_inference_data(fit)

<ipython-input-29-8ec449826632> in sel_observations(self, idx)
     20             print (ary)
     21             print (type(ary))
---> 22         y__i = [ary[~mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     23         y_ex = [ary[mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     24         data__i = (x[~mask], np.hstack(y__i)) # maybe vstack?

<ipython-input-29-8ec449826632> in <listcomp>(.0)
     20             print (ary)
     21             print (type(ary))
---> 22         y__i = [ary[~mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     23         y_ex = [ary[mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
     24         data__i = (x[~mask], np.hstack(y__i)) # maybe vstack?

TypeError: only integer scalar arrays can be converted to a scalar index

[OriolAbril]

These

y__i = [ary[~mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
y_ex = [ary[mask] for ary in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "

should be

y__i = [ary[~mask] for ary in y_list]
y_ex = [ary[mask] for ary in y_list]

we want to subset the actual arrays, not their names

[mservantufc]

Still an error :-(

IndexError                                Traceback (most recent call last)
<ipython-input-42-dcd11ffc8f30> in <module>
----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose)
    104                 _log.info("Refitting model excluding observation %d", idx)
    105             new_obs, excluded_obs = wrapper.sel_observations(idx)
--> 106             fit = wrapper.sample(new_obs)
    107             idata_idx = wrapper.get_inference_data(fit)
    108             log_like_idx = wrapper.log_likelihood__i(excluded_obs, idata_idx).values.flatten()

<ipython-input-40-ddd7af045286> in sample(self, observed_data)
      1 class PyMC3Model_1Wrapper(az.SamplingWrapper):
      2     def sample(self, observed_data):
----> 3         with self.model(*observed_data) as Model_1:
      4             idata = pm.sample(
      5                 **self.sample_kwargs,

<ipython-input-23-1c20af5a7205> in compile_Model_1(Age_years_indiv, params_indiv)
     31 
     32 
---> 33         y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_indiv[:, 0])
     34         y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_indiv[:, 1])
     35         y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_indiv[:, 2])

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

[OriolAbril]

Have you tried changing the hstack with vstack?

[mservantufc]

Yes, in this case I get a different error:

ValueError                                Traceback (most recent call last)
<ipython-input-50-dcd11ffc8f30> in <module>
----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose)
    104                 _log.info("Refitting model excluding observation %d", idx)
    105             new_obs, excluded_obs = wrapper.sel_observations(idx)
--> 106             fit = wrapper.sample(new_obs)
    107             idata_idx = wrapper.get_inference_data(fit)
    108             log_like_idx = wrapper.log_likelihood__i(excluded_obs, idata_idx).values.flatten()

<ipython-input-46-165a42955f70> in sample(self, observed_data)
      1 class PyMC3Model_1Wrapper(az.SamplingWrapper):
      2     def sample(self, observed_data):
----> 3         with self.model(*observed_data) as Model_1:
      4             idata = pm.sample(
      5                 **self.sample_kwargs,

<ipython-input-23-1c20af5a7205> in compile_Model_1(Age_years_indiv, params_indiv)
     31 
     32 
---> 33         y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_indiv[:, 0])
     34         y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_indiv[:, 1])
     35         y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_indiv[:, 2])

~\Anaconda3\lib\site-packages\pymc3\distributions\distribution.py in __new__(cls, name, *args, **kwargs)
     81         else:
     82             dist = cls.dist(*args, **kwargs)
---> 83         return model.Var(name, dist, data, total_size, dims=dims)
     84 
     85     def __getnewargs__(self):

~\Anaconda3\lib\site-packages\pymc3\model.py in Var(self, name, dist, data, total_size, dims)
   1115                     distribution=dist,
   1116                     total_size=total_size,
-> 1117                     model=self,
   1118                 )
   1119             self.observed_RVs.append(var)

~\Anaconda3\lib\site-packages\pymc3\model.py in __init__(self, type, owner, index, name, data, distribution, total_size, model)
   1738 
   1739             self.missing_values = data.missing_values
-> 1740             self.logp_elemwiset = distribution.logp(data)
   1741             # The logp might need scaling in minibatches.
   1742             # This is done in `Factor`.

~\Anaconda3\lib\site-packages\pymc3\distributions\continuous.py in logp(self, value)
    535         mu = self.mu
    536 
--> 537         return bound((-tau * (value - mu)**2 + tt.log(tau / np.pi / 2.)) / 2.,
    538                      sigma > 0)
    539 

~\Anaconda3\lib\site-packages\theano\tensor\var.py in __sub__(self, other)
    148         # and the return value in that case
    149         try:
--> 150             return theano.tensor.basic.sub(self, other)
    151         except (NotImplementedError, AsTensorError):
    152             return NotImplemented

~\Anaconda3\lib\site-packages\theano\gof\op.py in __call__(self, *inputs, **kwargs)
    672                 thunk.outputs = [storage_map[v] for v in node.outputs]
    673 
--> 674                 required = thunk()
    675                 assert not required  # We provided all inputs
    676 

~\Anaconda3\lib\site-packages\theano\gof\op.py in rval()
    860 
    861         def rval():
--> 862             thunk()
    863             for o in node.outputs:
    864                 compute_map[o][0] = True

~\Anaconda3\lib\site-packages\theano\gof\cc.py in __call__(self)
   1737                 print(self.error_storage, file=sys.stderr)
   1738                 raise
-> 1739             reraise(exc_type, exc_value, exc_trace)
   1740 
   1741 

~\Anaconda3\lib\site-packages\six.py in reraise(tp, value, tb)
    701             if value.__traceback__ is not tb:
    702                 raise value.with_traceback(tb)
--> 703             raise value
    704         finally:
    705             value = None

ValueError: Input dimension mis-match. (input[0].shape[0] = 5, input[1].shape[0] = 128)

[OriolAbril]

Stacking should convert five 128 arrays into a 2d array with shape 128, 5 array (same shape as params_indiv) but if the stacking could also generate a 1d array with shape 640 or even a 2d array with shape 5, 128 which won't work either.

I'd recommend playing around with five 1d aranges or arrays of ones, the content is not relevant and see how to stack them into the right shape so you get a better understanding of what is happening. See what does hstack do, what is the difference with vstack, how do they compare to stack and what effect has the axis argument of stack. As a last resort you can always transpose the array to get the 128, 5 from the 5, 128

[mservantufc]

Yes, I tried to transpose, printing the shape of yi, datai, and data_ex: def sel_observations(self, idx): x = self.idata_orig.constant_data["x"] y_list = [self.idata_orig.observed_data[var_name] for var_name in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add " mask = np.isin(np.arange(len(x)), idx) yi = [ary[~mask] for ary in y_list] y_ex = [ary[mask] for ary in y_list] yi = np.array(yi) yi = yi.T yi = yi.tolist() print (np.shape(yi)) datai = (x[~mask], yi) print (np.shape(datai)) data_ex = [x[mask], *y_ex] print (np.shape(data_ex)) return datai, data_ex

which returns:

(128, 5)
(2, 128)
(6, 1)
(128, 5)
(2, 128)
(6, 1)
C:\Users\mservant\Anaconda3\lib\site-packages\numpy\core\_asarray.py:83: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-130-dcd11ffc8f30> in <module>
----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose)
    104                 _log.info("Refitting model excluding observation %d", idx)
    105             new_obs, excluded_obs = wrapper.sel_observations(idx)
--> 106             fit = wrapper.sample(new_obs)
    107             idata_idx = wrapper.get_inference_data(fit)
    108             log_like_idx = wrapper.log_likelihood__i(excluded_obs, idata_idx).values.flatten()

<ipython-input-126-051efb0daa68> in sample(self, observed_data)
      1 class PyMC3Model_1Wrapper(az.SamplingWrapper):
      2     def sample(self, observed_data):
----> 3         with self.model(*observed_data) as Model_1:
      4             idata = pm.sample(
      5                 **self.sample_kwargs,

<ipython-input-23-1c20af5a7205> in compile_Model_1(Age_years_indiv, params_indiv)
     31 
     32 
---> 33         y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_indiv[:, 0])
     34         y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_indiv[:, 1])
     35         y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_indiv[:, 2])

TypeError: list indices must be integers or slices, not tuple

[OriolAbril]

yi = yi.tolist()

y__i must be a 2d numpy array, not a list.

You have to be careful if using np.shape because it gives the shape of the object after having passed it through np.array, which is why you are getting 2d shapes that can not really be indexed as 2d arrays because your object is a list and has no shape.

For example y.shape will work if y is a numpy array but not if it is a list. Therefore, if I have seen that y.shape is (128, 5) I know that y[:, 0] is a valid indexing operation. On the other hand, np.shape(y) will work with both, but as you can see, even if np.shape(y) is (128, 5), doing y[:, 0] does not neccessarily work

[mservantufc]

Ok, we are now very close. I solved the shape issue doing using np.stack(y__i,axis = 1):

class PyMC3Model_1Wrapper(az.SamplingWrapper):        
    def sample(self, observed_data):
        with self.model(*observed_data) as Model_1:
            idata = pm.sample(
                **self.sample_kwargs, 
                return_inferencedata=True, 
                idata_kwargs=self.idata_kwargs
            )
        return idata

    def get_inference_data(self, idata):
        return idata

    def sel_observations(self, idx):
        x = self.idata_orig.constant_data["x"]
        y_list = [self.idata_orig.observed_data[var_name] for var_name in ("y_v","y_a","y_max_ampl","y_tau","y_ter")] # add "
        mask = np.isin(np.arange(len(x)), idx)
        y__i = [ary[~mask] for ary in y_list]
        y_ex = [ary[mask] for ary in y_list]
        data__i = (x[~mask], np.stack(y__i,axis = 1)) 
        data_ex = [x[mask], *y_ex]
        return data__i, data_ex 

I then ran loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1) which gave:

TypeError                                 Traceback (most recent call last)
<ipython-input-68-dcd11ffc8f30> in <module>
----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1)

~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose)
    106             fit = wrapper.sample(new_obs)
    107             idata_idx = wrapper.get_inference_data(fit)
--> 108             log_like_idx = wrapper.log_likelihood__i(excluded_obs, idata_idx).values.flatten()
    109             loo_lppd_idx = scale_value * _logsumexp(log_like_idx, b_inv=len(log_like_idx))
    110             khats[idx] = 0

~\Anaconda3\lib\site-packages\arviz\wrappers\base.py in log_likelihood__i(self, excluded_obs, idata__i)
    181             *arys,
    182             kwargs=self.log_lik_kwargs,
--> 183             **self.apply_ufunc_kwargs,
    184         )
    185         return log_lik_idx

~\Anaconda3\lib\site-packages\xarray\core\computation.py in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, *args)
   1108             join=join,
   1109             exclude_dims=exclude_dims,
-> 1110             keep_attrs=keep_attrs,
   1111         )
   1112     # feed Variables directly through apply_variable_ufunc

~\Anaconda3\lib\site-packages\xarray\core\computation.py in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args)
    260 
    261     data_vars = [getattr(a, "variable", a) for a in args]
--> 262     result_var = func(*data_vars)
    263 
    264     if signature.num_outputs > 1:

~\Anaconda3\lib\site-packages\xarray\core\computation.py in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args)
    698             )
    699 
--> 700     result_data = func(*input_data)
    701 
    702     if signature.num_outputs == 1:

TypeError: calculate_log_lik() takes 22 positional arguments but 23 were given

It looks like my calculate_log_lik() function (below) is missing an argument, but I don't understand what is expected here.

def calculate_log_lik(x, y_v,y_a,y_max_ampl,y_tau,y_ter,  α_v, α_a,α_max_ampl, α_tau, α_ter,  β, γ_v,γ_a,γ_max_ampl,γ_tau,γ_ter,  σ_v,σ_a,σ_max_ampl,σ_tau,σ_ter):
    μ_v = -α_v*np.exp(-β*x) + γ_v
    μ_a = α_a*np.exp(-β*x) + γ_a      
    μ_max_ampl = α_max_ampl*np.exp(-β*x) + γ_max_ampl
    μ_tau = α_tau*np.exp(-β*x) + γ_tau 
    μ_ter = α_ter*np.exp(-β*x) + γ_ter         

    a = stats.norm(μ_v, σ_v).logpdf(y_v)
    b = stats.norm(μ_a, σ_a).logpdf(y_a)
    c = stats.norm(μ_max_ampl, σ_max_ampl).logpdf(y_max_ampl)
    d = stats.norm(μ_tau, σ_tau).logpdf(y_tau)
    e = stats.norm(μ_ter, σ_ter).logpdf(y_ter)
    f = a+b+c+d+e

    return f

[ahartikainen]

Maybe try with this

def calculate_log_lik(*args, *kwargs):
    print(args)
    print(kwargs)

And see if you notice something special (e.g. a model handle etc)

[mservantufc]

Are you suggesting to replace my calculate_log_lik with the one above? If so, I get this message: T``` ypeError Traceback (most recent call last)

in ----> 1 loo_relooed_Model_1 = az.reloo(pymc3_wrapper, loo_orig=loo_Model_1) ~\Anaconda3\lib\site-packages\arviz\stats\stats_refitting.py in reloo(wrapper, loo_orig, k_thresh, scale, verbose) 106 fit = wrapper.sample(new_obs) 107 idata_idx = wrapper.get_inference_data(fit) --> 108 log_like_idx = wrapper.log_likelihood__i(excluded_obs, idata_idx).values.flatten() 109 loo_lppd_idx = scale_value * _logsumexp(log_like_idx, b_inv=len(log_like_idx)) 110 khats[idx] = 0 ~\Anaconda3\lib\site-packages\arviz\wrappers\base.py in log_likelihood__i(self, excluded_obs, idata__i) 181 *arys, 182 kwargs=self.log_lik_kwargs, --> 183 **self.apply_ufunc_kwargs, 184 ) 185 return log_lik_idx ~\Anaconda3\lib\site-packages\xarray\core\computation.py in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, *args) 1108 join=join, 1109 exclude_dims=exclude_dims, -> 1110 keep_attrs=keep_attrs, 1111 ) 1112 # feed Variables directly through apply_variable_ufunc ~\Anaconda3\lib\site-packages\xarray\core\computation.py in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args) 260 261 data_vars = [getattr(a, "variable", a) for a in args] --> 262 result_var = func(*data_vars) 263 264 if signature.num_outputs > 1: ~\Anaconda3\lib\site-packages\xarray\core\computation.py in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args) 698 ) 699 --> 700 result_data = func(*input_data) 701 702 if signature.num_outputs == 1: TypeError: calculate_log_lik() takes 22 positional arguments but 23 were given ```

[OriolAbril]

The error message is telling you that the log lik function expects 22 arguments: 5 y + 5 alpha + 5 gamma + 5 sigma + x + beta but it it actually getting 23 inputs instead, there is one extra input, probably repeated.

In previous stages I remember seeing that data_ex is a list of length 6 (x + the 5 y) so this should be fine unless modified.

I would start double checking the posterior vars argument used when initializing the wrapper and them move to doble checking the length of data_ex until the problem is found.

As you can see in the notebook, the signature way that lok lik function is called is with *data_ex, *posterior_data_arrays so this should create a list of length 22, not 23

[mservantufc]

Alright guys, thank you! The results are interesting. The original arviz.loo for Model_1 and Model_2 gave the following results:

Model_1:

Computed from 8000 by 129 log-likelihood matrix

         Estimate       SE
elpd_loo   997.68    50.69
p_loo       32.74        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)      127   98.4%
 (0.5, 0.7]   (ok)          0    0.0%
   (0.7, 1]   (bad)         1    0.8%
   (1, Inf)   (very bad)    1    0.8%

Model_2:

Computed from 8000 by 129 log-likelihood matrix

         Estimate       SE
elpd_loo   998.04    51.06
p_loo       33.99        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)      125   96.9%
 (0.5, 0.7]   (ok)          3    2.3%
   (0.7, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    1    0.8%

The reloo function gave the following results: Model_1:

Computed from 8000 by 129 log-likelihood matrix

         Estimate       SE
elpd_loo   799.66   187.55
p_loo      230.76        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)      129  100.0%
 (0.5, 0.7]   (ok)          0    0.0%
   (0.7, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    0    0.0%

Model_2:

Computed from 8000 by 129 log-likelihood matrix

         Estimate       SE
elpd_loo   945.95    95.38
p_loo       86.08        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)      126   97.7%
 (0.5, 0.7]   (ok)          3    2.3%
   (0.7, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    0    0.0%

A couple of issues: 1) I still see a warning with reloo. Do you guys know why? 2) I don't understand the large p_loo for Model_1 (230.76). In addition, the large SE is somewhat strange. 3) Based on these findings, woiuld you guys conclude that Model_2 has better predictive accuracy than Model_1 ?

[OriolAbril]

I still see a warning with reloo. Do you guys know why?

The loo object stores if there has been a warning during the calculations as a boolean (i't not checked every time it's printed), my guess is that reloo does not set the warning field to false. reloo is still an experimental and untested feature.

I don't understand the large p_loo for Model_1 (230.76).

Ignore the p_loo parameter, I am quite sure it's wrong.

In addition, the large SE is somewhat strange.

This one is not so surprising but it is somewhat strange. plot_elpd may give more insights into this. And it would probably be interesting to ask about this on https://discourse.mc-stan.org/, reloo is an experimental feature here but it has been available in the stan R ecosystem for a while. They will have much more experience using reloo than me.

Based on these findings, woiuld you guys conclude that Model_2 has better predictive accuracy than Model_1 ?

After reloo it has a higher elpd but the large errors are troublesome. I would probably rerun reloo with a 0.5 threshold to be extra sure. The result should be exactly the same given that PSIS should be able to correctly estimate the elpd (same thing that happens in the example notebooks, PSIS worked so reloo gives the same result and confirms that reloo worked). I also wanted to point out that ArviZ has compare and plot_compare functions to ease interpretation of loo/waic results.