Version 3.4 Release 1 for Sun SPARC, SunOS 5.3 : 1996 (I think it
immaterial, but I have the SpatialStats module).
I have a simple function which calculates error weighted mean and
error (see below). I would like to call this function using
"apply", but can't figure out how to pass 2 parameters. In this
application the input values and errors are 3D arrays and I want to
find the error weighted mean and error over the 3rd dimension.
Input is typically 3 or 4 matrices each of dimension 360x180 and I
have to execute the procedure repeatedly, so calculation speed is
a consideration (though not yet sufficient for the jump to FORTRAN).
Example:
Given data matrices x1 with errors e1 and x2 with errors e2
> x1
[,1] [,2]
[1,] 1 1
[2,] 1 1
> x2
[,1] [,2]
[1,] 1 2
[2,] 3 4
> e1
[,1] [,2]
[1,] 0.2 0.2
[2,] 0.2 0.2
> e2
[,1] [,2]
[1,] 1 1
[2,] 1 1
I'd like to do something like:
apply(list(X,E),c(1,2),errmean)
where X is x1 stacked on x2 and E is e1 stacked on e2 to create 3D
arrays. Note that whenever a value (e.g. x1[i,j] is NA the associated
error (e1[i,j]) is also NA so potential misalignment is not a problem
within function errmean.
--------------------------------------------
> errmean
function(x, err)
{
#
# x = vector of values
# s = vector of errors
# error weighted average and error weighted error
# returned
erave <- sum(x/err^2,na.rm=T)/sum(1/err^2,na.rm=T)
erstd <- sqrt(1/sum(1/err^2,na.rm=T))
normave <- aveerr(x)
list(erave = erave, erstd = erstd, ave = normave[1], std = normave[2])
}
----------------------------------------------
#"aveerr" is another simple function that calculates average and
standard deviation with NAs removed (using "mean" and "var")
-----------------------------------------------
Note that for simple matrix means I use the following function (matmean)
which was developed with help from this group
matmean_function(inlist)
{
#
# Calculate the element by element mean of matrices
# NA's are ignored
# Usage:
# matmean(list of matrices)
#
if(!is.list(inlist))
mats <- list(...)
else
mats_inlist
A <- array(unlist(mats), dim = c(dim(mats[[1]]), length(mats)))
apply(A, c(1,2), mean, na.rm = T)
}
--------------------------------------------------
Hopefully this all makes sense. If the only way to proceed is via
"for" looping, no help required.
sincerely,
bob
|