Implement deterministic nontrivial function of random variables

Questions
1

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

2

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!

3

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…

4

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.*?

5

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...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [yc, xc, alpha, B0, d, f, N0, R, a]
Sampling 4 chains: 100%|██████████| 42000/42000 [05:55<00:00, 118.21draws/s]

in which data.values is a (53, 113) real matrix of solar magnetic field observations. It compiles now, although the traces for all the parameters doesn’t seem to converge (see example for parameter ‘R’ below).
R_trace

But this is of course another problem.

Thank you for the feedback!
Francisco