[R] SARIMA in R
Hi, Is R's implementation of Seasonal ARIMA in the arima() function a multiplicative or an additive model? e.g., is an ARIMA(0,1,1)(0,1,1)[12] from arima() the same as Box et al's ARIMA(0,1,1)x(0,1,1)[12] (from Time Series Analysis 1994, p.333). From another post http://tolstoy.newcastle.edu.au/R/help/04/07/0117.html I suspect it's additive but I'm not sure. Thanks, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Off topic: S.E. for cross validation
Hi, I'm performing (blocked) 10-fold cross-validation of a several time series forecasting methods, measuring their mean squared error (MSE). I know that the MSE_cv is the average over the 10 MSEs. Is there a way to calculate the standard error as well? The usual SD/sqrt(n) formula probably doesn't apply here as the 10 observations aren't independent. Thanks, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Allocating shelf space
A: Make the efficient use of space B: Minimise the spatial disclocation of related books (it is acceptable to separate large books from small books on the same subject, for the sake of efficient packing). Some comments, hope they make sense: Let f(x) be a function that maps from a specific book arrangement to a certain amount of space wastage. You're also trying to minimise some function g() of the books' location. You can't minimise two functions at once, unless you minimise some function of both: h(f(x), g(x)). It up to you to determine what h() is. For example, you could use a linear function, deciding that saving space is 10 times more important than keeping books close together. Then your objective function could be: minimise: h = f(x) + g(x) subject to: f(x) = 10g(x) f(x) = 0, g(x) = 0 (plus some nontrivial constraints on x) (You should also set a lower bound on the solution values, otherwise f will always be minimised at the expense of g, since f is worth more). Although I've stated the problem in terms of Linear Programming, it's really cheating. The much bigger issue is the combinatorial optimisation problem underneath --- different arrangements of x result in different values of h. This is much harder than LP, for anything but a small number of objects to arrange. I'd be tempted to set up a toy version, with small number of possible x values and simple constraints, and run some heuristic-driven optimisation method such as simulated annealing, Ant Colony Optimisation, Genetic Algorithms, etc. Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Not showing dvi with Hmisc latex()
Hi, I'm using latex() from Frank Harrell's Hmisc library to produce LaTeX files. By default, it calls xdvi and displays the dvi. How can I make xdvi not show? I couldn't find a clue in the extensive documentation. Thanks, Gad ps: Hmisc 3.3-1 on R 2.5.0 for Linux. -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Not showing dvi with Hmisc latex()
Frank E Harrell Jr wrote: Duncan Murdoch wrote: On 4/26/2007 9:20 PM, Gad Abraham wrote: Hi, I'm using latex() from Frank Harrell's Hmisc library to produce LaTeX files. By default, it calls xdvi and displays the dvi. How can I make xdvi not show? I couldn't find a clue in the extensive documentation. Unclass the result so it doesn't print as a latex object. For example, unclass(latex(1, file=test.tex)) $file [1] test.tex $style character(0) Alternatively, if you just assign the result you can print it later. It's when you print that the latex'ing happens. Duncan Murdoch Just say w - latex(object, file='foo.tex') So simple... thanks. Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] extracting the mode of a vector
Benoît Lété wrote: Hello, I have an elementary question (for which I couldn't find the answer on the web or the help): how can I extract the mode (modal score) of a vector? Assuming that your vector contains only integers: v - sample(1:5, size=20, replace=T) v [1] 1 1 1 1 2 3 5 1 1 5 2 4 1 3 1 1 5 4 1 5 vt - table(v) as.numeric(names(vt[vt == max(vt)])) [1] 1 Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Fractals with R
kone wrote: Hi everybody, I put some R-code to a web page for drawing fractals. See http://fractalswithr.blogspot.com/ If you have some R-code for fractal images, I'm willing to include them to the page. That's really nice, but what's with the annoying popups, courtesy of Webstats4U? Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ [EMAIL PROTECTED] 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] Fractals with R
Atte Tenkanen wrote: Hi, That is of counter for web page. Do you get some pop-up windows? Atte Hi Atte, Yes, annoying stuff from http://ilead.itrack.it. Your webcounter provider using their scripts to do this, without your intervention. How nice of them... Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ [EMAIL PROTECTED] 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] truehist bug?
Hi, Is this a bug in truehist()? library(MASS) x - rep(1, 10) truehist(x) Error in pretty(data, nbins) : invalid 'n' value Thanks, Gad R.version platform i486-pc-linux-gnu arch i486 os linux-gnu system i486, linux-gnu status major 2 minor 4.1 year 2006 month 12 day18 svn rev40228 language R version.string R version 2.4.1 (2006-12-18) sessionInfo() R version 2.4.1 (2006-12-18) i486-pc-linux-gnu locale: LC_CTYPE=en_AU.UTF-8;LC_NUMERIC=C;LC_TIME=en_AU.UTF-8; LC_COLLATE=en_AU.UTF-8;LC_MONETARY=en_AU.UTF-8; LC_MESSAGES=en_AU.UTF-8;LC_PAPER=en_AU.UTF-8;LC_NAME=C; LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_AU.UTF-8; LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods [7] base other attached packages: MASS 7.2-32 -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] ERROR: 'latex' needed but missing on your system.
Joel J. Adamson wrote:
After successfully building R on Slackware Linux v11.0 I went to make
the documentation; the texi files went fine and then I hopefully issued
make dvi
after having gotten the warning to the effect of You cannot build the
DVI or PDF manuals during compilation. And, as expected I got the
error
ERROR: 'latex' needed but missing on your system.
The problem is that latex is on my system and is in root's path:
/usr/src/R-2.4.1 Super-User echo $PATH
/usr/share/texmf/bin/:/opt/kde/bin/:/uCsr/local/stata{sic}:/usr/local/sbin:/usr/local/bin:/sbin:/usr/sbin:/bin:/usr/bin
I can issue latex from the command line as root (su'd to root, that
is) and it will run successfully. Also, whereis latex turns up
empty.
It's a bit strange, because by default files like latex should be
readable by all users. Did you install latex from source?
Try this:
As root, do 'which latex' to see where it's installed. Make sure that
the file and directories on its path are readable by your non-root user,
and that the directory is in the non-root user's path. The file 'latex'
might also be a symlink to some other file (as is in Ubuntu), so that
one will also need to be readable.
--
Gad Abraham
Department of Mathematics and Statistics
The University of Melbourne
Parkville 3010, Victoria, Australia
email: [EMAIL PROTECTED]
web: http://www.ms.unimelb.edu.au/~gabraham
__
R-help@stat.math.ethz.ch 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] expm() within the Matrix package
Gad If you convert to numeric, you can then assign it to Loglik: Loglik[1] - as.numeric(log(p %*% expm(Q * y[i]) %*% q)) Loglik[1] Gad [1] 134.5565 Hmm, I don't think that's Laura's problem (and actually I don't know what her problem is) : Assignment of a 1 x 1 matrix to a vector is not a problem, and hence the as.numeric(.) above really has no effect : ll - 1:2 (m - matrix(pi, 1,1)) [,1] [1,] 3.141593 ll[1] - m ll [1] 3.141593 2.00 Ah but you're using 'matrix' while she's using 'Matrix' (AFAIK there is no expm for matrix): library(Matrix) Loading required package: lattice m - Matrix(1, nrow=1, ncol=1) m 1 x 1 diagonal matrix of class ddiMatrix [,1] [1,]1 Warning message: Ambiguous method selection for diag, target ddiMatrix (the first of the signatures shown will be used) diagonalMatrix ddenseMatrix in: .findInheritedMethods(classes, fdef, mtable) a - vector(numeric, 0) a[1] - m Error in a[1] - m : incompatible types (from S4 to double) in subassignment type fix a[1] - as.numeric(m) a [1] 1 sessionInfo() R version 2.4.1 (2006-12-18) i486-pc-linux-gnu locale: LC_CTYPE=en_AU.UTF-8;LC_NUMERIC=C;LC_TIME=en_AU.UTF-8;LC_COLLATE=en_AU.UTF-8;LC_MONETARY=en_AU.UTF-8;LC_MESSAGES=en_AU.UTF-8;LC_PAPER=en_AU.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_AU.UTF-8;LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods [7] base other attached packages: Matrix lattice 0.9975-11 0.14-16 -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] expm() within the Matrix package
Laura Hill wrote:
Hi
Could anybody give me a bit of advice on some code I'm having trouble with?
I've been trying to calculate the loglikelihood of a function iterated over
a data of time values and I seem to be experiencing difficulty when I use
the function expm(). Here's an example of what I am trying to do
y-c(5,10) #vector of 2 survival times
p-Matrix(c(1,0),1,2) #1x2 matrix
Q-Matrix(c(1,2,3,4),2,2) #2x2 matrix
q-Matrix(c(1,2),2,1) #2x1 matrix
Loglik-rep(0,length(y))#creating empty vector of same length as y
for(i in 1:length(y)){
Loglik[i]-log((p %*% (expm(Q*y[i]))) %*% q) #calculating
# Loglikelihood for each y[i]
}
The code works perfectly if I use exp(Q*y[i]) but not for expm()
You need to do a type conversion here, because you get the loglik as a
1x1 Matrix, instead of a scalar:
(after running your code)
log(p %*% expm(Q * y[i]) %*% q)
1 x 1 Matrix of class dgeMatrix
[,1]
[1,] 134.5565
If you convert to numeric, you can then assign it to Loglik:
Loglik[1] - as.numeric(log(p %*% expm(Q * y[i]) %*% q))
Loglik[1]
[1] 134.5565
--
Gad Abraham
Department of Mathematics and Statistics
The University of Melbourne
Parkville 3010, Victoria, Australia
email: [EMAIL PROTECTED]
web: http://www.ms.unimelb.edu.au/~gabraham
__
R-help@stat.math.ethz.ch 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 standard error
Hi, Can anyone explain how the standard error in arima() is calculated? Also, how can I extract it from the Arima object? I don't see it in there. x - rnorm(1000) a - arima(x, order = c(4, 0, 0)) a Call: arima(x = x, order = c(4, 0, 0)) Coefficients: ar1 ar2 ar3 ar4 intercept -0.0451 0.0448 0.0139 -0.0688 0.0010 s.e. 0.0316 0.0316 0.0317 0.0316 0.0296 sigma^2 estimated as 0.9775: log likelihood = -1407.56, aic = 2827.12 names(a) [1] coef sigma2var.coef mask loglikaic arma residuals call seriescode n.condmodel Thanks, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] ARIMA standard error
Andrew Robinson wrote: Hi Gad, try: class(a) [1] Arima getAnywhere(print.Arima) Thanks Andrew. For the record, the standard error is the square root of the diagonal of the covariance matrix a$var.coef (itself obtained through some magic): ses[x$mask] - round(sqrt(diag(x$var.coef)), digits = digits) Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] multiplying matrix by vector of times
Thanks for the tip about the empty vector, I ever knew you could do that. I just have one problem, Lets say Q is a 2x2 matrix p is a 1x2 matrix q is a 2x1 matrix y is vector of times, say y = c(5, 10) How do I multiply Q by each time y[i]? I would like to get the answer to the equation loglik[i] - log(p %*% expm(Q * y[i]) %*% q) Where first y=5 and then y=10 so that the answers to loglik for each i are put into the empty vector. I'm sure that I am missing something fairly obvious here but can't put my finger on it. That's just what the code in my last message is doing: iterating over the y vector, multiplying Q by each y[i], and storing the result in a different cell of loglik. So you're not multiplying Q by a vector y, but Q by a scalar y[i]. Have a look at the R Introduction at http://cran.r-project.org/doc/manuals/R-intro.html, especially section 2. -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] AR(1) and gls
Kate Stark wrote: Hi there, I am using gls from the nlme library to fit an AR(1) regression model. I am wondering if (and how) I can separate the auto-correlated and random components of the residuals? Id like to be able to plot the fitted values + the autocorrelated error (i.e. phi * resid(t-1)), to compare with the observed values. I am also wondering how I might go about calculating confidence (or prediction) intervals around these new fitted values (i.e. fitted new = fitted + autocorrelated error component)? Since no one else has answered, I'll put in my 2c. Why not use the arima() function? Once you've fitted the AR(1) to it, you can extract the coefficients and the residuals, and look at the ACF and PACF of the residuals to see whether they are autocorrelated or not. tsdiag() is also useful. If an AR(1) captures all of the autocorrelation in the data, then the residuals will be iid. arima() also gives you standard errors for each estimated coefficient. To get a prediction from the model, you can use arima.sim(). One simple way to get prediction intervals about the point prediction is to sample by calling arima.sim() repeatedly and then calculate quantiles. -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Estimating parameters of 2 phase Coxian using optim
Laura Hill wrote:
On 7/3/07 00:15, Gad Abraham [EMAIL PROTECTED] wrote:
Andy Fugard wrote:
Hi There,
Perhaps the problem is the line
loglik-log(p %*% expm(Q * y(i)) %*% q)
You mention that y is a vector but here you're treating it as a
function. Maybe try
loglik-log(p %*% expm(Q * y[i]) %*% q)
?
Don't have a clue about the correctness of the contents of cox2.lik...
Andy
On 6 Mar 2007, at 08:54, Laura Hill wrote:
Hi,
My name is Laura. I'm a PhD student at Queen's University Belfast
and have
just started learning R. I was wondering if somebody could help me
to see
where I am going wrong in my code for estimating the parameters
[mu1, mu2,
lambda1] of a 2-phase Coxian Distribution.
cox2.lik-function(theta, y){
mu1-theta[1]
mu2-theta[2]
lambda1-theta[3]
p-Matrix(c(1, 0), nrow=1, ncol=2)
Q-Matrix(c(-(lambda1 + mu1), 0, lambda1, -mu2), nrow=2, ncol=2)
q-Matrix(c(mu1, mu2), nrow=2, ncol=1)
for (i in 1:length(y)){
loglik-log(p %*% expm(Q * y(i)) %*% q)
return(loglik)}
sumloglik-sum(loglik)
return(-sumloglik)
}
Just to add my 2 AU cents regarding the for loop:
You're trying to create a vector of log likelihoods to sum up later, but
that's not what's happening there. Instead, assign an empty vector of
same length as y, then assign the loglik from each iteration to a
different cell.
Lastly, there's no need to return anything from a for loop, it's not a
function.
HTH,
Gad
Hi Gad,
Yes that's exactly hat I am trying to do. If I gave you a simple example,
could you perhaps tell me how I could create a vector of log likelihoods.
Lets say I have 1x1 matrices:
p=[1]
Q=[0.05] i.e. [mu1]
q=[-0.05] i.e. [-mu1]
Where mu1 is the parameter that I would like to estimate and I have chosen
the initial value mu1=0.05
Loglik-p %*% expm(Q*y) %*% q
Where y=(5 10)
I want to sum the log likelihoods that I get for y=5 and y=10 using
Sumloglik-sum(allloglik)
Where allloglik = vector of log likelihoods
Any help would be greatly appreciated.
Thanks in advance
Laura
Hi Laura,
Make an empty vector of required length, then assign the loglik to each
of its cells, and don't return() anything:
loglik - rep(0, length(y))
for(i in 1:length(y)){
loglik[i] - log(p %*% expm(Q * y[i]) %*% q)
}
Then you can sum(loglik) like you did before.
Cheers,
Gad
--
Gad Abraham
Department of Mathematics and Statistics
The University of Melbourne
Parkville 3010, Victoria, Australia
email: [EMAIL PROTECTED]
web: http://www.ms.unimelb.edu.au/~gabraham
__
R-help@stat.math.ethz.ch 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] good procedure to estimate ARMA(p, q)?
Michael wrote: Hi all, I have some residuals from regression, and i suspect they have correlations in them... I am willing to cast the correlation into a ARMA(p, q) framework, what's the best way to identify the most suitable p, and q, and fit ARMA(p, q) model and then correct for the correlations in regression? I know there are functions in R, I have used them before, but I just want to see if I can do the whole procedure myself, just to improve my understanding ... Please give me some pointers! Thanks a lot I'm assuming the data is a time series, otherwise ARIMA models might not be applicable here. I think identifying the order of ARIMA models is something of an art, because most real world models aren't as clean and simple as textbook examples. When you have several similar models, each with its own strengths and weaknesses, which one is best? In short, you want to make sure your series is stationary, look at its ACF and PACF, then try different values of p and q based on that, and finally look at the residuals (autocorrelation, distribution, etc). This is basically the Box-Jenkins methodology. The most accessible descriptions I've seen are in Forecasting: Methods and Applications by Makridakis, Wheelwright and Hyndman (chapter 7), and Forecasting with Univariate Box-Jenkins Models by Pankratz. Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Estimating parameters of 2 phase Coxian using optim
Andy Fugard wrote:
Hi There,
Perhaps the problem is the line
loglik-log(p %*% expm(Q * y(i)) %*% q)
You mention that y is a vector but here you're treating it as a
function. Maybe try
loglik-log(p %*% expm(Q * y[i]) %*% q)
?
Don't have a clue about the correctness of the contents of cox2.lik...
Andy
On 6 Mar 2007, at 08:54, Laura Hill wrote:
Hi,
My name is Laura. I'm a PhD student at Queen's University Belfast
and have
just started learning R. I was wondering if somebody could help me
to see
where I am going wrong in my code for estimating the parameters
[mu1, mu2,
lambda1] of a 2-phase Coxian Distribution.
cox2.lik-function(theta, y){
mu1-theta[1]
mu2-theta[2]
lambda1-theta[3]
p-Matrix(c(1, 0), nrow=1, ncol=2)
Q-Matrix(c(-(lambda1 + mu1), 0, lambda1, -mu2), nrow=2, ncol=2)
q-Matrix(c(mu1, mu2), nrow=2, ncol=1)
for (i in 1:length(y)){
loglik-log(p %*% expm(Q * y(i)) %*% q)
return(loglik)}
sumloglik-sum(loglik)
return(-sumloglik)
}
Just to add my 2 AU cents regarding the for loop:
You're trying to create a vector of log likelihoods to sum up later, but
that's not what's happening there. Instead, assign an empty vector of
same length as y, then assign the loglik from each iteration to a
different cell.
Lastly, there's no need to return anything from a for loop, it's not a
function.
HTH,
Gad
--
Gad Abraham
Department of Mathematics and Statistics
The University of Melbourne
Parkville 3010, Victoria, Australia
email: [EMAIL PROTECTED]
web: http://www.ms.unimelb.edu.au/~gabraham
__
R-help@stat.math.ethz.ch 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] set.seed() and .Random.number
Taka Matzmoto wrote: Hi R-users I have two conditions. For each condition, 100 sets of 10 random numbers from N(0,1) need to be generated. Here is my question. At the begining I specify a seed number. I want to make the 100th set of the first condition and 1st set of the second conditon the same. What do I need to do ? After generating 99th set of 10 random numbers and then saving .Random.seed then using .Random.seed for generating 100th set of the first condition and 1st set of the second condtion. Is this right? Yes, that's right. .Random.seed is a vector with 626 numbers, but set.seed() only accepts a integer number. What do I need to do for saving .Random.seed and using the saved .Random.seed for setting set.seed() ? As Prof Ripley pointed to me in a post a while ago out, the help ?.Random.seed specifies that .Random.seed is the PRNG state, not the seed itself. So assign .Random.seed to a variable, and then assign the variable to it later: runif(3) [1] 0.919502 0.361622 0.847657 r - .Random.seed runif(3) [1] 0.3898361 0.4727398 0.5698874 .Random.seed - r runif(3) [1] 0.3898361 0.4727398 0.5698874 Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics The University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Initialising Mersenne-Twister with one integer
Hi, It seems to me that the Mersenne-Twister PRNG can be initialised using one integer instead of 624 integers, since inside RNG.c code there's a function defined as MT_sgenrand(Int32). How do I actually set this seed within R? I've tried: .Random.seed - c(3, 1) runif(1) Error in runif(1) : .Random.seed has wrong length In addition, is '3' actually the correct rng.kind for the Mersenne-Twister? I'm using R version 2.2.1, 2005-12-20 on Ubuntu Dapper Linux 686. Thanks, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch 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] Time series labeling with Zoo
Gabor Grothendieck wrote: When the axis labelling does not work well you will have to do it yourself like this. The plot statement is instructed not to plot the axis and then we extract into tt all the dates which are day of the month 1. Then we manually draw the axis using those. library(zoo) set.seed(1) z - zoo(runif(10), as.Date(2005-06-01) + 0:380) plot(z, xaxt = n) tt - time(z)[as.POSIXlt(time(z))$mday == 1] axis(1, tt, format(tt, %d%b%y)) Thanks Gabor, That works nicely. Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Time series labeling with Zoo
Hi, I'm using zoo because it can automatically label the months of a time series composed of daily observations. This works well for certain time series lengths, but not for others, e.g.: While: library(zoo) plot(zoo(runif(10), as.Date(2005-06-01) + 0:50)) Shows up the months and day of month, plot(zoo(runif(10), as.Date(2005-06-01) + 0:380)) results in a single label 2006 in the middle of the x axis, with no months labelled, although there's plenty of room for labels. How do I get zoo to label the months too, without having to manually work out which day was the first day of each month? I'm using zoo 1.0-3 with R 2.2.1. Thanks, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Parkville 3010, Victoria, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Multiple lag.plots per page
Prof Brian Ripley wrote: On Wed, 14 Jun 2006, Gad Abraham wrote: Prof Brian Ripley wrote: On Tue, 13 Jun 2006, Gad Abraham wrote: Prof Brian Ripley wrote: On Tue, 13 Jun 2006, Gad Abraham wrote: Hi, I'm trying to plot several lag.plots on a page, however the second plot replaces the first one (although it only takes up the upper half as it should): par(mfrow=c(2,1)) a-sin(1:100) b-cos(1:100) lag.plot(a) lag.plot(b) What's the trick to this? lag.plot itself calls par(mfrow). The trick is to get one call to do the plots you want: lag.plot(cbind(a,b)) Thanks, that works great for multiple lag.plots. Is it possible to have a lag.plot and another type of plot on the same page? The second plot() always replaces the lag.plot for me. Yes, if the other plot is second, e.g par(mfrow=c(2,1)) a-sin(1:100) lag.plot(a) par(mfg=c(2,1)) # move to second plot plot(1:10) Following from my previous questions, lag.plot doesn't recognise some of the standard plot variables, e.g. xaxt=n doesn't remove the x-axis, and setting xlab causes an error: lag.plot(sin(1:100), xlab=foo) Error in plotts(x = x, y = y, plot.type = plot.type, xy.labels = xy.labels, : formal argument xlab matched by multiple actual arguments Is this a bug or a feature? feature. Note that the help page says ...: Further arguments to 'plot.ts'. and not `graphical parameters'. Also, how can I make lag.plot behave nicely when plotted with other plots on the same page? it takes up more room than it's allocated by par(mfrow). Really you are not using it for its intended purpose, multiple plots at different lags. (Notice the plurals in the title and the description on the help page.) Why not use plot.ts directly? If you want to pursue lag.plot, try the version in R-devel which works better for single-plot displays. OK, I've experimented with plot.ts and it does what I need it to. Thanks for your help, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Updating R on an old Linux installation (was: Where is package gridBase?)
Jonathan Dushoff wrote: I am running R 2.2.1 on a University-supported linux installation based on Redhat EL3. I am sorry that it did not occur to me to mention this before; I updated R very recently, with the most recent version available for EL3 at http://cran.cnr.berkeley.edu/bin/linux/redhat/el3/. Looking at the gridBase documentation, I find that 2.2.1 is not in fact the most recent version. I have now spent hours trying to install 2.3. No available binary rpms seem to work. If I try the source code, configure crashes with a complaint that I do not have X, which is false (see below). Do I need a new version of Linux? Any help is appreciated. Jonathan -- ./configure (lots of positive stuff) checking for X... no configure: error: --with-x=yes (default) and X11 headers/libs are not available You might have X but probably not the X development libraries/headers, which are in a separate package(s). If RH EL3 still uses XFree86 and not X.org, then I think the main package is called XFree86-devel, and you'll probably need a few others, going by http://rpmfind.net//linux/RPM/redhat/8.0/i386/XFree86-devel-4.2.0-72.i386.html. Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Multiple lag.plots per page
No No wrote: Does this help? pdf(file=lag.pdf) lag.plot(a) lag.plot(b) dev.off() After that you can open each page of the lag.pdf file with GIMP for further manipulation. It gives each plot on a different page, but no plot is replaced. That would be a last resort, but obviously it's not very practical if you have more than a handful of plots. Thanks, Gad 2006/6/13, Gad Abraham [EMAIL PROTECTED] mailto:[EMAIL PROTECTED]: Hi, I'm trying to plot several lag.plots on a page, however the second plot replaces the first one (although it only takes up the upper half as it should): par(mfrow=c(2,1)) a-sin(1:100) b-cos(1:100) lag.plot(a) lag.plot(b) What's the trick to this? I'm using R 2.2.1 (2005-12-20 r36812) on Ubuntu Linux. -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Multiple lag.plots per page
Prof Brian Ripley wrote: On Tue, 13 Jun 2006, Gad Abraham wrote: Hi, I'm trying to plot several lag.plots on a page, however the second plot replaces the first one (although it only takes up the upper half as it should): par(mfrow=c(2,1)) a-sin(1:100) b-cos(1:100) lag.plot(a) lag.plot(b) What's the trick to this? lag.plot itself calls par(mfrow). The trick is to get one call to do the plots you want: lag.plot(cbind(a,b)) Thanks, that works great for multiple lag.plots. Is it possible to have a lag.plot and another type of plot on the same page? The second plot() always replaces the lag.plot for me. Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Multiple lag.plots per page
Prof Brian Ripley wrote: On Tue, 13 Jun 2006, Gad Abraham wrote: Prof Brian Ripley wrote: On Tue, 13 Jun 2006, Gad Abraham wrote: Hi, I'm trying to plot several lag.plots on a page, however the second plot replaces the first one (although it only takes up the upper half as it should): par(mfrow=c(2,1)) a-sin(1:100) b-cos(1:100) lag.plot(a) lag.plot(b) What's the trick to this? lag.plot itself calls par(mfrow). The trick is to get one call to do the plots you want: lag.plot(cbind(a,b)) Thanks, that works great for multiple lag.plots. Is it possible to have a lag.plot and another type of plot on the same page? The second plot() always replaces the lag.plot for me. Yes, if the other plot is second, e.g par(mfrow=c(2,1)) a-sin(1:100) lag.plot(a) par(mfg=c(2,1)) # move to second plot plot(1:10) Thanks, works great. Cheers, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Multiple lag.plots per page
Prof Brian Ripley wrote: On Tue, 13 Jun 2006, Gad Abraham wrote: Prof Brian Ripley wrote: On Tue, 13 Jun 2006, Gad Abraham wrote: Hi, I'm trying to plot several lag.plots on a page, however the second plot replaces the first one (although it only takes up the upper half as it should): par(mfrow=c(2,1)) a-sin(1:100) b-cos(1:100) lag.plot(a) lag.plot(b) What's the trick to this? lag.plot itself calls par(mfrow). The trick is to get one call to do the plots you want: lag.plot(cbind(a,b)) Thanks, that works great for multiple lag.plots. Is it possible to have a lag.plot and another type of plot on the same page? The second plot() always replaces the lag.plot for me. Yes, if the other plot is second, e.g par(mfrow=c(2,1)) a-sin(1:100) lag.plot(a) par(mfg=c(2,1)) # move to second plot plot(1:10) Following from my previous questions, lag.plot doesn't recognise some of the standard plot variables, e.g. xaxt=n doesn't remove the x-axis, and setting xlab causes an error: lag.plot(sin(1:100), xlab=foo) Error in plotts(x = x, y = y, plot.type = plot.type, xy.labels = xy.labels, : formal argument xlab matched by multiple actual arguments Is this a bug or a feature? Also, how can I make lag.plot behave nicely when plotted with other plots on the same page? it takes up more room than it's allocated by par(mfrow). Thanks for your help, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Multiple lag.plots per page
Hi, I'm trying to plot several lag.plots on a page, however the second plot replaces the first one (although it only takes up the upper half as it should): par(mfrow=c(2,1)) a-sin(1:100) b-cos(1:100) lag.plot(a) lag.plot(b) What's the trick to this? I'm using R 2.2.1 (2005-12-20 r36812) on Ubuntu Linux. Thanks, Gad -- Gad Abraham Department of Mathematics and Statistics University of Melbourne Victoria 3010, Australia email: [EMAIL PROTECTED] web: http://www.ms.unimelb.edu.au/~gabraham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
