Scott - I agree with Spencer Graves that there's a scoping issue here: Where does function Dk() pick up the values for n0 and w, and does it get them from the SAME place when it's called from inside FindLikelihood() as from outside ?
But more important is this one: All arithmetic on vectors or matrices is done element by element; every matrix or array is treated as a vector (no 'dim' attribute) during this process, and "the elements of shorter vectors are recycled as necessary" (quoting from help("Arithmetic")). Therefore, Dk(seq(-k,k), 0.2, 0.2) should return a vector of length (2 * k + 1), and Nk * log(Dk()) # (I omit the arguments to Dk() here.) should produce a vector of length max(length(Nk), 2 * k + 1) in which element 1 of Nk is paired with xk = -k, element 2 of Nk is paired with xk = (-k + 1), et cetera. This product then has a vector of length (2 * k + 1) subtracted from it and the resulting vector is summed. Now, maybe you have promised to only call FindLikelihood() with an argument Nk of length 15 = (2 * k + 1), in which case all the lengths match and element i from Nk is always paired with the value (i - 8) in the first argument of Dk(), but there's certainly a lack of defensive programming here. An alternate way to calculate the grid of likelihood values which seems to be your intention is to explicitly build four four-dimensional arrays named A, B, xk and Nk, all with the same dimensions, and with the values changing along only one dimension in each array. Then do whatever arithmetic you want with these four arrays (such as the expressions inside Dk() and FindLikelihood() and collapse the result by summing over rows or slices or whatever at the end. The functions array(), aperm(), matrix() and '%*%' are useful in this process. This business of four or five-dimensional arrays is one I use routinely. The result is equivalent to as.vector(outer( ...)), but it forces you to think carefully about the various dimensions. HTH - tom blackwell - u michigan medical school - ann arbor - On Mon, 27 Oct 2003, Scott Norton wrote: > I'm pulling my hair (and there's not much left!) on this one. Basically I'm > not getting the same result t when I "step" through the program and evaluate > each element separately than when I use the outer() function in the > FindLikelihood() function below. > > > > Here's the functions: > > > > Dk<- function(xk,A,B) > > { > > n0 *(A*exp(-0.5*(xk/w)^2) + B) > > } > > > > FindLikelihood <- function(Nk) > > { > > A <- seq(0.2,3,by=0.2) > > B <- seq(0.2,3,by=0.2) > > k <-7 > > L <- outer(A, B, function(A,B) sum( (Nk*log(Dk(seq(-k,k),A,B))) - > Dk(seq(-k,k),A,B) )) > > return(L) > > } > > > > > > where Nk <- c(70 , 67 , 75 , 77 , 74 ,102, 75, 104 , 94 , 74 , 78 , 79 , 83 > , 73 , 76) > > > > > > Here's an excerpt from my debug session.. > > > > > Nk > > [1] 70 67 75 77 74 102 75 104 94 74 78 79 83 73 76 > > > debug(FindLikelihood) > > > L<-FindLikelihood(Nk) > > debugging in: FindLikelihood(Nk) > > debug: { > > A <- seq(0.2, 3, by = 0.2) > > B <- seq(0.2, 3, by = 0.2) > > k <- 7 > > L <- outer(A, B, function(A, B) sum((Nk * log(Dk(seq(-k, > > k), A, B))) - Dk(seq(-k, k), A, B))) > > return(L) > > } > > Browse[1]> n > > debug: A <- seq(0.2, 3, by = 0.2) > > Browse[1]> n > > debug: B <- seq(0.2, 3, by = 0.2) > > Browse[1]> n > > debug: k <- 7 > > Browse[1]> n > > debug: L <- outer(A, B, function(A, B) sum((Nk * log(Dk(seq(-k, k), > > A, B))) - Dk(seq(-k, k), A, B))) > > Browse[1]> sum((Nk * log(Dk(seq(-k, k),0.2,0.2))) - Dk(seq(-k, k), 0.2, > 0.2)) # WHY DOES THIS LINE GIVE ME THE CORRECT RESULT WHEN I SUBSTITUTE > 0.2, 0.2 FOR A AND B > > [1] 2495.242 > > Browse[1]> outer(A, B, function(A, B) sum((Nk * log(Dk(seq(-k, k), > > + A, B))) - Dk(seq(-k, k), A, B))) > > [,1] [,2] [,3] [,4] [,5] [,6] [,7] > [,8] > > [1,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 # BUT ELEMENT (1,1) WHICH SHOULD ALSO BE (A,B) = (0.2, 0.2), > GIVES THE INCORRECT RESULT???? > > [2,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [3,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [4,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [5,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [6,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [7,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [8,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [9,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [10,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [11,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [12,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [13,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [14,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [15,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > 58389.48 > > [,9] [,10] [,11] [,12] [,13] [,14] [,15] > > [1,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [2,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [3,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [4,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [5,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [6,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [7,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [8,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [9,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [10,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [11,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [12,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [13,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [14,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > [15,] 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 58389.48 > > Browse[1]> > > > > As "commented" above, when I evaluate a single A,B element (i.e. A=0.2, > B=0.2) I get a different result than when I use OUTER() which should also be > evaluating at A=0.2, B=0.2?? > > > > Any help appreciated. I know I'm probably doing something overlooking > something simple, but can anyone point it out??? > > > > Thanks! > > -Scott > > > > Scott Norton, Ph.D. > > Engineering Manager > > Nanoplex Technologies, Inc. > > 2375 Garcia Ave. > > Mountain View, CA 94043 > > www.nanoplextech.com > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > [EMAIL PROTECTED] mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help > ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help