Hi

I have a problem that I am hoping someone can help me with.


I have been using the raster package and some loops to use a normal
probability density function to assess likelihoods that a sample in a
data set could have come from a location on a map (example script
below). 

library(raster)
##### example map
map <- raster(nrows=20, ncols=20)
map[] <-rnorm(ncell(map), -89, 15)
xmin(map)<- -5
xmax(map)<- 45
ymin(map)<- 25
ymax(map)<- 75

plot(map)

###### values (data) to be assess likelihoods for 
x <- rnorm(10,mean=-74, sd= 15)

################ Define the variation 
SD <-11.67114

######establish a set of lists to hold the results of subsequent loops
r <- list()
############ Apply Normal Probability Density Function to Assess
Likelihoods ####################
for (i in 1:length(x)){
    r[i] = (1/(sqrt(2*pi)*SD)*exp(-(1/(2*(SD)^2))*(x[i]-map)^2))
}

############ Normalize by Sum of Values....yielding map of probabilities
that sums to 1 ####################
for (i in 1:length(x)){
r[i] <- r[[i]] / cellStats(r[[i]],'sum')
}

##example plot
plot(r[[1]])


This method is working beautifully to assess the likelihoods, but I now
want to come up with a way of accurately selecting and reclassifying the
resulting probability maps cells into binary likely (1) versus unlikely
(0) maps based on selecting all cells in the map above a threshold
probability value, such that the sum of the probabilities in the
"likely" category is equal to a specified cumulative probability (e.g.
67%, 95%). I have previously been using some course approximations, but
now want a highly precise., accurate and efficient method to do this,
but appear to be stuck.

If this were a data frame, I could just sort based on the estimated
probabilities and calculate a new cumulative probability variable(s). I
could probably then use a spline function to estimate the minimum
probability at which the selected cumulative probability of interest was
reached (similar to the following post
(http://finzi.psych.upenn.edu/R/Rhelp01/archive/13524.html)). However, I
am not certain how I can do this for raster objects. 

Ideally I would like to apply a spline function to find the desired
threshold value and directly incorporate the result into a loop that
would reclass each r[i] into a binary map, 
e.g. r[[1]] <-reclass(r[[1]], c(-Inf, 0.000694, 0, 0.000694, Inf, 1))
where 0.000694 was pasted from the lookup spline function 

Then I want to sum the "likely" areas across all samples in the data set
E.g.
s<-stack(r)
i<-maxValue
rc <- overlay(s, fun=function(i){return(sum(i))} ) 



I would really appreciate it if anyone can provide advice on 
A) how to estimate the threshold value for each sample
B) how I could add a loop and paste the result from A into a list of
reclass definitions (a separate one for each sample) 

I am hoping to find the most efficient method possible to accomplish
this, especially since x for some of my applications can be very long
(e.g. 1000 samples). 

-Steve 






	[[alternative HTML version deleted]]

