Vectorize an anonymous function

I was wondering whether someone would be able to help/point me in the right direction for vectorizing a anonymous function?

I am trying to fit a multivariate gaussian distribution to data. Here is my code:

for n = 1:length(uniqPosOne)
    gaussElpt = @(param) (param(1) + param(2).*(exp(-1/2 .* (uniqPosOne(n, :) - [param(3), param(4)]) * ...
        ([param(5), param(6); param(6), param(7)]).^(-1) * (uniqPosOne(n, :) - [param(3), param(4)]).')));
    
    a = gaussElpt(startPoint);
    tempList(n) = a;
end

I am trying to optimize the parameters that allow for me to fit my data (not shown), I use lsqnonlin.

uniqPosOne is a n x 2 array, params(1) and (2) are offset/gain, params(3) and (4) are vector means, and params(5), (6), (7) are part of covaraince matrix. The params are the parameters used by lsqnonlin to find the parameters with the lowest error.

I want to be able to have a nx2 array (uniqPosOne) and have the function evaluate each pair separately in a vectorized manner, currently I have been unable to do so as attempting vectorization results in a nxn, while I want a nx1.

The main reason is due to time to compute is rather long for a small subset of data. Any help is appreciated.

[Matt J]

If this is supposed to be a Gaussian distirbution, with covariance [param(5), param(6); param(6), param(7)] then this,

([param(5), param(6); param(6), param(7)]).^(-1) 

should reallly be this,

inv([param(5), param(6); param(6), param(7)])

[AES]

Thank you!

[Walter Roberson]

Or, better yet,

[param(5), param(6); param(6), param(7)] \ (uniqPosOne(n, :) - [param(3), param(4)]).'

Answer 1

[Matt J]

Just replace all the uniqPosOne(n,:) with uniqPosOne:

gaussElpt = @(param) (param(1) + param(2).*(exp(-1/2 .* (uniqPosOne - [param(3), param(4)]) * ...
        ([param(5), param(6); param(6), param(7)]).^(-1) * (uniqPosOne - [param(3), param(4)]).')));
    

[AES]

I have tried that previously and it produces an nxn, when I am interested in a nx1. Is there a way to vectorize it so that I can get a nx1?

[Walter Roberson]

@Ezra Sutter

Because of the matrix multiplications, you would have to recode it using pagemtimes . That would probably require that you permute() some of your matrices so that data was in the 3rd dimension.

[Matt J]

n=10;

uniqPosOne=rand(n,1);

gaussElpt=@(param) gauss2dEval(param,uniqPosOne);

gaussElpt(rand(1,7))

ans = 10×1

1.6660 1.6798 1.6847 1.6562 1.6129 1.6887 1.6103 1.5560 1.4912 1.4979

function out=gauss2dEval(param,uniqPosOne)

 param=param(:).'; %ensure row

 Sigma=[param(5), param(6); param(6), param(7)];
 dev=(uniqPosOne-param(3:4))';
 
 out=param(1)+param(2)*exp(-0.5*sum((Sigma\dev).*dev) ).' ;

end

[AES]

Thank you, this is what I was looking to do. If you have time could you breifly explain, or point me towards some resources that may explain, this operation: sum((Sigma\dev).*dev)?

[Matt J]

Sigma\dev is the more recommended way of doing inv(Sigma)*dev.

See mldivide, \

See also Walter's comment above.

Answer 2

[Matt J]

This FEX submission already appears to do what you want,

https://www.mathworks.com/matlabcentral/fileexchange/69116-gaussfitn?s_tid=srchtitle