Perfect random numbers - How can I refine an initial dataset from randn to be perfrectly normal?

Hi, I would like to create "perfect" random numbers. By perfect I mean a set of random values that has

a mean of 0.0000, a std of 1.0000, a skewness of 0.0000, and a kurtosis of 3.0000

while each of numbers created with randn comes very close, it is not perfect in that sense. How can I further refine the initial dataset from randn to be "perfectly" normal?

Thanks

[Matt Kindig]

I'm not you can do this with a finite number of random values.

Answer 1

[Teja Muppirala]

If you want to generate a set of numbers that "looks random", has 0 mean, 1 variance, 0 skewness, 3 kurtosis, you can start with an initial dataset using RANDN and then kind of iteratively perturb it until you get what you want.

rng(0)
R = randn(1000,1);

R0 = R;
R_prev= 0;
opts = optimset('TolFun',0,'Display','Off');
smallNumber = 1e-12;

while norm(R-R_prev)/norm(R_prev) > smallNumber;
    R_prev = R;

      R = R-mean(R); %Fix Mean
      R = R/std(R); %Fix Std. Dev

      % Fix Skewness
      x = fzero( @(x) skewness( R .* exp(x*R) ) ,0 ,opts);
      R = R .* exp(x*R);

      % Fix Kurtosis
      x = fzero( @(x) kurtosis(sign(R) .* (abs(R).^x))-3 ,1 ,opts);
      R = sign(R) .* (abs(R).^x);

end

figure, plot(R0,'b');
hold on;
plot(R,'r');
legend({'Original', 'Slightly Modified'});

format long
disp('Original:')
[mean(R0) std(R0) skewness(R0) kurtosis(R0)]

disp('Modified:')
[mean(R) std(R) skewness(R) kurtosis(R)]

[Sven]

Teja, thanks - tested your code and it does what I was trying to achieve. As your code is using a couple of functions that I have not used before, it will take a bit to figure out what exactly is going on. That said, could you please point me in the right direction as to how I would be able to redefine my target distribution. If I wanted to have a target settings for std, skew and kurtosis, where in the code would I put them in? So for instance if I want random numbers that have exactly a std of 1.2, a skew of .5 and a kurtosis .4 - how can set those in your code? Thanks for guiding me in the right direction.

Everyone else: thanks for your comments. I guess I should have been a bit more precise in my question regarding being perfect: I was referring to almost perfect only...errors in the 3rd or 4th decimal point as by using this code is absolutely fine for what I am after.

[Teja Muppirala]

To generate a set of points with standard deviation of 1.2, a skew of 0.5, and an (excess) kurtosis of 0.4

You would change these three lines:

R = R/std(R); %Fix Std. Dev
x = fzero( @(x) skewness( R .* exp(x*R) ) ,0 ,opts);
x = fzero( @(x) kurtosis(sign(R) .* (abs(R).^x))-3 ,1 ,opts);

to

R = R/std(R)*1.2; %Fix Std. Dev      % Fix Skewness
x = fzero( @(x) skewness( R .* exp(x*R) )-0.5 ,0 ,opts);
x = fzero( @(x) kurtosis(sign(R) .* (abs(R).^x))-3.4 ,1 ,opts);

[Sven]

thank you Teja, this works fine for me....

[Sven]

Teja, I have changed the above to be part of a custom formula that reads like follows:

function [vRand ] = fPerfectRandomFigures(iNumber,dStd,dSkew,dKurt)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

rng('shuffle')
R = randn(iNumber,1);
R_prev= 0;
opts = optimset('TolFun',0,'Display','Off');
smallNumber = 1e-16;

while norm(R-R_prev)/norm(R_prev) > smallNumber;
        R_prev = R;
      R = R-mean(R); %Fix Mean
      R = R/std(R) * dStd; %Fix Std. Dev      
      x = fzero( @(x) skewness( R .* exp(x*R) )-dSkew ,0 ,opts); % Fix Skewness
      R = R .* exp(x*R);      % Fix Kurtosis
      x = fzero( @(x) kurtosis(sign(R) .* (abs(R).^x))-(3+dKurt) ,1 ,opts);
      R = sign(R) .* (abs(R).^x);
end

vRand=R;

end

I now would like to extend this a bit further to also allow to control for the serial correlation of the random numbers (1st order serial correlation). I would like to add a new parameter that can take any value between -1 and 1 and the vector of random numbers would be adjusted to have a serial correlation of that value. How would I do that?

Thanks for any hint...

Sven

[süleymand demir]

Hi Teja,

I wanted to use this code to genarate some datas. But I could not generate data for some conditions.

For example: a mean of 0.0000, a std of 1.0000, a skewness of 1.0000, and a kurtosis of 4.0000.

Could you help me please.

Answer 2

[Walter Roberson]

You cannot do this without constraining two of the values to be non-random (values that fix the mean and std.)

The Normal Distribution is defined by an infinite process; if you use a finite number of values, you can only approximate the Normal Distribution.

[Jan]

I even do not think, that any set of random numbers can have exactly the mean 0.0, most of all if the std is fixed to exactly 1.0 also. Therefore solving this problem requires to redefine "random" exactly at first.

[Walter Roberson]

Roger has a FEX contribution for a "random" set of numbers with a fixed sum. That technique can be used with the fixed sum being 0. Call that T. Then T/std(T) will have a sum of 0 and a std of 1.

I have difficulties accepting that the numbers are truly "random" if they have a fixed sum, but if Roger finds it meaningful then I consider that perhaps there is something lacking in my understanding.

[Jan]

I think, Roger applies a specific definition of "random" here: The maximal achievable indepedence between the values might be a solid base.

Considering the total interconnectedness of the properties in the universe, random is only vaguely defined at all. Einstein, Podolsky and Rosen hesitated to accept, that the spins of two conjugated photons can be oriented randomly although the sum of the spins will vanish after the measurement.

Answer 3

[Iain]

It isn't possible.

Just looking at the mean & std, if you do manage to get the mean to be exactly 0 and the std to be exactly 1, how can you possibly add another sample without either changing the mean (and retaining the std), or changing the std (and retaining the mean). Skewness and kurtosis just make it even harder.

Answer 4

[Sven]

To all that commented: thanks for your comments. I guess I should have been a bit more precise in my question regarding being perfect: I was referring to almost perfect only...errors in the 3rd or 4th decimal point as by using this code is absolutely fine for what I am after.

[Matt J]

In that case, why do you say randn isn't close enough? Using a sufficiently large sample set, the sample moments I get are well within 3 or 4 decimal places,

     x=randn(1,1e7);

    moment1=mean(x),
    moment2=std(x),
    moment3=mean(x.^3),
    moment4=mean(x.^4),

moment1 =

    -2.2094e-04

moment2 =

      0.9998

moment3 =

     -0.0012

moment4 =

      2.9995

[Sven]

while it is certainly possible to extend the number of draws to a very high number, I am trying to keep the size of the sample rather small.

[Matt J]

Whether the number of draws is large or small is irrelevant. The statistical mean and other moments are what they are independently of how you sample them. My test with large draws was only to measure the statistical moments better and show you that they have the values they should.

If you're saying that you want the empirical moments to deviate from the statistical moments by a limited amount for a certain sample size, then you have no control over that, I'm afraid. The empirical moments are random variables and will have a variance that depends on your sample size and the statistical variance. You cannot change the empirical moments' variances without first changing the statistical distribution of the things you're sampling.

Answer 5

[Iain]

Hmmm:

Take the inverse cdf, at regular intervals.

 rnds = icdf('norm',(1/samples):(1/samples):(1-1/samples), 0 , 1);

These numbers are not random, but they, as a set, have a mean of 0, std of close to 1, a skewness of 0, and a kurtosis of "close to" 3.

Generate a random order:

 [~, I] = sort(randn(size(rnds)));
 randoms = rnds(I);

[Sven]

Thanks lain, the first part works fine for generating the initial set of random numbers, but I am not sure I follow how the second part of your suggested answer would enforce a prespecified autocorrelation on the vector (see my comment from June 12th above), or was that not part of your answer? Sven

[Iain]

That was not part of my answer.