[R] Nonparametric Estimation of Tail Dependence Coefficients
Dear all,
in which package can I find an implementation of the nonparametric
estimation of tail dependence coefficients:
lambda_L = lim_{u\to 0} P[F_1(X_1)u|F_2(X_2)>u],
where (X_1,X_2) has marginal distribution functions F_1 and F_2?
(The nonparametric estimators typically require to choose a tuning
parameter, often called 'k', which can be automated by using a 'plateau
finding' algorithm, which I would like to be implemented as well).
Thank you very much for your help,
Gildas
__
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] Matrix package + mfcol/mfrow
Hi, When I load the "Matrix" package, I cannot get the par(mfrow=c(..,..)) to work, that is, I cannot get to display several images at a time. How can I fix this problem ? Thanks in advance, Gildas Mazo __ 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] Matrix package + mfcol/mfrow
Hi, When I load the "Matrix" package, I cannot get the par(mfrow=c(..,..)) to work, that is, I cannot get to display several images at a time. How can I fix this problem ? Thanks in advance, Gildas Mazo __ 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] function which saves an image of a dgtMatrix as png
Hi,
I'm getting crazy:
This does work:
library(Matrix)
a1<-b1<-c(1,2)
c1<-rnorm(2)
aDgt<-spMatrix(ncol=3,nrow=3,i=a1,j=b1,x=c1)
png("myImage.png")
image(aDgt)
dev.off()
But this doesn't !!!
f<-function(x){
png("myImage.png")
image(x)
dev.off()
}
f(aDgt)
My image is saved as a text file and contains nothing at all !!!
Thanks in advance,
Gildas Mazo
__
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.
Re: [R] function which saves an image of a dgtMatrix as png
Thanks so much
Douglas Bates a écrit :
> image applied to a sparseMatrix object uses lattice functions to
> create the image. As described in R FAQ 7.22 you must use
>
> print(image(x))
>
> or
>
> show(image(x))
>
> or even
>
> plot(image(x))
>
> when a lattice function is called from within another function.
> On Wed, Apr 28, 2010 at 1:20 PM, Gildas Mazo wrote:
>
>> Hi,
>>
>> I'm getting crazy:
>>
>> This does work:
>>
>> library(Matrix)
>> a1<-b1<-c(1,2)
>> c1<-rnorm(2)
>> aDgt<-spMatrix(ncol=3,nrow=3,i=a1,j=b1,x=c1)
>> png("myImage.png")
>> image(aDgt)
>> dev.off()
>>
>> But this doesn't !!!
>>
>> f<-function(x){
>> png("myImage.png")
>> image(x)
>> dev.off()
>> }
>> f(aDgt)
>>
>> My image is saved as a text file and contains nothing at all !!!
>> Thanks in advance,
>>
>> Gildas Mazo
>>
>> __
>> 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.
[R] optim and "the function should not" advice
Dear R users,
I'm using optim to optimize a pretty complicated function. This function takes
the parameter vector "theta" and within its body I use instructions like
sigma<-theta[a:b]; computations with sigma...
out<-c()
for (i in 1:d){
a<-theta[(3*d+i):c]
out[i]<-evaluation of an expression involving 'a' (I use symbolic
differentiation)
}
Unfortunately for certain problems 'optim' returns a parameter vector which
didn't move at all from the initial parameters, and the output says that
although the function has been evaluated a high number of times, the gradient
(which I fed the function with) has been evaluated only one time. I used the
BFGS method.
By chance I looked at the help and I read "The parameter vector passed to fn
has special semantics and may be shared between calls: the function should not
change or copy it" . Could the instructions above be the cause of the failure?
If so, how to deal with symbolic differentation?
Thanks in advance,
Gildas
[[alternative HTML version deleted]]
__
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] CRAN policies - citing others
Dear R users,
we try to submit a package to CRAN. One of our function pbcOptim.R uses
lbfgsb.cpp which is a slight modification of optim.c. In the documentation for
pbcOptim.R we said
"The code for pbcOptim is based on that of \code{\link{optim}}. In particular,
pbcOptim calls /src/lbfgsb.cpp which is a slight adaptation of lbfgsb.c and
optim.c part of the R software (http://www.r-project.org/). See
\code{\link{optim}} and ?optim."
However, CRAN responded that we were not the authors of lbfgsb.cpp and we did
not mention them. Do we have to spell their names explicitely? In the
documentation for pbcOptim.R? In the file DESCRIPTION? In the file lbfgsb.cpp?
Thanks you very much for your help
--
Gildas Mazo
PhD student
MISTIS team at INRIA
Grenoble, France
__
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] SQL/R
Dear R users,
I want to aggregate data in the following way:
###
X <- data.frame(u = c("T1","T1","T1","T2"), v=c("a","a","b","a"))
X
library(sqldf)
sqlOut <- sqldf("select count(distinct(v)) from X group by u")
sqlOut
###
Now I want to get the same result without using SQL. How can I achieve
that ?
Thanks for your help,
Gildas
__
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.
Re: [R] SQL/R
Thanks for your answers,
Best,
Gildas
Brian Diggs a écrit :
> On 7/22/2010 5:01 AM, Allan Engelhardt wrote:
>> There are so many ways Here is one:
>>
>> aggregate(v ~ u, data=X, function(...) length(unique(...)))
>> # u v
>> # 1 T1 2
>> # 2 T2 1
>>
>> Hope this helps
>
> Here is one other way, using the plyr package (which is very good for
> taking a data structure (data.frame, list, array), pulling it apart by
> some criteria, doing something on each of the parts, and putting the
> results back together):
>
> library("plyr")
> ddply(X, .(u), function(x) {length(unique(x$v))})
> # u V1
> #1 T1 2
> #2 T2 1
>
>
>> Allan.
>>
>> On 22/07/10 12:52, Gildas Mazo wrote:
>>> Dear R users,
>>>
>>> I want to aggregate data in the following way:
>>>
>>> ###
>>>
>>> X<- data.frame(u = c("T1","T1","T1","T2"), v=c("a","a","b","a"))
>>> X
>>> library(sqldf)
>>> sqlOut<- sqldf("select count(distinct(v)) from X group by u")
>>> sqlOut
>>>
>>> ###
>>>
>>> Now I want to get the same result without using SQL. How can I achieve
>>> that ?
>>>
>>> Thanks for your help,
>>>
>>> Gildas
>>>
>>> __
>>> 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.
[R] optimization subject to constraints
Dear R users, I'm looking for tools to perform optimization subject to constraints, both linear and non-linear. I don't mind which algorithm may be used, my primary aim is to get something general and easy-to-use to study simples examples. Thanks for helping, Gildas __ 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.
Re: [R] optimization subject to constraints
Thanks, but I still cannot get to solve my problem: consider this simple
example:
f <- function(x){
return(x[1]+x[2])
} # objective function
g <- function(x){
return(x[1]^2+x[2]^2)
} # constraint
#
I wanna Maximize f(x) subject to g(x) = 1. By hand the solution is
(1/sqrt(2), 1/sqrt(2), sqrt(2)). This is to maximizing a linear function
subject to a nonlinear equality constraint. I didn't find any suitable
function in the packages I went through.
Thanks in advance,
Gildas
Spencer Graves a écrit :
>To find every help page containing the term "constrained
> optimization", you can try the following:
>
>
> library(sos)
> co <- findFn('constrained optimization')
>
>
> "Printing" this "co" object opens a table in a web browser with
> all matches sorted first by the package with the most matches and with
> hot links to the documentation page.
>
>
> writeFindFn2xls(co)
>
>
> This writes an excel file, with the browser table as the second
> tab and the first being a summary of the packages. This summary table
> can be made more complete and useful using the "installPackages"
> function, as noted in the "sos" vignette.
>
>
> A shameless plug from the lead author of the "sos" package.
> Spencer Graves
>
>
> On 8/9/2010 10:01 AM, Ravi Varadhan wrote:
>> constrOptim can only handle linear inequality constraints. It cannot
>> handle
>> equality (linear or nonlinear) as well as nonlinear inequality
>> constraints.
>>
>> Ravi.
>>
>> -Original Message-
>> From: r-help-boun...@r-project.org
>> [mailto:r-help-boun...@r-project.org] On
>> Behalf Of Dwayne Blind
>> Sent: Monday, August 09, 2010 12:56 PM
>> To: Gildas Mazo
>> Cc: r-help@r-project.org
>> Subject: Re: [R] optimization subject to constraints
>>
>> Hi !
>>
>> Why not constrOptim ?
>>
>> Dwayne
>>
>> 2010/8/9 Gildas Mazo
>>
>>> Dear R users,
>>>
>>> I'm looking for tools to perform optimization subject to constraints,
>>> both linear and non-linear. I don't mind which algorithm may be
>>> used, my
>>> primary aim is to get something general and easy-to-use to study
>>> simples
>>> examples.
>>>
>>> Thanks for helping,
>>>
>>> Gildas
>>>
>>> __
>>> 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<http://www.r-project.org/posting
>>
>> -guide.html>
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>> [[alternative HTML version deleted]]
>>
>> __
>> 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.
>
__
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] [Fwd: Re: optimization subject to constraints]
--- Begin Message ---
Danke schön Matthias.
I had naively started with x0 = c(0,0) and I got a "Redundant
constraints were found" error. What's the problem with (0,0) ?
Matthias Gondan a écrit :
> try this (package Rsolnp)
>
> library(Rsolnp)
>
> g<- function(x)
> {
> return(x[1]^2+x[2]^2)
> } # constraint
>
> f<- function(x)
> {
> return(x[1]+x[2])
> } # objective function
>
> x0 = c(1, 1)
>
> solnp(x0, fun=f, eqfun=g, eqB=c(1))
>
>
>
> Am 10.08.2010 14:59, schrieb Gildas Mazo:
>> Thanks, but I still cannot get to solve my problem: consider this simple
>> example:
>>
>>
>>
>> f<- function(x){
>>return(x[1]+x[2])
>> } # objective function
>>
>> g<- function(x){
>>return(x[1]^2+x[2]^2)
>> } # constraint
>>
>> #
>>
>> I wanna Maximize f(x) subject to g(x) = 1. By hand the solution is
>> (1/sqrt(2), 1/sqrt(2), sqrt(2)). This is to maximizing a linear function
>> subject to a nonlinear equality constraint. I didn't find any suitable
>> function in the packages I went through.
>>
>> Thanks in advance,
>>
>> Gildas
>>
>>
>>
>>
>>
>> Spencer Graves a écrit :
>>> To find every help page containing the term "constrained
>>> optimization", you can try the following:
>>>
>>>
>>> library(sos)
>>> co<- findFn('constrained optimization')
>>>
>>>
>>>"Printing" this "co" object opens a table in a web browser with
>>> all matches sorted first by the package with the most matches and with
>>> hot links to the documentation page.
>>>
>>>
>>> writeFindFn2xls(co)
>>>
>>>
>>>This writes an excel file, with the browser table as the second
>>> tab and the first being a summary of the packages. This summary table
>>> can be made more complete and useful using the "installPackages"
>>> function, as noted in the "sos" vignette.
>>>
>>>
>>>A shameless plug from the lead author of the "sos" package.
>>>Spencer Graves
>>>
>>>
>>> On 8/9/2010 10:01 AM, Ravi Varadhan wrote:
>>>> constrOptim can only handle linear inequality constraints. It cannot
>>>> handle
>>>> equality (linear or nonlinear) as well as nonlinear inequality
>>>> constraints.
>>>>
>>>> Ravi.
>>>>
>>>> -Original Message-
>>>> From: r-help-boun...@r-project.org
>>>> [mailto:r-help-boun...@r-project.org] On
>>>> Behalf Of Dwayne Blind
>>>> Sent: Monday, August 09, 2010 12:56 PM
>>>> To: Gildas Mazo
>>>> Cc: r-help@r-project.org
>>>> Subject: Re: [R] optimization subject to constraints
>>>>
>>>> Hi !
>>>>
>>>> Why not constrOptim ?
>>>>
>>>> Dwayne
>>>>
>>>> 2010/8/9 Gildas Mazo
>>>>
>>>>> Dear R users,
>>>>>
>>>>> I'm looking for tools to perform optimization subject to constraints,
>>>>> both linear and non-linear. I don't mind which algorithm may be
>>>>> used, my
>>>>> primary aim is to get something general and easy-to-use to study
>>>>> simples
>>>>> examples.
>>>>>
>>>>> Thanks for helping,
>>>>>
>>>>> Gildas
>>>>>
>>>>> __
>>>>> 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<http://www.r-project.org/posting
>>>>
>>>>
>>>> -guide.html>
>>>>> and provide commented, minimal, self-contained, reproducible code.
>>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> __
>>>> 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.
>>
>
> __
> 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.
>
>
--- End Message ---
__
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.
Re: [R] [Fwd: Re: optimization subject to constraints]
Thanks Ravi.
Gildas
Ravi Varadhan a écrit :
> I think the problem is because the the Hessian of the augmented Lagrangian
> iis singular at c(0,0).
>
> Try this:
>
> require(alabama)
>
> heq <- function(x) {
> x[1]^2+x[2]^2 - 1
> }
>
>
>> constrOptim.nl(par=c(0,0), fn=f, heq=heq, control.outer=list(trace=FALSE))
>>
> $par
> [1] -0.7071067 -0.7071067
>
> $value
> [1] -1.414213
>
> $iterations
> [1] 10
>
> $lambda
> [1] -0.7068717
>
> $penalty
> [1] -6.496021e-08
>
> $counts
> function gradient
> 100 30
>
>
> Ravi.
>
>
>
> Ravi Varadhan, Ph.D.
> Assistant Professor,
> Division of Geriatric Medicine and Gerontology
> School of Medicine
> Johns Hopkins University
>
> Ph. (410) 502-2619
> email: rvarad...@jhmi.edu
>
>
> - Original Message -
> From: Gildas Mazo
> Date: Tuesday, August 10, 2010 10:11 am
> Subject: [R] [Fwd: Re: optimization subject to constraints]
> To: r-help@r-project.org
>
>
>
>> - Original Message -
>>
>
>
>> From Gildas Mazo
>>
>
>
>> Date Tue, 10 Aug 2010 15:49:19 +0200
>>
>
>
>> To Matthias Gondan
>>
> Subject Re: [R] optimization subject to constraints
>
>> Danke schön Matthias.
>>
>> I had naively started with x0 = c(0,0) and I got a "Redundant
>> constraints were found" error. What's the problem with (0,0) ?
>>
>>
>>
>>
>>
>>
>>
>> Matthias Gondan a écrit :
>> > try this (package Rsolnp)
>> >
>> > library(Rsolnp)
>> >
>> > g<- function(x)
>> > {
>> > return(x[1]^2+x[2]^2)
>> > } # constraint
>> >
>> > f<- function(x)
>> > {
>> > return(x[1]+x[2])
>> > } # objective function
>> >
>> > x0 = c(1, 1)
>> >
>> > solnp(x0, fun=f, eqfun=g, eqB=c(1))
>> >
>> >
>> >
>> > Am 10.08.2010 14:59, schrieb Gildas Mazo:
>> >> Thanks, but I still cannot get to solve my problem: consider this
>> simple
>> >> example:
>> >>
>> >>
>> >>
>> >> f<- function(x){
>> >>return(x[1]+x[2])
>> >> } # objective function
>> >>
>> >> g<- function(x){
>> >>return(x[1]^2+x[2]^2)
>> >> } # constraint
>> >>
>> >> #
>> >>
>> >> I wanna Maximize f(x) subject to g(x) = 1. By hand the solution is
>> >> (1/sqrt(2), 1/sqrt(2), sqrt(2)). This is to maximizing a linear function
>> >> subject to a nonlinear equality constraint. I didn't find any suitable
>> >> function in the packages I went through.
>> >>
>> >> Thanks in advance,
>> >>
>> >> Gildas
>> >>
>> >>
>> >>
>> >>
>> >>
>> >> Spencer Graves a écrit :
>> >>> To find every help page containing the term "constrained
>> >>> optimization", you can try the following:
>> >>>
>> >>>
>> >>> library(sos)
>> >>> co<- findFn('constrained optimization')
>> >>>
>> >>>
>> >>>"Printing" this "co" object opens a table in a web browser
>> with
>> >>> all matches sorted first by the package with the most matches and
>> with
>> >>> hot links to the documentation page.
>> >>>
>> >>>
>> >>> writeFindFn2xls(co)
>> >>>
>> >>>
>> >>>This writes an excel file, with the browser table as the second
>> >>> tab and the first being a summary of the packages. This summary
>> table
>> >>> can be made more complete and useful using the "installPackages"
>> >>> function, as noted in the "sos" vignette.
>> >>>
>> >>>
>> >>>A shameless plug from the lead author of the "sos" package.
>> >>>Spencer Graves
>> >>>
>> >>>
>> >>> On 8/9/2010 10:01 AM, Ravi Varadhan wrote:
>>
[R] Sweave & png
Hi, Is there a simple way to save my figures in png instead of pdf with Sweave ?? Thanks in advance, Gidas __ 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.
Re: [R] Sweave & png
Thank you Ben Bolker a écrit : > Gildas Mazo curie.fr> writes: > > >> Is there a simple way to save my figures in png instead of pdf with >> Sweave ?? >> > > > See > > http://sites.google.com/site/thibautjombart/r-packages > > (scroll to the bottom) > > __ > 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.
[R] Build Design Matrix with avoiding loops
Dear R users,
I'd like to build a simple design matrix in a efficient way. I naively
wrote the code below.
n = 15
k = 3
nbPerGrp = c(5,5,5)
xT <- list()
for (i in 1:k){
xT[[i]] <- rep(0, k)
xT[[i]][i] <- 1
}
X <- matrix(nrow = n, ncol = k) #design matrix
for (i in 1:nbPerGrp[1]){
X[i,] <- xT[[1]]
}
for (i in 1:k-1){
for (j in nbPerGrp[i]+1:nbPerGrp[i+1]){
X[j,] <- xT[[i]]
}}
for (i in 1:nbPerGrp[k]){
X[n - nbPerGrp[k] + i, ] <- xT[[k]]
}
X # That's I wanna get.
But as soon as n, k increase it takes too much time because of the loops.
Then my question is how can I get such a design matrix X without too
much loops ? Which function I should look at ?
Thanks in advance for responding me,
Gildas
__
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] message-passing algorithm
Dear R users, I'm searching a message-passing algorithm in R. Any help much appreciated. Thank you, Gildas __ 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] Numerical Inversion of Cumulative Distribution Function
Dear R users,
Given a user-defined cumulative distribution function F, I want to compute
F^{-1}(x). How is that possible with R?
Best Regards,
--
Gildas Mazo
PhD student
MISTIS team at INRIA
Grenoble, France
[[alternative HTML version deleted]]
__
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] Computing High Order Derivatives (Numerically)
Dear R users,
Let f be a function over d variables x1,..,xd. I want to compute the k^th-order
derivative with respect to x1,..,xk (k<=d). I have a by hand solution (see
below) using an iterating code using D. However, I expect d to be high and f to
be complicated. Then I want a vector x to be the input, instead of x1,..,xd.
How to avoid the x1 <- x[1]; x2 <- x[2], etc steps in the code below? Moreover,
D uses symbolic differentation and then eval evaluates the output to get a
numerical result. But is there a way to compute the desired derivatives
numerically directly (without using symbolic calculus at all)? Finally, what is
the most efficient and fast way to get a numerical result for such derivatives?
Thank you very much in advance,
Gildas
### Code ###
### dif takes a function f, an order k, and a vector x as input. f must be a
function of x1,..,xd with d >= k. The correspondance is done between xi and
x[i]. The expression for f must be at the last row of the body function.
dif <- function(f,k,x){
o <- list()
n <- length(body(f))
o[[1]] <- body(f)[[n]]
for (i in 1:k){
xi <- paste("x",i,sep="")
o[[i+1]] <- D(o[[i]],name=xi)
}
x1 <- x[1]
x2 <- x[2]
x3 <- x[3]
eval(o[[k+1]])
}
### Examples ###
## function to differentiate
f <- function(x){
x1 <- x[1]
x2 <- x[2]
x3 <- x[3]
0.5*x1*x2*x3^2
}
## derivative w.r.t. x1, x2 and x3 at the point (1,2,3).
dif(f,3,c(1,2,3))
### My Questions ###
## how to avoid to write by hand xi <- x[i] ??
## is there a way in R to compute such derivatives without using symbolic
calculation but numerical compuation instead.
__
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.
Re: [R] Computing High Order Derivatives (Numerically)
Dear Petr Savicky,
this helped indeed. Thank you very much.
Gildas
- Mail original -
> De: "Petr Savicky"
> À: r-help@r-project.org
> Envoyé: Vendredi 23 Mars 2012 09:39:37
> Objet: Re: [R] Computing High Order Derivatives (Numerically)
> On Fri, Mar 23, 2012 at 12:35:57AM +0100, Gildas Mazo wrote:
> > Dear R users,
> >
> > Let f be a function over d variables x1,..,xd. I want to compute the
> > k^th-order derivative with respect to x1,..,xk (k<=d). I have a by
> > hand solution (see below) using an iterating code using D. However,
> > I expect d to be high and f to be complicated. Then I want a vector
> > x to be the input, instead of x1,..,xd. How to avoid the x1 <- x[1];
> > x2 <- x[2], etc steps in the code below? Moreover, D uses symbolic
> > differentation and then eval evaluates the output to get a numerical
> > result. But is there a way to compute the desired derivatives
> > numerically directly (without using symbolic calculus at all)?
> > Finally, what is the most efficient and fast way to get a numerical
> > result for such derivatives?
> >
> > Thank you very much in advance,
> > Gildas
> >
> > ### Code ###
> > ### dif takes a function f, an order k, and a vector x as input. f
> > must be a function of x1,..,xd with d >= k. The correspondance is
> > done between xi and x[i]. The expression for f must be at the last
> > row of the body function.
> > dif <- function(f,k,x){
> > o <- list()
> > n <- length(body(f))
> > o[[1]] <- body(f)[[n]]
> > for (i in 1:k){
> > xi <- paste("x",i,sep="")
> > o[[i+1]] <- D(o[[i]],name=xi)
> > }
> > x1 <- x[1]
> > x2 <- x[2]
> > x3 <- x[3]
> > eval(o[[k+1]])
> > }
> >
> > ### Examples ###
> > ## function to differentiate
> > f <- function(x){
> > x1 <- x[1]
> > x2 <- x[2]
> > x3 <- x[3]
> > 0.5*x1*x2*x3^2
> > }
> > ## derivative w.r.t. x1, x2 and x3 at the point (1,2,3).
> > dif(f,3,c(1,2,3))
> >
> > ### My Questions ###
> > ## how to avoid to write by hand xi <- x[i] ??
> > ## is there a way in R to compute such derivatives without using
> > symbolic calculation but numerical compuation instead.
>
> Hi.
>
> For the first question, try the following
>
> dif <- function(f,k,x){
> o <- list()
> n <- length(body(f))
> o[[1]] <- body(f)[[n]]
> for (i in 1:k){
> xi <- paste("x",i,sep="")
> o[[i+1]] <- D(o[[i]],name=xi)
> assign(xi, x[i])
> }
> eval(o[[k+1]])
> }
>
> For the second question, try the following.
>
> x <- c(1, 2, 3)
> k <- length(x)
> grid <- as.matrix(expand.grid(rep(list(c(0, 1)), times=k)))
> signs <- 1 - 2*(rowSums(1 - grid) %% 2)
> for (eps in 2^-(5:20)) {
> xeps <- eps*grid + rep(x, each=nrow(grid))
> print(sum(signs*apply(xeps, 1, FUN=f))/eps^k)
> }
>
> [1] 3.015625
> [1] 3.007812
> [1] 3.003906
> [1] 3.001953
> [1] 3.000977
> [1] 3.000488
> [1] 3.000244
> [1] 3.000122
> [1] 3
> [1] 3
> [1] 3
> [1] 3
> [1] 4
> [1] 0
> [1] 0
> [1] 0
>
> If the above is computed in an exact arithmetic, then
> with "eps" converging to zero, the result converges to
> the required derivative. Since the numerical computations
> are done with a rounding error, too small "eps" yields
> a completely wrong result. The choice of a good "eps"
> depends on the function and on "k". For a high "k", there
> may even be no good "eps". See the considerations at
>
> http://en.wikipedia.org/wiki/Numerical_derivative
>
> where the choice of "eps" is discussed in the simplest
> case of a univariate function.
>
> Hope this helps.
>
> Petr Savicky.
>
> __
> 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.
--
Gildas Mazo
PhD student
MISTIS team at INRIA
Grenoble, France
__
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.
