How to code for marginal and conditional distribution involving nine variables?
3 views (last 30 days)
Show older comments
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229102/image.png)
where
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229107/image.png)
Using Bayes' Theorem, the joint PDf of
can be written as:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229112/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229117/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229122/image.png)
where Mis the length of available data,
. (Sorry for seemingly too much maths)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229127/image.png)
This implies that ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229132/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229132/image.png)
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):
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229137/image.png)
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
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229142/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229147/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1229152/image.png)
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?
0 Comments
Answers (1)
SANKALP DEV
on 12 Sep 2023
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.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!