From 4a4ddf0e9c2ca86152b7c877814d06c90845534f Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 7 Aug 2025 19:00:23 +0200 Subject: [PATCH] New script added to better understand cumulants. --- .../extractCustomCorrelation.m | 11 +- .../simulateDistribution.m | 39 ++- .../understandingCumulants.m | 223 ++++++++++++++++++ 3 files changed, 247 insertions(+), 26 deletions(-) create mode 100644 Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/understandingCumulants.m diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m index 61078de..dff7b7e 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m @@ -200,7 +200,7 @@ angle_threshold = 100; mean_max_g2_values = zeros(1, N_params); skew_max_g2_angle_values = zeros(1, N_params); var_max_g2_values = zeros(1, N_params); -kurt_max_g2_angle_values = zeros(1, N_params); +fourth_order_cumulant_max_g2_angle_values= zeros(1, N_params); fifth_order_cumulant_max_g2_angle_values = zeros(1, N_params); sixth_order_cumulant_max_g2_angle_values = zeros(1, N_params); @@ -265,10 +265,9 @@ for i = 1:N_params std_error_g2_values(i) = sqrt(kappa(2)) / sqrt(N_eff); skew_max_g2_angle_values(i) = kappa(3); - kurt_max_g2_angle_values(i) = kappa(4); + fourth_order_cumulant_max_g2_angle_values(i)= kappa(4); fifth_order_cumulant_max_g2_angle_values(i) = kappa(5); sixth_order_cumulant_max_g2_angle_values(i) = kappa(6); - end %% Plot PDF of order parameter @@ -384,10 +383,10 @@ ylabel('$\kappa_3$', 'Interpreter', 'latex', ... set(gca, 'FontSize', axis_fontsize); grid on; -% 4. Kurtosis +% 4. Binder Cumulant subplot(3,2,4); -plot(scan_vals, kurt_max_g2_angle_values, '^-', 'LineWidth', 1.5, 'MarkerSize', 6); -title('Kurtosis of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ... +plot(scan_vals, fourth_order_cumulant_max_g2_angle_values, '^-', 'LineWidth', 1.5, 'MarkerSize', 6); +title('Binder Cumulant of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ... 'Interpreter', 'latex', 'FontSize', title_fontsize); xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ... 'FontSize', label_fontsize); diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m index 68bd956..87d40a6 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m @@ -7,16 +7,16 @@ alpha_values = linspace(0, 45, N_params); all_data = cell(1, N_params); % Transition control -alpha_start = 5; % where transition begins -alpha_end = 40; % where transition ends +alpha_start = 5; % where transition begins +alpha_end = 40; % where transition ends -mu_start = 1.2; % high initial mean -mu_end = 0.8; % low final mean +mu_start = 1.2; % high initial mean +mu_end = 0.8; % low final mean -sigma_start = 0.2; % wide initial std -sigma_end = 0.07; % narrow final std +sigma_start = 0.2; % wide initial std +sigma_end = 0.07; % narrow final std -max_skew = 5; % peak skew strength +max_skew = 5; % peak skew strength % Loop through alpha for i = 1:N_params @@ -77,8 +77,7 @@ for i = 1:N_params elseif alpha <= (bimodal_start + transition_width) % Smooth transition window t = (alpha - (bimodal_start - transition_width)) / (2 * transition_width); - - + % Weights transition 1 -> 0.7 w1 = 1 - 0.3 * t; w2 = 1 - w1; @@ -104,8 +103,8 @@ for i = 1:N_params else % After transition: bimodal, but not strongly balanced - w1 = 0.7; - w2 = 0.3; + w1 = 0.5; + w2 = 0.5; mu1 = 1 - 0.3; mu2 = 1 + 0.3; @@ -123,14 +122,14 @@ for i = 1:N_params data = data(randperm(length(data))); end - all_data{i} = data; - kappa = computeCumulants(data,6); - mean_vals(i) = kappa(1); - var_vals(i) = kappa(2); - skew_vals(i) = kappa(3); - kurt_vals(i) = kappa(4); - kappa5_vals(i) = kappa(5); - kappa6_vals(i) = kappa(6); + all_data{i} = data; + kappa = computeCumulants(data,6); + mean_vals(i) = kappa(1); + var_vals(i) = kappa(2); + skew_vals(i) = kappa(3); + kurt_vals(i) = kappa(4); + kappa5_vals(i) = kappa(5); + kappa6_vals(i) = kappa(6); end %% Animate evolving distribution and cumulant value @@ -238,7 +237,7 @@ grid on; % 4. Kurtosis subplot(3,2,4); plot(scan_vals, kurt_vals, '^-', 'LineWidth', 1.5, 'MarkerSize', 6); -title('Kurtosis of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex'); +title('Fourth-order Cumulant of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex'); xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize); ylabel('$\kappa_4$', 'Interpreter', 'latex', 'FontSize', label_fontsize); set(gca, 'FontSize', axis_fontsize); diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/understandingCumulants.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/understandingCumulants.m new file mode 100644 index 0000000..80fc723 --- /dev/null +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/understandingCumulants.m @@ -0,0 +1,223 @@ +%% Main script: Sweep different parameter pairs +% Default parameters +defaults.mu1 = 0.5; +defaults.mu2 = 1.0; +defaults.sigma1 = 0.1; +defaults.sigma2 = 0.1; +defaults.weight1 = 0.5; + +% Parameter pair definitions +%{ +param_pairs = { + 'mu1', linspace(0.7, 1.0, 40), ... + 'mu2', linspace(1.0, 1.3, 40); + + 'mu1', linspace(0.7, 1.0, 40), ... + 'weight1', linspace(0.2, 0.8, 40); + + 'sigma1', linspace(0.05, 0.2, 40), ... + 'sigma2', linspace(0.05, 0.2, 40); + + 'mu1', linspace(0.7, 1.0, 40), ... + 'sigma1', linspace(0.05, 0.2, 40); + + 'mu2', linspace(1.0, 1.3, 40), ... + 'weight1', linspace(0.2, 0.8, 40); +}; +%} +param_pairs = { + 'mu1', linspace(0.1, 1.5, 40), ... + 'mu2', linspace(0.1, 1.5, 40); +}; + +% Cumulant index to visualize (2=variance, 3=skewness, 4=kurtosis) +cumulant_to_plot = 4; + +% Run sweep for each pair +for i = 1:size(param_pairs,1) + param1_name = param_pairs{i,1}; + param1_vals = param_pairs{i,2}; + param2_name = param_pairs{i,3}; + param2_vals = param_pairs{i,4}; + + fprintf('Sweeping %s and %s...\n', param1_name, param2_name); + + Z = sweepBimodalCumulants(param1_name, param1_vals, ... + param2_name, param2_vals, ... + defaults, cumulant_to_plot); +end + +%% +% Parameters +N_total = 10000; +mu1 = 0.5; +sigma1 = 0.1; +mu2 = 1.0; +sigma2 = 0.1; +weight1 = 0.7; + +% Generate data +data = generateBimodalDistribution(N_total, mu1, mu2, sigma1, sigma2, weight1); + +% Plot histogram +figure(3); +clf +set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +histogram(data, 'Normalization', 'pdf', 'EdgeColor', 'none', 'FaceAlpha', 0.5); +hold on; + +% Overlay smooth density estimate +[xi, f] = ksdensity(data); +plot(f, xi, 'r-', 'LineWidth', 2); + +% Labels and title +xlabel('Value'); +ylabel('Probability Density'); +legend('Histogram', 'Smoothed Density'); +grid on; + +%% +N = size(Z, 1); + +main_diag_values = diag(Z); +anti_diag_values = diag(flipud(Z)); + +param1_diag_main = param1_vals; +param2_diag_main = param2_vals; + +param1_diag_anti = param1_vals; +param2_diag_anti = flip(param2_vals); + +% For example, plot the main diagonal cumulants: +figure(4); +set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +plot(1:N, main_diag_values, '-o'); +xlabel('Index along diagonal'); +ylabel('$\kappa_4$', 'Interpreter', 'latex'); +title('$\kappa_4$ along anti-diagonal', 'Interpreter', 'latex'); + +% Plot anti-diagonal cumulants: +figure(5); +set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +plot(1:N, anti_diag_values, '-o'); +xlabel('Index along anti-diagonal'); +ylabel('$\kappa_4$', 'Interpreter', 'latex'); +title('$\kappa_4$ along anti-diagonal', 'Interpreter', 'latex'); + +%% === Helper: Bimodal Distribution === +function data = generateBimodalDistribution(N_total, mu1, mu2, sigma1, sigma2, weight1) +%GENERATEBIMODALDISTRIBUTION Generates a single bimodal distribution. +% +% data = generateBimodalDistribution(N_total, mu1, mu2, sigma1, sigma2, weight1) +% +% Inputs: +% N_total - total number of samples +% mu1, mu2 - means of the two modes +% sigma1, sigma2 - standard deviations of the two modes +% weight1 - fraction of samples from mode 1 (between 0 and 1) +% +% Output: +% data - shuffled samples from the bimodal distribution + + % Validate weight + weight1 = min(max(weight1, 0), 1); + weight2 = 1 - weight1; + + % Determine number of samples for each mode + N1 = round(N_total * weight1); + N2 = N_total - N1; + + % Generate samples + mode1_samples = normrnd(mu1, sigma1, [N1, 1]); + mode2_samples = normrnd(mu2, sigma2, [N2, 1]); + + % Combine and shuffle + data = [mode1_samples; mode2_samples]; + data = data(randperm(length(data))); +end + +%% === Helper: Cumulant Calculation === +function kappa = computeCumulants(data, max_order) + data = data(:); + mu = mean(data); + centered = data - mu; + + % Preallocate + c = zeros(1, max_order); + kappa = zeros(1, max_order); + + % Compute central moments up to max_order + for n = 1:max_order + c(n) = mean(centered.^n); + end + + % Assign cumulants based on available order + if max_order >= 1, kappa(1) = mu; end + if max_order >= 2, kappa(2) = c(2); end + if max_order >= 3, kappa(3) = c(3); end + if max_order >= 4, kappa(4) = c(4) - 3*c(2)^2; end + if max_order >= 5, kappa(5) = c(5) - 10*c(3)*c(2); end + if max_order >= 6 + kappa(6) = c(6) - 15*c(4)*c(2) - 10*c(3)^2 + 30*c(2)^3; + end +end + +%% === Helper: Cumulant Calculation === +function Z = sweepBimodalCumulants(param1_name, param1_vals, ... + param2_name, param2_vals, ... + fixed_params, ... + cumulant_index) +%SWEEPBIMODALCUMULANTS Sweep 2 parameters and return a chosen cumulant. +% +% Z = sweepBimodalCumulants(...) +% Returns a matrix Z of cumulant values at each grid point. + + % Setup grid + [P1, P2] = meshgrid(param1_vals, param2_vals); + Z = zeros(size(P1)); + + N_samples = 1000; + maxOrder = max(4, cumulant_index); + + for i = 1:numel(P1) + % Copy fixed parameters + params = fixed_params; + + % Override swept parameters + params.(param1_name) = P1(i); + params.(param2_name) = P2(i); + + % Generate and compute cumulants + data = generateBimodalDistribution(N_samples, ... + params.mu1, params.mu2, ... + params.sigma1, params.sigma2, ... + params.weight1); + + kappa = computeCumulants(data, maxOrder); + Z(i) = kappa(cumulant_index); + end + + % Plot full heatmap + figure; + set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); + imagesc(param1_vals, param2_vals, Z); + set(gca, 'YDir', 'normal'); + xlabel(param1_name); + ylabel(param2_name); + title(['Cumulant \kappa_', num2str(cumulant_index)]); + colorbar; + axis tight; + + % Optional binary colormap (red = ≥0, blue = <0) + figure; + set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); + imagesc(param1_vals, param2_vals, Z); + set(gca, 'YDir', 'normal'); + xlabel(param1_name); + ylabel(param2_name); + title(['Binary color split of \kappa_', num2str(cumulant_index)]); + clim([-1 1]); + colormap([0 0 1; 1 0 0]); % Blue (neg), Red (pos & zero) + colorbar; + axis tight; +end