Filter löschen
Filter löschen

Performing matrix operation without loop

3 Ansichten (letzte 30 Tage)
SChow
SChow am 17 Aug. 2020
Kommentiert: SChow am 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 Kommentare
Matt J
Matt J am 17 Aug. 2020
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 am 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

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Matt J
Matt J am 17 Aug. 2020
Bearbeitet: Matt J am 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 Kommentar
SChow
SChow am 17 Aug. 2020
Thanks Matt!!
That is a great solution!!
Many thanks for such perseverance

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Matt J
Matt J am 17 Aug. 2020
Bearbeitet: Matt J am 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 Kommentare
SChow
SChow am 17 Aug. 2020
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 am 17 Aug. 2020
Make sure the for-loop is running over the full range
for n=1:size(muC,3)

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by