Running into Theano issues when using DensityDist

Questions
1

Hello,

I am new to PyMC3 and Bayesian statistics in general. I am running into an issue trying to create my own likelihood function (following this example), but am running into some trouble.

Basically, my likelihood functions (I have one version that uses numpy and one theano.tensor) compare data to a forward model synthetic generated with parameters w and h, with uncertainty sigma. Here’s a code snippet:

import matplotlib.pyplot as plt
import pandas as pd
import pymc3 as pm
import theano
import theano.tensor as T

def loglike1_np(n,dat,syn,sigma):
    likelihood = np.zeros(n)
    likelihood = np.abs(np.divide(syn-dat,sigma)) - np.log(2.0*sigma)

    logp = np.sum(-1*likelihood)
    logp = pm.theanof.floatX(logp)
    return logp

def loglike1_tt(n, data,syn,sigma):
    
    likelihood = T.abs_(T.true_div(T.sub(syn-data),sigma)) - T.log(2*sigma)
    
    logp = -1*T.sum(likelihood)
    logp = pm.theanof.floatX(logp)
    return logp

def forward_model(w,h,...

    return synthetic

...

model = pm.Model()
with model:
    numiter = int(nburnin + nsample)
    data = pm.Data('data', amplitudes)
    H = pm.Uniform('h', lower=h1, upper=h2)
    W = pm.Uniform('w', lower=0.1, upper=100.0)
    sigma = pm.HalfCauchy('sigma',beta=2)

    synth = forward_model(W, H, ...)

    likelihood = pm.DensityDist('likelihood', loglike1_np, observed={'n':n, 'data':data, 'syn':synth, 'sigma':sigma})
#    likelihood = pm.DensityDist('likelihood', loglike1_tt, observed={'n':n, 'data':data, 'syn':synth, 'sigma':sigma})

    step = pm.Metropolis()

    trace = pm.sample(numiter, step=step)

    burned_trace = trace[int(nburnin):]

    pm.traceplot(burned_trace)

For both loglike1_ functions, I get this error within the DensityDist():

TypeError: float() argument must be a string or a number, not ‘TensorVariable’

I think this results from my usage of the data or synth variables within the loglike1_() functions, as they are both vectors containing elements of type >theano.tensor.var.TensorVariable

I am at a complete loss on how to resolve this, however. I’ve been beating my head against the table (read: searching for the solution here and on StackOverflow, Theano docs) for about 4 hours now. Any advice is greatly appreciated.

2

Maybe try something like:

def loglike1_tt(n, data,syn,sigma):
    
    likelihood = T.abs_((syn-data) / sigma) - T.log(2*sigma)
    
    logp = -1 * T.sum(likelihood)
    return logp
3

Thank you for the reply. This gives the same error as above. Here’s the full error output:

TypeError: float() argument must be a string or a number, not ‘TensorVariable’

The above exception was the direct cause of the following exception:

Traceback (most recent call last):

File “”, line 1, in
runfile(‘/home/user/pymc3_test/convert/please.py’, wdir=‘/home/user/pymc3_test/convert’)

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py”, line 827, in runfile
execfile(filename, namespace)

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py”, line 110, in execfile
exec(compile(f.read(), filename, ‘exec’), namespace)

File “/home/user/pymc3_test/convert/please.py”, line 340, in
likelihood = pm.DensityDist(“likelihood”, loglike1_tt, observed={‘n’:n,‘data’:amp[0:n],‘syn’:synth,‘sigma’:sigma})

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/pymc3/distributions/distribution.py”, line 83, in new
return model.Var(name, dist, data, total_size, dims=dims)

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/pymc3/model.py”, line 1102, in Var
model=self,

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/pymc3/model.py”, line 1791, in init
for name, data in data.items()

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/pymc3/model.py”, line 1791, in
for name, data in data.items()

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/pymc3/model.py”, line 1660, in as_tensor
data = pandas_to_array(data).astype(dtype)

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/pymc3/model.py”, line 1652, in pandas_to_array
return pm.floatX(ret)

File “/home/user/anaconda3/envs/test/lib/python3.7/site-packages/pymc3/theanof.py”, line 80, in floatX
return X.astype(theano.config.floatX)

ValueError: setting an array element with a sequence.

If I replace the variable synth with a static synthetic computed with random w and h (not PyMC3 RVs), the above code runs the MCMC. If I try to run loglike1_tt() outside of the DensityDist(), I get a different error:

AsTensorError: (‘Cannot convert [Elemwise{mul,no_inplace}.0 Elemwise{mul,no_inplace}.0\n Elemwise{mul,no_inplace}.0 Elemwise{mul,no_inplace}.0\n … repeat for length of vector … Elemwise{mul,no_inplace}.0\n Elemwise{mul,no_inplace}.0] to TensorType’, <class ‘numpy.ndarray’>)

So it seems to me that the error is occurring when trying to convert the TensorVariable synth into a theano.config.floatX, caused by an issue with reading the elements in the vector as Elemwise{mul, no_inplace}.0. Do you have any recommendations? Thanks!

4

Hard to say without knowing the type and shape of your output that goes into the logp function. Maybe you can put up a small notebook with simulated data?

5

Thanks again for the reply, and sorry for the delay. Here is the requested code with simulated data.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pymc3 as pm
import theano
import theano.tensor as T


""" functions """

def forward_model(nf,f,w,h):
    """ Return synthetic spectra given amplitudes """
        
    # constants for setup
    g = 9.80665
    n = 2.4
    re0 = 3280.
    aoa0 = 3.86 ; kval = 2.0 ; cfact = 1.0 ; mfact = 1.513
    vp = 5500.0 ; vs = 3175.0 ; rho = 2550.0
    const = 14.8

    # Intermediate Calculations
    mu = vs*vs*rho
    lamda = vp*vp*rho - 2.0*mu
    beta = (lamda + 2.0*mu)/(4.0*mu)
    nu = lamda/(2.0*(lamda+mu))

    rc = const * (w**0.29) * (h**-0.11) * mfact
    
    h_ft = h*3.28084 # meters to feet
    rho_gpcc = rho/1000.0
    re = re0 * (aoa0**(1.0/n)) * (w**(1.0/3.0)) / ((rho_gpcc*h_ft)**(1.0/n))

    p0s = 1.5*rho*g*h
    p0c = ((4.0*mu)/3.0)*(rc/re)**3.0*cfact   
    alpha = kval*vp/re
    
    omega0 = vp/re
    
    m0 = 4*np.pi*rho*(vp**2)*(re**3)*p0c/(4*mu)
    f0 = omega0/(2*np.pi)
    fc = np.sqrt(p0s/p0c)*vs/(np.pi/re)
    
    m0rp = (3.53e8*(1 + nu)**1.86*(((1-nu)**0.15)*(rho**0.13)*
                    (vp**1.7)*(w**0.87)*(h**-0.33)))
    fcrp = (0.029*(1+nu)**-0.93 *((1-nu)/((1-2*nu)**-0.75)* (rho**0.2267) * 
                   (vp**0.15) * (w**-0.2683) * (h**0.4567)))
    
    # Create the synthetic (spec)
    spec = []
    for i in range(0,nf):
        omega = 2*np.pi*f[i]
        iomega = complex(0,omega)
        p_hat = (p0s - p0c)/(alpha + iomega) + p0c/iomega
        rdp = p_hat*re/(4*mu)*vp**2/(omega**2 + omega0*iomega - beta*omega**2)
        rvp = iomega*rdp
        m0rate = 4*np.pi*rho*vp**2 * abs(rvp)
        spec.append(m0rate)
    spec = np.asarray(spec)
    return spec


def loglike1_np(n,dat,syn,sigma):

    likelihood = np.zeros(n)
    likelihood = np.abs(np.divide(syn-dat,sigma)) - np.log(2.0*sigma)

    logp = np.sum(-1*likelihood)
    logp = pm.theanof.floatX(logp)
    return logp

def loglike1_tt(n, data,syn,sigma):
    
    likelihood = T.abs_((syn - data) / sigma) - T.log(2*sigma)
    
    logp = -1*T.sum(likelihood)
    return logp

""" Begin script """

nburnin=0.5e3; nsample=3e3; # Small number for now
nchains = 2
cores = 1

# frequency list
f = np.ones(616)
f[0:26] = 0.5; f[27:81] = 0.7071; f[172:274] = 1.4142; f[275:369] = 2;
f[370:464] = 2.8284; f[465:547] = 4; f[548:613] = 5.6569; f[614:] = 8;

n = len(f) - 1
h1 = 1.; h2 = 1000.   # h prior lower & upper

# Simulated data in similar range to real data
amp = np.random.uniform(low=14.5, high=16.5, size=(n+1,))

# Begin PyMC3 model
model = pm.Model()

numiter = int(nburnin + nsample)
with model:
    data = pm.Data('data', amp[0:n])
    H = pm.Uniform("h",lower=h1, upper=h2)
    W = pm.Uniform('yield',lower=0.1, upper=100.0)
    sigma = pm.HalfCauchy('sigma',beta=2)
    
    synth = forward_model(n,f,W,H)
    
    likelihood = pm.DensityDist("likelihood", loglike1_tt, 
                                observed={'n':n,'data':data,
                                          'syn':synth,'sigma':sigma})
    step = pm.Metropolis()
    
    trace = pm.sample(numiter, step,chains=nchains, cores=cores)

The forward model is quite large, so I apologize for that. Everything else is as close as possible to original code without being excessive. The synthetic output by the forward model is a numpy array of size (n, ) and the data is of the same shape. n itself is an integer, and sigma is a PyMC3 TransformedRV

6

You probably need to rewrite your forward function, as spec = np.asarray(spec) does not do what you expected (ie. it does not turn a list into theano tensor), this is what causing problem down stream.
Ideal solution is to rewrite the for loop into theano.scan, i am skeptical of the complex(0,omega), so you might need to rewrite it with separating the real and complex part.
Other quick fix is to wrap the forward function with theano.as_op, but it will likely be pretty slow.

2 Likes
7

Thanks so much!