[R] modifying a built in function from the stats package (fixing arima)

2009-03-03 Thread Marc Vinyes
Dear members of the list,

I'm a beginner in R and I'm having some trouble with: "Error in
optim(init[mask], armafn, method = "BFGS", hessian = TRUE, control =
optim.control,  :
  non-finite finite-difference value [8]"

when running "arima".

I've seen that some people have come accross the same problem:
https://stat.ethz.ch/pipermail/r-help/2008-August/169660.html

So I'd like to modify the code of arima to change the optimization function
with another one that handles these problems automatically , however I don't
find the way to do it and
http://tolstoy.newcastle.edu.au/R/e6/help/09/01/2476.html points out a way
that doesn't work for me:

* If I type edit(arima) and I modify it, changes are not saved,
* If I copy the code and I save it like a different function, I get the hard
error: "Error in Delta %+% c(1, -1) : object "R_TSconv" not found"

Anybody can give me a hint? I miss matlab's easy way of doing this ("edit
function.m").

Thanks in advance

MarC (AleaSoft)

[[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.


Re: [R] scatter plot question

2009-03-04 Thread Marc Vinyes
plot(x,rho,pch=id)

-Mensaje original-
De: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]en
nombre de Dipankar Basu
Enviado el: 03 March 2009 20:11
Para: r-help@r-project.org
Asunto: [R] scatter plot question


Hi R Users,

I have a dataframe like this:

id  x   rho
A  1   0.1
B  20  0.5
C  2   0.9
...

I want to do a scatter plot of "x" versus "rho" but for each point on the
scatter plot I want the corresponding entry for "id" instead of points. In
STATA I can do so by

twoway (scatter x rho, mlabel(id))

How can I do the same in R? I am sure there is some simple way to do this.

Dipankar

[[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.


Re: [R] modifying a built in function from the stats package (fixing arima)

2009-03-04 Thread Marc Vinyes
Dear Carlos and Kjetil,

Thanks for your answer.

>I do not think that is the way to go. If you believe that your algorithm
>is better than the existing one, talk to the author of the package and
>discuss the improvement. The whole community will benefit.

I should be able to *easily* modify it and test it first!

>Copy the existing function into a new file, edit it and load it via
>source.

>3)  after downloading the source package (stats) containung arima,
> rename it (my.arima) and then do the changes.

I obviously saved it with a different name and I was expecting it to work
out of the box but I get an error that I don't know how to solve:
Error in Delta %+% c(1, rep(0, seasonal$period - 1), -1) :
  object "R_TSconv" not found

Other people have previously discussed this in this list with no success...
http://www.nabble.com/Foreign-function-call-td21836156.html

Any other hints or maybe help with the error that I'm getting?

-Mensaje original-
De: Carlos J. Gil Bellosta [mailto:c...@datanalytics.com]
Enviado el: 03 March 2009 21:30
Para: Marc Vinyes
CC: r-help@r-project.org
Asunto: Re: [R] modifying a built in function from the stats package
(fixing arima)


Hello,

I do not think that is the way to go. If you believe that your algorithm
is better than the existing one, talk to the author of the package and
discuss the improvement. The whole community will benefit.

If you want to tune the existing function and tailor it to your needs,
you have several ways to go, among them:

1) Copy the existing function into a new file, edit it and load it via
source.

2) Download the source package and modify it for your own purposes.

Best regards,

Carlos J. Gil Bellosta
http://www.datanalytics.com


On Tue, 2009-03-03 at 18:20 +0100, Marc Vinyes wrote:
> Dear members of the list,
>
> I'm a beginner in R and I'm having some trouble with: "Error in
> optim(init[mask], armafn, method = "BFGS", hessian = TRUE, control =
> optim.control,  :
>   non-finite finite-difference value [8]"
>
> when running "arima".
>
> I've seen that some people have come accross the same problem:
> https://stat.ethz.ch/pipermail/r-help/2008-August/169660.html
>
> So I'd like to modify the code of arima to change the optimization
function
> with another one that handles these problems automatically , however I
don't
> find the way to do it and
> http://tolstoy.newcastle.edu.au/R/e6/help/09/01/2476.html points out a way
> that doesn't work for me:
>
> * If I type edit(arima) and I modify it, changes are not saved,
> * If I copy the code and I save it like a different function, I get the
hard
> error: "Error in Delta %+% c(1, -1) : object "R_TSconv" not found"
>
> Anybody can give me a hint? I miss matlab's easy way of doing this ("edit
> function.m").
>
> Thanks in advance
>
> MarC (AleaSoft)
>
>   [[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.


Re: [R] modifying a built in function from the stats package (fixing arima)

2009-03-05 Thread Marc Vinyes
>If you ***look at the code*** for arima you will see that ``%+%'' is
>defined
>in terms of a call to ``.Call()'' which calls ``R_TSconv''.  So
>apparently
>R_TSconv is a C or Fortran function or subroutine in a ``shared
>object library''
>or dll upon which arima depends.  Hence to do anything with it you'll
>need to get
>that shared object library and dynamically load it.  (E.g. get the
>code, SHLIB it,
>and dynamically load the resulting shared object library.)

>The code is all available from the R source tarball.

>If this is a challenge for you then the best advice would be not to
>mess with it.

Hi Rolf,
It took me some time to come to the same conclusion (I didn't even know what
.Call() was) but I've found an easier way to modify the R file without
having to understand how to link dlls. I just downloaded the full R package,
Rtools and followed the instructions in
http://cran.r-project.org/doc/manuals/R-admin.html#Building-the-core-files
to build it. Then I can modify C:\R\src\library\stats\R\arima.R and run it.
It is quite exagerated that I have to build R in order to modify an R file
without messing with dlls, and I think it would be interesting to make this
process easier, but for now I'm happy to be productive again.

Thank you all for your help,

Best,
MarC

__
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] modifying a built in function from the stats package (fixing arima)

2009-03-05 Thread Marc Vinyes
>Just a quick note on your original question:
>if you use edit(arima), you have to remember that it returns the
>modified function, which then must be stored.

>I.e, use
>arima<-edit(arima)

>instead of just

>edit(arima)

>,and changes should be stored.

THIS IS IT.
IMHO, this should be written in BOLD LETTERS in the "Introduction to R"
manual
("edit" is only mentioned to edit DATA).

Best,
MarC

__
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] arima, xreg, and the armax model

2009-03-26 Thread Marc Vinyes
Hello all,
 
I'm having fun again with the arima function. This time I read in:
http://www.stat.pitt.edu/stoffer/tsa2/R_time_series_quick_fix.htm
 
<>
(by R.H. Shumway & D.S. Stoffer)
 
This is quite surprising... Does anybody know anything about it?
 
Marc Vinyes (AleaSoft)

[[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.


Re: [R] modifying a built in function from the stats package (fixing arima) (CONCLUSIONS)

2009-03-06 Thread MarC Vinyes Raso
Thanks a lot to everybody that helped me out with this.

Conclusions:

(1)
In order to edit arima in R:
>fix(arima)

or alternatively:
>arima<-edit(arima)

(2)
This is not contained in the "Introduction to R" manual.

(3)
A "productive" fix of arima is attached (arma coefficients printed out and
error catched so that it doesn't halt parent loops to search for candidate
coefficients):
Note 1: "productive" means I'm a beginner in R so there is probably a better
way to print the error message and fill the output arguments (I only return
NA in aic,var and sigma2).
Note 2: Changing BFGS to Nelder–Mead in "exitpoint 0" changes the
coefficients for which arima can't fit a model but results in terms of aic
and sigma2 also change significantly. By visual inspection I think that BFGS
works better.

function (x, order = c(0, 0, 0), seasonal = list(order = c(0,
0, 0), period = NA), xreg = NULL, include.mean = TRUE, transform.pars =
TRUE,
fixed = NULL, init = NULL, method = c("CSS-ML", "ML", "CSS"),
n.cond, optim.control = list(), kappa = 1e+06)
{
"%+%" <- function(a, b) .Call(R_TSconv, a, b)
upARIMA <- function(mod, phi, theta) {
p <- length(phi)
q <- length(theta)
mod$phi <- phi
mod$theta <- theta
r <- max(p, q + 1)
if (p > 0)
mod$T[1:p, 1] <- phi
if (r > 1)
mod$Pn[1:r, 1:r] <- .Call(R_getQ0, phi, theta)
else if (p > 0)
mod$Pn[1, 1] <- 1/(1 - phi^2)
else mod$Pn[1, 1] <- 1
mod$a[] <- 0
mod
}
arimaSS <- function(y, mod) {
.Call(R_ARIMA_Like, y, mod$phi, mod$theta, mod$Delta,
mod$a, mod$P, mod$Pn, as.integer(0), TRUE)
}
armafn <- function(p, trans) {
par <- coef
par[mask] <- p
trarma <- .Call(R_ARIMA_transPars, par, arma, trans)
Z <- upARIMA(mod, trarma[[1]], trarma[[2]])
if (ncxreg > 0)
x <- x - xreg %*% par[narma + (1:ncxreg)]
res <- .Call(R_ARIMA_Like, x, Z$phi, Z$theta, Z$Delta,
Z$a, Z$P, Z$Pn, as.integer(0), FALSE)
s2 <- res[1]/res[3]
0.5 * (log(s2) + res[2]/res[3])
}
armaCSS <- function(p) {
par <- as.double(fixed)
par[mask] <- p
trarma <- .Call(R_ARIMA_transPars, par, arma, FALSE)
if (ncxreg > 0)
x <- x - xreg %*% par[narma + (1:ncxreg)]
res <- .Call(R_ARIMA_CSS, x, arma, trarma[[1]], trarma[[2]],
as.integer(ncond), FALSE)
0.5 * log(res)
}
arCheck <- function(ar) {
p <- max(which(c(1, -ar) != 0)) - 1
if (!p)
return(TRUE)
all(Mod(polyroot(c(1, -ar[1:p]))) > 1)
}
maInvert <- function(ma) {
q <- length(ma)
q0 <- max(which(c(1, ma) != 0)) - 1
if (!q0)
return(ma)
roots <- polyroot(c(1, ma[1:q0]))
ind <- Mod(roots) < 1
if (all(!ind))
return(ma)
if (q0 == 1)
return(c(1/ma[1], rep(0, q - q0)))
roots[ind] <- 1/roots[ind]
x <- 1
for (r in roots) x <- c(x, 0) - c(0, x)/r
c(Re(x[-1]), rep(0, q - q0))
}
series <- deparse(substitute(x))
if (NCOL(x) > 1)
stop("only implemented for univariate time series")
method <- match.arg(method)
x <- as.ts(x)
if (!is.numeric(x))
stop("'x' must be numeric")
storage.mode(x) <- "double"
dim(x) <- NULL
n <- length(x)
if (!missing(order))
if (!is.numeric(order) || length(order) != 3 || any(order <
0))
stop("'order' must be a non-negative numeric vector of length
3")
if (!missing(seasonal))
if (is.list(seasonal)) {
if (is.null(seasonal$order))
stop("'seasonal' must be a list with component 'order'")
if (!is.numeric(seasonal$order) || length(seasonal$order) !=
3 || any(seasonal$order < 0))
stop("'seasonal$order' must be a non-negative numeric vector
of length 3")
}
else if (is.numeric(order)) {
if (length(order) == 3)
seasonal <- list(order = seasonal)
else ("'seasonal' is of the wrong length")
}
else stop("'seasonal' must be a list with component 'order'")
if (is.null(seasonal$period) || is.na(seasonal$period) ||
seasonal$period == 0)
seasonal$period <- frequency(x)
arma <- as.integer(c(order[-2], seasonal$order[-2], seasonal$period,
order[2], seasonal$order[2]))
narma <- sum(arma[1:4])
xtsp <- tsp(x)
tsp(x) <- NULL
Delta <- 1
for (i in seq_len(order[2])) Delta <- Delta %+% c(1, -1)
for (i in seq_len(seasonal$order[2])) Delta <- Delta %+%
c(1, rep(0, seasonal$period - 1), -1)
Delta <- -Delta[-1]
nd <- order[2] + seasonal$order[2]
n.used <- sum(!is.na(x)) - length(Delta)
if (is.null(xreg)) {
ncxreg <- 0
}
else {
nmxreg <- deparse(substitute(x