# Implement deterministic nontrivial function of random variables

Hi!

I am trying to implement a nontrivial function of several random variables. I used to write,

``````@T.compile.ops.as_op(itypes=[Tensor.dvector],otypes=[Tensor.dvector])
def fun_name(x = x):
complex function
``````

but I am having problems since now I am implementing a more complex function of several variables. I have also tried to implement the function using this sintaxis,

``````with pm.Model() as model:
det_funct = pm.Deterministic('function', function(x1, x2, x3))
``````

but I have problems also here, probably related to the fact the function is not in Theano sintaxis.

The question are:

1. Which is the recommended way to implement a deterministic function of complex variables in pymc3?
2. Are there similar examples to show? I dig into a similar unanswered question in this forum 6 months ago.

Any feedback is appreciated

Francisco

I think I should extend the answer with more specifics.

``````with pm.Model() as model:
a = pm.Uniform('a', 0, 60)
R = pm.Uniform('R', 0, 60)
N0 = pm.Uniform('N0', -2, 2)
f = pm.Uniform('f', -1, 1)
d = pm.Uniform('d', 0, 50)
B0 = pm.Uniform('B0', 0, 1000)
alpha = pm.Uniform('alpha', -np.pi/2, np.pi/2)
xc = pm.Uniform('xc', 0, 113)
yc = pm.Uniform('yc', 0, 53)
``````

#This is the problem I must solve, strategy #1 = using pm.Deterministic. Please note that xv and yv are constant matrices define outside the function. So the output mag is a matrix of the same dimension.

``````def torusmodel(input):
a, R, N0, f, d, B0, alpha, xc, yc = input

x=(xv-xc+0.5)*np.cos(alpha) + (yv-yc+0.5)*np.sin(alpha)
y=-(xv-xc+0.5)*np.sin(alpha) + (yv-yc+0.5)*np.cos(alpha)
xr=np.sqrt(x**2+d**2) - R
rho=np.sqrt(xr**2 + y**2)
u=np.sqrt(x**2 + d**2)
costh=xr/rho
Nt=N0*(1+f*(rho/a)**2)
mag = (x*(1+g*(costh)) - 2*Nt*d*y/u)*(-1)*B0*np.exp((-1)*(rho/a)**2)/u

return mag

with model:
torusmodel_pm = pm.Deterministic('torusmodel_pm', torusmodel(a,R,N0,f,g,d,B0,alpha,xc,yc))

with model:
observations = pm.Normal( "obs",  mu=torusmodel_pm, tau=1/0.1, observed=data.values)
``````

#########################

#This is strategy #2 = using decorators. Please note that xv and yv are constant matrices define outside the function. So the output mag is a matrix of the same dimension.

``````@T.compile.ops.as_op(itypes=[T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector],otypes=[T.dvector])
def torusmodel(input):
a, R, N0, f, d, B0, alpha, xc, yc = input

x=(xv-xc+0.5)*np.cos(alpha) + (yv-yc+0.5)*np.sin(alpha)
y=-(xv-xc+0.5)*np.sin(alpha) + (yv-yc+0.5)*np.cos(alpha)
xr=np.sqrt(x**2+d**2) - R
rho=np.sqrt(xr**2 + y**2)
u=np.sqrt(x**2 + d**2)
costh=xr/rho
Nt=N0*(1+f*(rho/a)**2)
mag = (x*(1+g*(costh)) - 2*Nt*d*y/u)*(-1)*B0*np.exp((-1)*(rho/a)**2)/u

return mag

with model:
``````

`observations = pm.Normal( "obs", mu=torusmodel, tau=1/0.1, observed=data.values)`

Both approaches fail, with obscure errors. At this point, I am not sure if this is related to: (1) the definition of the deterministic function or (2) the fact that I am building matrices of PyMC3 random variables. In any case, this simple approach do work.

``````# this dummy example that works fine
def torusmodel_dummy(a,R,N0,f,g,d,B0,alpha,xc,yc):
return a + R + N0 + f + g + d + B0 + alpha + xc + yc

with model:
torusmodel_pm = pm.Deterministic('torusmodel_pm', torusmodel_dummy(a,R,N0,f,g,d,B0,alpha,xc,yc))
``````

Any feedback is welcomed!

It would be helpful if you posted the exact code that results in an error, complete with (possibly randomly generated) values for xv and yv. Until then, nobody can actually run your code, making it harder to help…

Looking at the function `torusmodel` itself it doesnt look too difficult to rewrite it into theano - have you try `import theano.tensor as tt` and just replace `np.*` with `tt.*`?

Solved! When I was preparing the submission I found the error in the matrix dimensions.

The complete code is here,

``````import numpy as np
import pymc3 as pm
import theano.tensor as tt

with pm.Model() as model:
a = pm.Uniform('a', 0, 60)
R = pm.Uniform('R', 0, 60)
N0 = pm.Uniform('N0', -2, 2)
f = pm.Uniform('f', -1, 1)
d = pm.Uniform('d', 0, 50)
B0 = pm.Uniform('B0', 0, 1000)
alpha = pm.Uniform('alpha', -np.pi/2, np.pi/2)
xc = pm.Uniform('xc', 0, 113)
yc = pm.Uniform('yc', 0, 53)

sz1 = 112 #113
sz2 = 52 #53
nx, ny = (sz1, sz2)
x1 = np.linspace(0, sz1, sz1+1)
y1 = np.linspace(0, sz2, sz2+1)
xv, yv = np.meshgrid(x1, y1)

def torusmodel(inputParam):
a, R, N0, f, d, B0, alpha, xc, yc = inputParam

x=(xv-xc+0.5)*tt.cos(alpha) + (yv-yc+0.5)*tt.sin(alpha)
y=-(xv-xc+0.5)*tt.sin(alpha) + (yv-yc+0.5)*tt.cos(alpha)
xr=tt.sqrt(x**2+d**2) - R
rho=tt.sqrt(xr**2 + y**2)
u=tt.sqrt(x**2 + d**2)
costh=xr/rho
Nt=N0*(1+f*(rho/a)**2)
mag = (x*(1+g*(costh)) - 2*Nt*d*y/u)*(-1)*B0*tt.exp((-1)*(rho/a)**2)/u

return mag

inputParam = (10,30,0.5,1,0,15,1000,10,50,30)

with model:
torusmodel_pm = pm.Deterministic('torusmodel_pm', torusmodel(inputParam))

with model:
observations = pm.Normal( "obs",  mu=torusmodel_pm, tau=1/0.1, observed=data.values)

with model:
trace = pm.sample(10000)

Auto-assigning NUTS sampler... 