I'm trying to evaluate a multivariate normal CDF several times using theano's scan function, but I'm getting a ValueError.
Here's an example of the original function I'm trying to vectorize:
from scipy.stats.mvn import mvnun # library that calculates MVN CDF
low = [-1.96, 0 ] # lower bounds of integration
upp = [0 , 1.96] # upper bounds of integration
mean = [0 , 0 ] # means of the jointly distributed random variables
covs = [[1,0.25],[0.25,1]] # covariance matrix
print(mvnun(low,upp,mean,cov))
This produces the following output:
(0.19620339269649473, 0)
Simple and straightforward, right?
What I'm really trying to do is create 4 large input objects with 1500 elements each. That way, I get to evaluate the mvnun function 1500 times. The idea is that at each iteration, all inputs are different from the last, and no information from the previous iteration is necessary.
Here is my setup:
import theano
import numpy as np
lower = theano.tensor.dmatrix("lower") # lower bounds - dim: 1500 x 2
upper = theano.tensor.dmatrix("upper") # upper bounds - dim: 1500 x 2
means = theano.tensor.dmatrix("means") # means means - dim: 1500 x 2
covs = theano.tensor.dtensor3("covs") # cov matrices - dim: 1500 x 2 x 2
results, updates = theano.scan(fn=mvnun,
sequences=[lower,upper,means,covs])
f = theano.function(inputs=[lower, upper, means, covs],
outputs=results,
updates=updates)
However, when I try to run this block of code, I get an error on the line with the scan command. The error states: ValueError: setting an array element with a sequence.. The full traceback of the error is below:
Traceback (most recent call last):
File "", line 7, in sequences=[lower,upper,means,covs])
File "C:\Anaconda2\lib\site-packages\theano\scan_module\scan.py", line 745, in scan condition, outputs, updates = scan_utils.get_updates_and_outputs(fn(*args))
ValueError: setting an array element with a sequence.
I originally thought that the code wasn't working because the mvnun function returns a two-element tuple instead of a single value.
However, when I tried to vectorize a test function (that I created) that also returned a two-element tuple, things worked just fine. Here's the full example:
# Some weird crazy function that takes in three Nx1 vectors
# and an NxN matrix and spits out a tuple of scalars.
def test_func(low_i,upp_i,mean_i,cov_i):
r1 = low_i.sum() + upp_i.sum()
r2 = np.dot(mean_i,cov_i).sum()
test_func_out = (r1,r2)
return(test_func_out)
lower = theano.tensor.dmatrix("lower") # lower
upper = theano.tensor.dmatrix("upper") # upper
means = theano.tensor.dmatrix("means") # means
covs = theano.tensor.dtensor3("covs") # covs
results, updates = theano.scan(fn=test_func,
sequences=[lower,upper,means,covs])
f = theano.function(inputs=[lower, upper, means, covs],
outputs=results,
updates=updates)
np.random.seed(666)
obs = 1500 # number of elements in the dataset
dim = 2 # dimension of multivariate normal distribution
# Generating random values for the lower bounds, upper bounds and means
lower_vals = np.random.rand(obs,dim)
upper_vals = lower_vals + np.random.rand(obs,dim)
means_vals = np.random.rand(obs,dim)
# Creates a symmetric matrix - used for the random covariance matrices
def make_sym_matrix(dim,vals):
m = np.zeros([dim,dim])
xs,ys = np.triu_indices(dim,k=1)
m[xs,ys] = vals[:-dim]
m[ys,xs] = vals[:-dim]
m[ np.diag_indices(dim) ] = vals[-dim:]
return m
# Generating the random covariance matrices
covs_vals = []
for i in range(obs):
cov_vals = np.random.rand((dim^2 - dim)/2+dim)
cov_mtx = make_sym_matrix(dim,cov_vals)
covs_vals.append(cov_mtx)
covs_vals = np.array(covs_vals)
# Evaluating the test function on all 1500 elements
print(f(lower_vals,upper_vals,means_vals,covs_vals))
When I run this block of code, everything works out fine and the output I get is a list with 2 arrays, each containing 1500 elements:
[array([ 4.24700864, 3.80830129, 2.60806493, ..., 3.12995381, 4.41907055, 4.12880839]),
array([ 0.87814314, 1.01768617, 0.45072405, ..., 1.15788282, 0.15766754, 1.32393402])]
It's also worth noting that the order in which the vectorized function is getting elements from the sequences is perfect. I ran a sanity check with the first 3 numbers in the list:
for i in range(3):
print(test_func(lower_vals[i],upper_vals[i],means_vals[i],covs_vals[i]))
And the results are:
(4.2470086396797502, 0.87814313729162796)
(3.808301289302495, 1.017686166097616)
(2.6080649327828564, 0.45072405177076169)
These values are practically identical to the first 3 output values in the vectorized approach.
So back to the main problem: why can't I get the mvnun function to work when I use it in the scan statement? Why am I getting this odd ValueError?
Any kind of advice would be really helpful!!!
Thanks!!!
