% This script illustrates how to calculate noise corrected higher order moment from data.

% Load data to use as original distribution
load('Data_TH.mat')

%

% Get the example data and visualize distribution
Data_NoNoise = squeeze(nanmean(Data(10,:,:,15),2))';
figure
hist(Data_NoNoise)
title(['Kurtosis = ' num2str(kurtosis(Data_NoNoise))])

%

NoiseLev =  logspace(-5,-1,20) * std(Data_NoNoise);
clear Kurtosis
for nn = 1:length(NoiseLev)
    for iter = 1:100
        for tr = 1:15
            
            Data_Noise(:,tr) = Data_NoNoise + randn(1,140)*NoiseLev(nn);
            
        end
        [SigMom(nn,iter), NoiseMom,Kurtosis(nn,iter),Var(nn,iter)] = GetMomentsNoiseSignal(Data_Noise);
    end
end

%% Show for some example staitsitics that estimate tracks the true underlying value
figure
errorbar(NoiseLev,nanmean(Kurtosis'),std(Kurtosis'),'.-','MarkerSize',10)
hold on
line([min(NoiseLev) max(NoiseLev)],[1 1]*(kurtosis(Data_NoNoise)-3),'color','k')
set(gca,'XScale','log')
xlabel('Noise Level')
ylabel('Kurtosis')
legend('Estimate from noisy data','True value')

figure
TrueVariance = nanmean((Data_NoNoise - nanmean(Data_NoNoise)).^2);
errorbar(NoiseLev,nanmean(Var'),std(Var'),'.-','MarkerSize',10)
hold on
line([min(NoiseLev) max(NoiseLev)],[1 1]*(TrueVariance),'color','k')
set(gca,'XScale','log')
xlabel('Noise Level')
ylabel('Variance')
legend('Estimate from noisy data','True value')

figure
True3rdOrder= nanmean((Data_NoNoise).^3);
Est3rdOrder = reshape([SigMom.M3],nn,iter);
errorbar(NoiseLev,nanmean(Est3rdOrder'),std(Est3rdOrder'),'.-','MarkerSize',10)
hold on
line([min(NoiseLev) max(NoiseLev)],[1 1]*(True3rdOrder),'color','k')
set(gca,'XScale','log')
xlabel('Noise Level')
ylabel('3rd order - not centered')
legend('Estimate from noisy data','True value')