Re: [R] The quadprog package
G'day Thomas,
On Mon, 3 Sep 2007 10:08:08 +0200
[EMAIL PROTECTED] wrote:
What's wrong with my code?
Require(quadprog)
library(quadprog) ? :)
Dmat-diag(1,7,7)
# muss als quadratische Matrix eingegeben werden
Dmat
dvec-matrix(0,7,1) # muss als Spaltenvektor eingegeben werden
dvec
mu-0 # (in Mio. €)
bvec-c(1,mu,matrix(0,7,1)) # muss als Spaltenvektor eingegeben werden
bvec
mu_r-c(19.7,33.0,0.0,49.7, 82.5, 39.0,11.8)
Amat-matrix(c(matrix(1,1,7),7*mu_r,diag(1,7,7)),9,7,byrow=T)
# muss als Matrix angegeben werden, wie sie wirklich ist
Amat
meq-2
loesung-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
With the commands above, I get on my system:
loesung-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
Error in solve.QP(Dmat, dvec, Amat = t(Amat), bvec = bvec, meq = 2) :
constraints are inconsistent, no solution!
Which is a bit strange, since Dmat is the identity matrix, so there
should be little room for numerical problems.
OTOH, it is known that the Goldfarb-Idnani algorithm can have problem
in rare occasions if the problem is ill-scaled; see the work of
Powell.
Changing the definition of Amat to
Amat-matrix(c(matrix(1,1,7),mu_r,diag(1,7,7)),9,7,byrow=T)
leads to a successful call to solve.QP. And since mu=0, I really
wonder why you scaled up the vector mu_r by a factor of 7 when putting
it into Amat :)
loesung-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
Now, with the solution I obtained the rest of your commands show:
loesung$solution %*% mu_r
[,1]
[1,] 6.068172e-15
sum(loesung$solution)
[1] 1
for (i in 1:7){
a-loesung$solution[i]=0
print(a)
}
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] FALSE
[1] FALSE
[1] TRUE
But:
for (i in 1:7){
a-loesung$solution[i]=- .Machine$double.eps*1000
print(a)
}
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
for (i in 1:7) print(loesung$solution[i])
[1] 0
[1] 0
[1] 1
[1] 0
[1] -4.669169e-17
[1] -1.436586e-17
[1] 8.881784e-16
This is a consequence of finite precision arithmetic. Just as you
should not compare to numeric numbers directly for equality but rather
that the absolute value of their difference is smaller than an
appropriate chosen threshold, you should not check whether a number is
bigger or equal to zero, but whether it is bigger or equal to an
appropriately chosen negative threshold. More information are given in
FAQ 7.31.
HTH.
Cheers,
Berwin
=== Full address =
Berwin A TurlachTel.: +65 6515 4416 (secr)
Dept of Statistics and Applied Probability+65 6515 6650 (self)
Faculty of Science FAX : +65 6872 3919
National University of Singapore
6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED]
Singapore 117546http://www.stat.nus.edu.sg/~statba
__
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] question on using gl1ce from lasso2 package
G'day Li, On Wed, 25 Jul 2007 00:50:43 -0400 Li Li [EMAIL PROTECTED] wrote: I tried several settings by using the family=gaussian in gl1ce, but none of them works. [...] gl1ce(Petal.Width~Sepal.Length+Sepal.Width+Petal.Length, data=iris,family=gaussian()) Error in eval(expr, envir, enclos) : object etastart not found Does anyone have experience with this function? Any help will be appreciated, Actually, gl1ce(Petal.Width~Sepal.Length+Sepal.Width+Petal.Length, data=iris, + family=gaussian) should work. No need to say `family=gaussian()'. However, omitting the brackets leads to an even more obscure error message and using traceback() makes me believe that the function `family()' must have been changed since this code was ported from S-Plus to R. The version with brackets does not seem to work since the function that tries to determine initial values for the iterative process of fitting a GLM tries to access an object called `etastart', which exist in the list of formal parameters of glm() but is not a formal parameter of gl1ce(). I am not sure whether this problem always existed or is also new due to changes in glm() and accompanying functions. (Time to re-work the whole lasso2 package, I guess.) A way to solve your problem is to issue the following commands: etastart - NULL gl1ce(Petal.Width~Sepal.Length+Sepal.Width+Petal.Length, data=iris, + family=gaussian()) However, since you are using the gaussian family, why not use l1ce() directly? The following command leads to the same output: l1ce(Petal.Width~Sepal.Length+Sepal.Width+Petal.Length, data=iris, + absolute=TRUE) Hope this helps. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability+65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546http://www.stat.nus.edu.sg/~statba __ 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] Optimization
G'day Massimiliano, On Mon, 16 Jul 2007 22:49:32 +0200 massimiliano.talarico [EMAIL PROTECTED] wrote: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Note that given the constraint x1+x2+x3=1 and the positivity constraints on x1, x2 and x3, the conditions that all of them should be =1 are redundant. I would go about as follows: x1+x2+x3=1 means x1=1-x2-x3 plug this into the next constraint to obtain: a*(1-x2-x3)^2+b*x2^2+c*x3^2=0.04^2 for suitable a, b, c. Note that for any value of x3, you can solve this equation for x2. Hence, take a grid of x3 values between 0 and 1 and calculate the correpsonding x2 values. From x3 and x2 calculate x1. Discard all tuples that do not fulfill the positivity constraints and calculate the function values for the remaining ones. If the grid of x3 values is fine enough, that should give you an idea what the solution should be. Alternatively, you could write a function that takes values of x3 and does all the computations outline above (return -Inf if the value x3 is not feasible) and then pass the function to optimise() for numerical optimisation... Hope this helps. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability+65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546http://www.stat.nus.edu.sg/~statba __ 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] Optimization
G'day Moshe,
On Tue, 17 Jul 2007 17:32:52 -0700 (PDT)
Moshe Olshansky [EMAIL PROTECTED] wrote:
This is partially true since both the function to be
maximized and the constraint are non-linear.
I am not sure what your definition of non-linear is, but in my book,
and I believe by most mathematical/statistical definitions, the
objective function is linear.
The only non-linearity comes in through the second constraint.
One may substitute 1-x1-x2 for x3 and use (let say) Lagrange
multipliers to get two non-linear equations with 2
unknowns for which there should be a function solving
them.
Why would you want to use Lagrange multipliers? Isn't that a bit of an
overkill? Once you substitute 1-x1-x2 for x3 in the second constraint,
you have a quadratic equations in x1 and x2. So for any given value of
x1 you can solve for x2 (or for any given value of x2 you can solve for
x1). They still teach how to solve quadratic equations at school,
don't they? ;-)
Then you must find the points where the
constraint function intersects with the triangle
{x1=0,x2=0,x1+x2=1}, which is easier (for each of
the 3 edges you get a non-linear equation in one
variable).
Even easier. Take an x1 between 0 and 1. If for that x1 the quadratic
equation in x2 has no real solution, then x1 is not feasible.
Otherwise find the values of x2 that solve the equation. Use each of
these values together with x1 to calculate corresponding values of x3.
Then check these tuples for feasibility. If they are feasible,
evaluate the objective function and return the tuple with the larger
function value.
All the calculations outlined in the paragraph above are easily
implemented in R, e.g. the function polyroot() returns the roots of a
polynomial.
Cheers,
Berwin
=== Full address =
Berwin A TurlachTel.: +65 6515 4416 (secr)
Dept of Statistics and Applied Probability+65 6515 6650 (self)
Faculty of Science FAX : +65 6872 3919
National University of Singapore
6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED]
Singapore 117546http://www.stat.nus.edu.sg/~statba
__
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] The R Book by M. J. Crawley
G'day Uwe, On Tue, 03 Jul 2007 14:33:05 +0200 Uwe Ligges [EMAIL PROTECTED] wrote: Pietrzykowski, Matthew (GE, Research) wrote: I saw the new book, The R Book, by Michael J. Crawley and wanted to know what R users thoughts of it. The author seems to be an expert in (almost?) all available statistical programming languages I would have thought that honour would go to Brian S. Everitt. :-) M.J. Crawley only seems to write books using S-PLus and R. His book on Statistical Computing (using S-Plus) is about 750 pages and his book Statistics: An introduction using R is about 320 pages. I do not know and have not seen The R Book yet, so I cannot comment on it. The statistical material presented in the other two books is pretty sound and well explained. M.J. Crawley definitely has some strong opinions on how certain data should be analysed and how statistics should be used. The S-Plus code or R code he uses can, on occasions, be somewhat improved. Cheers, Berwin __ 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] [ANN] Static and dynamic graphics course, July 2007, Salt Lake City
G'day Hadley, On Wed, 30 May 2007 10:57:54 +0200 hadley wickham [EMAIL PROTECTED] wrote: We're pleased to announce a one day course covering static and dynamic graphics using R, ggplot and GGobi. The course will be held just before the JSM, on Saturday, 28 July 2007, in Salt Lake City. The course will be presented by Dianne Cook and Hadley Wickham. Where exactly is the course held? I guess Salt Lake City is big, so is it close to the place where the JSM is held or at the other end of town? I am considering of coming a day early and to attend the course. Thanks for your help. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability+65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546http://www.stat.nus.edu.sg/~statba __ 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] R is not a validated software package..
G'day Marc, On Fri, 08 Jun 2007 14:11:48 -0500 Marc Schwartz [EMAIL PROTECTED] wrote: On Fri, 2007-06-08 at 16:02 +0200, Giovanni Parrinello wrote: Dear All, discussing with a statistician of a pharmaceutical company I received this answer about the statistical package that I have planned to use: As R is not a validated software package, we would like to ask if it would rather be possible for you to use SAS, SPSS or another approved statistical software system. Could someone suggest me a 'polite' answer? TIA Giovanni The polite answer is that there is no such thing as 'FDA approved' software for conducting clinical trials. The FDA does not approve, validate or otherwise endorse software. I like this one. :) My polite answer would have been: Sure, can do. If you pay for the license of the finally agreed upon statistical software system and pay for the time it takes me to learn it so that I can do the analysis using that system instead of R. Most clients I know would withdraw such requests if they notice that it will probably double or triple their bill. ;-) Cheers, Berwin __ 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] Bad optimization solution
G'day Paul, On Mon, 7 May 2007 22:30:32 +0100 Paul Smith [EMAIL PROTECTED] wrote: I am trying to perform the below optimization problem, but getting (0.5,0.5) as optimal solution, which is wrong; Why? As far as I can tell you are trying to minimize |x1-x2| where x1 and x2 are between 0 and 1. The minimal value that the absolute function can take is zero and any point (x1,x2)=(x,1-x) where x is between 0 and 1 will achieve this value and also respect the constraints that you have imposed. Hence, any such point, including (0.5,0.5) is a solution to your problem. the correct solution should be (1,0) or (0,1). Why? Unless there are some additional constraint that you have not told optim() (and us) about, these are two possible solutions from an infinite set of solutions. As I said, any point of the form (x, 1-x) with x between 0 and 1 is a solution to your problem, unless I am missing something Cheers, Berwin __ 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] Bad optimization solution
G'day Paul,
On Mon, 7 May 2007 23:25:52 +0100
Paul Smith [EMAIL PROTECTED] wrote:
[...]
Furthermore, X^2 is everywhere differentiable and notwithstanding the
reported problem occurs with
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
(x1-x2)^2
}
Same argument as with abs(x1-x2) holds. (x1-x2)^2 is non-negative for
all x1, x2. All points of the form (x, 1-x) where x is between 0 and 1
satisfy the constraints and achieve a function value of 0. Hence, all
such points are solutions.
There is no problem. Except if there are further constraints that
reduce the set of possible solutions which you have not told us about.
Cheers,
Berwin
__
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] Bad optimization solution
G'day all, On Tue, 8 May 2007 12:10:25 +0800 Berwin A Turlach [EMAIL PROTECTED] wrote: I am trying to perform the below optimization problem, but getting (0.5,0.5) as optimal solution, which is wrong; Why? As far as I can tell you are trying to minimize |x1-x2| It was pointed out to me that you are actually trying to maximise the function, sorry, didn't read the command carefully enough. And I also noticed that my comments were anyhow wrong, for minimisation the points (x,x) and not (x, 1-x) would be the solutions; must have also misread the function as |x1+x2|... Looks as if I need more coffee to wake up and get my brain going... Please ignore my last two mails on this issue and apologise to the list Cheers, Berwin __ 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] ED50 from logistic model with interactions
On Wed, 02 May 2007 11:37:22 +1000 Kate Stark [EMAIL PROTECTED] wrote: [...] My model is: fit - glm(Mature ~ Season * Size - 1, family = binomial, data=dat) where Mature is a binary response, 0 for immature, 1 for mature. There are 3 Seasons. I would use: fit - glm(Mature ~ Season/Size - 1, family=binomial, data=dat) With this parameterisation you get the three intercepts and the three slopes directly (together with there standard errors from summary()). Makes life simpler for your calculations. In Faraway(2006) he has an example using the delta method to calculate the StdErr, but again without any interactions. I can apply this for the first Season, as there is just one intercept and one slope coefficient, but for the other 2 Seasons, the slope is a combination of the Size coefficient and the Size*Season coefficient, [...] Not with the above parameterisation, so life is easier. I don't have my copy of Faraway (2006) handy at the moment, so I cannot vouch that you can use the method the describes now directly. But I expect you can. :) I could divide the data and do 3 different logistic regressions, one for each season, but while the Mat50 (i.e. mean Size at 50% maturity) is the same as that calculated by the separate lines regression, Im not sure how this may change the StdErr? As far as I can tell, there should be no difference. Compare the estimates and their standard errors that you obtain from the 3 different logistic fits with the estimates and standard errors from the parameterisation that I suggested. They should be the same. Hope this helps. Cheers, Berwin __ 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] Wikibooks
G'day Kaltja, On Sat, 31 Mar 2007 10:19:00 +0300 (EEST) [EMAIL PROTECTED] wrote: [...] I did not even now there was R wiki. I couldn't find a link from the R organization pages to it or am I just blind? If you talk about CRAN (e.g. http://cran.r-project.org/) then no, but if you really talk about The R Project for Statistical Computing pages, i.e. http://www.r-project.org/ then yes. :) On www.r-project.org you have the following links under documentation: Manuals FAQs Newsletter Wiki Books Other Note that CRAN (and mirrors) have a link called R Homepage under About R. Cheers, Berwin __ 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] Strange integer result on Debian/amd64
G'day Dave, On Tue, 20 Mar 2007 14:59:38 + Dave Ewart [EMAIL PROTECTED] wrote: All help/suggestions/appreciated! The calculations are done in floating point arithmetic, not integer arithmetic. From the help page on `choose' one might already guess so much, but reading R-2.4.1/src/nmath/choose.c definitely confirms this fact. Thus, your problem is covered by FAQ 7.31: http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f HTH. Cheers, Berwin __ 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] order of values in vector
G'day Tord, On Mon, 19 Mar 2007 10:05:47 +0100 Tord Snäll [EMAIL PROTECTED] wrote: but I was hoping to get this: x [1,] 2 20 [2,] 3 30 [3,] 5 50 [4,] 4 40 [5,] 6 60 [6,] 1 10 You did try rank(), didn't you? x= c(20,30,50,40,60,10) cbind(rank(x), x) x [1,] 2 20 [2,] 3 30 [3,] 5 50 [4,] 4 40 [5,] 6 60 [6,] 1 10 HTH. Cheers, Berwin __ 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] Sweave question: prevent expansion of unevaluated reused code chunk
G'day Kevin, On Tue, 13 Mar 2007 20:18:11 -0500 Kevin R. Coombes [EMAIL PROTECTED] wrote: [1] In the example I presented, you're right; I can just use a verbatim environment. In the following more complicated example [...] I would still want to be able to see the outline of what is going on, and a verbatim environment wouldn't work. If you want to have a literate programming environment in R, in particular one where you can refer to code chunks that are only defined later in the document, then I would investigate the package relax by Peter Wolf. See: http://www.wiwi.uni-bielefeld.de/~wolf/software/relax/relax.html And, for an example output just like you what you want: http://www.wiwi.uni-bielefeld.de/~wolf/software/relax/hello-world.pdf HTH. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability+65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546http://www.stat.nus.edu.sg/~statba __ 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] problem with help.start and ?somefunction
G'day Jenny, On Wed, 28 Feb 2007 16:41:18 -0600 (CST) [EMAIL PROTECTED] wrote: I am going to be teaching a workshop next week using R and Bioconductor in one of our university's computer labs. They have recently installed R 2.4.1 for me, and I'm checking all my scripts. I just noticed that using the ?somefunction call to access the documentation for that function is not working. On my own PC, the ? call output changed between R 2.3 and 2.4; before it would open some sort of plain text file and now it opens a nice browser with all the functions of that package listed in a side frame and the function documentation listed in the main window. Yes, as far as I know, one of the changes on MS Windows from R 2.3.x to R 2.4.x was that the compiled html help system became the default. At the computer lab, calling ?somefunction opens the browser with the functions in the side frame, but the main window says The page cannot be displayed On MS Windows (which is not my preferred platform, so I might use some inaccurate terminology in the sequel) I always preferred the compiled html help system. So I usually asked our IT guys to make that the default via the etc/Rprofile.site configuration file, and ran into the same problem you have about a year ago. In a nutshell, yet another security issue was discovered if compiled html help files are run over a network. Hence MS released a patch that switched this possibility off. If I remember correctly, the discussion at: http://forums.techarena.in/showthread.php?t=322640goto=nextnewest proved to be very helpful to us and enabled us to allow again the use of compiled help files over the network. But I never really understood what the security implications are and what issues really existed. I trusted our IT guys to understand that and to decide whether it was o.k. to make the necessary registry changes to allow compiled html help again. :-) If I try help.start(), I get the warnings below, but it does open the html search page. I'm guessing the two are related, No, they are not related. :-) As the warnings message says, if you say help.start() R tries to update a file to which it does not have write permissions. As I understand it, help.start() runs by default with update=TRUE (you may want to try update=FALSE) on MS Windows. The update=TRUE options means that the functions make.packages.html() and make.search.html() are called to update some information (details are on the help pages of these two functions). The updated information is written to a file, in your case P:\xpapps\R\R-2.4.1/doc/html/search/index.txt, but you have no write access to this file. Hence the warning. This is a problem with how the permissions are set up on your lab/network. Not sure how to solve this one, our lab had a setup where users had write access to the location of that file, so we never experienced the problem (except once, when due to a server upgrade the permissions on the lab machines got messed up). Hope this helps. Cheers, Berwin __ 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] How to optimize this loop ? (correction)
G'day Martin and others,
On Thu, 18 Jan 2007 15:22:07 +0100
Martin Becker [EMAIL PROTECTED] wrote:
Sorry, I accidentaly lost one line of the function code (result -0),
see below...
But even with that line, the code doesn't work properly. :)
Since you changed the order of the for loop (from i in 1:(len-1) to
i in (len-1):1) you broke the logic of the code and the result
returned is wrong (should be something like result - len-i instead
of result - i-1).
I agree that your function is faster than the one based on rle(),
presumably because rle() does just too many additional computations
that are not necessary to solve the problem at hand. But I still don't
like to see for() loops in R code. :) I would suggest the following
solution (which is only marginally slower than yours on these example
inputs):
myfun2 - function(x){
len - length(x)
which.max(c(x[len] = rev(x[-len]), TRUE) ) - 1
}
The concatenation of TRUE at the end of the logical vector is necessary
for the case that all elements before the last are smaller than the
last one. It also has the side benefit that myfun(2) will return 0
instead of numeric(0) (or something similar). On my machine:
system.time(for (i in 1:1) erg-my_fun(c(3, 4, 10,14,8,3,4,6,9)))
[1] 2.860 0.016 2.950 0.000 0.000
system.time(for (i in 1:1) erg-myfun1(c(3, 4, 10,14,8,3,4,6,9)))
[1] 0.256 0.000 0.256 0.000 0.000
system.time(for (i in 1:1) erg-myfun2(c(3, 4, 10,14,8,3,4,6,9)))
[1] 0.280 0.000 0.283 0.000 0.000
Cheers,
Berwin
=== Full address =
Berwin A TurlachTel.: +65 6515 4416 (secr)
Dept of Statistics and Applied Probability+65 6515 6650 (self)
Faculty of Science FAX : +65 6872 3919
National University of Singapore
6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED]
Singapore 117546http://www.stat.nus.edu.sg/~statba
__
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] A faster way to calculate Trace?
G'day Yongwan, YC == YONGWAN CHUN [EMAIL PROTECTED] writes: YC I want to know how to get trace of product of matrices YC **faster** when the matrices are really big. Unfortunately the YC matrices are not symmetric. If anybody know how to get the YC trace of it, please help me. An example is as below. The first one is quite simple to speed up: n - 2500 a - matrix(rnorm(n*n),n,n) b - matrix(rnorm(n*n),n,n) sum(diag(a %*% b)) [1] 1890.638 tb - t(b) sum(a*tb) [1] 1890.638 For the second one, you may try: sum(diag(a %*% b %*% a %*% b)) [1] 10668786 cmat - a %*% b sum(cmat*t(cmat)) [1] 10668786 It gives somewhat a speedup, since you only have to multiply two huge matrices once instead of thrice, but I wonder whether further improvements are possible. Hope this helps. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] Sweave and the [ function
G'day all,
RD == Ross Darnell [EMAIL PROTECTED] writes:
RD Hi Vincent This would seem logical but in this case doesn't
RD work.
Well, he said that it was untested. :)
Same idea, but with a slightly different implementation (and this one
works, I have tested it):
=
str(women)
women$height
women[,1]
@
\begin{Sinput}
[(women,1)
\end{Sinput}
echo=FALSE, eval=TRUE=
[(women,1)
@
HTH
Cheers,
Berwin
== Full address
Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr)
School of Mathematics and Statistics+61 (8) 6488 3383 (self)
The University of Western Australia FAX : +61 (8) 6488 1028
35 Stirling Highway
Crawley WA 6009e-mail: [EMAIL PROTECTED]
Australiahttp://www.maths.uwa.edu.au/~berwin
__
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] From two plot() to two X11-vindows?
AT == kone [EMAIL PROTECTED] writes: AT but how to send the information from different plot- and AT lines-commands to a certain window? ?dev.set ?dev.next ?dev.prev ?dev.cur HTH. Cheers, Berwin __ 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] Solve non-linear equatuion, grid search... [Was] RE: minimization of a quadratic form with some coef fixed and someconstrained
G'day Yingfu, YX == Yingfu Xie [EMAIL PROTECTED] writes: YX Thanks a lot for answering! I prefer the first method, which YX seems easier to carry out. So, what I need now is some method YX to solve the non-linear equation [...] I searched the R YX archive somehow, but didn't find anything valuable. Is there YX any function in R available for this problem? I am grateful of YX any hints. You have to code the function f(lam) = c'(M_11-lam*I)^(-2)c-1 in R. Note that lam is univariate. If you could easily find lam1 and lam2 with f(lam1)0 and f(lam2)0, then you could use uniroot(). But optim() could be easier to use for finding lam s.t. f(lam)=0, just try to minimise f(lam)^2. YX When c=0, my problem reduces to the typical minimization of YX b'M_11b w.r.t b'b=1. The solution is the normalized YX eigenvector associated with the minimum eigenvalue of M_11, YX right? Depends on your M_11, if M_11 is the identity matrix, then you have infinitely many solutions. Intuitively, I would have thought that the solution is the normalised eigenvector (and its negation, so there are at least two solutions) of the maximum eigenvalue of M_11, but I might be drawing the wrong picture in my mind. It should be either the maximum or minimum. Note, if that eigenvalue has geometric multiplicity 1, then there will be infinitely many solutions. YX PS: There is a type error in the first condition for b: the YX '+' should write to '-'. I believe both versions are correct. We are enforcing an equality constraint, not an inequality constraint. So it doesn't matter whether we write the Lagrangian as objective function + lam *(b'b-1) or objective function - lam *(b'b-1) In the KKT conditions, we will only have the complimentary condition that lam*(b'b-1)=0, no condition that lam=0. When you enforce inequality constraints, then you have to take care with your signs and how you write the Lagrangian. Cheers, Berwin __ 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] minimization a quadratic form with some coef fixed and some constrained
YX == Yingfu Xie [EMAIL PROTECTED] writes:
YX Thanks for reply! But I think that solution is right without
YX the constrain b'b=1. With this constrain, the solution is not
YX so simple. :(
But simple enough. :)
Write down the Lagrange function for the problem. Say, 'lam' is the
Lagrange parameter for enforcing the constraint b'b=1. Then, using
Rolf's notation:
RT [...] Write M as
RT | M_11 c |
RT | c' m |
Then the system of equations that b and the Lagrange parameter have to
fulfill is:
b = (M_11 + lam*I)^{-1} c (with I being the identity matrix)
and lam = b' M_11 b - b'c
You can either use the first equation and do a (grid) search for the
value of 'lam' that gives you b'b=1 (could be negative!), or start
with lam=0 and then alternate between the two equations until
convergence.
At least I think that this will solve your problem. :) Thinking a bit
about the geometry of the problem, I actually believe that if c=0, you
might have an identifiability problem, i.e. there are at least two
solutions, or, depending on M_11, infinitely many.
Hope this helps.
Cheers,
Berwin
== Full address
Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr)
School of Mathematics and Statistics+61 (8) 6488 3383 (self)
The University of Western Australia FAX : +61 (8) 6488 1028
35 Stirling Highway
Crawley WA 6009e-mail: [EMAIL PROTECTED]
Australiahttp://www.maths.uwa.edu.au/~berwin
__
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] order() 'decreasing =' argument must be typed in full
G'day John, JT == Thaden, John J [EMAIL PROTECTED] writes: JT Isn't is standard practice to have R base functions accept JT truncated argument names as long as they are unambiguous? No. :) The help page of order states: Usage: order(..., na.last = TRUE, decreasing = FALSE) Standard practice in R is that arguments behind ... are only matched by exact matching, i.e., you have to specify the tag of that argument in full. Details on how arguments are matched are in the section Argument matching of The R Language Definition (and all though this is a draft, that part seems to be authoritative :) ). Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] How to pass eval.max from lme() to nlminb?
G'day Andrew,
AR == Andrew Robinson [EMAIL PROTECTED] writes:
AR I'm fitting a complex mixed-effects model that requires
AR numerous iterations and function evaluations. I note that
AR nlminb accepts a list of control parameters, including
AR eval.max. Is there a way to change the default eval.max value
AR for nlminb when it is being called from lme?
Looking at the code of lme.formula, I can only find this snippet:
[...]
optRes - if (controlvals$opt == nlminb) {
nlminb(c(coef(lmeSt)), function(lmePars) -logLik(lmeSt,
lmePars), control = list(iter.max = controlvals$msMaxIter,
trace = controlvals$msVerbose))
}
else {
optim(c(coef(lmeSt)), function(lmePars) -logLik(lmeSt,
lmePars), control = list(trace = controlvals$msVerbose,
maxit = controlvals$msMaxIter, reltol = if (numIter ==
0) controlvals$msTol else 100 * .Machine$double.eps),
method = controlvals$optimMethod)
}
[...]
this seems to indicate that you can only change the values for
'iter.max' and 'trace' in the call to 'nlminb()' by setting values for
'msMaxIter' and 'msVerbose', using 'lmeControl', when calling 'lme()'.
Cheers,
Berwin
== Full address
Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr)
School of Mathematics and Statistics+61 (8) 6488 3383 (self)
The University of Western Australia FAX : +61 (8) 6488 1028
35 Stirling Highway
Crawley WA 6009e-mail: [EMAIL PROTECTED]
Australiahttp://www.maths.uwa.edu.au/~berwin
__
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] Parameterization puzzle
G'day all,
BDR == Prof Brian Ripley [EMAIL PROTECTED] writes:
BDR R does not know that poly(age,2) and poly(age,1) are linearly
BDR dependent.
Indeed, I also thought that this is the reason of the problem.
BDR (And indeed they only are for some functions 'poly'.)
I am surprised about this. Should probably read the help page of
'poly' once more and more carefully.
BDR I cannot reproduce your example ('l' is missing), [...]
My guess is that 'l' is 'pyears'. At least, I worked under that
assumption.
Interestingly, on my machine (using R 2.3.1, 2.2.1 and 2.1.1) I cannot
fit any of the Poisson GLM that Murray tried. I always get the error
message:
Error: no valid set of coefficients has been found: please supply starting
values
But I have to investigate this further. I can fit binomial models
that give me similar answers.
BDR [...] but perhaps
BDR glm(deaths ~ poly(age,2) + poly(age,1)*Smoke + offset(l),
BDR poisson)
BDR was your intention?
In this parameterisation a 'poly(age,1)' term will appear among the
coefficients with an estimated value of NA since it is aliased with
'poly(age, 2)1'. So I don't believe that this was Murray's intention.
The only suggestion I can come up with is:
summary(glm(cbind(deaths, l-deaths) ~ age*Smoke+I(age^2), family=binomial))
[...]
Coefficients:
Estimate Std. Error z value Pr(|z|)
(Intercept) -10.798950.45149 -23.918 2e-16 ***
age2.378920.20877 11.395 2e-16 ***
SmokeYes 1.445730.37347 3.871 0.000108 ***
I(age^2) -0.197060.02749 -7.168 7.6e-13 ***
age:SmokeYes -0.308500.09756 -3.162 0.001566 **
[...]
Which doesn't use orthogonal polynomials anymore. But I don't see how
you can fit the model that Murray want to fit using orthogonal
polynomials given the way R's model language operates.
So I guess the Poisson GLM that Murray wants to fit is:
glm(deaths~ age*Smoke+I(age^2)+offset(l), family=poisson)
Cheers,
Berwin
== Full address
Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr)
School of Mathematics and Statistics+61 (8) 6488 3383 (self)
The University of Western Australia FAX : +61 (8) 6488 1028
35 Stirling Highway
Crawley WA 6009e-mail: [EMAIL PROTECTED]
Australiahttp://www.maths.uwa.edu.au/~berwin
__
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] Test for equality of coefficients in multivariate multipleregression
G'day John, JF == John Fox [EMAIL PROTECTED] writes: JF Simply stacking the problems and treating the resulting JF observations as independent will give you the correct JF coefficients, but incorrect coefficient variances Yes, after Andrew's (off-list) answer I realised this too. If I am not mistaken, all variances/covariances should be off by a factor of 1/2 or something like that. JF and artificially zero covariances. Well, I must admit that I misread Ulrich's code for most of the day. I hadn't realised that the variable `tmp' introduces a correlation between `y1' and `y2' in his code: DF-data.frame(x1=rep(c(0,1),each=50),x2=rep(c(0,1),50)) tmp-rnorm(100) DF$y1-tmp+DF$x1*.5+DF$x2*.3+rnorm(100,0,.5) DF$y2-tmp+DF$x1*.5+DF$x2*.7+rnorm(100,0,.5) for some reason, my brain kept parsing this as generate *one* random intercept for y1 and *one* random intercept for y2, not that each individual observation has a random intercept. Under the model that my brain kept parsing, one would have zero covariances. :) Now I understand why Andrew suggested the use of mixed models and would go down that way too. But I believe your approach is valid too. JF (BTW, it's a bit unclear to me how much of this exchange was JF on r-help, Easy, all those that have r-help either in the TO: or CC: field. Those were Ulrich's original message and the answer by you and Andrew, I kept all my mails so far off-list. JF but I'm copying to r-help since at least one of Ulrich's JF messages referring to alternative approaches appeared there. Yes, I noticed that and answered off-list. In that message, if I read it correctly, he had confused Andrew and me. JF I hope that's OK.) Sure, why not? :) Cheers, Berwin __ 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] Problems using quadprog for solving quadratic programming problem
G'day Fabian, FB == Fabian Barth [EMAIL PROTECTED] writes: FB I'm using the package quadprog to solve the following FB quadratic programming problem. FB I want to minimize the function FB (b_1-b_2)^2+(b_3-b_4)^2 FB by the following constraints b_i, i=1,...,4: FB b_1+b_3=1 FB b_2+b_4=1 FB 0.1=b_1=0.2 FB 0.2=b_2=0.4 FB 0.8=b_3=0.9 FB 0.6=b_4=0.8 FB In my opinion the solution should be b_1=b_2=0.2 und b_3=b_4=0.8. This could well be the correct solution. For these values for the b_i the function evaluates to zero and you definitely won't get anything smaller than that. :) FB Unfortunately R doesn't find this solution and what's FB surprising to me, evaluation the solution of solve.QP with my FB function doesn't lead the minimal value calculated by FB solve.QP FB I would be very happy, if anyone could help and tell me, FB where's my mistake. Mmh, the first code that you are sending seems to involve only three unknowns. Did you try to eliminate one by hand? If so, which one. Also, in that part of the code you do not tell solve.QP that there are equality constraints. So I am a bit confused what problem the first part is supposed to solve. :-) The second code, labelled `my non working sample' seems to implement the above problem directly. The problem with that code is that the matrix Dmat that you construct is not symmetric, there is the (apparently undocumented and unchecked) assumption that Dmat is symmetric. (Note if the matrix D in the quadratic form is not summetric, then you can always replace it by (D+D^T)/2 without changing the objective function. Thus, without loss of generality, the matrix D in the objective function is always assumed to be symmetric). So adding code like Dmat - (Dmat+t(Dmat))/2 in front of your `print(Dmat)' statement would be necessary. Two minor observations. First, note that the help page for solve.QP states that the function to be minimised is -d^T b + 1/2 b^T D b where D is passed in Dmat and d is passed in dvec. In your case it might not matter, but strictly speaking, you should multiply your Dmat by 2 to take care of the 1/2 factor in the definition of the objective function. Secondly, the FORTRAN code that is used by solve.QP only uses the upper triangular part of the matrix Dmat that is passed to solve.QP. Thus, essentially you were minimizing min 1/2 (b_1^2 + b_2^2 + b_3^2 + b_4^2) under all the conditions stated. And, for that problem solve.QP gave the correct answer. Once you correct the problem, and a faster way of coding your problem is probabably the following: Dmat1 - matrix(c(2, -2, -2, 2), 2, 2) Dmat2 - matrix(0, 2, 2) Dmat - rbind(cbind(Dmat1, Dmat2), cbind(Dmat2, Dmat1)) dvec - rep(0,4) Amat - cbind(c(1,0,1,0), c(0,1,0,1), diag(4), -diag(4)) bvec - c(1, 1, 0.1, 0.2, 0.8, 0.6, -0.2, -0.4, -0.9, -0.8) you will see that there is still a problem: solve.QP(Dmat,dvec,Amat,bvec=bvec, meq=2) Error in solve.QP(Dmat, dvec, Amat, bvec = bvec, meq = 2) : matrix D in quadratic function is not positive definite! The Goldfarb-Idnani algorithm starts off by calculating the unconstrained solution. Thus, it requires that the matrix D in the objective function is positive definite. In your problem the matrix is indefinite. Hope this helps. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] predict function does not provide SE estimates for multivariate timeseries VAR models?
Michael == Michael [EMAIL PROTECTED] writes: Michael What can I do? Implement this feature and send it to the person responsible for the code for inclusion? Michael Thanks a lot! Indeed, you contribution would be very much appreciated. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] How can you buy R?
G'day Duncan, DM == Duncan Murdoch [EMAIL PROTECTED] writes: DM On 5/22/2006 3:55 AM, Berwin A Turlach wrote: I agree with you on this. Probably I was to terse in my writing and produced misunderstandings. I never intended to say something about the rights that the user has with regards to P alone. My comments were directed towards the linked product P+Q. In particular, it is not clear to me whether one can execute such a product without violating copyright laws. DM The GPL is quite explicit on this: as Deepayan said, it DM confers rights to copy, modify and redistribute P. [..] Yes, so he did. But I refer above to copyright laws, not the GPL. It is not clear to me under which rules/licence/laws the combined product falls and whether you have the right to execute it. And in some places, including Australia, copyright laws are getting really strange and restrictive, so it seems. DM Now, I suppose you might argue that executing P+Q makes a copy DM of it in memory, [...] Yes, I read arguments along these lines on gnu.misc.disucss. Personally, I wouldn't use them because I know too little about these processes (and what happens under static linking, dynamic linking and so on) to make any statements on such issues. From following such discussion, I only got the impression that these points seem to be relevant in defining when (under copyright law) a derivative work is produced. DM but I think countries that have modernized their copyright DM laws recognize that this is something you have a right to do DM with a legally acquired copy. I doubt it, I don't think that the lawyers understand these technicalities either. ;-)) Cheers, Berwin __ 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] How can you buy R?
G'day Deepayan, DS == Deepayan Sarkar [EMAIL PROTECTED] writes: DS On 5/22/06, Berwin A Turlach [EMAIL PROTECTED] wrote: DS [...] [...] Should perhaps better be formulated as: My understanding was that in that moment a product was created that would have to be wholly under the GPL, so the person who did the linking was violating the GPL and it is not clear whether anyone is allowed to use the linked product. DS I think you are still missing the point. [...] Quite possible, as I said early on IANAL. And these discussion really starts to remind me too much of those that I read in gnu.misc.discuss. Since I never participated in them, I don't see why I should here. And that group is probably a better forum to discuss all these issues. If some of the guys who always tried to argue that they found a way to circumvent the GPL are still hanging around, I am sure they are happy if you come along and confirm that according to your understanding of the GPL eveything they are doing is o.k.. :) DS The act of creating a derivative work is NOT governed by the DS GPL, Yes, as a part of the GPL that I quoted earlier states. it is the (local?) copyright law which defines when a derivative work is created. The GPL just stipulates under which licence this derivative work has to be. DS so it cannot possibly by itself violate the GPL. Fair enough. There are probably several people in gnu.misc.discuss who would be happy to hear this. :) DS The question of violation only applies when the creator of DS this derivative work wishes to _distribute_ it. This is like DS me writing a book that no one else ever reads; it doesn't DS matter if I have plagiarized huge parts of it. This point is DS not as academic as you might think. Hey, I work in an academic environment, so it is hard to imagine that I would view any point as being too academic. :) DS It is well known that Google uses a customized version of DS Linux for their servers; however, they do not distribute this DS customized version, and hence are under no obligation to DS provide the changes (and they do not, in fact). This is NOT a DS violation of the GPL. I agree and would have never claimed anything different. You are stating the obvious here. [...] If one scenario is not on, I don't see how the other one could be acceptable either. Except that in the first scenario there is a clear intend of circumventing the GPL. [...] DS That's your choice, but the situations are not symmetric, and DS quite deliberately so. That's why I studied mathematics and not law. I readily accept that there is some logic in law, it is just that I never got it. For me, if I make someone else link a GPL product P with a non-GPL product Q, then this is the same, whether I was the provider of P or Q. DS The FSF's plan was not to produce a completely independent and DS fully functional 'GNU system' at once (which would be DS unrealistic), but rather produce replacements of UNIX tools DS one by one. It was entirely necessary to allow these new DS versions to operate within the older, proprietary system. Wasn't your argument above, in response to the scenario that I was describing, that it is not necessary to explicitly allow this because a user can never violate the GPL? As long as you operate on a proprietary system and not distributing anything, why would there all of a sudden be a problem? DS In fact, GCC was not the first piece of software released DS under the GPL, My guess is that the first piece of software released under the GPL was Emacs, but it is quite likely that I will be corrected on this point. DS and until then the only way to use GPL software was to compile DS them using a non-free compiler. I know, I have compiled a lot of GPL software with non GCC compilers; and I have been using GCC when it was still standing for GNU C Compiler. But thanks for the history lesson anyhow. :) Cheers, Berwin __ 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] Accessing The Output of NLS
LI == Lorenzo Isella [EMAIL PROTECTED] writes: LI Hi, Probably a trivial question: if you do type something LI like: out-nls( here goes the model ), then you can type out LI or summary(out) to see the fitted parameters, the quality of LI the fit etc... but what if you want to get the fitted LI parameters as a vector to re-use them straight away in the LI following computations? ?coef ?? Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] glmmADMB and the GPL -- formerly-- How to buy R.
is about 18 months old, and glmmADMB is at version 0.3. I also note that the person seeking advice was asking specifically whether what he wanted to do was possible to be done using R. In my opinion, he got a sensible and correct answer at that time. DF [...] The R code that Anders wrote is simply an interface DF which takes the R specification for the model and outputs a DF data file in the format the an ADMB program expects. The ADMB DF program is a stand alone exe. The R script then reads the ADMB DF output files and presents the results to the user in a more DF familiar R format. [...] Yes, by looking at the package and the R code that is there, it is quite obvious that this is what is done. DF [...] Now it appears at some revision someone put a GPL notice DF on this package although Anders states that he did not do so, DF and and he is certain that it was not originally included by DF him. [...] Yes, someone must have done this. If you cannot figure out who did it, you have to tighten the procedures in your company. A question which could be of interest here is whether the person who put the GPL notice into the package had the authority to do so? Of course, to answer that question you first have to find out who did it. DF [...] In any event the R script is easily extracted from the DF package by those who know how to do so [...] Yes, download glmmADMB_0.3.tar.gz from your site (it was still there yesterday after Hans Skaug's e-mail), uncompress it and untar it. DF [...] and we have no problem with making the ADMB-RE source to DF the exe (TPL file) available. In fact the original was on our DF web site but was modified as we made to program more robust to DF deal with difficult data sets. The compiled TPL file links DF with our proprietary libraries and we have no intention of DF providing the source for these, [...] Well, speak with your lawyers. They can probably advise you want your obligations under the GPL is. My understanding is that if you release the exe under GPL, then releasing the ADMB-RE source for the exe would be the minimum, possibly even the source to all the libraries that the exe links too. But, as I said several times in the last days, IANAL. DF [...] Prof. Ripley seems to feel that he is a qualified DF spokesman for the open source community. I have no idea what DF the community at large feels about this. Well, I have no idea how Brian feels or what the community at large feels but: 1) The R community, and you are posting to an R mailing list, is only a small part of a community that at large subscribe to the ideas of the GPL. 2) open source is a term that encompasses many things and, IIRC, Richard Stallman and the FSF are very critical about open source and that community. 3) If asked about who are spokesman for the open source community, other names would spring to my mind. And, see 2), these names would not include Richard Stallman or the FSF. DF What follows is Hans Skaug's post with Prof. Ripley's reply. Which seems to be a private e-mail and, hence, completely a matter between Hans and Brian. DF Hans' post was an attempt to reach some sort of consensus with DF the R community so that users who so wished could continue to DF use the glmmADMB software. So far this is the only response DF we have received. I was tempted to respond, if only to point out that glmmADMB_0.3.tar.gz was still available from your web-site. But then, since IANAL, I refrained since: a) I don't feel qualified on giving you advice on how to sort out your licencing issues; and b) I already pointed out earlier in a reply to Spencer Grave how the license of glmmADMB should probably be formulated to be above water and not in conflict with the GPL. Just read up that e-mail. DF I guess it is up to the R community to decided whether DF Prof. Ripley speaks for all of you. I might flog a dead horse here, but I hope that you will eventually get the message: It was a private, off-list e-mail. You can't seriously expect that opinions expressed in such an e-mail can be taken as speaking for the R community. Even if it had been an e-mail to the R mailing list, such e-mails nly express the opinions of the sender, not that of the members of the mailing list and not that of the R community (which could encompass more than the readers of a mailing list) at large. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin
Re: [R] How can you buy R?
G'day Deepayan, DS == Deepayan Sarkar [EMAIL PROTECTED] writes: DS let me first summarize this sub-discussion so far: [...] Sound like a perfect summary. :) DS As far as I can tell (and please correct me if I'm wrong), DS your contention is that by linking a GPL component P with a DS non-GPL component Q, a user may lose the rights granted to him DS by the GPL to the GPL-d part P. I don't think that I said this explicitly, but I can see how what I said can be interpreted in such a way. The point is rather that at the moment component P and Q are linked (and I perhaps carelessly assumed that the user was doing this) a product is produced that should be completely under the GPL. Obviously it is not. Hence, the status of this linked product, and whether it can be used by anybody, is an open question. And the answer is probably given by the copyright laws (and others?) of the country in the linking happens. DS Let's assume this is true. All that means is that the user has DS lost his rights to copy, modify and redistribute P. He does DS NOT lose the rights to use P. I agree with you on this. Probably I was to terse in my writing and produced misunderstandings. I never intended to say something about the rights that the user has with regards to P alone. My comments were directed towards the linked product P+Q. In particular, it is not clear to me whether one can execute such a product without violating copyright laws. Thus, the last sentence of mine that you quoted: My understanding was that in that moment a product was created that would have to be wholly under the GPL, so the user was violating the GPL and lost the write to use your package. Should perhaps better be formulated as: My understanding was that in that moment a product was created that would have to be wholly under the GPL, so the person who did the linking was violating the GPL and it is not clear whether anyone is allowed to use the linked product. A simple google search would have confirmed to you that the linux kernel is developed under the GPL. [...] DS Linux is under GPL2, and not GPL2 or later. [...] Oh, I wasn't aware that they did not use the typical(?) or later phrase. Thanks for pointing this out and I note that we both agree that the linux kernel is definitely not under LGPL. DS In any case, this is the complete opposite of the situation we DS were originally discussing: [...] [...] So I have to wonder to what you are referring to as the situation we were originally discussing. DS I was referring to your question (quoted above) about use of DS GPL'd code in S-PLUS, which is what I was replying to. As I DS was saying, that situation is the opposite of the one in your DS example. O.k., sorry, I used a different scale with the time point of origin at Spencer's e-mail and my answer to that mail. Now I am with you. Agreed, the situation is the opposite, but that was the example discussed in gnu.misc.discuss. From an abstract point of view the situations are the same. You make someone else link a GPL product with a non-GPL product creating a derived work, the derived work would have to be under the GPL but is not. Hence, the derived work has a legal status that is in limbo and it is not clear whether anyone has to right to use it. The discussions on gnu.misc.discuss were centred on cases were people provided non-GPL binaries, asked their users to download GPL software from elsewhere, compile and link everything together and then use the combined product. As you say it is the exact opposite (and hence mirror image) from the situation that I was worried about, where I provide GPL software and ask others to compile and link it with non-GPL binaries and then use the combined product. If one scenario is not on, I don't see how the other one could be acceptable either. Except that in the first scenario there is a clear intend of circumventing the GPL. But I was not sure whether such kind of intent makes any difference. Thus, to avoid all these problems I decided to rather use the LGPL since that licence definitely seemed to allow both. Hope this clarifies some of my comments. Cheers, Berwin __ 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] How can you buy R?
G'day Brian, BDR == Prof Brian Ripley [EMAIL PROTECTED] writes: BDR The issue in the glmmADMB example is not if they were BDR required to release it under GPL I should probably bow to your superior command of the English language and trust that you can interpret Spencer's questions much better than I, but I was addressing, amongst other things, the following comment: SG A boundary case is provided by the glmmADMB package. As I SG read the GPL, this package must operate under GPL. which to me seems to ask exactly about this issue. BDR (my reading from the GPL FAQ is that they probably were not, BDR given that communication is between processes and the R code BDR is interpreted). So, it seems we agree. :) (Though for different reasons) BDR Rather, it is stated to be under GPL Indeed, and I noted so. Furthermore, I thought it was rather pointless to confirm that under the licence of the package as it is stated at the moment they would actually be required to provide the source code of the binaries. My apologies for not doing a thorough discussion of all possible scenarios. I just pointed out that if the developers of this package do not want to provide the source code for these binaries, they should probably state another licence for them in the DESCRIPTION file and, since the binaries are not loaded dynamically, they would not be obliged to release the source code; a statement that you seem to agree to. BDR [...] As the executables are not for my normal OS and I BDR would like to exercise my freedom to try the GPLed code, I BDR have requested the sources from the package maintainer. Good luck. :) BDR Once again, the GPL FAQ and its references, BDR http://www.gnu.org/licenses/gpl-faq.html, are a more informed BDR source than mailing lists. If you think you understand it, BDR try the exam at BDR http://www.gnu.org/cgi-bin/license-quiz.cgi BDR (cheaper than testing in court). Well, if you read this material, you get the opinions of the FSF on these matters, other people might have other opinions/interpretations. If you read the material really careful, you will notice that there is one point, namely what exactly constitutes combining two parts into one program, for which even the FSF concedes that it is a legal question, which ultimately judges will decide. And it is exactly this point which often raises discussions about the GPL, e.g. (parts of) the current discussion. Luckily, the GPL is very well written and so far nobody with deep enough pockets was found who really wanted to have a definite answer to this question. Cheers, Berwin __ 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] How can you buy R?
G'day Deepayan, DS == Deepayan Sarkar [EMAIL PROTECTED] writes: DS A user can never violate the GPL. The GPL does not govern use, DS it governs distribution. Specifically, As I said, I stopped reading gnu.misc.discuss long time ago, but if I remember correctly sometimes in the (early?) '90s the following case was discussed. A company made a binary module available for download and gave the instructions go to the FSF site (or a mirror), download version X.Y.Z of program U and compile it with these options, then link our module and start the program, now you can use these features of ours and are in business. (Remember, these were the days when most free software was only available in .tar.gz form, people were used to compile their own software and slackware was the dominant (only?) linux distribution, no .rpm or .deb files.) Note also that they did not distribute any GPL code, they said go and get it. As far as I remember, they were told by the FSF that they cannot do this and had to stop. And, IIRC, the argument was that whether they did the linking or let the end user do the linking was the same and, hence, the GPL was violated. But it is quite possible that the argument was based on other sections of the GPL. And, obviously, there must be some mechanism in the GPL that prohibits the above procedure, otherwise it would be very easy to circumvent the GPL. (This idea of circumventing the GPL was regularly floated on gnu.misc.discuss while I followed it.) Cheers, Berwin __ 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] How can you buy R?
G'day Deepayan, DS == Deepayan Sarkar [EMAIL PROTECTED] writes: DS On 5/21/06, Berwin A Turlach [EMAIL PROTECTED] wrote: DS A user can never violate the GPL. The GPL does not govern use, DS it governs distribution. Specifically, As I said, I stopped reading gnu.misc.discuss long time ago, but if I remember correctly sometimes in the (early?) '90s the following case was discussed. A company made a binary module available for download and gave the instructions go to the FSF site (or a mirror), download version X.Y.Z of program U and compile it with these options, then link our module and start the program, now you can use these features of ours and are in business. (Remember, these were the days when most free software was only available in .tar.gz form, people were used to compile their own software and slackware was the dominant (only?) linux distribution, no .rpm or .deb files.) Note also that they did not distribute any GPL code, they said go and get it. As far as I remember, they were told by the FSF that they cannot do this and had to stop. And, IIRC, the argument was that whether they did the linking or let the end user do the linking was the same and, hence, the GPL was violated. DS I'll readily concede that my interpretation may be flawed, but DS this example doesn't seem to contradict anything I said. This DS binary module was clearly (in the opinion of the FSF) a DS derivative work of something that was GPL, and hence the DS company was violating the GPL by distributing the binary DS module under a license other than the GPL. Well, the question is whether the binary module is a derivative work. Note, the GPL makes four (4) references to derivatives and derivative work. The one most relevant in this discussion is probably the one in paragraph 0.: [...] The Program, below, refers to any such program or work, and a work based on the Program means either the Program or any derivative work under copyright law: that is to say, a work containing the Program or a portion of it, either verbatim or with modifications and/or translated into another language. [...] So they are actually referring to copyright law to define what a derivative work is, and this law might be different in different countries. Thus, whether you violate the GPL or not may well depend on where you are located. Another question is when does the work becomes a derived work? I believe few will disagree that at the moment the two components are dynamically linked, a derivative work is produced. Whether something is a derivative work before the linking, is a debatable question. But, as I wanted to illustrate with this example, even if you make the user do the linking, in that moment a product is created that violates the GPL and the user looses all rights to the GPL part of this new product (at least that is may understanding). Hence, a user can (be made to) violate the GPL. DS Also, as far as I can tell, your description applies to the DS situation with Nvidia's binary kernel drivers for the Linux DS kernel (which is GPL and not LGPL AFAIK), A simple google search would have confirmed to you that the linux kernel is developed under the GPL. There are actually reports that the developers are currently discussing to move to GPL 3 (a bit strange, since GPL 3 is, AFAIK, open to discussion but not yet released) and many of them wanting to stick with GPL 2. (Another thing I find strange, because given the standard clause, one can take GPL 2 code, modify it and then release the new version under GPL 3 [or later].) DS which is obviously tolerated, so there must have been some DS other nuances. It is, but there are problems. I have computers with an Nvidia card and run Debian. Quite often, when the kernel image pacakge is updated (even if the same version numbers are kept!), my system breaks and until I recompile the interface to the Nvidia binaries, I cannot use X. Interfacing proprietary binary modules to linux seems to be a perennial problem. Recently (0.5-1 year ago), there were some reports on relevant sites that somebody who was providing support for binary drivers of some digital camera got very upset with how it was made harder and harder for him to provide this support. And when some code was removed from the kernel which made it possible for him to provide the support (and the argument for the removal was that the code only serves this purpose), he really spit the dummy and asked that his code and any other code that he had provided to the linux kernel be removed. As far as I remember, Linus complied but a huge discussion ensued on whether he was actully allowed to withdraw all his code. I did not follow this incidence though to see how it was eventually resolved. DS In any case
Re: [R] intervals from cut() as numerics?
G'day Gavin, GS == Gavin Simpson [EMAIL PROTECTED] writes: GS The problem is getting the range/interval for each group from GS (4,4.3], so I can automate this. Most likely there is an easier way, but this seems to work: ## get the levels of groups: tmp - levels(groups) ## remove the opening ( and closing ] from the string: tmp1 - sapply(tmp, function(x) substr(x, 2, nchar(x)-1)) ## split into two character strings: tmp2 - strsplit(tmp1, ,) ## turn into results into two numbers: tmp3 - lapply(tmp2, as.numeric) ## Of course, we can do everything in one go: lapply(strsplit(sapply(levels(groups), function(x) substr(x, 2, nchar(x)-1)), ,), as.numeric) $(4.05,4.32] [1] 4.05 4.32 $(4.32,4.6] [1] 4.32 4.60 $(4.6,4.87] [1] 4.60 4.87 $(4.87,5.15] [1] 4.87 5.15 $(5.15,5.43] [1] 5.15 5.43 $(5.43,5.7] [1] 5.43 5.70 $(5.7,5.98] [1] 5.70 5.98 $(5.98,6.25] [1] 5.98 6.25 $(6.25,6.53] [1] 6.25 6.53 $(6.53,6.8] [1] 6.53 6.80 Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] How can you buy R?
definite answers to these licences questions by testing them in court. Good luck if you choose to do so. :-) Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] No output in sourced R program
G'day Sigbert, SK == Sigbert Klinke [EMAIL PROTECTED] writes: SK Hi, If I type it in the command line I get, as expected: 1:30 SK [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 SK 23 24 25 [26] 26 27 28 29 30 q() SK Save workspace image? [y/n/c]: n SK If I create a program SK 02e451444d6a46acf551996579092c911b90aa8e.R and run it I get SK R : Copyright 2006, The R Foundation for Statistical Computing SK Version 2.3.0 (2006-04-24) source(02e451444d6a46acf551996579092c911b90aa8e.R) SK Save workspace image? [y/n/c]: n SK mars:/srv/www/htdocs/mediawiki/teachwiki/Rfiles # more SK 02e451444d6a46acf551996579092c911b90aa8e.R SK rfiles-/srv/www/htdocs/mediawiki/teachwiki/Rfiles 1:30 q() SK No output from 1:30 to the screen. Any idea what I do wrong ? Autoprint of objects is only enabled on the command line, on all other levels you have to explictly call the print() command. So your script should be: print(1:30) q() BTW, interesting name for the file. How do you select that one? ;-) Cheers, Berwin PS: I heard rumours that Wiwi had HU-Berlin was using R instead of XploRe for some tasks. Slowly I start to believe it... :-)) == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] vector-factor operation
G'day Murray, MJ == Murray Jorgensen [EMAIL PROTECTED] writes: MJ I found myself wanting to average a vector [vec] within each MJ level of a factor [Fac], returning a vector of the same length MJ as vec. I presume that the vector that you want as result should not just have the same length as vec, should it? :-) MJ But there must be another way to do this, and it would be good MJ to be able to apply other functions than mean() in this way. Something like: Fac - sample(gl(2,4)) Fac [1] 2 1 1 2 2 1 1 2 Levels: 1 2 vec - rnorm(length(Fac)) tapply(vec, Fac, mean) 1 2 -0.6435816 -0.9267021 tapply(vec, Fac, mean)[Fac] 2 1 1 2 2 1 1 -0.9267021 -0.6435816 -0.6435816 -0.9267021 -0.9267021 -0.6435816 -0.6435816 2 -0.9267021 tapply(vec, Fac, sum)[Fac] 2 1 1 2 2 1 1 2 -3.706808 -2.574326 -2.574326 -3.706808 -3.706808 -2.574326 -2.574326 -3.706808 Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] Odd anova(lm()) order phenomenon, looking for an explanation
G'day Andrew, AR == Andrew Robinson [EMAIL PROTECTED] writes: AR I would have expected any term to explain less Sum Sq if AR listed second than if listed first. Is my intuition awry? Yes. :-) I would not expect that *any* term explains less Sum Sq if listed second, then life and (linear) modelling would be simple. The problem with multiple regression is that a covariate might look unimportant if used first (i.e. has small Sum Sq associated with it in the anova table), but if we first correct for other regressor, then this covariate becomes important all of a sudden (i.e. has large Sum Sq associated with it in the anova table). What surprised me, was that you observed this phenomenon with respect to both regressors. If only one had displayed this behaviour, I would have readily explained it as above, but that both display it, I found surprising too. AR Does anyone have any modelling insight to help me interpret AR what I'm seeing? Don't know if the following example, which shows the same behaviour, leads to any insight. n - 100 x1 - runif(n, -1,1) x2 - runif(n, -1,1) y - x1*x1*x2 + rnorm(n, sd=0.05) y - y - mean(y) anova(lm(y~x1+x2)) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(F) x1 1 0.0055 0.0055 0.1485 0.7008 x2 1 5.0071 5.0071 134.3499 2e-16 *** Residuals 97 3.6151 0.0373 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(lm(y~x2+x1)) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(F) x2 1 4.9930 4.9930 133.9723 2e-16 *** x1 1 0.0196 0.0196 0.5261 0.47 Residuals 97 3.6151 0.0373 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] Derivative of a splinefun function.
G'day Rolf, RT == Rolf Turner [EMAIL PROTECTED] writes: RT Is there a way of calculating the derivative of a function RT returned by splinefun()? Yes, of course. :) As somebody else once said: this is R, anything can be done, the question is just how easy it is to do so. It may depend on what kind of spline you are fitting with splinefun()... RT Such a function is a cubic spline, whence it has a calculable RT derivative, but is there a (simple) way of getting at it? Depends on your definition of simple. ;-)) RT One workaround that I have thought of is to take a fine grid RT of points, evaluate the function returned by splinefun() at RT these points, put an interpolating spline through these points RT using smooth.spline() with spar=0, and then calculate the RT derivative with predict(...,deriv=1). Why easy if one can also make it complicated? There is the function interpSpline() in the package splines for which, according to its help page, the predict method also works. So using smooth.spline with spar=0 does indeed look a bit kludgy. :) RT Is there a better way? splinefun() essentially returns a function that contains one object in its environment, namely z, which contains all the necessary information for calculating the spline and whose body is just a call to an internal C function. The components y, b, c and d of the object z contain the piecewise polynomial representation between knots of the spline if you choose the methods fmm or natural. If you use the method periodic then the component e of the object z is used too, but I don't know whether it is just used to hold intermediate results while calculating the spline or whether it is actually necessary for calculating the spline. (An analysis of the C code that is called to calculate the spline and to evaluate the spline should answer that question but until now I never bothered of doing this since I didn't have use for periodic splines yet.) Thus, since the spline is stored in a piecewise polynomial fashion, it shouldn't be hard to define a function that calculates the derivative of the function returned by splinefun(). If my calculus is correct, the following code should construct such a function: par(mfrow=c(1,2)) n - 9 x - 1:n y - rnorm(n) f - splinefun(x, y) ls(envir = environment(f)) [1] z splinecoef - eval(expression(z), envir = environment(f)) curve(f(x), 1, 10, col = green, lwd = 1.5) points(splinecoef, col = purple, cex = 2) ## ## construct a function that calculates the derivative of f() ## g - f environment(g) - new.env(parent=baseenv()) z - get(z, envir=environment(f)) z$y - z$b z$b - 2*z$c z$c - 3*z$d z$d - rep(0, z$n) assign(z, z, envir=environment(g)) curve(g(x), 1, 10, col = green, lwd = 1.5) As indicated above, I would only vouch for this method if splinefun() is called with method=fmm or method=natural. And the second caveat is, of course, that you are assuming that I get my calculus correct on a Saturday morning which might be a bit much to ask. ;-)) Hope that helps. Cheers, Berwin __ 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] Additional arguments in S3 method produces a warning
G'day Martin,
MM == Martin Maechler [EMAIL PROTECTED] writes:
MM [...]
MM The grain of truth in Henrik's statement is that for a generic
MM function, S3 or S4, one should carefully consider which
MM arguments should be shared by all methods and which not. But
MM having at least one `` non-... '' argument should be the rule,
MM and is even a necessity for S4. Hence I'd strongly discourage
MM defining generic functions with only a (...) argument list.
I was pondering this question recently too and was already wondering
about sending an e-mail to ask for clarification to r-help or r-devel.
What are the arguments against having only a (...) argument list for
the generic?
I can think of several arguments why using only a (...) argument list
should be o.k.
First, I noticed that there are examples in R base that use this
construct, e.g.:
seq
function (...)
UseMethod(seq)
environment: namespace:base
Secondly, if I plan to write several methods for a generic and for
these methods the canonical name for the first parameter differ, it
would make sense to me to use only a (...) argument list for the
generic. As an example from R, take the princomp() command:
princomp
function (x, ...)
UseMethod(princomp)
environment: namespace:stats
princomp.default
Error: object princomp.default not found
getAnywhere(princomp.default)
[...]
function (x, cor = FALSE, scores = TRUE, covmat = NULL, subset =
rep(TRUE,
nrow(as.matrix(x))), ...)
{
[...]
}
environment: namespace:stats
getAnywhere(princomp.formula)
[...]
function (formula, data = NULL, subset, na.action, ...)
{
[...]
}
environment: namespace:stats
I would argue that it is natural to use formula as the first
argument to princomp.formula and x as the first argument to
princomp.default. But as the generic is princomp(x, ...), the
princomp.formula method does not comply with two of the three
recommendations in the chapter on Generic functions and methods in
Writing R Extensions.
Thus, without more information what happens if these recommendations
are not followed, I would hesitate to use such code in my own
packages. I would either make the generic of princomp have only a
(...) argument, or call the first argument of the princomp.formula
method x (which would obfuscate the code of that function a bit).
Lastly, I would suggest if the generic for transform would have been:
transform
function (...)
UseMethod(transform)
environment: namespace:base
instead of the current
transform
function (x, ...)
UseMethod(transform)
environment: namespace:base
Then a certain R user would not have reported a 10 sec puzzle to
r-devel recently:
https://stat.ethz.ch/pipermail/r-devel/2006-February/036492.html
(And I guess for many other R users solving this puzzle would have
taken more than 10 sec.)
And it is hard to see what purpose the argument x is actually
serving in argument list of the generic. So why not use a (...)
argument list?
Thus, what speaks against a (...) argument list for a S3 generic?
Cheers,
Berwin
__
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] Help - lm, glm, aov results inconsistent with other stati stical package
BR == Ben Ridenhour [EMAIL PROTECTED] writes: BR Alright, I'll try to give some sample code. BR # create A with 2 levels - 2 and 4 O.k., let us add a set.seed(1) to make thinks reproducible. BR A-c(rep(2,times=30),rep(4,times=42)) BR # make A a factor BR A-as.factor(A) BR #generate 72 random x points BR x-rnorm(72) BR #create different slopes intercepts for A=2 and A=4 BR #add a random error term BR y-c(x[(1:30)]+1,2*x[(31:72)]-2)+rnorm(72,mean=0,sd=2) BR #use model y~A*x for lm (or glm) BR test.lm-lm(y~A*x) BR This essentially creates something like my data set and uses BR the same model. In response to (1), I was just using 0/1 BR because these are codings for the 2 levels, correct? Actually, no, the internal coding of the factor is with 1 and 2: str(A) Factor w/ 2 levels 2,4: 1 1 1 1 1 1 1 1 1 1 ... unique(as.numeric(A)) [1] 1 2 and how factors are internally coded is documented on the help page of `factor'. BR In response to (2), yes, I do want different slopes for the BR two categories (that is what I am interested in testing). Whether you get estimates of different slopes does not depend on how the factor is internally coded but how your design matrix is constructed given your model specification and how the factor is coded when the design matrix is created. R has plenty of options here. In your above example, you created data from 1+x for one group and -2+2*x for the other groups. Assuming that you have not changed the way R treats factors, some possibilities would be: lm(y~x*A) Call: lm(formula = y ~ x * A) Coefficients: (Intercept)x A4 x:A4 0.9515 1.2449 -3.1825 0.8744 In this case, the first and second values are estimates for the intercept and the slope, respectively, of the regression line in the first group. The other two values are the estimates for the *differences* in intercept and slope for the two groups. lm(y~A/x) Call: lm(formula = y ~ A/x) Coefficients: (Intercept) A4 A2:x A4:x 0.9515 -3.1825 1.2449 2.1193 In this case, the first value is the estimate for the intercept in the first group. The second value is the *difference* in intercept for the two regression lines. The last two values are estimates for the slopes in the two groups. lm(y~A/(x-1)) Call: lm(formula = y ~ A/(x - 1)) Coefficients: A2 A4 A2:x A4:x 0.9515 -2.2309 1.2449 2.1193 In this case, the first two values are the estimates for the intercepts for the two regression lines and the last two values give the two slopes. Quite reasonable estimates in my book. BR If I export the data created above to JMP and run what I BR believe to be the same model, I get a different answer for my BR equations :( Presumably because you are using different parameterisations and you interpret the estimates incorrectly. To interpret the output, you have to know what parameters your estimates are estimating. I do not know what JMP is doing or what kind of parameterisation it uses given a certain model specification. [...] Why can't I get the same answer from JMP as in R? This is very disturbing to me! Rest assured, no matter how disturbing this might be for you, it is much more disturbing for us to see people using statistical packages without knowing what they are doing. :) And I would find it particular disturbing if there is no statistician (or someone with sufficient statistical training) at Washigton State University who could explain all this to you. Best wishes, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] command line completion in R?
DR == David Ruau [EMAIL PROTECTED] writes: DR Hi, I was wondering if there were a option to have command DR line completion in R like in the Bash shell? Of course, via Emacs and ESS :-) For ESS see: http://ess.r-project.org/ For Emacs see: http://www.gnu.org/software/emacs/emacs.html In fact, I got so used to command line completion, that I am hopeless on any other interface to R since I cannot remeber the correct spelling of (some) commands ;-) HTH. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] List Conversion
LD == Liz Dem [EMAIL PROTECTED] writes: LD I have a list (mode and class are list) in R that is many elements long and of the form: length(list) LD [1] 5778 list[1:4] LD $ID1 LD [1] num1 LD $ID2 LD [1] num2 num3 LD $ID3 LD [1] num4 LD $ID4 LD [1] NA LD I'd like to convert the $ID2 value to be in one element rather LD than in two. It shows up as c(\num2\, \num3\) if I try to LD use paste(list[2], collapse=). You want list[[2]], not list[2]: tt - list(ID1=num1, ID2=c(num2, num3), ID3 = num4, ID4 =NA) tt $ID1 [1] num1 $ID2 [1] num2 num3 $ID3 [1] num4 $ID4 [1] NA paste(tt[2], collapse=) [1] c(\num2\, \num3\) paste(tt[[2]], collapse=) [1] num2num3 paste(tt[[2]], collapse= ) [1] num2 num3 LD I need to do this over the length of the entire list, however LD list2 - apply(list, 1, paste, collapse=) tells me 'Error in LD apply : dim(X) must have a positive length'. You want to use `lapply' not `apply': tt1 - lapply(tt, paste, collapse= ) tt1 $ID1 [1] num1 $ID2 [1] num2 num3 $ID3 [1] num4 $ID4 [1] NA LD What I want to get is: list[1:4] LD $ID1 LD [1] num1 LD $ID2 LD [1] num2 num3 LD $ID3 LD [1] num4 LD $ID4 LD [1] NA In that case, if you don't want NA's to turn into strings: tt2 - lapply(tt, function(x) if(is.na(x[1])) NA else paste(x, collapse= )) tt2 $ID1 [1] num1 $ID2 [1] num2 num3 $ID3 [1] num4 $ID4 [1] NA HTH. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] reading in a tricky computer program output
G'day Taka,
TM == Taka Matzmoto [EMAIL PROTECTED] writes:
TM and then assign the character vector to the numeric vector by
TM names-first.10
TM first.10 = numeric.vector
TM combined.one - cbind(names,first.10)
TM container - diag(10)
TM for (i in 1:(10*10))
I don't really understand this loop. If I reverse-engineer this code
correctly thenthe matrix `combined.one' is not a 2*100 matrix, so you
should get an error while exectuting this loop.
TM Is there any other neat way to do this?
Neat way to create those character vectors? Or a neat way to read in
the data from a file?
If the latter, I would use the following code:
matzmoto - function(file, diag=TRUE){
dat - scan(file)
if(diag){
nvar - sqrt(2*length(dat)+0.25) - 0.5
nn - nvar
}else{
nvar - sqrt(2*length(dat)+0.25) + 0.5
nn - nvar - 1
}
res - matrix(0,nvar,nvar)
ind - upper.tri(res, diag=diag)
rind - 1:5
while(nn 0){
if( nn 5 ){
rind - rind[1:nn]
tmp - matrix(0,nn,nvar)
}else{
tmp - matrix(0,5,nvar)
}
how.many - sum(ind[rind,])
tmp[ind[rind,]] - dat[1:how.many]
res[rind,] - tmp
dat - dat[-(1:how.many)]
rind - rind + 5
nn - nn - 5
}
t(res)
}
res - matzmoto(matzmoto1.dat, TRUE)
print(res)
res - matzmoto(matzmoto2.dat, FALSE)
print(res)
After storing the two examples that you posted into the files
matzmoto1.dat and matzmoto2.dat, respectively, and removing the part
that you said you have added, I get the following result on my machine
when sourcing the above code:
source(matzmoto.R)
Read 55 items
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0
[2,] 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0
[3,] -0.002 0.019 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0
[4,] 0.012 -0.004 -0.020 0.000 0.000 0.000 0.000 0.000 0.00 0
[5,] -0.015 0.003 0.011 0.008 0.000 0.000 0.000 0.000 0.00 0
[6,] 0.005 -0.008 -0.005 0.002 0.005 0.000 0.000 0.000 0.00 0
[7,] 0.008 -0.007 0.013 0.003 0.007 -0.037 0.000 0.000 0.00 0
[8,] -0.014 -0.011 -0.010 -0.025 0.002 0.010 0.027 0.000 0.00 0
[9,] 0.006 0.003 -0.010 0.002 -0.020 0.032 -0.004 0.008 0.00 0
[10,] 0.006 0.010 -0.006 0.005 0.008 -0.008 -0.011 0.015 -0.02 0
Read 45 items
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0
[2,] -0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0
[3,] 0.003 -0.053 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0
[4,] -0.026 0.010 0.045 0.000 0.000 0.000 0.000 0.000 0.000 0
[5,] 0.023 -0.008 -0.025 -0.016 0.000 0.000 0.000 0.000 0.000 0
[6,] -0.012 0.023 0.013 -0.005 -0.011 0.000 0.000 0.000 0.000 0
[7,] -0.031 0.031 -0.054 -0.013 -0.027 0.127 0.000 0.000 0.000 0
[8,] 0.040 0.042 0.031 0.075 -0.007 -0.035 -0.166 0.000 0.000 0
[9,] -0.012 -0.009 0.023 -0.005 0.037 -0.083 0.015 -0.027 0.000 0
[10,] -0.013 -0.027 0.014 -0.013 -0.020 0.021 0.047 -0.052 0.048 0
The documentation of the function is quite simple. Just pass the name
of the file in which the output is and whether it is a file that
includes the output of the diagonal or not. If you don't trust the
calculation of how many variables are involved, then you might want to
change the function so that this is another input paramter.
HTH.
Cheers,
Berwin
__
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] reading in a tricky computer program output
BAT == Berwin A Turlach [EMAIL PROTECTED] writes:
TM == Taka Matzmoto [EMAIL PROTECTED] writes:
TM and then assign the character vector to the numeric vector by
TM names-first.10
TM first.10 = numeric.vector
TM combined.one - cbind(names,first.10)
TM container - diag(10)
TM for (i in 1:(10*10))
BAT I don't really understand this loop. If I reverse-engineer this code
BAT correctly thenthe matrix `combined.one' is not a 2*100 matrix, so you
BAT should get an error while exectuting this loop.
This should have been: ... `combined.one' is not a 100*2 matrix ...
TM Is there any other neat way to do this?
BAT Neat way to create those character vectors? Or a neat way to read in
BAT the data from a file?
BAT If the latter, I would use the following code:
Before somebody else points out the obvious optimisation of my code,
here is a (slightly) improved version of the function
matzmoto - function(file, diag=TRUE){
dat - scan(file)
nn - nvar - sqrt(2*length(dat)+0.25) - 0.5
if(!diag){
nvar - nvar + 1
}
res - matrix(0,nvar,nvar)
ind - upper.tri(res, diag=diag)
rind - 1:5
while(nn 0){
if( nn 5 ){
rind - rind[1:nn]
}
how.many - sum(ind[rind,])
res[rind,][ind[rind,]] - dat[1:how.many]
dat - dat[-(1:how.many)]
rind - rind + 5
nn - nn - 5
}
t(res)
}
Cheers,
Berwin
__
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] How to convert decimals to fractions
G'day Muhammad, MS == Muhammad Subianto [EMAIL PROTECTED] writes: MS Are there any functions to convert decimals to fractions in R? MS I have the result: Something like: library(MASS) as.fractions(c(0, 0.0133, 0.04, 0.0667, 0.0933, 0.107, 0.133, 0.147, 0.16, 0.187,0.2, 0.227, 0.253, 0.267, 0.293, 0.32, 0.333, 0.347, 0.373)) [1] 0 1/75 1/25 1/15 7/75 8/75 2/15 11/75 4/25 14/75 1/5 17/75 [13] 19/75 4/15 22/75 8/25 1/3 26/75 28/75 ? Cheers, Berwin __ 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] D(dnorm...)?
G'day Spencer, SG == Spencer Graves [EMAIL PROTECTED] writes: SG I'm not qualified to make this suggestion since I'm incapable SG of turning it into reality, [...] This statement applies to me too, nevertheless I would like to point out the following GPL library: http://www.gnu.org/software/libmatheval/ I am wondering since some time how hard it would be to incorporate that library into R and let it provide symbolic differentiation capabilities for R. Cheers, Berwin __ 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] Powell's unconstrained derivative-free nonlinear least squares routine, VA05AD
G'day David, DK == David Kinniburgh [EMAIL PROTECTED] writes: DK I have used Mike Powell's optimization routine (VA05AD) from DK the Harwell Subroutine Library (HSL) for more than 20 years. DK [...] DK Now that I have converted to R, I will miss my trusted DK friend. I have started using nls() but have not accumulated DK enough experience to compare the two. It would be great if DK VA05AD could be an option in there (algorithm=Powell). To DK this end, I recently enquired of the custodians of the HSL DK whether it would be possible to make it freely available to DK the R community. The answer was basically 'yes'. [...] You better ask again. :-) I presume VA05AD is in what is now called the HSL archive? If it is part of HSL 2004, then I don't see a way how it could be incorporated into R given the information on their web site. Even for the HSL archive, it is stated at http://hsl.rl.ac.uk/archive/hslarchive.html that [...] HSL Archive codes are now available without charge to anyone, so long as they are not then incorporated into a commercial product; the latter is, of course, still permitted subject to a commercial licence. [...] which sounds as if it would be possible to interface (parts of) HSL archive to R. But then that URL goes on ans states: Access to the Archive is by means of a short-lived individual password-controlled account. Potential users are asked for brief details of the use they intend to make of the package(s) they aim to download. Users are also asked to accept a conditions-of-use form, and are not permitted to divulge their userid and password to anyone else, nor to distribute any codes they obtain to a third party. IANAL, but as I understand the GPL, the last part of that sentence will make it pretty impossible to incorporate into R (or distribute as an R package) (parts of) the HSL archive routines. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] how to specify dev.print target by a variable?
G'day Leif, LK == Leif Kirschenbaum [EMAIL PROTECTED] writes: LK How do I do this such that I can specify FOO to be one of LK several choices? (GDD, PNG, postscript, etc.) If I make FOO a LK character variable, then dev.print complains. Mmh, I am not sure what the complaint of R 2.2.0 on MS Windows is, but I guess it is the same as under linux: DEVw=500 DEVh=350 fname=my_plot plot(rnorm(300)) FOO - png dev.print(file=fname, device=FOO, width=DEVw, height=DEVh, bg=transparent) Error in dev.copy(device = png, file = fname, width = DEVw, height = DEVh, : 'device' should be a function Which is very informative. `device' is supposed to be a function, not a character variable, thus: FOO - png dev.print(file=fname, device=FOO, width=DEVw, height=DEVh, bg=transparent) X11 2 FOO - pdf dev.print(file=fname, device=FOO, width=DEVw, height=DEVh, bg=transparent) X11 2 all seem to work. HTH. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] How to simulate correlated data
G'day Lisa, LW == Lisa Wang [EMAIL PROTECTED] writes: LW I would like to simulate X --Normal (20, 5) Y-- Normal (40, LW 10) and the correlation between X and Y is 0.6. LW How do I do it in R? That depends on what you want the joint distribution to be. :) If you want the joint distribution to be normal, you could use the function mvrnorm() from the MASS package. HTH. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] R is GNU S, not C.... [was how to get or store .....]
vic == vincent [EMAIL PROTECTED] writes: vic ronggui a __ 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] apply() and dropped dimensions
G'day Robin, RH == Robin Hankin [EMAIL PROTECTED] writes: RH How do I rewrite jj() so that it consistently returns a RH matrix? How about explicitly returning a matrix with the desired dimensions? jj function(m1,m2,f,...) matrix(apply(m1, 1, function(y) apply(m2, 1, function(x) f(x, y, ...) )), nrow=nrow(m2), ncol=nrow(m1)) jj(matrix(1:20,4,5),matrix(1:10,2,5),f=sum) [,1] [,2] [,3] [,4] [1,] 70 75 80 85 [2,] 75 80 85 90 jj(matrix(1:20,4,5),matrix(1:5,1,5),f=sum) [,1] [,2] [,3] [,4] [1,] 60 65 70 75 Or, using the outer command: jj function(m1,m2,f,...) outer(1:nrow(m2), 1:nrow(m1), function(i,j) apply(cbind(i,j), 1, function(ii) f(m2[ii[1],], m1[ii[2],], ...))) jj(matrix(1:20,4,5),matrix(1:10,2,5),f=sum) [,1] [,2] [,3] [,4] [1,] 70 75 80 85 [2,] 75 80 85 90 jj(matrix(1:20,4,5),matrix(1:5,1,5),f=sum) [,1] [,2] [,3] [,4] [1,] 60 65 70 75 Cheers, Berwin __ 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] extracting rows of a dataframe
RH == Robin Hankin [EMAIL PROTECTED] writes: RH Three questions: RH (1) why is as.vector(a[2,-1]) not a vector? Did you read the help file of as.vector() ? a - data.frame(male=c(T,T,F),mass=c(1,4,3),height=c(4,3,2)) a male mass height 1 TRUE1 4 2 TRUE4 3 3 FALSE3 2 x - as.vector(a[2,-1]) mode(x) [1] list str(x) `data.frame': 1 obs. of 2 variables: $ mass : num 4 $ height: num 3 Clearly you want: x - as.vector(a[2,-1], mode=numeric) str(x) num [1:2] 4 3 x [1] 4 3 RH (2) How come setting names to NULL gives me bad weirdness? names(x) - NULL str(x) `data.frame': 1 obs. of 2 variables: $ : num 4 $ : num 3 x structure(4, class = AsIs) structure(3, class = AsIs) 2 4 3 With your names(x)-NULL command you are wiping out the names of the variables in your data frame. Looking at print.data.frame(), the answer to the so-called weirdness can presumably be found in format.data.frame(). RH (3) Can I structure my dataframe in a better way, so that RH this problem does not occur? Not really, it's more a question of the appropriate use of as.vector(). :-)) Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] floor()
G'day Werner, WB == Werner Bier [EMAIL PROTECTED] writes: floor((5.05-floor(5))*100) WB [1] 4 WB I would expect 5, or am I wrong? You are wrong. :) Consider: (5.05-floor(5))*100 [1] 5 (5.05-floor(5))*100 - 5 [1] -1.776357e-14 and read FAQ 7.31 Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] force apply() to return a list
G'day Robin,
RH == Robin Hankin [EMAIL PROTECTED] writes:
RH Still, it's a little disconcerting that apply() as generally
RH used can return a list 99.9% of the time and a matrix the
RH other 0.1%.
It probably depends on how it is used. :-)
I usually use apply in situations where I know that the lengths of the
result is always the same and want a matrix returned. Thus, in my
applications in 99.99% of the time a matrix is returned and the other
0.01% a list (and I forget that apply could return a list).
(BTW, 67.8% of all statistics are made up on the spot, just as the
last three presumably. :) )
RH It took me a long time to track this down when debugging my
RH code the other night.
It is called defensive programming (and it is great if it works). :-)
If your code expects a list, make sure that it gets a list and spits a
dummy if not, then these kind of problems are easy to find.
RH Does anyone else agree? Would it be possible to add an
RH argument such as force.list (with default FALSE) to apply()
RH that makes apply() return a list under all circumstances?
RH Also, adding Dimitris's excellent workaround to apply()'s
RH manpage would be helpful.
If you supply appropriate patches against the SVN source, I am sure
that R core will consider it.
And although Dimitris solution is quite nice, I find that it obfuscate
the code a bit. What is wrong with checking whether the result
returned by apply is a matrix, and if so turn the matrix into a list.
Probably the easiest way to turn a matrix into a list is to use
as.data.frame(), if you want that the printed version of the object
looks like a list, then use as.list() too:
a2 - cbind(1:3,3:1)
f - function(x){1:max(x[1],4)}
apply(a2,1,f)
[,1] [,2] [,3]
[1,]111
[2,]222
[3,]333
[4,]444
res - apply(a2,1,f)
if(is.matrix(res)) res - as.list(as.data.frame(res))
res
$V1
[1] 1 2 3 4
$V2
[1] 1 2 3 4
$V3
[1] 1 2 3 4
Cheers,
Berwin
__
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] A 'sweave' strange problem !!!
G'day Stéphane, SD == Stéphane Dray [EMAIL PROTECTED] writes: SD I have found a problem (bug?) with Sweave. [...] The value SD of 'a' has change between the two Schunks. It seems that the SD problem only appear when there are plot (fig=T) in the first SD one. Without plot, there are no problems: a remains SD unchanged. Is it a bug or have I misundertsood something ? You have misunderstood something :) But it is quite subtle and it took me some time to realise what was going on. Note, that if an Sweave chunk produces a figure, then by default a PDF and an EPS version of the figure is produced. But R can produce only one figure at a time, thus the chunk will be executed at least twice if you set 'fig=TRUE' (and do not change any of the other arguments). I noticed that even if I add 'eps=FALSE' or 'pdf=FALSE', the value of 'a' changed, if 'fig=TRUE'. Only with 'fig=TRUE' and 'eps=FALSE' and 'pdf=FALSE' (don't ask me why I tried it), did the value of 'a' not change. Hence, my guess is that chunks that have 'fig=TRUE' are executed once to produce the output for the .tex file and then they are executed again to produce EPS and/or PDF output. Thus such a chunk is executed once, twice or thrice; depending on the settings of 'eps'and 'pdf'. Thus, it is not a good idea to have statements in such chunks that produce (pseudo-)random results. SD Thanks a lot ! My pleasure. HTH. SD R version 2.1.0 (Debian). Well, I guess the standard on this mailing list is to point out that this is quite an old version of R and that the current one is R 2.2.0 (but that one has the same behaviour). :) Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] Optim with two constraints
AD == Alexis Diamond [EMAIL PROTECTED] writes: AD I have a follow-up from Jens's question and Professor Ripley's AD response. AD Jens wants to do quadratic optimization with 2 constraints: # I need two constraints: # 1. each element in par needs to be between 0 and 1 # 2. sum(par)=1, i.e. the elements in par need to sum to 1 AD how does one set both constraints in quadprog, per AD Prof. Ripley's suggestion? i know how to get quadprog to AD handle the second constraint, The first is actually not one constraint but 2*k constraints, where k is the number of elements in par. But there is quite some redundancy in this set of equation. It suffices to constrain each element to be bigger or equal to 0 and that they should sum to 1. If these constraints are fulfilled, then each element must be less or equal to one. AD but not BOTH, since quadprog only takes as inputs the AD constraint matrix A and constraint vector b-- So what stops you from coding the k constraints from 1.) in the form that quadprog requires them? From memory, i.e. untested: m - length(par) Amat - cbind(rep(1,m), diag(m)) bvec - c(1,rep(0,m)) meq - 1 solve.QP(Dmat, dvec, Amat, bvec, meq) AD unlike in ipop (kernlab), there is no additional option for AD box constraints. Well, the problems that I had (and still have) usually don't involve box constraints, but I see that other people use them again and again. So probably it would be a good idea to implement them... But, more importantly, would be to implement Powell's modifications of the Goldfarb-Idnani algorithm to make it numerically more robust... Oh, yeah, and a warm start option from a feasible point would be nice too Probably all in a future version which should be released sometime before Xmas 20xx. :) AD apologies if i am not seeing something obvious here. Apologies accepted. :) Cheers, Berwin __ 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] Doubt about Sweave
G'day Ronaldo, RRJ == Ronaldo Reis-Jr [EMAIL PROTECTED] writes: RRJ I reading about Sweave and it make a good output. But all RRJ example is made with R commands mades in a file. Is possible RRJ to make an output with sweave interactively in R and after RRJ the analysis end export the latex code to a file? As far as I know, no. But you may want to look at Peter Wolf's package 'relax'; see: http://www.wiwi.uni-bielefeld.de/~wolf/software/relax/relax.html That package may have the functionality that you are looking for. Cheers, Berwin __ 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] R binaries
G'day Brian, BDR == Prof Brian Ripley [EMAIL PROTECTED] writes: BDR On Wed, 31 Aug 2005, Berwin A Turlach wrote: available.packages() does not seem to have a type argument according to its documentation, so I guess that even if it is run under Windows that it returns a list of all source packages available in contrib. BDR Depends on what contriburl is set to, but the default under BDR Windows is binary packages. See my article in the current BDR R-Newsletter. BDR The default argument is contrb.url(getOption(repos)), and BDR contrib.url does have a type argument (and its default is BDR getOption(pkgType) ). Thanks for pointing this out; seems as if I didn't read the documentation in sufficient details. Thus, since Nam-Ky works on a Linux box (private e-mail) the commands should be options(repos=c(CRAN=http://cran.au.r-project.org/;), + pkgType=source) download.packages(available.packages(), destdir=.) to download all the source files; followed by options(pkgType=win.binary) download.packages(available.packages(), destdir=.) to download all contributed binary packages for Windows. Cheers, Berwin __ 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] R binaries
G'day Nam-Ky, NKN == Nam-Ky Nguyen [EMAIL PROTECTED] writes: NKN I intend to burn some R CDs to colleagues in Vietnam. I want NKN to put all binary files for base as well as contributed NKN packages (for both Windows and Linux). I don't believe that precompiled binary files for Linux exist, since they would depend very much on which flavour of Linux you are running. NKN It is very time consuming if I download files by files. Is NKN there a place a can buy a CD with all the mentioned files? Not that I am aware off. But I would do the following for easy download of all contributed files: Issue from an R session the following commands (mostly untested): options(repos=c(CRAN=http://cran.au.r-project.org/;)) tmp - available.packages() # Or just use CRAN.packages() ?? ## To download all source files download.packages(tmp, destdir=.) ## To download all windows binary download.packages(tmp, destdir=., type=win.binary) If you do it from a Windows machine, you may have to adapt the destdir argument in those commands. Also, if you sit behind a proxy, then you might have to first issue an appropriate 'Sys.putenv(http_proxy=put your proxy here)' command. available.packages() does not seem to have a type argument according to its documentation, so I guess that even if it is run under Windows that it returns a list of all source packages available in contrib. If some of those are not available as a precompiled windows binary, the download.packages() command will probably fail with an error, so you have to do the download in bits and pieces, removing offending entries from `tmp' -- but that should be still faster than downloading file by file (presumably by clicking in a web browser?). (The development version of R should now be more robust (Thanks again, Brian) when downloading and just continue downloading the remaining files if it encounters an error. But I don't know where you run (or what to run) R-devel ;-) ). Hope this helps. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] lm.ridge
G'day Daniel,
DR == daniel [EMAIL PROTECTED] writes:
DR First: I think coefficients from lm(Employed~.,data=longley)
DR should be equal coefficients from
DR lm.ridge(Employed~.,data=longley, lambda=0) why it does not
DR happen?
Which version of R and which version of MASS are you using?
lm(Employed~.,data=longley)
Call:
lm(formula = Employed ~ ., data = longley)
Coefficients:
(Intercept) GNP.deflator GNPUnemployed Armed.Forces
-3.482e+03 1.506e-02-3.582e-02-2.020e-02-1.033e-02
Population Year
-5.110e-02 1.829e+00
lm.ridge(Employed~.,data=longley, lambda=0)
GNP.deflator GNPUnemployed Armed.Forces
-3.482259e+03 1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
Population Year
-5.110411e-02 1.829151e+00
These coefficients look pretty identical to me, except that they are
printed to different numbers of significant digits.
In fact, the following shows that they are identical (upto numerical
precision):
fm1 - lm(Employed~.,data=longley)
fm2 - lm.ridge(Employed~.,data=longley, lambda=0)
coef2 - print(fm2)
GNP.deflator GNPUnemployed Armed.Forces
-3.482259e+03 1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
Population Year
-5.110411e-02 1.829151e+00
max(abs(coef(fm1)-coef2))
[1] 7.275958e-12
DR Second: if I have for example Ridge-lm.ridge(Employed~.,
DR data=longley, lambda = seq(0,0.1,0.001)), I suppose intercept
DR coefficient is defined implicit,
Yes.
DR why it does not appear in Ridge$coef?
If you look at the code of lm.ridge, you will see that, if an
intercept is included in the model, all non-constant regressors are
centered (i.e. made orthogonal to the intercept term) and scaled to
have the same variance. Further more, the intercept term is typically
*not* penalised. The components in Ridge$coef are the coefficients on
this transformed scale. No need of including the intercept here,
since it is the same for all values of lambda. If you print the
model, then the ridge coefficients on the original scale are
calculated, see:
getAnywhere(print.ridgelm)
A single object matching 'print.ridgelm' was found
It was found in the following places
registered S3 method for print from namespace MASS
namespace:MASS
with value
function (x, ...)
{
scaledcoef - t(as.matrix(x$coef/x$scales))
if (x$Inter) {
inter - x$ym - scaledcoef %*% x$xm
scaledcoef - cbind(Intercept = inter, scaledcoef)
}
print(drop(scaledcoef), ...)
}
environment: namespace:MASS
DR Third: I suppose that if I define
DR 1) y-longley$Employed
DR 2) X-as.matrix(cbind(1,Longley[,1:6])
DR 3) I = identity matrix the
DR following should be true: Coef=(X'X+kI)^(-1) X'y
No, as noted above, the intercept term is usually not penalised.
DR and if a take k=Ridge$kHKV, Coef should be approx equal to
DR Ridge$Coef[near value of kHKV]
No, as noted above the estimates in the coef component of an object
returned by lm.ridge are the coefficients on a different scale.
DR and it does not seem to happen, why?
Because the intercept is not penalised by lm.ridge and the
non-constant columns of the design matrix are rescaled; hence the
returned coefficients are on another scale.
DR Any help, suggestion or orientation?
HTH.
Cheers,
Berwin
__
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] Kronecker matrix product
RH == Robin Hankin [EMAIL PROTECTED] writes: RH I want to write a little function that takes a matrix X of RH size m-by-n, and a list L of length m, whose elements are RH matrices all of which have the same number of columns but RH possibly a different number of rows. RH I then want to get a sort of dumbed-down kronecker product in which RH X[i,j] is replaced by X[i,j]*L[[j]] RH where L[[j]] is the j-th of the m matrices. For example, if RH X = matrix(c(1,5,0,2),2,2) RH and RH L[[1]] = matrix(1:4,2,2) RH L[[2]] = matrix(c(1,1,1,1,1,10),ncol=2) RH I want RH [,1] [,2] [,3] [,4] RH [1,]1300 RH [2,]2400 RH [3,]5522 RH [4,]5522 RH [5,]5 502 20 tmp - sapply(1:length(L), function(j, mat, list) kronecker(X[j,,drop=FALSE], L[[j]]), mat=X, list=L) do.call(rbind, tmp) [,1] [,2] [,3] [,4] [1,]1300 [2,]2400 [3,]5522 [4,]5522 [5,]5 502 20 HTH. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] Kronecker matrix product
BT == Berwin A Turlach [EMAIL PROTECTED] writes: tmp - sapply(1:length(L), function(j, mat, list) kronecker(X[j,,drop=FALSE], L[[j]]), mat=X, list=L) Uups, should proof read more carefully before hitting the send button. This should be, of course: tmp - sapply(1:length(L), function(j, mat, list) kronecker(mat[j,,drop=FALSE], list[[j]]), mat=X, list=L) Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ 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] nls(): Levenberg-Marquardt, Gauss-Newton, plinear - PI curve fitting
G'day Chris,
CK == Christfried Kunath [EMAIL PROTECTED] writes:
CK With the nls()-function i want to fit following formula
CK whereas a,b, and c are variables: y~1/(a*x^2+b*x+c)
CK [...]
CK The algorithm plinear give me following error:
The algorithm plinear is inappropriate for your data and your model
since none of the parameters are linear. You are actually trying to
fit the model
y ~ d/(a*x^2+b*x+c)
where `d' would be the linear parameter.
CK phi function(x,y) {
CK
k.nls-nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.0005,b=0.02,c=1.5),alg=plinear)
CK coef(k.nls) }
CK [...]
CK The commercial software Origin V.6.1 solved this problem
CK with the Levenberg-Marquardt algorithm how i want. The
CK reference results are: a = 9.6899E-6, b = 0.00689, c = 2.72982
CK What are the right way or algorithm for me to solve this
CK problem and what means this error with alg=plinear?
The error means that at some point along the way a matrix was
calculated that needed to be inverted but was for all practical
purposes singular. This can happen in numerical optimisation
problems, in particular if derivatives have to be calculated
numerically.
How to solve this problem:
1) Don't use the algorithm plinear since it is inappropriate for
your model.
2) You may want to specify the gradient of the function that you are
minimising to make life easier for nls(), see Venables Ripley
(2002, page 215) for an example.
3) You can call nls() directly without specifying the plinear option:
(I renamed the variables in the data frame to x and y for
simplicity)
nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.0005,b=0.02,c=1.5),data=k)
Nonlinear regression model
model: y ~ 1/(a * (x^2) + b * x + c)
data: k
a b c
-7.326117e-05 4.770514e-02 2.490643e+00
residual sum-of-squares: 0.1120086
But the results seem to be highly depended on your starting values:
nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.5,b=0.002,c=2.5),data=k)
Nonlinear regression model
model: y ~ 1/(a * (x^2) + b * x + c)
data: k
abc
9.690204e-06 6.885570e-03 2.729825e+00
residual sum-of-squares: 0.000547369
Which is of some concern.
4) If you really want to fit the above model, you may also consider to
just use the glm() command and fit it within a generalised linear
model framework:
glm(y ~ I(x^2) + x, data=k, family=gaussian(link=inverse))
Call: glm(formula = y ~ I(x^2) + x, family = gaussian(link = inverse),
data = k)
Coefficients:
(Intercept) I(x^2)x
2.730e+009.690e-066.886e-03
Degrees of Freedom: 8 Total (i.e. Null); 6 Residual
Null Deviance:0.09894
Residual Deviance: 0.0005474 AIC: -53.83
HTH.
Cheers,
Berwin
== Full address
Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr)
School of Mathematics and Statistics+61 (8) 6488 3383 (self)
The University of Western Australia FAX : +61 (8) 6488 1028
35 Stirling Highway
Crawley WA 6009e-mail: [EMAIL PROTECTED]
Australiahttp://www.maths.uwa.edu.au/~berwin
__
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] problem with model.matrix
BN == Balasubramanian Narasimhan [EMAIL PROTECTED] writes: BN Trevor's problem is reproducible on R 1.8.1 (which is what he BN was using) on both Linux RHEL 3.0 and Solaris and R 1.8.0 on BN Windows. Yep, on my Linux box (details below) I have the same problem with R 1.8.1 but no problem under R 1.9.0 or R 1.9.1. Time to upgrade I would say. :) Cheers, Berwin --please do not edit the information below-- Version: platform = i686-pc-linux-gnu arch = i686 os = linux-gnu system = i686, linux-gnu status = major = 1 minor = 8.1 year = 2003 month = 11 day = 21 language = R Search Path: .GlobalEnv, package:methods, package:ctest, package:mva, package:modreg, package:nls, package:ts, Autoloads, package:base __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] problem with update.packages()
G'day Stefano, SC == Stefano Calza [EMAIL PROTECTED] writes: SC Hi. I'm experiencing a problem with updating packages on R SC 1.8.1 (2003-11-21) on Debian testing. I get the following SC message when updating for example Design: SC ... ... /usr/bin/ld: cannot find -lg2c-pic ... [...] SC Any suggestion? I ran into a similar problem some time ago. I had no problems with my original g77-3.3 from testing. But then g77-3.3 was upgraded and when I wanted to compile a package that contained Fortran code, I ran into this problem. A bit of investigation showed that the current g77-3.3 version in testing has (on my machine) in /usr/lib/gcc-lib/i486-linux/3.3.3 symbolic links for libg2c-pic.a and libg2c-pic.so to libg2c.a and libg2c.so, respectively, as if the latter libraries are in the same directory. But those two libraries are actually installed into /usr/lib/. So the symbolic links point into nirvana and ld can't find the libraries. My solution was to upgrade to g77-3.3 in unstable (by changing my preference file in /etc/apt (attached below). No problems since then. Hope that helps. Cheers, Berwin Package: * Pin: release a=stable Pin-Priority: 500 Package: * Pin: release a=testing Pin-Priority: 600 Package: * Pin: release a=unstable Pin-Priority: 50 Package: g77-3.3 Pin: release a=unstable Pin-Priority: 999 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Odd behaviour of step (and stepAIC)?
JO == Jari Oksanen [EMAIL PROTECTED] writes: JO There seems to be some obscure features in step() when you JO have interaction terms: their interpretation is order JO sensitive. [...] Just because R internally decides to order JO terms differently than in the scope (this may happen even when JO you have produced the scope by first fitting the maximal JO model, and then extracting that using scope=formula()). I once JO was diligent enough to trace this to a certain C function JO where this ordering was done, but I was not clever enough to JO instantly know how to change the behaviour. It is well possible that this may be occasionally a problem. (I remember that a colleague of mine once showed me a problem were the behaviour of step depended on how you specified the code. However, I believe more recent version of R don't show that behaviour anymore.) But I don't believe that this is the problem in Bob's case: I'm getting the following from a stepwise selection (with both step and stepAIC): step(lm(sqrt(Grids)~ SE + Edge + NH), scope=~ (Edge + SE + NH)^2) Start: AIC= 593.56 sqrt(Grids) ~ SE + Edge + NH Df Sum of Sq RSSAIC none 2147.0 593.6 + Edge:NH 1 3.0 2143.9 595.1 + SE:NH4 23.2 2123.8 598.4 - NH 1 75.8 .8 601.6 - Edge 1 448.7 2595.7 646.4 - SE 41033.7 3180.6 699.1 My problem is that the SE:Edge term is not added. If you look at the output, adding the SE:Edge term decreases RSS but increases AIC. The increase in parameters is not worth the decrease in RSS, according to the AIC criterion. step tries to find the best AIC model, i.e., the current model is the best. Cheers, Berwin == Full address Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics+61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009e-mail: [EMAIL PROTECTED] Australiahttp://www.maths.uwa.edu.au/~berwin __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Probit predictions outside (0,1) interval
AM == Arnab mukherji [EMAIL PROTECTED] writes: AM Any advice would be really helpful. Read the documentation of predict and then the one of predict.glm? ;-) I guess you actually wanted to do one of the following: yhat-predict(g, dat, type=response) range(yhat) [1] 0.2760238 0.9229622 or, yhat-fitted(g) range(yhat) [1] 0.2760238 0.9229622 Cheers, Berwin __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
