Understanding the following line of code regarding Gaussian distribution.

Hello all, I am trying to plot PDF of Gaussian distribution and for that I came across the following code:

In this code, I understand every line except the one which is commented as "Not understood".

I also observed that if we omit this sentence then also code works properly.

Any help in this regard will be highly appreciated.

clc;
clear all;
close all;
% Number of samples
num_samples = 1000;
% Generate random samples from Gaussian distribution
mu_gaussian = 0; % Mean
sigma_gaussian = 1; % Standard deviation
samples_gaussian = mu_gaussian + sigma_gaussian * randn(num_samples, 1); % Not understood
x = linspace(-5, 5, 1000);
pdf_gaussian = (1 / (sigma_gaussian * sqrt(2*pi))) * exp(-(x - mu_gaussian).^2 / (2 * sigma_gaussian^2));
cdf_gaussian = 0.5 * (1 + erf((x - mu_gaussian) / (sigma_gaussian * sqrt(2))));

% Plotting PDFs
figure;
subplot(3, 1, 1);
plot(x, pdf_gaussian, 'b', 'LineWidth', 2);
title('Gaussian Distribution (PDF)');
xlabel('x');
ylabel('Probability Density');

[Voss]

See the definitions of standard normal and general normal distributions here:

https://en.wikipedia.org/wiki/Normal_distribution#Definitions

[Bruno Luong]

In Voss's wikipedia link, read especially the section

Operations on a single normal variable

first bullet aX + b ...

Answer 1

[Sam Chak]

Hi @charu shree

This variable samples_gaussian is unused in plotting the pdf and cdf. The code runs without error.

% Number of samples
num_samples      = 1000;

% Generate random samples from Gaussian distribution
mu_gaussian      = 0; % Mean
sigma_gaussian   = 1; % Standard deviation
samples_gaussian = mu_gaussian + sigma_gaussian*randn(num_samples, 1);    % Not understood

x = linspace(-5, 5, 1000);
pdf_gaussian = (1/(sigma_gaussian*sqrt(2*pi)))*exp(- (x - mu_gaussian).^2/(2*sigma_gaussian^2));
cdf_gaussian = 0.5*(1 + erf((x - mu_gaussian)/(sigma_gaussian*sqrt(2))));

% Plotting
figure;

subplot(2, 1, 1);
plot(x, pdf_gaussian, 'b', 'LineWidth', 2); grid on
xlabel('x');
ylabel('PDF');
title('Probability Density');

subplot(2, 1, 2);
plot(x, cdf_gaussian, 'b', 'LineWidth', 2); grid on
xlabel('x');
ylabel('CDF');
title('Cumulative Density')

[charu shree]

Thank u all for ur reply...

My doubt is why the code is written like: samples_gaussian = mu_gaussian + sigma_gaussian*randn(num_samples, 1); % Not understood

Cant we write like: samples_gaussian = sigma_gaussian + mu_gaussian*randn(num_samples, 1); % Not understood

[Bruno Luong]

"Cant we write like: samples_gaussian = sigma_gaussian + mu_gaussian*randn(num_samples, 1); % Not understood"

NO, yours is wrong. You must multiply by the desired standard deviation sigma, then you shift by the desire mean mu, as in the original code.

Often when you try to explain to other people precisely the reason of your doubt you'll figure out automatically the mistake.

[dpb]

@charu shree, look at <generating random normal with specified mean, variance>

Simply trying it both ways and plotting results will quickly demonstrate which way is correct, too, of course.

mn_g=5; std_g=2;
v=mn_g+std_g*randn(1000,1);
disp('Mean + Std*N(0,1)')

Mean + Std*N(0,1)

disp([mean(v) std(v)])

    4.9016    2.0720

v=std_g+mn_g*randn(1000,1);
disp('Std + Mean*N(0,1)')

Std + Mean*N(0,1)

disp([mean(v) std(v)])

    1.9326    5.1066

You see the first gives the expected mean, std, while the second has reversed the two values numerically, the mean is about what expect of the std while the std is where the mean should be.

[charu shree]

@Bruno Luong, @dpb Thank u so much sir for ur answer...

[Sam Chak]

Hi @charu shree

The "Not understood" line is actually applying a linear transformation to the random variable.

For more info:

https://stattrek.com/random-variable/transformation

https://bookdown.org/kevin_davisross/probsim-book/linear-rescaling.html

https://proofwiki.org/wiki/Expectation_of_Linear_Transformation_of_Random_Variable

You can visualize the generated random samples in this script:

% Number of samples
num_samples   = 1000;
x = linspace(-5, 5, num_samples);

% Normal Distribution with Specific Mean and Variance
desired_mu    = 0; % Mean
desired_sigma = 1; % Standard deviation
desired_pdf   = (1/(desired_sigma*sqrt(2*pi)))*exp(- (x - desired_mu).^2/(2*desired_sigma^2));
plot(x, desired_pdf, 'LineWidth', 2); hold on

% Generate random samples according to the Linear Transformation of Random Variable
for j = 1:5
    norm_samples    = desired_mu + desired_sigma*randn(num_samples, 1);
    actual_mu(j)    = mean(norm_samples);
    actual_sigma(j) = std(norm_samples);
    
    pdf_gaussian = (1/(actual_sigma(j)*sqrt(2*pi)))*exp(- (x - actual_mu(j)).^2/(2*actual_sigma(j)^2));
    plot(x, pdf_gaussian, '--')
    hold on
end
grid on
xlabel('x');
ylabel('PDF');
title('Probability Density');

[charu shree]

@Sam Chak Thank u sir for your prompt answer...

Answer 2

[dpb]

It does precisely what the comment ahead of it says it does -- generates a set of psuedo-random normally distributed values with mean mu_gaussian, std deviation sigma_gaussian. Since those are set to 0,1 in the given code, it's the same as what randn returns, but this lets one change the location and dispersion generally.

The code does the same thing whether it is/is not there because while the array was generated, it's never used for anything in the given code; the editor will point that out for you...unless it is used later, of course, but in that case then you would discover that if you were to comment it out and the code tried to reference it. So, one presumes it was "just one of those things" that the original coder had an idea of using/needing a sampled array and then changed mind/direction and never removed the superfluous definition. Or, maybe they never got done...we have no way to know.