% Maximum number of samples nmax (recommended nmax=10, 100, ... 1e6)
function moments(nmax,prm)

% Parameter prm in (0:1)
%prm = 0.5;

% Create random variables
n = 1:nmax;
x = genrandom(nmax,prm);

% Sample sizes
% (magic code that creates up to 1000 points between 10 and nmax,
%  that determine the sample size to be evaluated)
sample_sizes = unique(ceil(logspace(0,log10(nmax),1000)))';
sample_sizes = sample_sizes((sample_sizes >= 10) & (sample_sizes <= nmax));

M = zeros(size(sample_sizes));        % mean <v>
V = zeros(size(sample_sizes));        % variance <v^2>
S = zeros(size(sample_sizes));        % skewness <v^3>/<v^2>^1.5
F = zeros(size(sample_sizes));        % flatness/kurtosis <v^4>/<v^2>^2

% Calculate moments of fixed sample size 
for i=1:size(sample_sizes,1)
    sample_size = sample_sizes(i);
    [M(i), V(i), S(i), F(i)] = calc_moments(x(1:sample_size));
end

% Compare to results from probability density function
% (this function is usually not known)
x2 = linspace(-100,100,1e7);
p2 = genpdf(x2,prm)*(x2(2)-x2(1));
M2 = sum(x2.*p2);
V2 = sum((x2-M2).^2.*p2);
S2 = sum((x2-M2).^3.*p2)/(V2^1.5);
F2 = sum((x2-M2).^4.*p2)/(V2^2);
ons = ones(size(sample_sizes));

% Figure 1: Plot random samples
figure(1)
clf;
plot(n,x)
xlabel('n');
ylabel('x');
suptitle('Random signal x');

% Figure 2: Histogram of the distribution of the x-value
figure(2)
clf;

[counts,centers]=hist(x,50);                    % make histogram
w=(centers(2)-centers(1));                      % bin width
norm_factor = sum(counts)*w;                    % factor for normalization

subplot(1,2,1)
bar(centers,counts/norm_factor,1.0,'b');        % blue: bar diagram of sampled data
hold on;
plot(...%centers,genpdf(centers,prm),'r',...    % red: analytic curve of sampled distribution
    centers,genpdf2(centers,M2,sqrt(V2)),'g')   % green: gaussian distribution
xlabel('x')
ylabel('p(x)')
xlim([0 15])
ylim([0 0.5])

subplot(1,2,2)
semilogy(centers,counts/norm_factor,'b',...     % blue: sampled data
    ...%centers,genpdf(centers,prm),'r',...     % red: analytic curve of sampled distribution
    centers,genpdf2(centers,M2,sqrt(V2)),'g');  % green: gaussian distribution
xlabel('x')
ylabel('p(x)')
xlim([0 15])
ylim([0 0.5])
suptitle('Probability distribution of x');


% Figure 3: Convergence of the moments to the final value
figure(3)
clf;
subplot(2,2,1)
semilogx(...%sample_sizes, M2*ons,'r',...       % red: analytic value
    sample_sizes,M,'b')                         % blue: sampled data
xlabel('sample size')
ylabel('Mean M')

subplot(2,2,2)
semilogx(...%sample_sizes, V2*ons,'r',...       % red: analytic value
    sample_sizes,V,'b')                         % blue: sampled data
xlabel('sample size')
ylabel('Variance V')

subplot(2,2,3)
semilogx(...%sample_sizes, S2*ons,'r',...       % red: analytic value
    sample_sizes,S,'b')                         % blue: sampled data
xlabel('sample size')
ylabel('Skewness S')

subplot(2,2,4)
semilogx(...%sample_sizes, F2*ons,'r',...       % red: analytic value
    sample_sizes,F,'b')                         % blue: sampled data
xlabel('sample size')
ylabel('Flatness F')

suptitle('Dependence on Sample Size');
end

% Calculate moments
% Input: Time series x of size (n x 1)
function [meanx, variance, skewness, flatness] = calc_moments(x)

    % Calculate mean
    meanx = mean(x);
    
    % Make x centered
    x = x - meanx;
    
    % Calculate 2nd, 3rd, 4th moment
    v2 = mean(x.^2);
    v3 = mean(x.^3);
    v4 = mean(x.^4);
    
    % Calculate variance, skewness, flatness
    variance = v2;
    skewness = v3 / (v2^1.5);
    flatness = v4 / (v2^2);
end

% Create (N x 1) random variables of a certain distribution
% (this is usually a measurement from a chaotic system)
function x = genrandom(N,prm)
    r = double(rand(N,1)<prm);
    x = (1-r).*genrandom2(5,1,N,1)+r.*genrandom2(6,2,N,1);
end

function r = genrandom2(mu, sigma, N, ~)
    r1 = rand(N,1);
    r2 = rand(N,1);
    r = sigma * sqrt(-2*log(r1)) .* cos(2*pi*r2) + mu;
end

% Probability density function for genrandom() variables
% (this function is usually not known)
function p = genpdf(x,prm)
    p = (1-prm)*genpdf2(x,5,1)+prm*genpdf2(x,6,2);
end

function p = genpdf2(x, mu, sigma)
    p = 1/sqrt(2*pi*sigma^2)*exp(-(x-mu).^2/(2*sigma^2));
end