Performing matrix operation without loop

19 views (last 30 days)
SChow
SChow on 17 Aug 2020
Commented: SChow on 17 Aug 2020
Hi,
I have two matrices (A and B), each of size, 4000x8000x3 (i,j,k).
The 'k th' dimension of both matrix contains median and 95% CI values.
To be robust by combining the 95%CI of both matrices and calculate the product of these matrices, I distribute each ' ith & jth' value of a matrix across the 'k th' values (95%CI) 1000 times for A and B. and obtain a 1000x1000 matrix for each ith and jth element. From which I calculate the mean and 95% CI for each ith and jth element.
However the below code takes > 4.5 hours to complete, I would request an easier and quicker way to do it. Thank you for the help
tic
for i=1:4000
for j=1:8000;
%%%distribute A lognormally
A_err= (log(A(i,j,3)) - log(A(i,j,1)))/1.96;
A_dist= lognrnd(log(A(i,j,1)), A_err,[1,1000]);
A_dist_m=repmat(A_dist,1000,1);
%%%%distribute B
B_err= (log(B(i,j,3)) - log(B(i,j,1)))/1.96;
B_dist= lognrnd(log(B(i,j,1)), B_err,[1,1000]);
% histogram (B_dist)
B_dist=B_dist';
B_dist_m=repmat(B_dist,1,1000);
C=B_dist_m.*A_dist_m;
range = prctile(C(:),[2.5 97.5]);
middle = mean(mean(C,1),2);
pr=[middle,range];
C1 (i,j,:)=C;
end;
end
toc
  3 Comments
SChow
SChow on 17 Aug 2020
Hi Matt, I am trying to combine the uncertainties (95%CI) of A and B> Thought this would be a good way to progress, any suggestions are welcome

Sign in to comment.

Accepted Answer

Matt J
Matt J on 17 Aug 2020
Edited: Matt J on 17 Aug 2020
Seems to me you could also just use the known quantile function for a log-normal distribution,
That wouldn't require any loops at all and you could compute any confidence interval that you want...
muA=log(A(:,:,1));
muB=log(B(:,:,1));
muC=muA+muB;
sigA=(log(A(:,:,3))-log(A(:,1)) )/1.96;
sigB=( log(B(:,:,3))-log(B(:,1)) )/1.96;
sigC=sigA+sigB;
Quant=@(p) exp( muC + sqrt(2).*erfinv(2*p-1).*sigC ); %quantile function
  1 Comment
SChow
SChow on 17 Aug 2020
Thanks Matt!!
That is a great solution!!
Many thanks for such perseverance

Sign in to comment.

More Answers (1)

Matt J
Matt J on 17 Aug 2020
Edited: Matt J on 17 Aug 2020
I think this should be faster. Obviously, it will be even faster if you can accept fewer than nSamples=1000 randomized samples for the calculations. Maybe 100 would be enough?
nSamples=1e3;
[M,N,~]=size(A);
e=ones(1,1,nSamples,'single');
res=@(z) reshape(single(z),M,10,[]);
muA=log(A(:,:,1));
muB=log(B(:,:,1));
muC=res(muA+muB);
sigA=(log(A(:,:,3))-log(A(:,1)) )/1.96;
sigB=( log(B(:,:,3))-log(B(:,1)) )/1.96;
sigC=res(sigA+sigB);
[middle,rangeLOW,rangeHIGH]=deal(nan(M,10));
tic
for n=1:size(muC,3)
C_dist=lognrnd(muC(:,:,n).*e, sigC(:,:,n).*e);
middle(:,:,n)=mean(C_dist,3);
range = prctile(C_dist,[2.5 97.5],3);
rangeLOW(:,:,n)=range(:,:,1);
rangeHIGH(:,:,n)=range(:,:,2);
end
toc
C1=cat(3,middle,rangeLOW,rangeHIGH);
C1=reshape(C1,M,[],3);
  4 Comments
Matt J
Matt J on 17 Aug 2020
Make sure the for-loop is running over the full range
for n=1:size(muC,3)

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by