Performing matrix operation without loop
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
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?
采纳的回答
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,
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 个)
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);
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!