Compare commits

...

2 Commits

Author SHA1 Message Date
Renken baac9e2032 Added the first attempt at CCRE 2020-02-20 18:14:55 +01:00
Renken 2c7a9f90f0 Updated hist command to histcount in order to comply with matlab development 2020-02-19 09:59:14 +01:00
2 changed files with 26 additions and 12 deletions

View File

@ -76,11 +76,13 @@ classdef CREClass < handle
if isempty(obj.nBin) if isempty(obj.nBin)
obj.nBin=floor(numel(obj.Data)./10); obj.nBin=floor(numel(obj.Data)./10);
end end
[N,bin]=hist(obj.Data(:),obj.nBin); %update part below to prefered matlab commands
N=transpose(N); % [N,bin]=hist(obj.Data(:),obj.nBin);
P=N./sum(N); % N=transpose(N);
CP=cumsum(P); % P=N./sum(N);
FC=1-CP; % CP=cumsum(P);
[CP,bin]=histcounts(obj.Data(:),obj.nBin,'Normalization','cdf');
FC=1-transpose(CP);
FC(FC<0)=0; FC(FC<0)=0;
dl=ones(size(FC))./numel(FC); dl=ones(size(FC))./numel(FC);
else else
@ -90,7 +92,7 @@ classdef CREClass < handle
for k=1:numel(q) for k=1:numel(q)
FC(k)=sum(obj.Data(:)>=q(k)); FC(k)=sum(obj.Data(:)>=q(k));
if k==1 if k==1
dl(k)=(q(k+1)-q(k))/2; %should I divide by 2? dl(k)=(q(k+1)-q(k))/2;
elseif k==numel(q) elseif k==numel(q)
dl(k)=(q(k)-q(k-1))/2; dl(k)=(q(k)-q(k-1))/2;
else else
@ -101,8 +103,14 @@ classdef CREClass < handle
end end
LogFC=log(FC); LogFC=log(FC);
LogFC(~isfinite(LogFC))=0;%the log of 0 is -inf. however in Fc.*logFc it should end up as 0. to avoid conflicts removing the -inf LogFC(~isfinite(LogFC))=0;%the log of 0 is -inf. however in Fc.*logFc it should end up as 0. to avoid conflicts removing the -inf
if any(isnan(FC));keyboard;end if any(isnan(FC));error('something went wrong');end %catch a posible error.
obj.CRE=-transpose(dl)*(FC.*LogFC)./(transpose(dl)*FC); %CRE zografos if transpose(dl)*FC==0;
%if only one value in the most left bin in the distribution I may get a 0 divided by 0
%as the CRE of a delta function is 0, enforce this outcome
obj.CRE=0;
else
obj.CRE=-transpose(dl)*(FC.*LogFC)./(transpose(dl)*FC); %CRE zografos
end
%% get the RB %% get the RB
dl(FC>0.5)=0; %set the weight for all in the histogram to the left side (i.e. <p50) to 0 dl(FC>0.5)=0; %set the weight for all in the histogram to the left side (i.e. <p50) to 0
CRE_Med_Inf=-transpose(dl)*(FC.*LogFC)./(transpose(dl)*FC); %CRE zografos CRE_Med_Inf=-transpose(dl)*(FC.*LogFC)./(transpose(dl)*FC); %CRE zografos

View File

@ -1,5 +1,5 @@
S=CREClass; S=CREClass;
S.Data=randn(100000,1); S.Data=randn(100,1);
S.nBin=(1000); S.nBin=(1000);
%% use histogram %% use histogram
S.UseHistProxy=true; S.UseHistProxy=true;
@ -13,10 +13,16 @@ std(S.Data)
%% multiply S %% multiply S
S.Data=100*S.Data; S.Data=100*S.Data;
S.UseHistProxy=true; S.UseHistProxy=true;
S.Calc S.Calc;
S S
std(S.Data) std(S.Data)
%% ones again the "new" way %% ones again the "new" way
S.UseHistProxy=false; S.UseHistProxy=false;
S.Calc S.Calc;
S S
%% Test CCRE class
CS=CCREClass;
CS.Data=randn(100,1);
CS.nBin=10;
CS.DataRef=randn(100,1);
CS.nBinRef=10;