Robert,
That worked.
 
Beautiful - thank you so much for your help!!

Cheers,
Diann
________________________________


The warning message still appears (but the matrix is written, and output 
look ok).

> x <- calc(s, fun)
Warning messages:
1: In runif(nnfin, min = pinf, max = psup) : NAs produced
2: In function (n, min = -1, mode = 0, max = 1)  : NaN in rtriang

Because the runif and rtriang components of the warning message made me 
concerned as to whether the mcstoc subfunctions were calling the correct 
parameters as opposed to the defaults printed in the warning, I did a 
number of diagnostics on the function to determine that indeed the correct 
parameters were being used. All looks good.




From:
"Robert J. Hijmans" 
To:
Diann J Prosser
Date:
03/13/2012 12:08 PM
Subject:
Re: [R-sig-Geo] apply monte carlo simulation for each cell in a matrix 
(originally raster)



Diann,
I think you can fix that by using:

quantile(risk[], na.rm=TRUE)

Robert

On Tue, Mar 13, 2012 at 12:43 AM, Diann Prosser  wrote:

Hi Robert, 

Thank you - the upgrade did fix the S4 issue. 

Kind regards, 
Diann 


FYI, now with new NA issues not seen before. 

I still can't get quantile(risk[])  to work: 

> # Function that works for a vector of three numbers 
> ndunc(1001) 
[1] 1001 
> fun <- function(x) { 
+ dPmn <- mcdata(x[1], type="0") 
+ dPsd <- mcdata(x[2], type="0") 
+ dWab <- mcdata(x[3], type="0") 
+ stPte <- mcstoc(rnorm, type="U", dPmn, dPsd, rtrunc=TRUE, linf=0) 
+ stWab <- mcstoc(rtriang, type="U", min=0, mode=dWab, max=75) 
+ risk <- stPte * stWab 
+ quantile(risk[]) 
+ } 
> x <- calc(s, fun) 
Error in .local(x, values) : 
  length(values) is not equal to ncell(x), or to 1 
In addition: Warning messages: 
1: In runif(nnfin, min = pinf, max = psup) : NAs produced 
2: In function (n, min = -1, mode = 0, max = 1)  : NaN in rtriang 

But c(min(risk), mean(risk), max(risk))  seems to. 
... 
+ c(min(risk), mean(risk), max(risk)) 
+ } 
> x <- calc(s, fun) 
Warning messages: 
1: In runif(nnfin, min = pinf, max = psup) : NAs produced 
2: In function (n, min = -1, mode = 0, max = 1)  : NaN in rtriang 
> x <- writeRaster(x, 'mc.tif', overwrite=TRUE)  # for ArcGIS 
> as.matrix(x) 
          mc.1      mc.2     mc.3 
[1,] 11.720980 140.53212 271.2712 
[2,]  5.760174 125.81305 224.5858 
[3,]  3.631572  74.98318 164.4021 
[4,]  3.637625 128.85443 239.7214 
> 

___________________________________________________________________ 
R version 2.14.2 (2012-02-29) 
Platform: i386-pc-mingw32/i386 (32-bit) 

locale: 
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United 
States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                   
        
[5] LC_TIME=English_United States.1252     

attached base packages: 
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages: 
[1] mc2d_0.1-12      mvtnorm_0.9-9992 rgdal_0.7-8      raster_1.9-70   
[5] sp_0.9-97       

loaded via a namespace (and not attached): 
[1] grid_2.14.2    lattice_0.20-6 tools_2.14.2  


From: 
"Robert J. Hijmans"
To: 
Diann Prosser 
Date: 
03/13/2012 01:52 AM 
Subject: 
Re: [R-sig-Geo] apply monte carlo simulation for each cell in a matrix 
(originally raster)

Dianne, see inline comments. Best, Robert

On Mon, Mar 12, 2012 at 10:30 PM, Diann Prosser <dprosser@usgs.gov> wrote: 

Robert, thank you very much! I believe this sets me on the right path.

2 issues I'm looking into:

1) I do receive an S4 error after the as.matrix(x) step:
> as.matrix(x)
/Error in as.vector(data) :
 no method for coercing this S4 class to a vector/ 


That is strange. Are you perhaps using an old version of 'raster' ?  Can 
you send your sessionInfo() 



2) Perhaps related, I'm not sure if I'm losing the variance measures 
around
the estimate (mean risk) or if I'm not using the correct function to
retrieve them. For example, "summary(node)" should give me the median, 
mean,
2.5% and 97.5% variance measures for each cell. "summary(x)" gives me:

> summary(x)
Cells:  4
NAs  :  0
Min.     77.36
1st Qu. 114.10
Median  126.70
Mean    118.30
3rd Qu. 130.90
Max.    142.30


which is information summarized across the matrix.
            [,1]     [,2]
[1,] 142.13677 125.7188
[2,]  75.59303 128.8319



Instead of having the function return 
mean(risk) 
you can do things like: 
quantile(risk[]) 
or  
c(min(risk), mean(risk), max(risk)) 


as in the below: 


fun <- function(x) { 
 dPmn <- mcdata(x[1], type="0") 
 dPsd <- mcdata(x[2], type="0") 
 dWab <- mcdata(x[3], type="0") 
 stPte <- mcstoc(rnorm, type="U", dPmn, dPsd, rtrunc=TRUE, linf=0) 
 stWab <- mcstoc(rtriang, type="U", min=0, mode=dWab, max=75) 
 risk <- stPte * stWab 
 # now what output do you want?  
 # there are many options, e.g.: 
 # quantile(risk[]) 
      # or simply: 
 #mean(risk) 
 c(min(risk), mean(risk), max(risk)) 
} 
  
  
x <- calc(s, fun) 
layerNames(x) = c('min', 'mean', 'max') 
plot(x) 
  

  
Ultimately, I want matrix/raster outputs for mean risk, plus the upper and
lower bounds.
Will keep digging. Thank you for your time!



--
View this message in context: 
http://r-sig-geo.2731867.n2.nabble.com/apply-monte-carlo-simulation-for-each-cell-in-a-matrix-originally-raster-tp7367094p7367555.html 

Sent from the R-sig-geo mailing list archive at Nabble.com.

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo 






	[[alternative HTML version deleted]]

