Random perturbation of a value in MATLAB

67 Views Asked by At

Let's say I have a value 100 and I want to perturb this value randomly by 1% based on a normal distribution. How could I do this in MATLAB? I seem to have some confusion with randn etc.

I have used randn so far but I do not seem to be inside this 1% threshold all the time.

1

There are 1 best solutions below

0
Wolfie On

The output of the normal distribution randn is unbounded, it does not give values in the range [0,1] like rand does.

This paper offers a good explanation of the "tool" you'll need, which is the truncated normal distribution: https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf

Informally, the truncated normal probability density function is defined in two steps. We choose a general normal PDF by specifying parameters µ and σ, and then a truncation range (a, b). The PDF associated with the general normal distribution is modified by setting values outside the range to zero, and uniformly scaling the values inside the range so that the total integral is 1.

This explicitly sets the function output to 0 outside your bounded range, and scales it by the difference of the cumulative distribution evaluated at b and a. Note that the standard deviation of this new distribution will not be the same as the original standard deviation, although if we choose symmetrical a and b around the mean then that should remain unchanged.

I'm going to skip the maths and implement one suggested sampling method from the linked literature which simply "imposes the truncation interval by rejection" with the following condition:

repeat:
    p = rand()
    x = Ф⁻¹(µ, σ; p)
until (a ≤ x ≤ b)

Where Ф⁻¹(µ, σ; p) is the inverse CDF of your normal distribution with mean µ and standard deviation σ, and rand() generates a uniform random number. MATLAB has an in-built function for sampling the inverse CDF if you have the Statistics toolbox: icdf('normal',p,mu,sigma)

So we can write a function like

function r = normrndt( mu, sigma, a, b, varargin )
% Output values in the truncated normal distribution, equivalent syntax to
% normrnd( mu, sigma, sz1, sz2, ... szN )
% But constrained to interval [a,b]

    % Start with uniform random values
    r = rand(varargin{:});
    % Feed them into the inverse CDF
    r = icdf('normal',r,mu,sigma);

    % Impose truncation interval by rejection
    % Repeating the above step per value if outside range
    for ii = 1:numel(r)
        while r(ii) < a || r(ii) > b
            r(ii) = icdf('normal',rand(),mu,sigma);
        end
    end
end

Testing this for different values of sigma, and a symmetrical distribution around mean 0 and interval [-1,1]. For your example, you could add these values to a base value of 100:

sigma = 0.4; % standard deviation of normal dist
N = 1e5;     % number of points to sample
n = normrnd(0,sigma,N,1);        % Generate "proper" normally distributed values
nt = normrndt(0,sigma,-1,1,N,1); % Generate truncated normally distributed values

% Plot the results for comparison
edges = linspace(-sigma*5, sigma*5, 100); % bin edges for histo
figure; histogram( n, edges ); hold on; histogram( nt, edges ) % plot overlay of outputs

You can see that for sigma values where the truncation interval lies several standard deviations away from the mean, there is hardly any change to the distribution. However as the standard deviation gets closer to the boundary we see the distribution become less and less "normal" looking.

comparison for various sigma