Hi:
I am trying to compare the results of "evir" and "fExtreme" packages. I could
not figure out how to save the "evir" package results. Also, how to pass the
results to "fExtreme" function "gevrlevelPlot" and "evir" function "rlevel.gev"
to get the return levels. I just need the values and not graphs.
#Example
library(fExtremes)
library(evir)
require(plyr)
y<- data.frame(rgev(300, xi = 0.1, mu = .5, sigma = 1.6))
colnames(y) <- c("y")
n <- data.frame(n = rep(c(1,2,3), each = 100))
z <- cbind(n,y)
y <- z$y
# Model for grouped data
##fExtremes package
z1 <- split(z,z$n)
rgf <- lapply(z1, function(x){
m <- as.numeric(x$y)
gevFit(m, block = 2, type = "mle")
})
#Save results
resf<- ldply(rgf, function(x) x@fit$par.ests)
#Qs: How to transfer rge object to "gevrlevelPlot" function to get the values?
##evir package
rge<- by(z,z[,"n"], function(x) gev(x$y, 2, method = "BFGS", control
=list(maxit
= 10000)))
#Qs:How to save par.ests and par.ses?
#Qs:How to transfer rge object to "rlevel.gev" function to get the
values?
#Model for single vector
rlf <- gevFit(y, block = 2, type = "mle")
rlfp <- gevrlevelPlot(rlf, kBlocks = 2)
rlfp
rle <- gev(y, 2, method = "BFGS", control =list(maxit = 10000))
rlep<- rlevel.gev(rle, k.blocks = 2, add = FALSE)
rlep
----- Original Message ----
From: Dennis Murphy <djmu...@gmail.com>
To: Peter Maclean <pmaclean2...@yahoo.com>
Sent: Thu, June 30, 2011 3:03:05 PM
Subject: Re: [R] Saving fExtremes estimates and k-block return level with
confidence intervals.
Hi:
The plyr package can help you out as well. As far as the estimates go,
library(plyr)
> ldply(res2, function(x) x@fit$par.ests)
.id xi mu beta
1 1 0.1033614 2.5389580 0.9092611
2 2 0.3401922 0.5192882 1.5290615
3 3 0.5130798 0.5668308 1.2105666
You could also look into the l_ply() function, which takes a list as
input and outputs nothing; however, it is usually called if one want
to make a series of similar plots. You need to write a function that
takes a generic list component and outputs a plot. Here's a fairly
trivial example:
testdf <- data.frame(gp = rep(LETTERS[1:3], each = 10), x = 1:10,
y = rnorm(30))
l_ply(split(testdf, testdf$gp), function(df) plot(y ~ x, data = df))
Alternatively,
# Observe that a data frame is used as input
pfun <- function(df) plot(y ~ x, data = df)
l_ply(split(testdf, testdf$gp), pfun)
In this case, l_ply() plots y vs. x in each of the three subgroups of
testdf separately. It appears you want to do a similar thing, but with
the list of model outputs instead; in your case, the function you
write should take a list as its sole argument. Any slots or components
that are accessed are then relative to the input list object - see the
anonymous function I wrote for the output data frame above as an
illustration.
Since it appears that gevrlevelPlot() also outputs a vector, you may
want to write a function that returns the output vector as a data
frame and call ldply() instead. I don't know enough about the
fExtremes package to write this myself, but it should be possible to
write a function for a generic model that generates both a plot and an
output data frame. The plyr package is pretty good at this sort of
thing.
HTH,
Dennis
On Wed, Jun 29, 2011 at 9:16 PM, Peter Maclean <pmaclean2...@yahoo.com> wrote:
> I am estimating a large model by groups. How do you save the results
>and returns
> the associated quantiles?
> For this example I need a data frame
> n xi mu beta
> 1 0.1033614 2.5389580 0.9092611
> 2 0.3401922 0.5192882 1.5290615
> 3 0.5130798 0.5668308 1.2105666
> I also want to apply gevrlevelPlot() for each "n" or group.
>
> #Example
> n <- c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,3)
> y <- c(2,3,2,3,4,5,6,1,0,0,0,6, 2, 1, 0, 0,9,3)
> z <- as.data.frame(cbind(n,y))
> colnames(z) <- c("n","y")
> library(fExtremes)
> z <- split(z, z$n)
> res2 <-lapply(z, function(x){
> m <- as.numeric(x$y)
> gevFit(m, block = 1, type = c("pwm"))
> })
>> res2
> $`1`
> Title:
> GEV Parameter Estimation
> Call:
> gevFit(x = m, block = 1, type = c("pwm"))
> Estimation Type:
> gev pwm
> Estimated Parameters:
> xi mu beta
> 0.1033614 2.5389580 0.9092611
> Description
> Wed Jun 29 23:07:48 2011
>
> $`2`
> Title:
> GEV Parameter Estimation
> Call:
> gevFit(x = m, block = 1, type = c("pwm"))
> Estimation Type:
> gev pwm
> Estimated Parameters:
> xi mu beta
> 0.3401922 0.5192882 1.5290615
> Description
> Wed Jun 29 23:07:48 2011
>
> $`3`
> Title:
> GEV Parameter Estimation
> Call:
> gevFit(x = m, block = 1, type = c("pwm"))
> Estimation Type:
> gev pwm
> Estimated Parameters:
> xi mu beta
> 0.5130798 0.5668308 1.2105666
> Description
> Wed Jun 29 23:07:48 2011
>
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
