Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
a random number obeying Gaussian + within a range

Subject: a random number obeying Gaussian + within a range

From: hailiang shen

Date: 18 Apr, 2009 23:18:02

Message: 1 of 16

Dear friends,
I am trying to generate 50 random numbers obeying Gaussian distribution (0,1), and within a range (-1, 1).

I tried to set the random number out of the range to 0 (i.e. the mean), but not sure whether the generated series also obeying Gaussian distribution. Is there theory behind this? How Matlab is dealing with this issue?

Thanks for any information,
Hailiang

Subject: a random number obeying Gaussian + within a range

From: John D'Errico

Date: 18 Apr, 2009 23:30:04

Message: 2 of 16

"hailiang shen" <hlshen2005@gmail.com> wrote in message <gsdn3a$g84$1@fred.mathworks.com>...
> Dear friends,
> I am trying to generate 50 random numbers obeying Gaussian distribution (0,1), and within a range (-1, 1).


Not possible. If it lives in a finite interval, then it is
not gaussian.

> I tried to set the random number out of the range to 0 (i.e. the mean), but not sure whether the generated series also obeying Gaussian distribution. Is there theory behind this?

Absolutely not, because your goal is not possible.
A normal (gaussian) distribution has infinite tails.

Certainly, setting anything outside of [-1,+1] to
zero will make the result strongly non-gaussian.

Do you wish to form a truncated Gaussian?


> How Matlab is dealing with this issue?

How does matlab deal with what issue? I don't
perceive any issue at all.

John

Subject: a random number obeying Gaussian + within a range

From: TideMan

Date: 18 Apr, 2009 23:35:29

Message: 3 of 16

On Apr 19, 11:18=A0am, "hailiang shen" <hlshen2...@gmail.com> wrote:
> Dear friends,
> I am trying to generate 50 random numbers obeying Gaussian distribution (=
0,1), and within a range (-1, 1).
>
> I tried to set the random number out of the range to 0 (i.e. the mean), b=
ut not sure whether the generated series also obeying Gaussian distribution=
. Is there theory behind this? How Matlab is dealing with this issue?
>
> Thanks for any information,
> Hailiang

I think you need to either go back to your textbook or Google
Gaussian.
If the numbers are restricted to being between -1 and 1, then they
will not be Gaussian distributed.
Indeed, the theory says that the expected maximum of a Gaussian
sequence is sqrt(2*log(n))*SD, so for 50 numbers from a (0,1) Gaussian
distribution you'd expect the range to be -2.8 to +2.8.
Perhaps you're looking for uniformly distributed random numbers?

Subject: a random number obeying Gaussian + within a range

From: hailiang shen

Date: 18 Apr, 2009 23:44:02

Message: 4 of 16

Thanks John. My actual problem is this: my physical parameter is meaningless for negative values, however, from literature, the distribution of this parameter is Gaussian distribution.
Does Matlab deals with ‘truncated Gaussian’?

Subject: a random number obeying Gaussian + within a range

From: hailiang shen

Date: 18 Apr, 2009 23:49:02

Message: 5 of 16

Thanks for all your replies.
For example, hydraulic conductivity can be described by normal or lognormal distribution, so if I get one negative number, I can not accept it.

My case is for nodal demand in water distribution system, which is describing how much water is taken in a specific location. It is normally distributed under uncertainty. Again I can not accept negative numbers.

Subject: a random number obeying Gaussian + within a range

From: Roger Stafford

Date: 19 Apr, 2009 01:03:03

Message: 6 of 16

"hailiang shen" <hlshen2005@gmail.com> wrote in message <gsdn3a$g84$1@fred.mathworks.com>...
> Dear friends,
> I am trying to generate 50 random numbers obeying Gaussian distribution (0,1), and within a range (-1, 1).
>
> I tried to set the random number out of the range to 0 (i.e. the mean), but not sure whether the generated series also obeying Gaussian distribution. Is there theory behind this? How Matlab is dealing with this issue?
>
> Thanks for any information,
> Hailiang

  As others in this thread have told you, a strictly Gaussian distribution, by definition, cannot be confined to a finite range. It must extend to plus and minus infinity. However you can generate a "truncated" Gaussian. Remember though, it is no longer Gaussian, only truncated Gaussian.

  If you have the Statistics Toolbox, do this:

 P1 = normcdf(-1,0,1);
 P2 = normcdf(+1,0,1);
 x = norminv(P1+(P2-P1)*rand(m,n));

The m by n matrix x should have values lying between -1 and +1 in accordance with the standard Gaussian distribution but truncated to lie between -1 and +1. If you don't have that toolbox, you can do the above in terms of the erfc and erfcinv functions which are contained in the basic matlab package.

  Beware; its standard deviation will no longer be 1 but necessarily something smaller than that.

Roger Stafford

Subject: a random number obeying Gaussian + within a range

From: Roger Stafford

Date: 19 Apr, 2009 01:30:04

Message: 7 of 16

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsdt87$ifj$1@fred.mathworks.com>...
> .....
> P1 = normcdf(-1,0,1);
> P2 = normcdf(+1,0,1);
> x = norminv(P1+(P2-P1)*rand(m,n));
> .....

  I left the necessary statistical parameters out of the third line. It should read:

 x = norminv(P1+(P2-P1)*rand(m,n),0,1);

Roger Stafford

Subject: a random number obeying Gaussian + within a range

From: Bruno Luong

Date: 19 Apr, 2009 09:32:01

Message: 8 of 16

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsdt87$ifj$1@fred.mathworks.com>...
> If you don't have that toolbox, you can do the above in terms of the erfc and erfcinv functions which are contained in the basic matlab package.
>

Like this. Beware that the standard deviation of truncated is no longer sigma, as Roger pointed out. I can be adjusted thought.

sigma=1;
range=3;

c=1/(sqrt(2)*sigma);
cdfinv = @(y) erfinv(y*2*erf(c*range) + erf(-c*range))/(c)
x=cdfinv(rand(1,10000))

std(x)
mean(x)
hist(x,64)

% Bruno

Subject: a random number obeying Gaussian + within a range

From: Matt

Date: 19 Apr, 2009 13:15:03

Message: 9 of 16

"hailiang shen" <hlshen2005@gmail.com> wrote in message <gsdok2$m0i$1@fred.mathworks.com>...
> Thanks John. My actual problem is this: my physical parameter is meaningless for negative values, however, from literature, the distribution of this parameter is Gaussian distribution.
------

This is a contradiction.

If the literature says the parameter is Gaussian and the literature is to be believed, then negative values for that parameter cannot be meaningless. By telling you that the parameter is Gaussian, the literature is implying that the parameter takes negative values with some probability (because all Gaussian variables do).

Subject: a random number obeying Gaussian + within a range

From: Bruno Luong

Date: 19 Apr, 2009 13:30:04

Message: 10 of 16

Here is a function generating truncated Gaussin with correct standard deviation:

% Save this in TruncatedGaussian.m file (included 2 nested function)

function x = TruncatedGaussian(sigma, range, n)
% x = TruncatedGaussian(sigma, range, n)
% Generate a random vector of length n; truncated Gaussian Distribution X
% in [-range,range]; mean(X) = 0 and std(X)=sigma

% Search for "sigmac" such that the truncated Gaussian having sigmac in
% the formula in the pdf procides a standard deviation equal to sigma
[sigmac res flag] = fzero(@scz,1.5*sigma,[],sigma,range);

if flag<0 % Someting is wrong
    fprintf('warning=range and sigma incompatible\n');
    sigmac = sigma;
else
    fprintf('sigmac=%g\n',sigmac)
    sigtarget=stdtrunc(range, sigmac);
    fprintf('sigma=%g\n',sigtarget)
end

% Inverse pg the cdf
c=1/(sqrt(2)*sigmac);
cdfinv = @(y) erfinv(y*2*erf(c*range) + erf(-c*range))/(c);
% Generate random
x=cdfinv(rand(1,n));

end

function stdt=stdtrunc(upper, sigma)
% Compute the standard deviation of a trunctated gaussian distribution
aa=upper./sigma;
it=erf(aa/sqrt(2));
corr=erf(aa/sqrt(2))-(2*aa)/sqrt(2*pi).*exp(-aa.^2/2);
stdt = sigma.*sqrt(corr./it);
end

function res=scz(sc, sigma, upper)
% Gateway for fzero
res = stdtrunc(upper, sc) - sigma;
end

% End of TruncatedGaussian.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Test this on command line:

sigma=2;
range=4;

x=TruncatedGaussian(sigma, range, 1e5);

mean(x)
std(x)
hist(x,64)

% Bruno

Subject: a random number obeying Gaussian + within a range

From: Bruno Luong

Date: 19 Apr, 2009 13:44:01

Message: 11 of 16

Obviously this distribution is "popular" enough to be documented in Wikipedia.
http://en.wikipedia.org/wiki/Truncated_normal_distribution

Bruno

Subject: a random number obeying Gaussian + within a range

From: Bruno Luong

Date: 19 Apr, 2009 15:11:01

Message: 12 of 16

I add few modifications to make the function more stable:

function x = TruncatedGaussian(sigma, range, varargin)
% x = TruncatedGaussian(sigma, range, n)
%
% Generate a random vector of size n; truncated Gaussian Distribution X
% in [-range,range]; mean(X) = 0 and std(X)=sigma

range=abs(range);
sigma=abs(sigma);

if isempty(varargin)
    n={};
else
    n=varargin;
end

if range^2<3*sigma^2
    warning('TruncatedGaussian: range and sigma incompatible\n');
    sigmac = Inf;
else
    % Search for "sigmac" such that the truncated Gaussian having sigmac in
    % the formula in the pdf procides a standard deviation equal to sigma
    [sigmac res flag] = fzero(@scz,1.5*sigma,[],sigma,range);
    sigmac=abs(sigmac);
    if flag<0 % Someting is wrong
        sigmac = Inf;
    end
end

% Debug
% fprintf('sigmac=%g\n',sigmac)
% sigtarget=stdtrunc(range, sigmac);
% fprintf('sigma=%g\n',sigtarget);

% Inverse of the cdf
if isinf(sigmac)
    % Uniform distribution to maximize the standard deviation
    cdfinv = @(y) range*(2*y-1);
else
    c=sqrt(2)*sigmac;
    cdfinv = @(y) c*erfinv(erf(range/c)*(2*y-1));
end

% Generate random
x = cdfinv(rand(n{:}));
% Prevent some nasty numerical issue with of erfinv function
x(x==inf)=range;
x(x==-inf)=-range;

end %

function c=corrtrunc(upper, s)
% Compute the standard deviation of a trunctated gaussian distribution
if isinf(s)
    c = upper^2/3;
else
    aa=(upper/sqrt(2))./s;
    corr = 1-(2/sqrt(pi))*aa.*exp(-aa.^2)./erf(aa);
    c = s.^2.*corr;
end
end

function stdt=stdtrunc(upper, s)
stdt = sqrt(corrtrunc(upper, s));
end

function res=scz(sc, sigma, upper)
% Gateway for fzero
res = corrtrunc(upper, sc) - sigma^2;
end

% End of TruncatedGaussian.m

%%%%

sigma=2
range=4

x=TruncatedGaussian(sigma, range, [1 1e5]);

disp(mean(x))
disp(std(x))
hist(x,64)

% Bruno

Subject: a random number obeying Gaussian + within a range

From: hailiang shen

Date: 20 Apr, 2009 18:50:17

Message: 13 of 16

Thanks for your entire replay, I get the new term 'truncated Gaussian distribution', and excellent code to implement them.

Subject: a random number obeying Gaussian + within a range

From: hailiang shen

Date: 21 Apr, 2009 18:55:02

Message: 14 of 16

I am confused by one term: a truncated normal distribution with a mean 0.075 and standard deviation 1. The (0.075, 1) may be of two cases:
i) of the normal distribution where truncated from,
ii) of the truncated normal distribution.

I think the 'norminv' function can support the case i). Is it necessary for case ii)? if so, how to generate a random number from case ii).

Thanks,
Hailiang

Subject: a random number obeying Gaussian + within a range

From: Bruno Luong

Date: 21 Apr, 2009 19:33:01

Message: 15 of 16

"hailiang shen" <hlshen2005@gmail.com> wrote in message <gsl4q6$g3n$1@fred.mathworks.com>...
> I am confused by one term: a truncated normal distribution with a mean 0.075 and standard deviation 1. The (0.075, 1) may be of two cases:
> i) of the normal distribution where truncated from,
> ii) of the truncated normal distribution.
>
> I think the 'norminv' function can support the case i). Is it necessary for case ii)? if so, how to generate a random number from case ii).
>

Necessary or not depends on what user needs. Nobody can answer at his/her place.

I have submitted this code on FEX with few additional improvements: select case (i) or (ii); possibility to adjust range with difference left and right bound, including infinity:

http://www.mathworks.com/matlabcentral/fileexchange/23832

Bruno

Subject: a random number obeying Gaussian + within a range

From: Zack

Date: 16 Aug, 2011 11:16:13

Message: 16 of 16

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gsfeu5$7oa$1@fred.mathworks.com>...
> I add few modifications to make the function more stable:
>
> function x = TruncatedGaussian(sigma, range, varargin)
> % x = TruncatedGaussian(sigma, range, n)
> %
> % Generate a random vector of size n; truncated Gaussian Distribution X
> % in [-range,range]; mean(X) = 0 and std(X)=sigma
>
> range=abs(range);
> sigma=abs(sigma);
>
> if isempty(varargin)
> n={};
> else
> n=varargin;
> end
>
> if range^2<3*sigma^2
> warning('TruncatedGaussian: range and sigma incompatible\n');
> sigmac = Inf;
> else
> % Search for "sigmac" such that the truncated Gaussian having sigmac in
> % the formula in the pdf procides a standard deviation equal to sigma
> [sigmac res flag] = fzero(@scz,1.5*sigma,[],sigma,range);
> sigmac=abs(sigmac);
> if flag<0 % Someting is wrong
> sigmac = Inf;
> end
> end
>
> % Debug
> % fprintf('sigmac=%g\n',sigmac)
> % sigtarget=stdtrunc(range, sigmac);
> % fprintf('sigma=%g\n',sigtarget);
>
> % Inverse of the cdf
> if isinf(sigmac)
> % Uniform distribution to maximize the standard deviation
> cdfinv = @(y) range*(2*y-1);
> else
> c=sqrt(2)*sigmac;
> cdfinv = @(y) c*erfinv(erf(range/c)*(2*y-1));
> end
>
> % Generate random
> x = cdfinv(rand(n{:}));
> % Prevent some nasty numerical issue with of erfinv function
> x(x==inf)=range;
> x(x==-inf)=-range;
>
> end %
>
> function c=corrtrunc(upper, s)
> % Compute the standard deviation of a trunctated gaussian distribution
> if isinf(s)
> c = upper^2/3;
> else
> aa=(upper/sqrt(2))./s;
> corr = 1-(2/sqrt(pi))*aa.*exp(-aa.^2)./erf(aa);
> c = s.^2.*corr;
> end
> end
>
> function stdt=stdtrunc(upper, s)
> stdt = sqrt(corrtrunc(upper, s));
> end
>
> function res=scz(sc, sigma, upper)
> % Gateway for fzero
> res = corrtrunc(upper, sc) - sigma^2;
> end
>
> % End of TruncatedGaussian.m
>
> %%%%
>
> sigma=2
> range=4
>
> x=TruncatedGaussian(sigma, range, [1 1e5]);
>
> disp(mean(x))
> disp(std(x))
> hist(x,64)
>
> % Bruno

Hi Bruno

Is there a way to have this fonction for a general truncated gaussian with mean mu, standard deviation sigma and a rang [a,b] perhapes a and b equal to inf

Thanks a lot for your help

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us