conditional generation of random numbers using matlab

Question

I have a function that generates normal random number matrix having normal distribution using normrnd.

             values(vvvv)= normrnd(0,0.2);

The output is from round1 is: ans =

0.0210    0.1445    0.5171   -0.1334    0.0375   -0.0165       Inf   -0.3866   -0.0878   -0.3589

The output from round 2 is: ans =

0.0667    0.0783    0.0903   -0.0261    0.0367   -0.0952    0.1724   -0.2723       Inf       Inf

The output from round 3 is: ans =

0.4047   -0.4517    0.4459    0.0675    0.2000   -0.3328   -0.1180   -0.0556    0.0845       Inf

the function will be repeated 20 times.

It is obvious that the function is completely random. What I seek is to add a condition.

What I need is: if any entry has a value between 0.2 and 0.3. that value will be fixed in the next rounds. Only the remaining entries will be subjected to change using the function rand.

I have found the rng(sd) which seeds the random number generator using the nonnegative integer sd so that rand, randi, and randn produce a predictable sequence of numbers.

How to set custom seed for pseudo-random number generator

but how to make several entries of the matrix only effected!!

Another problem: seems that rng is not available for matlab r2009 How to get something similar without entering in the complication of probability & statistics

Answer 1

You can do this more directly than actually generating all these matrices, and it's pretty easy to do so, by thinking about the distribution of the final output.

The probability of a random variable distributed by N(0, .2) lying between .2 and .3 is p ~= .092.

Call the random variable of the final output of your matrix X, where you do this n (20) times. Then either (a) X lies between .2 and .3 and you stopped early, or (b) you didn't draw a number between .2 and .3 in the first n-1 draws and so you went with whatever you got on the nth draw.

The probability of (b) happening is just b=(1-p)^(n-1): the independent events of drawing outside [.2, .3], which have probability 1-p, happend n-1 times. Therefore the probability of (a) is 1-b.

If (b) happened, you just draw a number from normrnd. If (a) happened, you need the value of a normal variable, conditional on its being between .2 and .3. One way to do this is to find the cdf values for .2 and .3, draw uniformly from the range between there, and then use the inverse cdf to get back the original number.

Code that does this:

mu = 0;
sigma = .2;
upper = .3;
lower = .2;
n = 20;
sz = 15;

cdf_upper = normcdf(upper, mu, sigma);
cdf_lower = normcdf(lower, mu, sigma);
p = cdf_upper - cdf_lower;
b = (1-p) ^ (n - 1);

results = zeros(sz, sz);
mask = rand(sz, sz) > b; % mask value 1 means case (a), 0 means case (b)
num_a = sum(mask(:));

cdf_vals = rand(num_a, 1) * p + cdf_lower;
results(mask) = norminv(cdf_vals, mu, sigma);

results(~mask) = normrnd(mu, sigma, sz^2 - num_a, 1);

---

If you want to simulate this directly for some reason (which is going to involve a lot of wasted effort, but apparently you don't like "the complications of statistics" -- by the way, this is probability, not statistics), you can generate the first matrix and then replace only the elements that don't fall in your desired range. For example:

mu = 0;
sigma = .2;
n = 10;
m = 10;
num_runs = 20;
lower = .2;
upper = .3;

result = normrnd(mu, sigma, n, m);
for i = 1 : (num_runs - 1)
    to_replace = (result < lower) | (result > upper);
    result(to_replace) = normrnd(mu, sigma, sum(to_replace(:)), 1);
end

---

To demonstrate that these are the same, here's a plots of the empirical CDFs of doing this for 1x1 matrices 100,000 times. (That is, I ran both functions 100k times and saved the results, then used cdfplot to plot values on the x axis vs portion of the obtained values that are less than that on the y axis.)

They're identical. (Indeed, a K-S test for identity of distribution gives a p-value of .71.) But the direct way was a bunch faster to run.