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.