A question about the distribution of random numbers

[spunky]

No it actually does mean division in this case. But the fact is that you don't get 0. You get

fail on my part for not reading the order of operations properly. CDF FIRST, then INVERSE.

darn it, but thanks for the correction!

[BigBugBuzzz]

I believe your original question can be thought of equivalently as the maximum of A,B.

Therefore take the product of the two CDFs (in one variable, say, X) and then take the derivative of the resulting product with respect to X and that should be your PDF.

My skill at statistics could be much better, and I am working on it, but right now, I don't know how to arrive at the PDF. I bought a licence for Mathematica, and I am working on learning that skill. If anyone on here would be kind enough to write the Mathematica code for the solution give by either Dragan above or BGM (or both), I would be extremely appreciative. That way, I know I haven't made an error, and I will be able to study the code and learn exactly what is going on.

And thank you for all the help so far!

[Dragan]

My skill at statistics could be much better, and I am working on it, but right now, I don't know how to arrive at the PDF. I bought a licence for Mathematica, and I am working on learning that skill. If anyone on here would be kind enough to write the Mathematica code for the solution give by either Dragan above or BGM (or both), I would be extremely appreciative. That way, I know I haven't made an error, and I will be able to study the code and learn exactly what is going on.

And thank you for all the help so far!

This is what I mean in Mathematica. A hypothetical with MuA=10, Sigma=3; and MuB=9, Sigma=2

phiA=CDF[NormalDistribution[10,3], X]
phiB=CDF[NormalDistribution[9,2], X]
u=phiA*PhiB
v=D[u,X]

v is the derivative of u (the product of the two CDFs) and should be you probability density function.

[BigBugBuzzz]

Thank you Dragan. I have implemented the code in Mathematica and plotted the PDF:

phiA = CDF[NormalDistribution[10, 1], X]
phiB = CDF[NormalDistribution[10, 5], X]
u = phiA*phiB
v = D[u, X]

Plot[(E^((12*(-10 + X)^2)/25)*Erfc[-((-10 + X)/Sqrt[2])] +
5*Erfc[-((-10 + X)/(5*Sqrt[2]))])/(10*E^((-10 + X)^2/2)*
Sqrt[2*Pi]), {X, 4.6, 15}]

The result is the "Dragan" plot attached.

Comparing that to the distribution derived from simulation, I think there is something wrong with the proposed solution. This distribution is also attached, and called "Simulation".

Is BGM correct? What would his Mathematica code look like? I realize this is asking a lot...

[BigBugBuzzz]

I would very much appreciate it if someone would post the Mathematica code for this...

[Dragan]

phiA = CDF[NormalDistribution[10, 1], X]
phiB = CDF[NormalDistribution[10, 5], X]

This distribution is also attached, and called "Simulation".

Is BGM correct? What would his Mathematica code look like? I realize this is asking a lot...

I'm not getting the empirical (or simulation) result that your getting...How are you performing your data generation and selection process?

Note: The approach I gave for the MAX(A,B) is correct - the simulation matches the theoretical.

[BigBugBuzzz]

I'm not getting the empirical (or simulation) result that your getting...How are you performing your data generation and selection process?

I generate two vectors of random numbers, one columns from the Normal Distribution A and one columns from the Normal Distribution B.
Then I identify all those instances where A > B (pair-wise).
I finally import these to SPSS and draw the histogram shown.

Example:

A B
2 1*
3 2*
1 2

So, 2 and 3 form part of the distribution examined...

[Dason]

@Dragan - I don't think max(A, B) is really the same thing as what they're doing. I didn't check all of the steps but BGMs derivation of the pdf looked reasonable to me.

[Dragan]

I would very much appreciate it if someone would post the Mathematica code for this...

Okay, I just looked at BGM's post, which I think is correct. What I suggested is a special case of what BGM proposed when both distributions are the same.

Anyway, here's the Mathematica code - it should work.

ma=10;
mb=10;
sa=1;
sb=5;

con=(ma-mb)/(Sqrt[sa^2 + sb^2];

cmp=CDF[NormalDistribution[0,1],con];

FB=CDF[NormalDistribution[mb,sb],x];

fA=PDF[NormalDistribution[ma,sa],x];

fx=(cmp^-1)*(1/sa)*FB*fA

Plot[fx,{x,6,14}]

[BigBugBuzzz]

Anyway, here's the Mathematica code - it should work.

I really appreciate you taking the time Dragan. Unfortunately, while the code works, this does not match the simulated distribution either. I guess BGM's solution was not correct after all.

I have attached the image of the distribution for the means and variances values you set , which match the means and variances used to generate the simulated distribution. As you will see, the distributions are not the same.

Perhaps this nut is too hard to crack...

Here is the reason why I think the simulated distribution is correct: The negative skew in the simulated distribution makes sense when the variance of A > B and means are equal, because larger values of A are allowed to "pass" at greater rate - even if smaller than muA - compared to the situation where the distributions are equal. The skewness of the distribution is intuitively determined by the relation between means and variance:

muA = muB, sA = sB --> No Skew
muA = muB, sA > sB --> Negative Skew

I am yet to examine:

muA > muB, sA = sB
muA > muB, sA > sB

I think the problem is interesting and I wish I had the capabilities to formalize, rather than just simulate! :-)

[BGM]

I think we may misunderstand your original problem set up.

Here is what I try with R:

mA<-10
mB<-10
sA<-1
sB<-5

con<-(mA-mB)/sqrt(sA^2+sB^2)
cmp<-pnorm(con)
FB<-function(x) pnorm(x,mB,sB)
fA<-function(x) dnorm(x,mA,sA)
f<-function(x) FB(x)*fA(x)/cmp

n<-100000
d<-1:n
for (i in 1:n) {
	A<-0
	B<-1
	while (A < B) {
		A<-rnorm(1,mA,sA)
		B<-rnorm(1,mB,sB)
	}
	d[i]<-A
}

hist(d,freq=F)
curve(f,add=T)

And the resulting empirical distribution seems match with the proposed density.

[Dragan]

I think we may misunderstand your original problem set up.

Yes, BGM, I suspected (yesterday) that there's a problem "somewhere" when I could not match BigBugBuzz's simulation result. I cannot say for sure where the problem is.

[BigBugBuzzz]

Perhaps I am not so clear in my formulation. But consider that B is just a line below the mean of A. Then C would surely not be symmetrical but would be cut sharply at B. Now consider that B is almost a line, but is in fact a very thin probability distribution, then C would still not be symmetrical, but would not be cut so sharply at the left, and so on. Does that make sense?

Attached is a visual example of the selection of values for C, in case that makes things more clear.

[BigBugBuzzz]

HOLD THE PRESS!! You guys appear to be right after all! Note that if sA > sB we get asymmetry. If sB > sA symmetry. Beautiful! Thank you so much for the help! :-)

I used Dragan's code and just switched sA and sB....

Dragan's solutions appears to be correct. Do you agree BGM?

[BGM]

From the pdf derived we can see that it is asymmetric - the (normal) CDF is a monotonic increasing function.