How to vectorize a for loop that computes the expected value with multidimensional arrays?

I am working on a dynamic programming problem and after profiling my code I noticed that most of the computing time is spent calculating expected values. The problem seems to be a loop that multiplies each element of a multidimensional array by each element of a transition probability array and then add the results. I tried the different ways of multiplication (.* and times) and the sum commands (e.g. sum(sum(A)) and sum(A(:)) but there is not significant improvement. I have not been able to find an adequate way to vectorize the loop or convert it to a parfor loop. Any suggestions would be very appreciated. Please let me know if the description/code is not clear. Thank you.

% 
u_s = omega_rev.*revenue';   % Realization of the stochastic component
for ii = 1:1:no_observations                
  for lu = 1:1:4               
      u_j_i(ii,lu) = u_s(lu) + matrix_k(ii,lu);
      EV(ii,lu) = sum(sum(sum(Vinit(:,:,:,next_age(lu),lu,ii).* P_array(:,:,:))));                                    
      exp_v_j_i(ii,lu) = exp(u_j_i(ii,lu) + disc_fact*EV(ii,lu));
      % Alternative way to compute the EV
%       EV_tmp = times(Vinit(:,:,:,next_age(lu),lu,ii), P_array(:,:,:));
%       EV(ii,lu) = sum(EV_tmp(:));
  end
end

Answer 1

[Matt J]

    EV=zeros(no_observations,4);
    for lu=1:4

      V = reshape(Vinit(:,:,:,next_age(lu),lu,:),[],no_observations);

      EV(:,lu)=(P_array(:).' * V).';

    end

    u_j_i=bsxfun(@plus, u_s(:).', matrix_k);

    exp_v_j_i = exp(u_j_i + disc_fact*EV);

[Raymundo Marcos-Martinez]

Thank you very much Matt. Your faster version of the code motivates me to continue learning the Matlab language. Thanks again.