How to define model depending on the output from the external program

Questions
1

Hello!

I would like to make a model like that:

with pm.Model() as linear_model:
   tx = pm.Uniform('tx', upper=70, lower=0)
   ty = pm.Uniform('ty', upper=70, lower=0)
   mu_sim = pm.MvNormal('mu_sim', 
                         mu=func0(tx,ty),
                         conv=R,
                         observed=mu)
    p = pm.sample()

It is simple but the problem that function func0 is calling some external program. (R and mu are numpy arrays) I have to transform variables tx and ty from stochastic to deterlministic to make it work, right? But I do not completelly understand how exactly it should work.
Thanks

2

You need to make func0 a theano op. What does it do exactly? If you post the code perhaps someone can help convert it.

3

Thanks for your reply. I googled about theano.Op and here I was trying to adapt one of the examples:

class Case(): 
.....
    def H(self, xs, ys, zs, qs):
        H = np.zeros([self.d, 1])
        for i, scalar in enumerate(self.adj_scalar_names):
            H[i] = qs * self.adj_solution.probe_real(scalar, xs, ys, zs);
        return H
....

case_2681829 = Case(mu_name = "probes_conc.csv",
                   mu_coord = "probes_coords.csv",
                   adj_sol_RESULT = 'RESULTS.case')

class theano_wrapper_H(tt.Op):
    itypes = [tt.dscalar]
    otypes = [tt.dscalar]

def perform(self, node, inputs, outputs):
    xs = inputs[0]
    ys = inputs[1]
    zs = inputs[2]
    qs = inputs[3]
    mu = case_2681829.H(xs, ys, zs, qs) 
    #mu should be a numpy array of length d
    outputs[0][0] = np.array(mu)

with pm.Model() as model:
    xs = pm.Uniform('theta_x', upper=70, lower=0)
    ys = pm.Uniform('theta_y', upper=105, lower=0)
    zs = pm.Uniform('theta_z', upper=7,  lower=0)
    qs = pm.Uniform('theta_q', upper=200, lower=0)

mu_sim = pm.Deterministic('mu_sim',theano_wrapper_H(xs,ys,zs,qs));
res = pm.MvNormal('res',
                  mu=mu_sim,
                  cov=case_2681829.R,
                  observed=case_2681829.mu);
step = pm.Metropolis()
trace = pm.sample(10000, step)

My core function H require that arguments xs,ys,zs,qs are scalars, not vectors. Is this code should be ok? Thanks

4

Besides the indentations, this looks reasonable.

5

Thanks a lot for your help. This case works now. But I have the other issue: it is very slow it shows me 20 hours to calculate 4 chains for 1000 samples. The dimension of MvNormal is 48. Could you give me a possible ways to decreasing time of PyMC3 calculations?

6

Probably the slowdown is the time it takes to evaluate H. If so, there’s nothing you can do except re-write H yourself, possibly in terms of native theano elements, to be more efficient.