Inquiry about how to define a function as pm.Uniform in pymc3

development
theano

#1

Hi
I have a question when I using the with pm.model(). I defined a function, and I want to use the return value of ff function on the pm.model(). But it shows an error ‘we expected 3 inputs but got 5’

@theano.compile.ops.as_op(itypes=[tt.lscalar, itypes=[tt.lscalar], itypes=[tt.lscalar], otypes=[tt.dvector])
def ff(vx,vy,sx,sy,ro):
    x = 0.1
    y= 0.1
    return 1/2/np.pi/sx/sy/tt.sqrt(1-ro**2)*tt.exp(-1/2/(1-ro)**2*((x-vx)**2/sx**2+(y-vy)**2/sy**2-2*ro*(x-vx)*(y-vy)/sx/sy))

with pm.Model() as converge_model: 
    vx_min = -5
    vx_max = 5
    vy_min = -5
    vy_max = 5
    Dx_min = 0.01
    Dx_max = 0.89
    Dy_min = 0.01
    Dy_max = 0.89
    ro_min = -0.999
    ro_max = 0.999

    vx = pm.Uniform('vx', lower=vx_min, upper=vx_max)
    Dx = pm.Uniform('Dx', lower=Dx_min, upper=Dx_max)
    vy = pm.Uniform('vy', lower=vx_min, upper=vx_max)
    Dy = pm.Uniform('Dy', lower=Dy_min, upper=Dy_max)
    ro = pm.Uniform('ro', lower=ro_min, upper=ro_max)  
    mu=tt.stack(vx,vy)
    cov = tt.stack([[Dx**Dx,ro*Dx*Dy],[ro*Dx*Dy,Dy**2]])
    x = pm.MvNormal('x',mu=mu,cov=cov,shape=(2,2),observed=xobser)

    ff = pm.Deterministic('ff',ff(vx,vy,Dx,Dy,ro))#,1/2/np.pi/Dx/Dy/tt.sqrt(1-ro**2)*tt.exp(-1/2/(1-ro)**2*((x[0]-vx)**2/Dx**2+(x[1]-vy)**2/Dy**2-2*ro*(x[0]-vx)*(x[1]-vy)/Dx/Dy)))
    C = pm.Exponential('C',lam=1/ff,observed=Cobser)  

How can I cooperated the function of ff into the C and ff?

Thanks for your help!!


#2

You dont need to wrap the function in as_op, this is a theano function.


#3

Hi Sir
Thanks for your reply.
Now it is working when I put x = 0.1 y=0.1 for the whole function.
But in realistic, the x and y are an array from the real field data. For example, x=[1,4,7,3,9] y=[2,4,8,9,0]. So how can I transfer the value of x,y to the ff and C function?
For the line:

 ff = pm.Deterministic('ff',1/2/np.pi/Dx/Dy/tt.sqrt(1-ro**2)*tt.exp(-1/2/(1-ro)**2*((0.1-vx)**2/Dx**2+(0.1-vy)**2/Dy**2-2*ro*(0.1-vx)*(0.1-vy)/Dx/Dy))).

When I put the real value 0.1=x, 0.1=y, the whole code works.
But when I wrote

 ff = pm.Deterministic('ff',1/2/np.pi/Dx/Dy/tt.sqrt(1-ro**2)*tt.exp(-1/2/(1-ro)**2*((x[0]-vx)**2/Dx**2+(x[1]-vy)**2/Dy**2-2*ro*(x[0]-vx)*(x[1]-vy)/Dx/Dy)))

It didn’t work.
I want to get a result for ff is:

(1/2/np.pi/Dx/Dy/tt.sqrt(1-ro**2)*tt.exp(-1/2/(1-ro)**2*((x[0]-vx)**2/Dx**2+(x[0]-vy)**2/Dy**2-2*ro*(y[0]-vx)*(y[0]-vy)/Dx/Dy)*(1/2/np.pi/Dx/Dy/tt.sqrt(1-ro**2)*tt.exp(-1/2/(1-ro)**2*((x[1]-vx)**2/Dx**2+(x[1]-vy)**2/Dy**2-2*ro*(y[1]-vx)*(y[1]-vy)/Dx/Dy)

How can I realize it? Can I write a loop on with models? Or can I put the value of x,y as a parameter as vc, Dx?

Thanks!!
My new code is below:

with pm.Model() as converge_model:

    vx_min = -5
    vx_max = 5
    vy_min = -5
    vy_max = 5
    Dx_min = 0.01
    Dx_max = 0.89
    Dy_min = 0.01
    Dy_max = 0.89
    ro_min = -0.999
    ro_max = 0.999

    vx = pm.Uniform('vx', lower=vx_min, upper=vx_max)
    Dx = pm.Uniform('Dx', lower=Dx_min, upper=Dx_max)
    vy = pm.Uniform('vy', lower=vx_min, upper=vx_max)
    Dy = pm.Uniform('Dy', lower=Dy_min, upper=Dy_max)
    ro = pm.Uniform('ro', lower=ro_min, upper=ro_max)  
    mu=tt.stack(vx,vy)
    cov = tt.stack([[Dx**Dx,ro*Dx*Dy],[ro*Dx*Dy,Dy**2]])
    x = pm.MvNormal('x',mu=mu,cov=cov,shape=(2,2),observed=xobser)
    ff = pm.Deterministic('ff',1/2/np.pi/Dx/Dy/tt.sqrt(1-ro**2)*tt.exp(-1/2/(1-ro)**2*((x[0]-vx)**2/Dx**2+(x[1]-vy)**2/Dy**2-2*ro*(x[0]-vx)*(x[1]-vy)/Dx/Dy)))
    C = pm.Exponential('C',lam=1/ff,observed=Cobser)

#4

What do you mean? The function seems to run, and I dont get any error running your code (without the Xobser). What kind of result you are expecting?


#5

Okaii. I think it works. Thanks.