How to separate into unimodal or bimodal

Hello, I have some histograms. How can I separate them through Matlab code and tell which one is unimodal or bimodal ?

thanks in advance

Answer 1

[Image Analyst]

Maybe use kstest(), lillietest(), or chi2gof().

Or fit the data to two Gaussians and see if the peak locations are separated by enough distance. See attached demos.

Also see links on the right, such as

https://www.mathworks.com/matlabcentral/answers/474774-kstest-for-unimodal-and-bimodal-models?s_tid=answers_rc1-2_p2_MLT

[Jórdan Venâncio Leite]

Thanks for the answer Image Analyst!

[Jórdan Venâncio Leite]

I'm trying to use the kstest, but i'm not understanding how I'm going to use a standard normal distribution to achieve my goal of differentiating unimodal from bimodal pattern

[Image Analyst]

Did you try fitting the data to two Gaussians using the other demo?

[Jórdan Venâncio Leite]

Yes, however, I still don't know what this part of the code here below does. I know before this is the part where I enter my noisy data. However, after that I got a little confused.

% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X(:), Y(:));
% Define the model as Y = a + b*x + c*exp(-(x-d)^2/e) + d * exp(-(x-f)^2/g)
% Note how this "x" of modelfun is related to big X and big Y.
% x(:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) + b(2) * x(:, 1) + b(3) * exp(-(x(:, 1) - b(4)).^2/b(5)) + b(6) * exp(-(x(:, 1) - b(7)).^2/b(8));  
beta0 = [6, 0.1, 35, 10, 3, 25, 50, 9]; % Guess values to start with.  Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!

% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}

% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,
% so that we'll get a much smoother curve.
X = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.

% Create smoothed/regressed data using the model:
yFitted = coefficients(1) + coefficients(2) * X + coefficients(3) * exp(-(X - coefficients(4)).^2 / coefficients(5)) + ...
	coefficients(6) * exp(-(X - coefficients(7)).^2 / coefficients(8));
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'northeast');
legendHandle.FontSize = 25;

[Jórdan Venâncio Leite]

My code is like this

load('image66.mat');
img=PerfilCinzasVertical;
N = histogram(img,12);

using the archive attached

[Jórdan Venâncio Leite]

Image Analyst, i can keep using this model here in my case which is a histogram? MODEL: Y = a + b*x + c*exp(-(x-d)^2/e) + d * exp(-(x-f)^2/g). How do I convert the histogram into a table to be used in the line 10 of my code knowing that the variable of this histogram in the workspace is a graph and not values?

clc
clear
close all

load('image66.mat');
img=PerfilCinzasVertical;
N = histogram(img,12);

% Convert the histogram into a table
tbl = table(X(:), Y(:));

% Define the model as Y = a + b*x + c*exp(-(x-d)^2/e) + d * exp(-(x-f)^2/g)
% Note how this "x" of modelfun is related to big X and big Y.
% x(:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.

modelfun = @(b,x) b(1) + b(2) * x(:, 1) + b(3) * exp(-(x(:, 1) - b(4)).^2/b(5)) + b(6) * exp(-(x(:, 1) - b(7)).^2/b(8));  
beta0 = [6, 0.1, 35, 10, 3, 25, 50, 9]; % Guess values to start with.  Just make your best guess.

% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!

% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}

% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,
% so that we'll get a much smoother curve.
X = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.

% Create smoothed/regressed data using the model:
yFitted = coefficients(1) + coefficients(2) * X + coefficients(3) * exp(-(X - coefficients(4)).^2 / coefficients(5)) + ...
	coefficients(6) * exp(-(X - coefficients(7)).^2 / coefficients(8));

% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'northeast');
legendHandle.FontSize = 25;

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

[Image Analyst]

X would be the bin values, and Y would be the bin counts.

[Jórdan Venâncio Leite]

Image Analyst, which model you suggest me to use with the figure 2 (second histogram)?

Thanks in advance

[Image Analyst]

Looks more of a log normal or exponential decay for the left hump and then possibly a gaussian for the right hump.

[Jórdan Venâncio Leite]

Image Analyst. I'm having doubts about which model to use to represent the gaussian. Can you give me some tips ?

clc;
clear
close all

load('image66.mat');
img=PerfilCinzasVertical;
N = histogram(img,12);

% Get the bin values into their own separate variable.
x = N.Values

% Get the bin counts into their own separate variable.
y = N.BinEdges(1:12)

% Convert the histogram into a table
tbl = table(x(:), y(:));

% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!

% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}

% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,
% so that we'll get a much smoother curve.
X = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.

% Create smoothed/regressed data using the model:
yFitted = coefficients(1) + coefficients(2) * X + coefficients(3) * exp(-(X - coefficients(4)).^2 / coefficients(5)) + ...
	coefficients(6) * exp(-(X - coefficients(7)).^2 / coefficients(8));

% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'northeast');
legendHandle.FontSize = 25;

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

[Image Analyst]

Try this:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
s = load('image66.mat')
subplot(2, 1, 1);
img = s.PerfilCinzasVertical;
plot(img, 'b-');
grid on;
title('Original Data', 'FontSize', fontSize);
subplot(2, 1, 2);
histogramObject = histogram(img, 12)
xlabel('Value', 'FontSize', fontSize);
ylabel('Count', 'FontSize', fontSize);
grid on;
% Convert the histogram into a table
X = histogramObject.BinEdges(1:end-1);
Y = histogramObject.BinCounts;
tbl = table(X(:), Y(:))
% Define the model as Y = a + b*x + c*exp(-(x-d)^2/e) + d * exp(-(x-f)^2/g)
% Note how this "x" of modelfun is related to big X and big Y.
% x(:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) * exp(-(x(:, 1) - b(2)).^2/b(3)) + b(4) * exp(-(x(:, 1) - b(5)).^2/b(6));  
beta0 = [100, 190, 200, 150, 235, 200]; % Guess values to start with.  Just make your best guess.
% Plot guess
xl = xlim; % Get range of x axis.
xGuess = linspace(xl(1), xl(2), 500)';
yGuess = modelfun(beta0, xGuess)
hold on;
plot(xGuess, yGuess, 'r-');
hold off;
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,
% so that we'll get a much smoother curve.
X = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) * exp(-(X - coefficients(2)).^2 / coefficients(3)) + ...
    coefficients(4) * exp(-(X - coefficients(5)).^2 / coefficients(6));
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Guess', 'Fitted Y', 'Location', 'northwest');
legendHandle.FontSize = 25;
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

[Jórdan Venâncio Leite]

Thank You Very Much for your answer ImageAnalyst !!! Do you think I can have automatic processing of my histograms? That is, for each different histogram, do you believe I can automatically determine if it is unimodal, bimodal or multimodal? or do you believe that i would need to manually configure a model every time i have a different histogram? thanks in advance!

[Image Analyst]

Yes, I'd probably first fit a single mode gaussian, and then a two Gaussian model. Then whichever one has the lowest MSE I'd say is the one it is.

[Image Analyst]

Mean Square Error. You can use immse() or just use the formula

mse = mean((trueValue - estimatedValue) .^ 2)

[Jórdan Venâncio Leite]

Image Analyst, are you saying that it is not possible to make this distinction in the histograms with respect to whether it is unimodal, bimodal or multimodal automatically? txs in advance.

[Image Analyst]

To be clear, I am not a statistician. I don't know that much about kstest, lillietest, anova, and other stats things. I suggest you ask your local statistician.

Why do you need to know anyway? What would you do with your signal if you knew the histogram indicated your signal came from two populations or one population?

[Jórdan Venâncio Leite]

I just need to know if the histogram is unimodal, bimodal or multimodal. This is my only problem for now.

[Image Analyst]

Strange. So you don't know why you need to know? Is it your homework then, and it's just something the professor asked? Did you discuss some method in class? Or if it's not your homework, why don't we just forget about that and figure out what to do with your original signal. Of course any signal that has even a little bit of noise could be considered the composition of two Gaussians. So let's just assume it is. OK, now what? What do you now need to do with your original signal knowing that its histogram is two Gaussians?

Well sorry I couldn't help more (I'm not an expert statistician) but good luck.

[Jórdan Venâncio Leite]

THANKS for your answer IMAGE ANALYST!! See if you can help me with this https://www.mathworks.com/matlabcentral/answers/1693350-how-to-define-if-a-histogram-is-unimodal-or-bimodal. It's almost finished already .