Performing matrix operation without loop

2 次查看（过去 30 天）
SChow
SChow 2020-8-17
评论： SChow 2020-8-17
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 个评论
显示 1更早的评论
Matt J
Matt J 2020-8-17
I find it peculiar that you are calculating ensemble estimates of mean and percentiles for data whose statistical distribution you already know? Is there an ulterior motive to this?
SChow
SChow 2020-8-17
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2020-8-17
编辑：Matt J 2020-8-17
Seems to me you could also just use the known quantile function for a log-normal distribution,
https://en.wikipedia.org/wiki/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 个评论
SChow
SChow 2020-8-17
Thanks Matt!!
That is a great solution!!
Many thanks for such perseverance

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2020-8-17
编辑：Matt J 2020-8-17
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 个评论
显示 2更早的评论
SChow
SChow 2020-8-17
Thanks Matt, I redid with your latest version.
Sorry but the resultant output is of size 4000x10x3. It should ideally be required to be 4000x8000. I tried without your reshape function, but Matlab shows out of memory error
Matt J
Matt J 2020-8-17
Make sure the for-loop is running over the full range
for n=1:size(muC,3)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by