Peter Maclean
2011-Jun-30 04:16 UTC
[R] Saving fExtremes estimates and k-block return level with confidence intervals.
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?
Jorge Ivan Velez
2011-Jun-30 07:37 UTC
[R] Saving fExtremes estimates and k-block return level with confidence intervals.
Hi Peter, Try data.frame(n = names(res2), t(sapply(res2, function(l) l@fit$par.ests))) for the first part. HTH, Jorge On Thu, Jun 30, 2011 at 12:16 AM, Peter Maclean <> 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. >[[alternative HTML version deleted]]
Peter Maclean
2011-Jul-06 08:29 UTC
[R] Saving fExtremes estimates and k-block return level with confidence intervals.
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 at 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 <djmuser at gmail.com>
To: Peter Maclean <pmaclean2011 at 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 at 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 <pmaclean2011 at 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 at 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.
>