You can do at least part of it in vectorized form.
The following deliberately creates random numbers in the same order as if you had run your function over multiple iterations
N = length(x);
rands = rand(3,3,2,N);
mu_rc = reshape(rands(:,:,1,:), 3, []); %arrange them as 3 x 3 beside each other
sigs = reshape(rands(:,:,2,:), 3, []); %arrange them as 3 x 3 beside each other
ct = cos(t(:));
st = sin(t(:));
g_r1 = reshape([ct, -st, x(:)].', 1, []);
g_r2 = reshape([st, ct, y(:)].', 1, []);
g_r3 = repmat([0 0 1], 1, N);
g_r = [g_r1; g_r2; g_r3];
mu_c = mat2cell(mu_rc, 3, 3*ones(1,N));
g_c = mat2cell(g_r, 3, 3*ones(1,N));
mu_blk = blkdiag(mu_c{:});
g_blk = blkdiag(g_c{:});
mug_blk = mu_blk \ g_blk;
mug_c = mat2cell(mug_blk, 3*ones(1,N), 3*ones(1,N));
y_mc = cellfun(@logm, mug_c(1:N+1:end), 'uniform', 0);
At this point, y_mc is a cell array of the y_m results. What remains is to do the vectorized calculation of the remaining two steps, namely
y_v = [y_m(1,3);y_m(2,3);y_m(2,1)];
Prob = exp(-0.5*y_v'/sig*y_v);
I am not going to do this at the moment; I am tired and need to do some other things.
