Getting an "Error using * " message when mixing scalars and vectors in an anonymous function

I'm trying to solve a multivariate normal with a positive semi-definite correlation matrix, but mvncdf only works for positive definite. Therefore, I am trying to solve the anonymous function for a 3 variable multivariate standard normal with a positive semi-definite correlation matrix.

SIGMA = [1 1 1; 1 1 1; 1 1 1] ;
my_tri = @(x,y,z) (2*pi).^(-3./2).*exp(-0.5.*([x y z] * SIGMA * [x y z]')) ;
my_int = integral3(my_tri,-Inf,-2.3263,-Inf,-2.0537,-Inf,-1.8808) ;

However, my "my_tri" function gives me the error message: Error using * Inner matrix dimensions must agree.

If I try the following though: (2.*pi).^(-3./2).*exp(-0.5.*([1 2 3] * ones(3) * [1 2 3]')) ans =

   9.6701e-10

How do I solve this error which is not an error (clearly my programming error)? Thanks.

Note: if I try to use the trick [x' y' z'] (or combine it with my_tri = ()') I get the: Error using integralCalc/finalInputChecks (line 515) Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

[dpb]

The anonymous function definition is ok for the input as specified as a triplet--

>> tri = @(x,y,z) (2*pi)^(-3/2)*exp(-0.5.*([x y z]*SIGMA*[x y z]'))
tri = 
  @(x,y,z)(2*pi)^(-3/2)*exp(-0.5.*([x,y,z]*SIGMA*[x,y,z]'))
>> tri(1,2,3)
ans =
 9.6701e-10
>>

So, the issue has to be that the function isn't called that way inside integral3. I don't have time at the moment; your mission, should you choose to accept it, is to set a breakpoint inside integral3 and see how the arguments to the function are actually built internally and passed to the function and recast it in that form.

Answer 1

[Matthew]

The issue is that integral3 calls the function with a row array. You can fix your function to work with integral3 by either stacking the rows, or by converting the rows into columns. I've only tested option 1, but both options should work with a bit of tweaking.

%Option 1
my_tri = @(x,y,z) (2*pi).^(-3./2).*exp(-0.5.*([x; y; z]' * SIGMA * [x; y; z])) ;

%Option 2
my_tri = @(x,y,z) (2*pi).^(-3./2).*exp(-0.5.*([x' y' z'] * SIGMA * [x' y' z']')) ;

Note that using the below makes diagnosing this kind of thing pretty easy.

dbstop if caught error

Also, make sure you are calling my_tri - your example has 'my_try' being called instead, which may give you different behavior depending on how it is defined.