How to do linear fitting with pymc3 by using pm.DensityDist

Questions
1

Hi

I’m trying to do linear fitting by using pymc3. I want to use pm.DensityDist function to make custom likelihood in the future. The linear fitting by using pm.DensityDist is a kind of practice for that.

my code is

import numpy as np
import matplotlib.pylab as plt

points = 40
x = np.random.uniform(-2,2,points)
x = np.sort(x)

mu, sigma = 0, 1.0
noise = np.random.normal(mu,sigma,x.size)

beta = np.array([1.5, 0.5])
x1 = np.array([x**1, x**0])
y = beta.dot(x1)+noise

import pymc3 as pm
import theano.tensor as tt

def likelihood(alpha,beta,sigma,x,y):
    x1 = np.array([x**1,x**0])
    const = tt.log(2*np.pi*sigma**2)*(-x.size/2)
    w = tt._shared(np.array([beta, alpha]))
    log_exp = -1/(2*sigma**2)*tt.dot((tt._shared(y)-tt.dot(tt._shared(x1.T),w) ),(tt._shared(y)-tt.dot(tt._shared(x1.T),w)))
    return const+log_exp

with pm.Model() as model:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10, shape = 2)
    sigma = pm.HalfNormal('sigma',sd = 10)

    like = pm.DensityDist('like', likelihood, observed=dict(alpha=alpha,beta=beta,
                                                            sigma = sigma,
                                                            x=x,
                                                            y=y))
    start = pm.find_MAP(model = model)
    step = pm.Metropolis()
    trace = pm.sample(10000,step=step,start=start)

and there is error like

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~/anaconda3/lib/python3.6/site-packages/theano/tensor/type.py in dtype_specs(self)
    268                 'complex64': (complex, 'theano_complex64', 'NPY_COMPLEX64')
--> 269             }[self.dtype]
    270         except KeyError:

KeyError: 'object'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
<ipython-input-9-fba3eee669e0> in <module>()
      7                                                             sigma = sigma,
      8                                                             x=x,
----> 9                                                             y=y))
     10     start = pm.find_MAP(model = model)
     11     step = pm.Metropolis()

~/anaconda3/lib/python3.6/site-packages/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
     35             total_size = kwargs.pop('total_size', None)
     36             dist = cls.dist(*args, **kwargs)
---> 37             return model.Var(name, dist, data, total_size)
     38         else:
     39             raise TypeError("Name needs to be a string but got: {}".format(name))

~/anaconda3/lib/python3.6/site-packages/pymc3/model.py in Var(self, name, dist, data, total_size)
    819             with self:
    820                 var = MultiObservedRV(name=name, data=data, distribution=dist,
--> 821                                       total_size=total_size, model=self)
    822             self.observed_RVs.append(var)
    823             if var.missing_values:

~/anaconda3/lib/python3.6/site-packages/pymc3/model.py in __init__(self, name, data, distribution, total_size, model)
   1341         self.missing_values = [datum.missing_values for datum in self.data.values()
   1342                                if datum.missing_values is not None]
-> 1343         self.logp_elemwiset = distribution.logp(**self.data)
   1344         # The logp might need scaling in minibatches.
   1345         # This is done in `Factor`.

<ipython-input-7-ee4be7b8b0a0> in likelihood(alpha, beta, sigma, x, y)
      2     x1 = np.array([x**1,x**0])
      3     const = tt.log(2*np.pi*sigma**2)*(-x.size/2)
----> 4     w = tt._shared(np.array([beta, alpha]))
      5     log_exp = -1/(2*sigma**2)*tt.dot((tt._shared(y)-tt.dot(tt._shared(x1.T),w) ),(tt._shared(y)-tt.dot(tt._shared(x1.T),w)))
      6     return const+log_exp

~/anaconda3/lib/python3.6/site-packages/theano/tensor/sharedvar.py in tensor_constructor(value, name, strict, allow_downcast, borrow, broadcastable, target)
     50     if broadcastable is None:
     51         broadcastable = (False,) * len(value.shape)
---> 52     type = TensorType(value.dtype, broadcastable=broadcastable)
     53     return TensorSharedVariable(type=type,
     54                                 value=np.array(value, copy=(not borrow)),

~/anaconda3/lib/python3.6/site-packages/theano/tensor/type.py in __init__(self, dtype, broadcastable, name, sparse_grad)
     49         # True or False
     50         self.broadcastable = tuple(bool(b) for b in broadcastable)
---> 51         self.dtype_specs()  # error checking is done there
     52         self.name = name
     53         self.numpy_dtype = np.dtype(self.dtype)

~/anaconda3/lib/python3.6/site-packages/theano/tensor/type.py in dtype_specs(self)
    270         except KeyError:
    271             raise TypeError("Unsupported dtype for %s: %s"
--> 272                             % (self.__class__.__name__, self.dtype))
    273 
    274     def to_scalar_type(self):

TypeError: Unsupported dtype for TensorType: object

I really want to do that with DensityDist
Can you help me?

2

beta and alpha are already theano tensors, so there is no need to transform them back and forth between numpy and theano.
Try: w = tt.stack([beta, alpha])

3

Also, FYI, we do not recommend using pm.find_MAP and Metropolis sampler. You should replace

    start = pm.find_MAP(model = model)
    step = pm.Metropolis()
    trace = pm.sample(10000,step=step,start=start)

with

trace = pm.sample(1000, tune=1000)
4

I did that as you said but there is still some error. like this

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-8-762927ba1f28> in <module>()
      7                                                             sigma = sigma,
      8                                                             x=x,
----> 9                                                             y=y))
     10     trace = pm.sample(1000, tune=1000)

~/anaconda3/lib/python3.6/site-packages/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
     35             total_size = kwargs.pop('total_size', None)
     36             dist = cls.dist(*args, **kwargs)
---> 37             return model.Var(name, dist, data, total_size)
     38         else:
     39             raise TypeError("Name needs to be a string but got: {}".format(name))

~/anaconda3/lib/python3.6/site-packages/pymc3/model.py in Var(self, name, dist, data, total_size)
    819             with self:
    820                 var = MultiObservedRV(name=name, data=data, distribution=dist,
--> 821                                       total_size=total_size, model=self)
    822             self.observed_RVs.append(var)
    823             if var.missing_values:

~/anaconda3/lib/python3.6/site-packages/pymc3/model.py in __init__(self, name, data, distribution, total_size, model)
   1341         self.missing_values = [datum.missing_values for datum in self.data.values()
   1342                                if datum.missing_values is not None]
-> 1343         self.logp_elemwiset = distribution.logp(**self.data)
   1344         # The logp might need scaling in minibatches.
   1345         # This is done in `Factor`.

<ipython-input-7-45c36df08df6> in likelihood(alpha, beta, sigma, x, y)
      3     const = tt.log(2*np.pi*sigma**2)*(-x.size/2)
      4     w = tt.stack([beta, alpha])
----> 5     log_exp = -1/(2*sigma**2)*tt.dot((tt._shared(y)-tt.dot(tt._shared(x1.T),w) ),(tt._shared(y)-tt.dot(tt._shared(x1.T),w)))
      6     return const+log_exp
      7 

~/anaconda3/lib/python3.6/site-packages/theano/tensor/sharedvar.py in tensor_constructor(value, name, strict, allow_downcast, borrow, broadcastable, target)
     43 
     44     if not isinstance(value, np.ndarray):
---> 45         raise TypeError()
     46 
     47     # if no broadcastable is given, then the default is to assume that

TypeError:

and Thank you for your fast and kind reply. Can you help me more?

5

The likelihood function which I want to make is

log( p(y|X,w)) = log(2\pi\sigma^2) \times (-n/2) - 1/(2\sigma^2) \times |\vec{y}-X^{T}*\vec{w}|^2
6

The error occur as you dont need to cast numpy array input into theano array of a DensityDist function, as they are already cast into theano by default.

I did not check your logp function in detail, but just fix the type error and shape error:

def likelihood(alpha, beta, sigma, x1, y):
    const = tt.log(2*np.pi*sigma**2)*(-x.size/2)
    w = tt.concatenate([beta, alpha], axis=-1)
    log_exp = -1/(2*sigma**2)*tt.dot((y-w.dot(x1)),
                                     (y-w.dot(x1)))
    return const+log_exp


with pm.Model() as model:
    alpha = pm.Normal('alpha', mu=0, sd=10,  shape=1)
    beta = pm.Normal('beta', mu=0, sd=10, shape=1)
    sigma = pm.HalfNormal('sigma', sd=10)

    like = pm.DensityDist('like', likelihood, observed=dict(alpha=alpha,
                                                            beta=beta,
                                                            sigma=sigma,
                                                            x1=x1,
                                                            y=y))

    trace = pm.sample()
7

Dear junpenglao

PyMC3 is very wonderful package for mcmc but It is still hard to use for me.

Thank you so much for your help. It is very helpful to me.

1 Like