Hi Michelle, Please cc the list on your replies so others (far more knowledgeable than I) can answer.
On Mon, Sep 24, 2012 at 12:12 AM, Michelle Cipullo <mlcipull at ncsu.edu> wrote:
I am doing some extreme value analysis on model output (WRF) which have the following dimensions: speed(time,lat,lon) I am trying to fit the GPD (gpd.fit) to each point (time,lat,lon) to get a return level plot with values at each grid point. (Map with return level by location.)
Finally, you should see if the function gpd.fit() is (in R speak) vectorized: if so, you will be able to drop the explicit loops and make things much much faster.
This function is from the ismev package: http://www.maths.lth.se/matstat/staff/nader/stint/R_Manuals/ismev.pdf It doesn't say explicitly (that I saw) if it is or isn't vectorized. When I tried doing the indexing the correct way (thanks!), I get the following error: Error in optim(init, gpd.lik, hessian = TRUE, method = method, control = list(maxit = maxit, : non-finite value supplied by optim
Untested in the absence of a reproducible example: http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example If I understand the model you want to give correctly, you are fitting multiple GPDs at each spatial point where data is collected over time. If that's the case, you likely should place the entire time "slice" to gpd.fit() all at once. Cheers, Michael
I've tried googling this to find a solution, but most resources i have found don't list a solution. The problem I have with not looping it over individual points is that the value in the fit should vary due to the fact that at each location, the terrain and factors driving the wind speed will vary greatly (some stuff over water, some stuff over mountains, etc.) so the plan is to do a gpd.fit for each individual point so that i can plot a return level map similar to what can be seen on page 8 of this paper (http://www.stat.duke.edu/~es112/Research/ims-MannshardtShamseldinSmithSain.pdf ). This paper takes into account spatial statistics, which is beyond the level of what I am trying to do, which is simply generating a return level at each individual point based on the data points of that individual place, not taking into account the neighbors of that point. Let me know if this doesn't make sense.