On 19-Feb-11 14:48:53, David Winsemius wrote:
> On Feb 19, 2011, at 5:47 AM, danielepippo wrote:
>> Hi,
>> I have two vector with the marginal distribution like this:
>>> a
>> [1] -0.419 -0.364 -0.159 -0.046 -0.010 -0.002 0.000 0.000 0.000
>>> b
>> [1] 0.125 0.260 0.270 0.187 0.097 0.041 0.014 0.004 0.001
>>
>> How can I calculate the joint distribution with R?
>
> Marginal distributions do not uniquely determine the joint
> distribution, Furthermore, the first one doesn't even look like a
> distribution or even a density. Densities are positive and CDFs range
> from 0 to 1. (The second on could be a discrete density for some set
> of values.) So you need to explain where those numbers come from and
> why you think we should apply sort of "distributional" assumptions
> about them.
> --
> David Winsemius, MD
> West Hartford, CT
Following up on David's reply (I was in the midst of thinking about
this when he responded):
I strongly suspect that the values of 'a' have had their signs
changed. If we simply make them all positive:
a <- c(0.419,0.364,0.159,0.046,0.010,0.002,0.000,0.000,0.000)
then now sum(a) = 1 and these can be probabilities. Also, with
the values given, sum(b) = 0.999; so (in the relatively least
disturbing way) I shall increase the largest value by 0.001:
b <- c(0.125,0.260,0.271,0.187,0.097,0.041,0.014,0.004,0.001)
(0.270 -> 0.271). Now also sum(b) = 1.
Now, as a somewhat trivial example to illustrate the non-uniqueness
that David points out, I shall invent a very arbitrary set of
probabilities P = {P[i,j]} such that they have the given marginal
distributions.
Since we have the marginals given to 3 decimal places, now consider
allocating 1000 pairs of random numbers (x,y). Then we pick out the
smallest 419 of the x's to go into the bottom class of 'a'
(proportion = 419/1000 = 0.419), then the 364 next smallest for the
second class of 'a', and so on. And similarly for the numbers 'y'
to go into the classes of 'b'.
Then the inequalities which define these marginal divisions can be
applied jointly to produce the joint divisions. In the first instance
the results will be computed as counts
Therefore:
X0 <- runif(1000) ; X<- sort(X0)
Y0 <- runif(1000) ; Y<- sort(Y0)
A <- cumsum(c(419,364,159, 46, 10, 2, 0, 0, 0))
B <- cumsum(c(125,260,271,187, 97, 41, 14, 4, 1))
Xdiv <- numeric(8)
Ydiv <- numeric(8)
for(i in (1:8)){
Xdiv[i] <- 0.5*(max(X[1:A[i]])+min(X[(A[i]+1):A[i+1]]))
Ydiv[i] <- 0.5*(max(Y[1:B[i]])+min(Y[(B[i]+1):B[i+1]]))
Xdiv[is.na(Xdiv)] <- 1.0
Ydiv[is.na(Ydiv)] <- 1.0
}
Xdiv <-c(0,Xdiv,1)
Ydiv <- c(0,Ydiv,1)
N <- matrix(0,nrow=9,ncol=9)
for(i in (1:9)){
for(j in (1:9)){
N[i,j] <- sum((Xdiv[i]<X0)&(X0<=Xdiv[i+1])&
(Ydiv[j]<Y0)&(Y0<=Ydiv[j+1]))
}
}
N
rowSums(N)
colSums(N)
Here is one typical result (depending on the random numbers) for N:
N
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,] 54 113 112 76 36 20 7 1 0
# [2,] 44 90 101 75 35 13 3 2 1
# [3,] 20 40 42 28 21 6 2 0 0
# [4,] 5 13 13 6 5 2 1 1 0
# [5,] 2 3 3 1 0 0 1 0 0
# [6,] 0 1 0 1 0 0 0 0 0
# [7,] 0 0 0 0 0 0 0 0 0
# [8,] 0 0 0 0 0 0 0 0 0
# [9,] 0 0 0 0 0 0 0 0 0
rowSums(N)
# [1] 419 364 159 46 10 2 0 0 0 ## Agreeing with (A)
colSums(N)
# [1] 125 260 271 187 97 41 14 4 1 ## Agreeing with (B)
And here is another:
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,] 52 109 113 89 35 15 3 3 0
# [2,] 45 95 108 50 40 16 9 0 1
# [3,] 22 44 33 33 16 9 1 1 0
# [4,] 4 9 14 13 4 1 1 0 0
# [5,] 2 3 3 1 1 0 0 0 0
# [6,] 0 0 0 1 1 0 0 0 0
# [7,] 0 0 0 0 0 0 0 0 0
# [8,] 0 0 0 0 0 0 0 0 0
# [9,] 0 0 0 0 0 0 0 0 0
rowSums(N)
# [1] 419 364 159 46 10 2 0 0 0
colSums(N)
# [1] 125 260 271 187 97 41 14 4 1
Then N/1000 gives the joint probability distribution, and
rowSums(N) and colSums(N) give the marginal distributions.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.hard...@wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 19-Feb-11 Time: 16:30:08
------------------------------ XFMail ------------------------------
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
