How to code for marginal and conditional distribution involving nine variables?

Hi Everybody. I am trying to obtain a marginal and conditional distribution for the a reservoir control problem. Just to give you the background, following is the maths involved:

where

Using Bayes' Theorem, the joint PDf of can be written as:

where Mis the length of available data,. (Sorry for seemingly too much maths)

This implies that

because α and τ are independent and every member of these sets is also independent of the others.

The marginal and conditional PDFs can be obtained using the following equations (last two equations, I promise):

Here pstr = posterior

Now, finally to the coding. Here's what I have tried to do:

%Noise
mu = 0; % mean
sigma = 1e-2; % standard deviation

% Defining coefficients distribution

mu_1 = 5e-6; % mean
sigma_1 = 1e-6; % standard deviation

mu_2 = 8e-6; % mean
sigma_2 = 1e-6; % standard deviation

mu_3 = 2e-6; % mean
sigma_3 = 1e-6; % standard deviation

mu_4 = 8e-6; % mean
sigma_4 = 1e-6; % standard deviation

mu_5 = 5e-6; % mean
sigma_5 = 1e-6; % standard deviation

mu_6 = 5.5e-6; % mean
sigma_6 = 1e-6; % standard deviation


min_x = min([mu_1, mu_2, mu_3, mu_4, mu_5, mu_6]); % start of query points
max_x = max([mu_1, mu_2, mu_3, mu_4, mu_5, mu_6]); % end of query points
max_std = max([sigma_1, sigma_2, sigma_3, sigma_4, sigma_5, sigma_6]); % maximum standard deviation

x = min_x-4*max_std:1e-9:max_x+4*max_std; % query points for the PDFs

pdf_1 = getPDF(x, mu_1, sigma_1);%getPDF returns for the query points x from the PDF defined by mu and sigma
pdf_2 = getPDF(x, mu_2, sigma_2);
pdf_3 = getPDF(x, mu_3, sigma_3);
pdf_4 = getPDF(x, mu_4, sigma_4);
pdf_5 = getPDF(x, mu_5, sigma_5);
pdf_6 = getPDF(x, mu_6, sigma_6);


% Defining delays' distributions
t1_l = 6; % minimum 1st delay
t1_u = 12; % maximum 1st delay
P_t1 = 1/(t1_u - t1_l + 1); % Prior PMF of tau_1

t2_l = 6; % minimum 2nd delay
t2_u = 12; % maximum 2nd delay
P_t2 = 1/(t2_u - t2_l + 1); % Prior PMF of tau_2

t3_l = 1; % minimum 3rd delay
t3_u = 3; % maximum 3rd delay
P_t3 = 1/(t3_u - t3_l + 1); % Prior PMF of tau_3

% Creating map objects for access of PDF/PMF values using the values of
% corresponding alpha instead of linear indexing (i.e. 1, 2, 3 ...)
f1 = containers.Map(x, pdf_1); % PDF of alpha_1
f2 = containers.Map(x, pdf_2); % PDF of alpha_2
f3 = containers.Map(x, pdf_3); % PDF of alpha_3
f4 = containers.Map(x, pdf_4); % PDF of alpha_4
f5 = containers.Map(x, pdf_5); % PDF of alpha_5
f6 = containers.Map(x, pdf_6); % PDF of alpha_6

p1 = containers.Map(t1_l:t1_u, repmat(P_t1, 1, t1_u-t1_l+1));%PMF of tau_1
p2 = containers.Map(t2_l:t2_u, repmat(P_t2, 1, t2_u-t2_l+1));%PMF of tau_2
p3 = containers.Map(t3_l:t3_u, repmat(P_t3, 1, t3_u-t3_l+1));%PMF of tau_3

max_tau = max([t1_u, t2_u, t3_u]);


%% Posterior Distribution of tau = [tau_1, tau_2, tau_3]

k = 1;
pdf_y = [];
for t1 = t1_l:t1_u
    for t2 = t2_l:t2_u
        for t3 = t3_l:t3_u 
            for a1 = x % x represents the query points
                for a2 = x
                    for a3 = x
                        for a4 = x 
                            for a5 = x 
                                for a6 = x 
                                    temp_y = [];
                                    for i = max_tau+1:length(y)%y represents data; max to avoid index error
                                        temp = (1/sqrt(2*pi*(sigma^2)))*...
                                            exp((-0.5*(y(i)-y(i-1)...
                                            -ah*q1(i-t1)-a2*q2(i-t2)...
                                -a3*q3(i-t3)-a4*q4(i)-a5*q5(i)-a6*q6(i)...
                                                     - mu)^2)/(sigma^2));
                                                 temp_y = [temp_y, temp];
                                    end
                                    pdf_y(k, 1) = t1;
                                    pdf_y(k, 2) = t2;
                                    pdf_y(k, 3) = t3;
                                    pdf_y(k, 4) = a1;
                                    pdf_y(k, 5) = a2;
                                    pdf_y(k, 6) = a3;
                                    pdf_y(k, 7) = a4;
                                    pdf_y(k, 8) = a5;
                                    pdf_y(k, 9) = a6;
                                    pdf_y(k, 10) = prod(temp_y)*...
                                        f1(a1)*f2(a2)*f3(a3)*f4(a4)*f5(a5)*...
                                        f6(a6)*p1(t1)*p2(t2)*p3(t3);
                                    k = k + 1;
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

This is where I am stuck because I haven't been able to get pdf_y. Because as expected, it takes a lot of time. I've waited for upto 4 hours but in vain. My idea is that after getting pdf_y, I can sum up the values corresponding to an alpha to simulate the integral (e.g. summing up values corresponding to a1 will mean integrating w.r.t ). And following the same procedure for the rest of alphas, I can obtain.

But the first step is to actually obtain the pdf_y.

Is there a faster way to do what I am trying to do with nine loops?

Answer 1

[SANKALP DEV]

Hey Muhammad,

I understand you want assistance in reducing the computational time of the code.

It appears that the primary issue in your code is the computational complexity due to the nested loops, which is causing significant runtime. To address this, I suggest some improvements by vectorizing operations and precomputing ‘PDF’ and ‘PMF’ values. This approach should help you calculate ‘pdf_y’ more efficiently.

Please refer to the following code.

% Precompute PDF values for alpha variables 

pdf_alpha = zeros(length(x), 6); 

pdf_alpha(:, 1) = getPDF(x, mu_1, sigma_1);  

pdf_alpha(:, 2) = getPDF(x, mu_2, sigma_2); 

pdf_alpha(:, 3) = getPDF(x, mu_3, sigma_3); 

pdf_alpha(:, 4) = getPDF(x, mu_4, sigma_4); 

pdf_alpha(:, 5) = getPDF(x, mu_5, sigma_5); 

pdf_alpha(:, 6) = getPDF(x, mu_6, sigma_6); 

 

% Precompute PMF values for tau variables 

P_tau = zeros(1, max_tau); 

P_tau(t1_1:t1_u) = 1 / (t1_u - t1_l + 1); 

P_tau(t2_1:t2_u) = 1 / (t2_u - t2_l + 1); 

P_tau(t3_1:t3_u) = 1 / (t3_u - t3_l + 1); 

 

% Vectorized computation of temp_y 

temp_y = (1/sqrt(2*pi*(sigma^2))) * exp((-0.5*(y(2:end)-y(1:end-1)... 

    - bsxfun(@times, [pdf_alpha(a1, 1), pdf_alpha(a2, 2), pdf_alpha(a3, 3), pdf_alpha(a4, 4), pdf_alpha(a5, 5), pdf_alpha(a6, 6)], [q1(t1), q2(t2), q3(t3), q4, q5, q6])... 

    - mu).^2) / (sigma^2)); 

 

% Calculate pdf_y using vectorized operations 

pdf_y = zeros(numel(x), numel(x), numel(x), numel(x), numel(x), numel(x), max_tau, max_tau, max_tau); 

for t1 = t1_l:t1_u 

    for t2 = t2_l:t2_u 

        for t3 = t3_l:t3_u 

            pdf_y(:, :, :, :, :, :, t1, t2, t3) = bsxfun(@times, bsxfun(@times, bsxfun(@times, bsxfun(@times, bsxfun(@times, bsxfun(@times, temp_y, pdf_alpha), P_tau(t1)), P_tau(t2)), P_tau(t3)), P_tau(t4)), P_tau(t5)), P_tau(t6)); 

        end 

    end 

end 

To know more about vectorization, please refer to following MATLAB documentation- https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html?s_tid=srchtitle_site_search_1_vectorization

Hope this helps.