Azzera filtri
Azzera filtri

Performing matrix operation without loop

2 visualizzazioni (ultimi 30 giorni)
SChow
SChow il 17 Ago 2020
Commentato: SChow il 17 Ago 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 Commenti
Matt J
Matt J il 17 Ago 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 il 17 Ago 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

Accedi per commentare.

Risposta accettata

Matt J
Matt J il 17 Ago 2020
Modificato: Matt J il 17 Ago 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 Commento
SChow
SChow il 17 Ago 2020
Thanks Matt!!
That is a great solution!!
Many thanks for such perseverance

Accedi per commentare.

Più risposte (1)

Matt J
Matt J il 17 Ago 2020
Modificato: Matt J il 17 Ago 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 Commenti
SChow
SChow il 17 Ago 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 il 17 Ago 2020
Make sure the for-loop is running over the full range
for n=1:size(muC,3)

Accedi per commentare.

Community Treasure Hunt

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

Start Hunting!

Translated by