How to draw random number from a conditional distribution

Hi, dear friends,

I am working a drawing 10000 random numbers for variable, X, based on the following two conditions:

condition 1: X follows a normal distribution with mean=12.21 and standard deviation = 0.63;

condition 2: given another variable Y, the probability of X >= Y is 0.97: prob (X(n)>=Y(n) =0.97).

I know how to draw a random variable with specific mean and standard deviation, but I don't know how to plug in the second condition to my code.

Can anyone give me some help?

I do really appreciate.

Thank you again

Best,

Yanting

[the cyclist]

I don't understand condition 2. What does n represent? What is X(n)?

Regardless, your condition 1 fully specifies the distribution. There is only one distribution that can obey all those conditions. So, I don't see how condition 2 can also be satisfied, in general.

Or maybe I am misunderstanding something.

Edit: Maybe Adam's answer is actually what you meant, in which case you can generate values using rejection sampling.

[yanting wu]

Hi,

Thank you for your quick reply.

The second condition is that: I have another variable Y, given the value of each nth Y, my nth variable X should follow a normal distribution and should have 97% of chance greater or equal than Y.

I did the following matlab loop and it takes me almost infinite time to run, so I want to add the second condition into each n.

The matlab code I used is:

Y=randn(10000,1) * 0.61+12.01;
count=zeros(10000,1);
control =1;
while control==1
X=randn(10000,1) * 0.63+12.21;  
for n1=1:10000
if X>=Y
    count(n1)=1;
else
    count(n1)=0;
end
s=sum(count);
pn=s/10000;
if pn<0.97
    control=1
else 
    control=control+1
end
end   
end

Answer 1

[Adam Danz]

The first condition is easy,

X = randn(1,100000) * .63 + 12.21;
figure
histogram(X)
mean(X)  % <--- ~12.21
std(X)   % <--- ~0.63

To find the cuttoff of the distribution that defines 97% of the data on the right and 3% of the data on the left, you need to find the 3rd percentile of the data. Any value in X has a 97% chance of being greater than 'pt'.

pt = prctile(X,3);  %stats toolbox

Now add a vertical reference line to the plot and convince yourself that 97% of the data are greater than this value.

xline(pt, 'k-', '3rd percentile')
mean(X > pt)  %  <--- ~97%

[yanting wu]

Hi, Adam,

Thank you for your quick reply.

The second condition is that: I have another variable Y, given the value of each nth Y, my nth variable X should follow a normal distribution and should have 97% of chance greater or equal than Y.

I did the following matlab loop and it takes me almost infinite time to run, so I want to add the second condition into each n.

The matlab code I used is:

Y=randn(10000,1) * 0.61+12.01;
count=zeros(10000,1);
control =1;
while control==1
X=randn(10000,1) * 0.63+12.21;  
for n1=1:10000
if X>=Y
    count(n1)=1;
else
    count(n1)=0;
end
s=sum(count);
pn=s/10000;
if pn<0.97
    control=1
else 
    control=control+1
end
end   
end

[Adam Danz]

Let's get this straight before we think about code.

Tell me if I got this right or where I get it wrong.

1. You have a vector X.
2. X are values from a normal distribution with a known mean (12.21) and known std (0.63).
3. You have a vector Y the same length as X.
4. Given Y(n), X(n) should be a member of the distribution described in [1,2] but X(n) needs to have a 97% chance of being greater or equal to Y(n).

If that's what you're trying to do, unless Y is a vector of the same value which is equal to the 3rd percentile of X, your goal is impossible. The distribution X only has 1 value where 97% of the data are greater than that value.

Perhaps a visualization would help. The blue distribution below is your vector X (100k values) with mean 12.21 and std 0.63. The black line is the 3rd percentile (11.02) where 97% of the distribution is on the right and 3% is on the left. 11.02 is the only value that satisfies that requirement for this distribution.

[yanting wu]

Thank you, Adam,

Your response is very clear.

The original code I code is to generate 10000 observations for two normal distributed variables X (mean= 12.21, SD= 0.63) and Y (mean =12.01, SD=0.61). And in the created 10000 X and Y, I need at least 9700 X is greater or equal to Y, otherwise the matlab need to repeat the process and draw new 10000 observation of X and Y.

It takes almost infinite time to do the creating, so I want to say if I can add the condition 2 in each 10000 observation generating process: the probability of X is greater or equal to Y is at least 97%.

Thank you

Best,

Yanting

[yanting wu]

Thank you Adam,

I think I understand your point. Based on the given distribution of X and Y, it is impossible for X to have 97% of time to be greater or equal to Y if we draw them randomly.

Thank you

[Adam Danz]

" I need at least 9700 X is greater or equal to Y"

How are you defining that? When you draw 'n' values from both distributions, are you just comparing X(n) to Y(n)? If so, is there a reason for that comparison?

If you sort the values first, you'll definitely get more than 97% of X>Y.

mean(X >= Y)*100             % unsorted percentage of X>Y  ~59%
mean(sort(x) >= sort(y))*100 % sorted percentage of X>Y    ~99%

Here's a demo.

% Define distributions
mu = [12.21, 12.01];
sig = [.63, .61];

% run 1000 repetitions of the test
ntests = 1000;
p = zeros(1,ntests); 
ps = p; 
wbh = waitbar(0,sprintf('%d repetitions',ntests)); 
for i = 1:ntests
    d1 = randn(1,10000) * sig(1) + mu(1);
    d2 = randn(1,10000) * sig(2) + mu(2);
    p(i) = mean(d1>=d2);               % unsorted
    ps(i) = mean(sort(d1)>=sort(d2));  % sorted
    waitbar(i/ntests,wbh)
end
delete(wbh)

% Plot results
figure()
plot(p*100, 'b-', 'DisplayName', 'unsorted')
hold on
plot(ps*100, 'r-', 'DisplayName', 'sorted')
xlabel('repetition number')
ylabel('X>Y (%)')
ylim([-inf, 105])
legend()

[Adam Danz]

@yanting wu , sorry for changing my comment while you were responding. I had a 2nd thought about your goals which is described above.

[yanting wu]

Hi, Adam,

Yes, I need to compare X(n) with Y(n). In total I need 10000 pairs of observations (X , Y), and at least 9700 pairs of them with X>= Y. Because Y are created in the previous step with normal distribution, I can consider Y(n) as a constant number for generating X.

Thank you

Yanting

[Adam Danz]

Yeah, that will be impossible as the figure shows in an earlier comment of mine (blue line). You'll never get 97% values in X>=Y if x and y are numbers from normal distributions with the means and std's you specified.

Answer 2

[Bruno Luong]

Your problem is not well defined, there are many Y that can meet such requirement (if you guys think it is impossible beside a scalar you are wrong). One these is given by this code:

mux = 12.21;
sigma = 0.63;

muy = fzero(@(muy) integral2(@(x,t) pdffun(x,t,muy), -Inf, Inf, 0, Inf)-0.97, 0);
muy = mux + sigma*muy; % 10.534298385905991 you can enter this value without using FZERO

X = mux + sigma*randn(1,1e6);
Y = muy + sigma*randn(1,1e6);

% Check how many X >= Y, it should be about 0.97
sum(X>=Y) / length(X)

function p = pdffun(x,t,muy)
y = x - t - muy;
p = 1/(2*pi)*exp(-x.^2/2).*exp(-y.^2/2);
end

[Adam Danz]

"Your problem is not well defined, there are many Y that can meet such requirement (if you guys think it is impossible beside a scalar you are wrong)."

Agreed that the problem is not well defined but some valid interpretations of the problem (including the correct interpretation, according to OP) would result in impossibilities as I've shown above.

Interpretation 1

The distribution X only has 1 value where 97% of the data are greater than that value (the 3rd percentile). So it would be impossible to have a normal distribution with a known mean and known std and have multiple y values that define that boundary.

Interpretation 2

If 'n' random Ys are drawn from distribution y and 'n' random Xs are drawn from a distribution x, 97% of the Y(1:n) should be greater than X(1:n). Given the means and std's the OP defined, I've shown that this is indeed impossible and plotted the simulated results across 1000 repetitions where there were never more than ~60% Y values greater than their X pair.

Note that in your solution, the two distributions do not have the means and std's given by OP: "...generate 10000 observations for two normal distributed variables X (mean= 12.21, SD= 0.63) and Y (mean =12.01, SD=0.61)."-OP