A wrapper function that allows the value of a TransformedRV on a given realization would be wonderful

[steveehlert]

I need to run some MCMC analysis on a rather complex model initially written in Fortran. Because of the way the initial Fortran code was structured, the best course of action right now is to have Python make a parameters txt file that gets input into the Fortran, have the Fortran do its things and write its output txt files to disk, and then have the Python read those outputs and then calculate the likelihood. A rough example is something like this:

with pymc3.Model() as MCMCModel:
    with open("PARAMS.TXT",'w') as file:
        rinf = pymc3.Lognormal('rinf',mu = np.log(self.rinf), sd = 0.5*np.log(10))  # 
        rinf_send = T.cast(rinf,dtype='float64').eval()  #not correct, but roughly what I think I am 
        #supposed to do
        towrite="""                                                                                                                                                                                                                                                                         
            %0.3f          RINF  """ %(rinf_send)
        file.write(towrite)

    #Run Fortran
    #Open and Load model outputs 
    obs = pymc3.Normal('obs', mu=modelmags, sd=0.5, observed=self.datadict['measurements'])

with MCMCModel: 
    trace = pymc3.sample(1000,njobs=4)

However, I will need to get the value of the TransformedRV rinf for each realization, which I am finding surprisingly challenging to do with the Theano tensors. I have done a lot of MCMC using the older version of pymc, but this is my first rodeo with pymc3. I have attempted some theano as_ops wrappers, but as the code above suggests the entire model run is wrapped in an object that may make those unwieldy.

Two questions/issues I therefore have:

1. What is the easiest way to get the value of rinf for each realization? Are there wrapper functions for it? I have had no such luck finding them in the code thus far.
2. Can you add that functionality to pymc3? I suspect that other people may find this useful, since I am finding Theano tensors very abstract and unwieldy to deal with for my problems.

Thank you for your help.

Steven Ehlert

Versions and main components

• PyMC3 Version: PyMC3 3.2
• Theano Version: 0.9
• Python Version: 2.7.3
• Operating system: Debian 3.2.88-1 x86_64 GNU/Linux
• How did you install PyMC3: pip

[fonnesbeck]

Can you not write an f2py wrapper for your Fortran code that turns it into a Python function? Then, you are more likely to have success with the Theano as_op. I assume that it just returns something that is treated as a Deterministic?

[steveehlert]

I have tried/thought about doing this, and while it would be ideal the structure of the Fortran makes this rather unfeasible right now. Most of the variables I need are arrays in general but written out to txt files one line at a time so I would have to make major and unpleasant changes to the memory structure of the Fortran to make that work. Since this is an operational model those changes are very strongly discouraged.

Since I wrote this, I did find a way to get some of the numbers working, but the sampler does not appear to be writing out new parameter files for each iteration. Any help/advice as to how to get that to work would be most appreciated.

[fonnesbeck]

There is a related discussion in #1925 which recommends using SMC.