How to compute unit vector based on RVs?

I’m analyzing time series data with a physics model and one of the parameters is a 2D vector, and I have priors on each of the two components, so I set up the model to sample them independently, and then create the vector from those:

import numpy as np
import pymc3 as pm

with pm.Model():
a0 = pm.Normal(‘a0’,mu=0,sd=1)
a1 = pm.Normal(‘a1’,mu=0,sd=1)
vector = np.array([a0,a1])
unit_vector = vector/np.linalg.norm(vector)

This throws a convoluted error message (below) that I’ve been able to find in other posts, but those don’t seem to have a solution that makes sense for me, and I’m too new to Theano to know how to do this properly. I’m hoping for an easy fix so I can avoid refactoring a lot of code wrapped up in model class functions.


KeyError Traceback (most recent call last)
/anaconda3/lib/python3.7/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)
/anaconda3/lib/python3.7/site-packages/theano/tensor/basic.py in constant(x, name, ndim, dtype)
245 try:
–> 246 ttype = TensorType(dtype=x_.dtype, broadcastable=bcastable)
247 if not constant.enable:

/anaconda3/lib/python3.7/site-packages/theano/tensor/type.py in init(self, dtype, broadcastable, name, sparse_grad)
50 self.broadcastable = tuple(bool(b) for b in broadcastable)
—> 51 self.dtype_specs() # error checking is done there
52 self.name = name

/anaconda3/lib/python3.7/site-packages/theano/tensor/type.py in dtype_specs(self)
271 raise TypeError(“Unsupported dtype for %s: %s”
–> 272 % (self.class.name, self.dtype))
273

TypeError: Unsupported dtype for TensorType: object

During handling of the above exception, another exception occurred:

TypeError Traceback (most recent call last)
/anaconda3/lib/python3.7/site-packages/theano/tensor/basic.py in as_tensor_variable(x, name, ndim)
193 try:
–> 194 return constant(x, name=name, ndim=ndim)
195 except TypeError:

/anaconda3/lib/python3.7/site-packages/theano/tensor/basic.py in constant(x, name, ndim, dtype)
265 except Exception:
–> 266 raise TypeError(“Could not convert %s to TensorType” % x, type(x))
267

TypeError: (‘Could not convert [a0 a1] to TensorType’, <class ‘numpy.ndarray’>)

During handling of the above exception, another exception occurred:

AsTensorError Traceback (most recent call last)
in
3 a1 = pm.Normal(‘a1’,mu=0,sd=1)
4 vector = np.array([a0,a1])
----> 5 unit_vector = vector/np.linalg.norm(vector)

/anaconda3/lib/python3.7/site-packages/theano/tensor/var.py in rtruediv(self, other)
201
202 def rtruediv(self, other):
–> 203 return theano.tensor.basic.true_div(other, self)
204
205 def rfloordiv(self, other):

/anaconda3/lib/python3.7/site-packages/theano/gof/op.py in call(self, *inputs, **kwargs)
613 “”"
614 return_list = kwargs.pop(‘return_list’, False)
–> 615 node = self.make_node(*inputs, **kwargs)
616
617 if config.compute_test_value != ‘off’:

/anaconda3/lib/python3.7/site-packages/theano/tensor/elemwise.py in make_node(self, *inputs)
478 using DimShuffle.
479 “”"
–> 480 inputs = list(map(as_tensor_variable, inputs))
481 out_dtypes, out_broadcastables, inputs = self.get_output_info(
482 DimShuffle, *inputs)

/anaconda3/lib/python3.7/site-packages/theano/tensor/basic.py in as_tensor_variable(x, name, ndim)
198 except Exception:
199 str_x = repr(x)
–> 200 raise AsTensorError(“Cannot convert %s to TensorType” % str_x, type(x))
201
202 # this has a different name, because _as_tensor_variable is the

AsTensorError: (‘Cannot convert [a0 a1] to TensorType’, <class ‘numpy.ndarray’>)

You can define a in single line to avoid the concatenation.
Then, use use pm.Deterministic to define a function of RVs. Something like this (untested):

from theano import tensor as T

with pm.Model():
	a = pm.Normal('a', mu=0, sd=1, shape=2)
	unit_vector = pm.Deterministic('unit_vector', a/T.norm(a))

I think T.sqrt(T.dot(a,a)) is a working denominator there. I’m finding that I’ll probably have to refactor the rest of my code to handle this anyhow, but this helps get me there, so thank you.

I often find that unit-normed parameters are hard to sample because the norm of the input vector is unconstrained. In this specific case, you may want to parametrize with a von-Mises distribution

theta = pm.VonMises('theta')
unit_vector = pm.Deterministic('unit_vector', T.cos(theta), T.sin(theta))
1 Like

+1 to @chartl’s comment. You might also find this discussion useful: How to model sinusoids in white gaussian noise