Re: [R] How to convert c:\a\b to c:/a/b
On Wed, 29 Jun 2005, David Duffy wrote: I couldn't resist adding a more literal answer This can only work for escapes which are preserved. The parser maps \n to a character (LF) and the deparser maps it back to \n. This happens to be true of \a \b \f \n \r \t \v \\ but no others. For example, \s is mapped to s, and there is no difference between \s and s in the parsed input. unback - function(x) { chars - unlist(strsplit(deparse(x),)) chars - chars[-c(1,length(chars))] paste(gsub(,/,chars),collapse=) } unback(\n) unback(\s) [1] s Spencer Graves keeps on insisting there is a better way, but all the solutions are to avoid sending the string to the parser, and hence avoiding having the string directly in an R script. This is common in shell scripts, which use 'here' documents to avoid 'quoting hell'. We can do that in R too. Here are two variants I have not seen in the thread test1.R: scan(, , allowEscapes=FALSE, n=1, quiet=TRUE) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt catIx, \n, sep=) R --slave --vanilla test1.R D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt (This one does not allow quoted strings.) test2.R: x - readLines(stdin(), n=1) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt x - gsub('^(.*)$', \\1, x) cat(x, \n) R --slave --vanilla test2.R D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt (This one allows surrounding quotes or not.) -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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] where can i download the metrics package?
There are several packages within rmetrics such as fSeries, fBasics, fExtremes, and so on. You can download those in the usual way. From: Spencer Graves [EMAIL PROTECTED] To: ronggui [EMAIL PROTECTED] CC: R [EMAIL PROTECTED] Subject: Re: [R] where can i download the metrics package? Date: Tue, 28 Jun 2005 08:27:47 -0700 How about Rmetrics, listed under misc: related projects on www.r-project.org. spencer graves ronggui wrote: i learn it from metrics: Towards a package for doing econometrics in Rbut i can not find it in cran. -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] quick way to construct formula
Here is a way: it uses 'paste' but I dont think it is a problem anyway to use it. Nevertheless, it is surely a bad idea to fit any model with more than 25 terms... main_effects = paste(nam,collapse=+) inter - outer(nam,nam,paste,sep=:) inter - paste(inter[upper.tri(inter)],collapse=+) log_effects - paste(log(,nam,),sep=,collapse=+) as.formula(paste(~,main_effects,+,inter,+,log_effects,sep=)) ~x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x1:x2 + x1:x3 + x2:x3 + x1:x4 + x2:x4 + x3:x4 + x1:x5 + x2:x5 + x3:x5 + x4:x5 + x1:x6 + x2:x6 + x3:x6 + x4:x6 + x5:x6 + x1:x7 + x2:x7 + x3:x7 + x4:x7 + x5:x7 + x6:x7 + x1:x8 + x2:x8 + x3:x8 + x4:x8 + x5:x8 + x6:x8 + x7:x8 + x1:x9 + x2:x9 + x3:x9 + x4:x9 + x5:x9 + x6:x9 + x7:x9 + x8:x9 + x1:x10 + x2:x10 + x3:x10 + x4:x10 + x5:x10 + x6:x10 + x7:x10 + x8:x10 + x9:x10 + log(x1) + log(x2) + log(x3) + log(x4) + log(x5) + log(x6) + log(x7) + log(x8) + log(x9) + log(x10) HTH, Eric Eric Lecoutre UCL / Institut de Statistique Voie du Roman Pays, 20 1348 Louvain-la-Neuve Belgium tel: (+32)(0)10473050 [EMAIL PROTECTED] http://www.stat.ucl.ac.be/ISpersonnel/lecoutre If the statistics are boring, then you've got the wrong numbers. -Edward Tufte -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Luke Sent: mercredi 29 juin 2005 4:49 To: R-help@stat.math.ethz.ch Subject: [R] quick way to construct formula Dear R users, I have a data with 1000 variables named x1, x2, ..., x1000, and I want to construct a formula like this format: ~x1+x2+...+x1000+x1:x2+x1:x3+x999:x1000+log(x1)+...+log(x1000) That is: the base variables followed by all interaction terms and all base feature log-transformations. I know I can use several paste functions to construct it. But is there any other handy way to do it? -Luke __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] enhanced MDS
Try also ?sphpca (library(psy) -- new version 0.7) Best Bruno Bruno Falissard INSERM U669, PSIGIAM Paris Sud Innovation Group in Adolescent Mental Health Maison de Solenn 97 Boulevard de Port Royal 75679 Paris cedex 14, France tel : (+33) 6 81 82 70 76 fax : (+33) 1 45 59 34 18 web site : http://perso.wanadoo.fr/bruno.falissard/ -Message d'origine- De : [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] De la part de Karen Kotschy Envoyé : mardi 28 juin 2005 13:39 À : r-help@stat.math.ethz.ch Objet : [R] enhanced MDS Hi again Sorry, in looking again at sammon and isoMDS I see that they seem to do exactly what I want, except that they are non-metric, which means, as I understand it, that they relate the rank orders of the variables rather than the actual distances. Could I use these non-metric MDS packages even if my distances are metric? Thanks Karen -- Karen Kotschy Centre for Water in the Environment University of the Witwatersrand Johannesburg P/Bag X3, Wits, 2050 Tel: +2711 717-6425 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] quick way to construct formula
On Tue, 28 Jun 2005, Luke wrote: Dear R users, I have a data with 1000 variables named x1, x2, ..., x1000, and I want to construct a formula like this format: ~x1+x2+...+x1000+x1:x2+x1:x3+x999:x1000+log(x1)+...+log(x1000) That is: the base variables followed by all interaction terms and all base feature log-transformations. I know I can use several paste functions to construct it. But is there any other handy way to do it? Do you really want a formula with over half a million terms in? I think it is likely that R will hit recursion limits in handing such a formula, quite possibly C-level stack overflow. For a more modest example: dd - data.frame(x1=1, x2=2, x3=3) terms(~ . + .^2, data=dd, simplify=TRUE) will get you all terms and their interactions. Using paste() to get the log terms is the easiest way I know. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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] enhanced multidimensional scaling?
Karen == Karen Kotschy [EMAIL PROTECTED] on Tue, 28 Jun 2005 13:13:31 +0200 writes: Karen Dear R list Would anyone be able to tell me whether Karen it is possible to do enhanced multidimensional Karen scaling (enhanced MDS) in R? In other words, Karen something that goes beyond cmdscale by iteratively Karen improving the fit between observed dissimilarities Karen and inter-object distances, using the KYST algorithm Karen (Kruskal, 1964). Since you know about cmdscale(), do read it's help page more carefully, in particular the section See Also: {which, in an ESS help buffer, is reachable by simple s s The help page of the first function mentioned there is entitled `` Kruskal's Non-metric Multidimensional Scaling '' Karen I have found several implementations of non-metric Karen MDS in various packages but nothing like what I have Karen described above. well, you've mentioned Kruskal(1964), but not described much more, hence I hope the above hint does help. Regards, Martin Maechler, ETH Zurich __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] comparison of packages for Unit Root test
Dear R Users, Could somebody please compare the packages of unit root test ( Uroot, Ucra, tseries and fseries ) regarding the type of test ( without constant and trend, with constant , and with constant and trend ) ? Regards, Amir Safari - Rekindle the Rivalries. Sign up for Fantasy Football [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Running SVM {e1071}
Dear David, Dear Friends, After any running svm I receive different results of Error estimation of 'svm' using 10-fold cross validation. What is the reason ? It is caused by the algorithm, libsvm , e1071 or something els? Which value can be optimal one ? How much run can reach to the optimality.And finally, what is difference between Error estimation of svm using 10-fold cross validation and MSE ( Mean Square Error ) ? So many thanks in advance for your help. Best Regards, Amir __ [[alternative HTML version deleted]] __ 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] comparison of packages for Unit Root test
On Wed, 29 Jun 2005 03:08:23 -0700 (PDT) Amir Safari wrote: Dear R Users, Could somebody please compare the packages of unit root test ( Uroot, Ucra, tseries and fseries ) regarding the type of test ( without constant and trend, with constant , and with constant and trend ) ? IMHO, the tseries implementation is still the reference implementation for the significance tests which also have a slim interface. urca (sic!) offers some more tools in an integrated way. If you are not just interested in one test or the other, urca certainly offers much tools in a unit root analysis. If you use the Rmetrics suite of packages, then fSeries (sic!) is probably most useful. uroot (sic!) offers a GUI which is a bit cumbersome to set up and the output is a bit crude, e.g., doesn't offer p-values. my 2 cent, Z Regards, Amir Safari - Rekindle the Rivalries. Sign up for Fantasy Football [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] Running SVM {e1071}
Amir Safari wrote: Dear David, Dear Friends, After any running svm I receive different results of Error estimation of 'svm' using 10-fold cross validation. What is the reason ? It is caused by the algorithm, libsvm , e1071 or something els? Which value can be optimal one ? How much run can reach to the optimality.And finally, what is difference between Error estimation of svm using 10-fold cross validation and MSE ( Mean Square Error ) ? So many thanks in advance for your help. 10-fold cross validation chooses the traing/test sets randomly, hence you won't get the same results each time you run it ... Uwe Ligges Best Regards, Amir __ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] svm and scaling input
[EMAIL PROTECTED] wrote: Dear All, I've a question about scaling the input variables for an analysis with svm (package e1071). Most of my variables are factors with 4 to 6 levels but there are also some numeric variables. I'm not familiar with the math behind svms, so my assumtions maybe completely wrong ... or obvious. Will the svm automatically expand the factors into a binary matrix? yes. If I add numeric variables outside the range of 0 to 1 do I have to scale them to have 0 to 1 range? svm() will scale your data by default. Cheers, David -- Dr. David Meyer Department of Information Systems and Operations Vienna University of Economics and Business Administration Augasse 2-6, A-1090 Wien, Austria, Europe Fax: +43-1-313 36x746 Tel: +43-1-313 36x4393 HP: http://wi.wu-wien.ac.at/~meyer/ __ 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] Running SVM {e1071}
Dear David, Dear Friends, After any running svm I receive different results of Error estimation of 'svm' using 10-fold cross validation. using tune.svm(), or the `cross' parameter of svm()? What is the reason ? It is caused by the algorithm, libsvm , e1071 or something els? The splits are chosen randomly. Which value can be optimal one? The Bayes Error. How much run can reach to the optimality. What do you mean by `How much run'? And finally, what is difference between Error estimation of svm using 10-fold cross validation and MSE ( Mean Square Error ) ? the former is an error estimation _procedure_, the latter is an error _measure. Cheers, David __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Memory Management under Linux: Problems to allocate large amounts of data
Dear Group I'm still trying to bring many data into R (see older postings). After solving some troubles with the database I do most of the work in MySQL. But still I could be nice to work on some data using R. Therefore I can use a dedicated Server with Gentoo Linux as OS hosting only R. This Server is a nice machine with two CPU and 4GB RAM which should do the job: Dual Intel XEON 3.06 GHz 4 x 1 GB RAM PC2100 CL2 HP Proliant DL380-G3 I read the R-Online help on memory issues and the article on garbage collection from the R-News 01-2001 (Luke Tierney). Also the FAQ and some newsgroup postings were very helpful on understanding memory issues using R. Now I try to read data from a database. The data I wanted to read consists of 158902553 rows and one field (column) and is of type bigint(20) in the database. I received the message that R could not allocate the 2048000 Kb (almost 2GB) sized vector. As I have 4BG of RAM I could not imagine why this happened. In my understanding R under Linux (32bit) should be able to use the full RAM. As there is not much space used by OS and R as such (free shows the use of app. 670 MB after dbSendQuery and fetch) there are 3GB to be occupied by R. Is that correct? After that I started R by setting n/vsize explicitly R --min-vsize=10M --max-vsize=3G --min-nsize=500k --max-nsize=100M mem.limits() nsize vsize 104857600NA and received the same message. A garbage collection delivered the following information: gc() used (Mb) gc trigger (Mb) limit (Mb) max used (Mb) Ncells 217234 5.9 50 13.4 280050 13.4 Vcells 87472 0.7 157650064 1202.8 3072 196695437 1500.7 Now I'm at a loss. Maybe anyone could give me a hint where I should read further or which Information can take me any further Dubravko Dolic Statistical Analyst Tel: +49 (0)89-55 27 44 - 4630 Fax: +49 (0)89-55 27 44 - 2463 Email: [EMAIL PROTECTED] Komdat GmbH Nymphenburger Straße 86 80636 München - ONLINE MARKETING THAT WORKS - This electronic message contains information from Komdat Gmb...{{dropped}} __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] coxph, survfit and Brier score
Hello and apologies for a very long question. I thought it better to be verbose and clear than short and imprecise. I am trying to compute the brier score comparing true surv object of test data to predictions from train data (using sbrier in ipred package). I am having trouble getting the right format for my prediction object I think. I have split the DLBCL dataset and fit a coxph on a training set then I am trying to get out a prediction using the survfit function with newdata=test. My code is below. library(survival) library(ipred) data(DLBCL) DLBCL[complete.cases(DLBCL),]-DLBCL train-DLBCL[1:19,] test-DLBCL[20:38,] train.surv-Surv(train$time, train$cens) test.surv-Surv(test$time, test$cens) train.mod-coxph(train.surv~IPI, data=train) pred- survfit(train.mod, newdata=test) class(pred) [1] survfit.cox survfit pred Call: survfit.coxph(object = train.mod, newdata = test) n events median 0.95LCL 0.95UCL [1,] 19 10Inf27.1 Inf [2,] 19 10Inf Inf Inf [3,] 19 10Inf Inf Inf [4,] 19 10 23.7 4.1 Inf [5,] 19 10 23.7 4.1 Inf [6,] 19 10Inf27.1 Inf [7,] 19 10Inf Inf Inf [8,] 19 10Inf Inf Inf [9,] 19 10Inf27.1 Inf [10,] 19 104.1 2.9 Inf [11,] 19 10Inf27.1 Inf [12,] 19 104.1 2.9 Inf [13,] 19 10Inf Inf Inf [14,] 19 10 23.7 4.1 Inf [15,] 19 10Inf27.1 Inf [16,] 19 10Inf27.1 Inf [17,] 19 10Inf Inf Inf [18,] 19 10Inf Inf Inf [19,] 19 10Inf27.1 Inf sbrier(test.surv, pred) Error in switch(ptype, survfit = { : switch: EXPR must return a length 1 vector The sbrier function clearly does not recognize the format of this survfit object pred. Indeed it seems to have the same variables but does not look like the object you receive from a normal survfit call e.g. KM - survfit(train.surv) KM Call: survfit(formula = train.surv) n events median 0.95LCL 0.95UCL 19.010.071.315.5 Inf sbrier(train.surv, KM) integrated Brier score 0.2220228 attr(,time) [1] 2.4 102.4 My question is can I either force survfit(etc..., newdata=x) to return a useful survfit object for use with sbrier, or alternatively coerce the pred object above into a survfit object? Or am I missing something ..is there a reason survfit returns the predicted object in a different format? PS I realize the fit is not good but just want to get the code to work. Thank You Stephen Henderson ** This email and any files transmitted with it are confidentia...{{dropped}} __ 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] Producing character given i.e. | with plotmath
Hello! Does someone know how to produce L(y|mu) with plotmath? Some code with unsuccessfull results: plot(dnorm(x = seq(from = -4, to = 4, by = 0.1)), type = l) ## Not what I want legend(legend = c(expression(L(y:mu))), x = topright) ## Strange, is this a bug? legend(legend = c(expression(L(y|mu))), x = top) No, | is a logical Operator that can be rewritten in its original function form as follows: |(FALSE, TRUE) Hence the result is expected. Yes, that makes sense. Thanks! ## Group produces an error legend(legend = c(expression(group(L(y, |, mu, x = topleft) You have not specified any delimiter. Ok, I got this point wrong from ?plotmath . ## Paste keeps commas in expression legend(legend = c(expression(paste(L(y, |, mu, x = bottomleft) correct ## This one is OK, but braces are not as they should be legend(legend = c(expression(paste(L(y, |, mu, x = bottom) What's wrong with the braces?` They are not the same as in previous cases. They are not bold. What you really want is: legend(legend = c(expression(L(group(, y, |) * mu))), x = center) Yes, that is what I want. Thanks! I got additionally response from Roger D. Peng, which stated: How about legend(legend = expression(L(y | * mu * )), x = topleft) That's also OK, but note above comment about bold braces. Regards, Gregor __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] x*x*x*... vs x^n
Hi I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [so f1() and f2() and f3() are algebraically identical] a - rnorm(100) system.time(ignore - f1(a)) [1] 0.50 0.17 2.88 0.00 0.00 system.time(ignore - f2(a)) [1] 0.31 0.03 1.40 0.00 0.00 system.time(ignore - f3(a)) [1] 0.10 0.07 0.18 0.00 0.00 [these figures show little variance from trial to trial] I was expecting f2() and f3() to be about the same. I was not expecting a factor of 3 there! anyone got any comments? -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] x*x*x*... vs x^n
Something like this is exploited very nicely in the mtx.exp for matrix powers in the Malmig package, actually. Hi I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [so f1() and f2() and f3() are algebraically identical] a - rnorm(100) system.time(ignore - f1(a)) [1] 0.50 0.17 2.88 0.00 0.00 system.time(ignore - f2(a)) [1] 0.31 0.03 1.40 0.00 0.00 system.time(ignore - f3(a)) [1] 0.10 0.07 0.18 0.00 0.00 [these figures show little variance from trial to trial] I was expecting f2() and f3() to be about the same. I was not expecting a factor of 3 there! anyone got any comments? -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 Ken Knoblauch Inserm U371, Cerveau et Vision Department of Cognitive Neurosciences 18 avenue du Doyen Lepine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: 06 84 10 64 10 http://www.lyon.inserm.fr/371/ __ 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] x*x*x*... vs x^n
I tried your code and got different results: system.time(ignore - f1(a)) [1] 0.83 0.09 1.08 NA NA system.time(ignore - f2(a)) [1] 0.38 0.01 0.41 NA NA system.time(ignore - f3(a)) [1] 0.32 0.04 0.43 NA NA So I tried it again but with a loop and got: for(i in 1:10) cat(system.time(ignore - f2(a)), \n) 0.36 0.04 0.44 NA NA 0.32 0.01 0.34 NA NA 0.28 0.03 0.32 NA NA 0.29 0.03 0.35 NA NA 0.3 0.02 0.32 NA NA 0.28 0.03 0.32 NA NA 0.3 0.02 0.32 NA NA 0.29 0.02 0.34 NA NA 0.23 0.03 0.32 NA NA 0.42 0 0.45 NA NA for(i in 1:10) cat(system.time(ignore - f3(a)), \n) 0.19 0.04 0.25 NA NA 0.17 0.04 0.25 NA NA 0.21 0.02 0.25 NA NA 0.21 0.02 0.23 NA NA 0.18 0.04 0.23 NA NA 0.18 0.05 0.23 NA NA 0.18 0.04 0.25 NA NA 0.17 0.06 0.23 NA NA 0.2 0.02 0.23 NA NA 0.14 0.06 0.25 NA NA It seems to me that f3 is 50% slower than f2 not 300%. My System is: - R version: R 2.1.1 - Operating System: Win XP - Compiler: mingw32-gcc-3.4.2 Jarek \=== Jarek Tuszynski, PhD. o / \ Science Applications International Corporation \__,| (703) 676-4192 \ [EMAIL PROTECTED] `\ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Robin Hankin Sent: Wednesday, June 29, 2005 7:32 AM To: r-help Subject: [R] x*x*x*... vs x^n Hi I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [so f1() and f2() and f3() are algebraically identical] a - rnorm(100) system.time(ignore - f1(a)) [1] 0.50 0.17 2.88 0.00 0.00 system.time(ignore - f2(a)) [1] 0.31 0.03 1.40 0.00 0.00 system.time(ignore - f3(a)) [1] 0.10 0.07 0.18 0.00 0.00 [these figures show little variance from trial to trial] I was expecting f2() and f3() to be about the same. I was not expecting a factor of 3 there! anyone got any comments? -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] x*x*x*... vs x^n
On 6/29/2005 7:32 AM, Robin Hankin wrote: Hi I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [so f1() and f2() and f3() are algebraically identical] a - rnorm(100) system.time(ignore - f1(a)) [1] 0.50 0.17 2.88 0.00 0.00 system.time(ignore - f2(a)) [1] 0.31 0.03 1.40 0.00 0.00 system.time(ignore - f3(a)) [1] 0.10 0.07 0.18 0.00 0.00 [these figures show little variance from trial to trial] I was expecting f2() and f3() to be about the same. I was not expecting a factor of 3 there! anyone got any comments? If you look in src/main/arithmetic.c, you'll see that R does a general real-valued power (using C's pow() function) whenever either one of the args is real (except for a few special cases, e.g. non-numbers, or powers of 2 or 0.5). There is an internal R function equivalent to your f3, but it is not used in the situation of real^integer (and in any case, x^8 is real^real). I suppose if you wanted to submit a patch someone would take a look, but the question is whether there is really any calculation whose execution time would be materially affected by this. Most computations are not dominated by integer power calculations, so is this really worth the trouble? Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] range of input data
Hi, I've written this function: g = function(test,p1,p2) { test=sort(test) merke=0 for (z in 1:length(test)) { F1=((2*z-1)/length(test)) F21=log(plnorm(test[z],p1,p2)) F22=log(1-plnorm(test[length(test)+1-z],p1,p2)) F2= F21+F22 merke=merke+F1*F2 } return(-length(test)*-merke) } When I call it for e.g. test = 1:10 it only returns ONE value. Instead what I want to have is 10 values (like it makes the dlnorm for instance). I could make that with a loop, but maybe there is there a more simple trick ? Carsten [[alternative HTML version deleted]] __ 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] Memory Management under Linux: Problems to allocate large amounts of data
Let's assume this is a 32-bit Xeon and a 32-bit OS (there are 64-bit-capable Xeons). Then a user process like R gets a 4GB address space, 1GB of which is reserved for the kernel. So R has a 3GB address space, and it is trying to allocate a 2GB contigous chunk. Because of memory fragmentation that is quite unlikely to succeed. We run 64-bit OSes on all our machines with 2GB or more RAM, for this reason. On Wed, 29 Jun 2005, Dubravko Dolic wrote: Dear Group I'm still trying to bring many data into R (see older postings). After solving some troubles with the database I do most of the work in MySQL. But still I could be nice to work on some data using R. Therefore I can use a dedicated Server with Gentoo Linux as OS hosting only R. This Server is a nice machine with two CPU and 4GB RAM which should do the job: Dual Intel XEON 3.06 GHz 4 x 1 GB RAM PC2100 CL2 HP Proliant DL380-G3 I read the R-Online help on memory issues and the article on garbage collection from the R-News 01-2001 (Luke Tierney). Also the FAQ and some newsgroup postings were very helpful on understanding memory issues using R. Now I try to read data from a database. The data I wanted to read consists of 158902553 rows and one field (column) and is of type bigint(20) in the database. I received the message that R could not allocate the 2048000 Kb (almost 2GB) sized vector. As I have 4BG of RAM I could not imagine why this happened. In my understanding R under Linux (32bit) should be able to use the full RAM. As there is not much space used by OS and R as such (free shows the use of app. 670 MB after dbSendQuery and fetch) there are 3GB to be occupied by R. Is that correct? Not really. The R executable code and the Ncells are already in the address space, and this is a virtual memory OS, so the amount of RAM is not relevant (it would still be a 3GB limit with 12GB of RAM). After that I started R by setting n/vsize explicitly R --min-vsize=10M --max-vsize=3G --min-nsize=500k --max-nsize=100M mem.limits() nsize vsize 104857600NA and received the same message. A garbage collection delivered the following information: gc() used (Mb) gc trigger (Mb) limit (Mb) max used (Mb) Ncells 217234 5.9 50 13.4 280050 13.4 Vcells 87472 0.7 157650064 1202.8 3072 196695437 1500.7 Now I'm at a loss. Maybe anyone could give me a hint where I should read further or which Information can take me any further -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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] x*x*x*... vs x^n
Hi Duncan On Jun 29, 2005, at 02:04 pm, Duncan Murdoch wrote: On 6/29/2005 7:32 AM, Robin Hankin wrote: Hi I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [snip] If you look in src/main/arithmetic.c, you'll see that R does a general real-valued power (using C's pow() function) whenever either one of the args is real (except for a few special cases, e.g. non-numbers, or powers of 2 or 0.5). There is an internal R function equivalent to your f3, but it is not used in the situation of real^integer (and in any case, x^8 is real^real). I suppose if you wanted to submit a patch someone would take a look, but the question is whether there is really any calculation whose execution time would be materially affected by this. Most computations are not dominated by integer power calculations, so is this really worth the trouble? Duncan Murdoch well, the Gnu Scientific Library has the pow_int() functions, which are a generalization of f3(), so someone thinks so. I did a speed test of them but they were much slower than R (for any of f1(), f2(), f3()): library(gsl) system.time(ignore - pow_int(a,8)) [1] 1.07 1.11 3.08 0.00 0.00 why the slow execution time? But I'm far more interested in the philosophy behind your comments. I would say that it definitely *is* worth the trouble because someone, somewhere, will want fast integer powers, and possibly use R for nothing else. Ken's point about matrix exponentiation is relevant here too. This is a stated design consideration in Mathematica, I think. (Of course, I'm not suggesting that other programming tasks be suspended! All I'm pointing out is that there may exist a user to whom fast integer powers are very very important) very best wishes rksh -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ 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] x*x*x*... vs x^n
On Wed, 29 Jun 2005, Duncan Murdoch wrote: On 6/29/2005 7:32 AM, Robin Hankin wrote: I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [so f1() and f2() and f3() are algebraically identical] a - rnorm(100) system.time(ignore - f1(a)) [1] 0.50 0.17 2.88 0.00 0.00 system.time(ignore - f2(a)) [1] 0.31 0.03 1.40 0.00 0.00 system.time(ignore - f3(a)) [1] 0.10 0.07 0.18 0.00 0.00 [these figures show little variance from trial to trial] I was expecting f2() and f3() to be about the same. I was not expecting a factor of 3 there! anyone got any comments? If you look in src/main/arithmetic.c, you'll see that R does a general real-valued power (using C's pow() function) whenever either one of the args is real (except for a few special cases, e.g. non-numbers, or powers of 2 or 0.5). There is an internal R function equivalent to your f3, but it is not used in the situation of real^integer (and in any case, x^8 is real^real). I suppose if you wanted to submit a patch someone would take a look, but the question is whether there is really any calculation whose execution time would be materially affected by this. Most computations are not dominated by integer power calculations, so is this really worth the trouble? As Luke Tierney frequently points out, selecting special cases can take more time than you save. The assembler code used by modern OSes will have worked out that compromise pretty effectively: even for real^real it spots simple cases of the exponent. Also, it depends on the CPU: I get on Athlon 2600 system.time(ignore - f1(a), gcFirst=T) [1] 0.08 0.05 0.14 0.00 0.00 system.time(ignore - f2(a), gcFirst=T) [1] 0.20 0.01 0.20 0.00 0.00 system.time(ignore - f3(a), gcFirst=T) [1] 0.03 0.02 0.05 0.00 0.00 and Opteron 248 system.time(ignore - f1(a), gcFirst=T) [1] 0.08 0.06 0.14 0.00 0.00 system.time(ignore - f2(a), gcFirst=T) [1] 0.19 0.01 0.20 0.00 0.00 system.time(ignore - f3(a), gcFirst=T) [1] 0.04 0.01 0.05 0.00 0.00 Note 1) the use of gcFirst=T 2) these need to be run several times to tune the gc() behaviour. After which f1(a) caused a couple of gc()s and the others none. 3) the Opteron is in general a much faster machine so these are all surprisingly slow. Occasionally the C-level pow() is slow: that happened in one MinGW update and for R Windows users we replaced it. Even there I was unsure that it would make enough difference on a real problem, but eventually found one where it did (MDS with Minkowski distances). It was because of that real problem that I added the special case for 0.5, after careful timing (and hearing Luke's comment alluded to above). -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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] x*x*x*... vs x^n
On 6/29/2005 9:31 AM, Robin Hankin wrote: Hi Duncan On Jun 29, 2005, at 02:04 pm, Duncan Murdoch wrote: On 6/29/2005 7:32 AM, Robin Hankin wrote: Hi I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [snip] If you look in src/main/arithmetic.c, you'll see that R does a general real-valued power (using C's pow() function) whenever either one of the args is real (except for a few special cases, e.g. non-numbers, or powers of 2 or 0.5). There is an internal R function equivalent to your f3, but it is not used in the situation of real^integer (and in any case, x^8 is real^real). I suppose if you wanted to submit a patch someone would take a look, but the question is whether there is really any calculation whose execution time would be materially affected by this. Most computations are not dominated by integer power calculations, so is this really worth the trouble? Duncan Murdoch well, the Gnu Scientific Library has the pow_int() functions, which are a generalization of f3(), so someone thinks so. I did a speed test of them but they were much slower than R (for any of f1(), f2(), f3()): library(gsl) system.time(ignore - pow_int(a,8)) [1] 1.07 1.11 3.08 0.00 0.00 why the slow execution time? Shouldn't you ask the gsl maintainer that? :-) But I'm far more interested in the philosophy behind your comments. I would say that it definitely *is* worth the trouble because someone, somewhere, will want fast integer powers, and possibly use R for nothing else. Ken's point about matrix exponentiation is relevant here too. This is a stated design consideration in Mathematica, I think. (Of course, I'm not suggesting that other programming tasks be suspended! All I'm pointing out is that there may exist a user to whom fast integer powers are very very important) Then that user should submit the patch, not you. But whoever does it should include an argument to convince an R core member that the change is worth looking at, and do it well enough that the patch is accepted. Changes to the way R does arithmetic affect everyone, so they had better be done right, and checking them takes time. Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] all connections are in use error during lazyload stage of packa ge installation
Hi, I suddenly started getting strange errors while working on my caTools package: RCMD install C:/programs/R/rw2011/src/library/caTools .. preparing package caTools for lazy loading Error in file(file, r, encoding = encoding) : all connections are in use Execution halted make: *** [lazyload] Error 1 *** Installation of caTools failed *** I searched through R website and mailing lists but did not found many clues. I am not running R environment so showConnections and CloseAllConnections functions are not helpful. Did I run into some OS specific limit on number of files allowed? I closed all other programs and rebooted my machine but it did not help. Any hints would be appreciated. My System is: - R version: R 2.1.1 - Operating System: Win XP - Compiler: mingw32-gcc-3.4.2 Jarek \=== Jarek Tuszynski, PhD. o / \ Science Applications International Corporation \__,| (703) 676-4192 \ [EMAIL PROTECTED] `\ __ 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] all connections are in use error during lazyload stage of packa ge installation
On 6/29/05, Tuszynski, Jaroslaw W. [EMAIL PROTECTED] wrote: Hi, I suddenly started getting strange errors while working on my caTools package: RCMD install C:/programs/R/rw2011/src/library/caTools .. preparing package caTools for lazy loading Error in file(file, r, encoding = encoding) : all connections are in use Execution halted make: *** [lazyload] Error 1 *** Installation of caTools failed *** I searched through R website and mailing lists but did not found many clues. I am not running R environment so showConnections and CloseAllConnections functions are not helpful. Did I run into some OS specific limit on number of files allowed? I closed all other programs and rebooted my machine but it did not help. I don't know what the problem is but you could see if adding LazyLoad: no to the DESCRIPTION file causes it to go away. Even if it does it would probably be best to track down what the cause is and in the past when I have had problems running R CMD on a package of mine I repeatedly removed half my package and reran R CMD until I found the offending portion. __ 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] x*x*x*... vs x^n
On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote: On 6/29/2005 9:31 AM, Robin Hankin wrote: Hi Duncan library(gsl) system.time(ignore - pow_int(a,8)) [1] 1.07 1.11 3.08 0.00 0.00 why the slow execution time? Shouldn't you ask the gsl maintainer that? :-) well I did ask myself, but in this case the gsl maintainer told me to ask the gsl co-author, who is no doubt much better informed in these matters ;-) (Of course, I'm not suggesting that other programming tasks be suspended! All I'm pointing out is that there may exist a user to whom fast integer powers are very very important) Then that user should submit the patch, not you. But whoever does it should include an argument to convince an R core member that the change is worth looking at, and do it well enough that the patch is accepted. Changes to the way R does arithmetic affect everyone, so they had better be done right, and checking them takes time. yes, that's a fair point. But including a native R command pow.int(), say, wouldn't affect anyone, would it? One could even use the (tested) GSL code, as it is GPL'ed. This would just be a new function that users could use at their discretion, and no existing code would break. I assume that such a function would not suffer whatever performance disadvantage that the GSL package approach had, so it may well be quite a significant improvement over the method used by R_pow_di() in arithmetic.c especially for large n. best wishes rksh Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] moving correlation coef ?
Hello, R gives us the correlation functions cor(). (Many thanks ;-)) Does it also exist a moving correlation coefficient ? (like the moving average). If not, could someone give me some infos or link on how to practically implement such a function in R. (I did search for moving correlation on the R homepage but didn't find anything.) Thank you. Vincent __ 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] moving correlation coef ?
vincent wrote: Hello, R gives us the correlation functions cor(). (Many thanks ;-)) Does it also exist a moving correlation coefficient ? (like the moving average). If not, could someone give me some infos or link on how to practically implement such a function in R. (I did search for moving correlation on the R homepage but didn't find anything.) Thank you. Vincent See ?running in the gtools package: library(gtools) X - rnorm(100); Y - rnorm(100) running(X, Y, cor) HTH, --sundar __ 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] all connections are in use error during lazyload stage of packa ge installation
I found the problem, by doing comparison of directories and files of working and not working versions, and applying changes one by one until one caused install to fail. It was a case of having a call to a source function somewhere in my code, that I forgot about. I was definitely doing something silly , but it was not a meaningful error message. Jarek \=== Jarek Tuszynski, PhD. o / \ Science Applications International Corporation \__,| (703) 676-4192 \ [EMAIL PROTECTED] `\ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Sent: Wednesday, June 29, 2005 10:07 AM To: Tuszynski, Jaroslaw W. Cc: r help Subject: Re: [R] all connections are in use error during lazyload stage of packa ge installation On 6/29/05, Tuszynski, Jaroslaw W. [EMAIL PROTECTED] wrote: Hi, I suddenly started getting strange errors while working on my caTools package: RCMD install C:/programs/R/rw2011/src/library/caTools .. preparing package caTools for lazy loading Error in file(file, r, encoding = encoding) : all connections are in use Execution halted make: *** [lazyload] Error 1 *** Installation of caTools failed *** I searched through R website and mailing lists but did not found many clues. I am not running R environment so showConnections and CloseAllConnections functions are not helpful. Did I run into some OS specific limit on number of files allowed? I closed all other programs and rebooted my machine but it did not help. I don't know what the problem is but you could see if adding LazyLoad: no to the DESCRIPTION file causes it to go away. Even if it does it would probably be best to track down what the cause is and in the past when I have had problems running R CMD on a package of mine I repeatedly removed half my package and reran R CMD until I found the offending portion. __ 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] x*x*x*... vs x^n
I ran 100 repetitions of the 3 multiplications that Robin had compared. Here are the summaries of system times (I only took the first component of system.time) that I obtained. It is clear that f1() is nearly twice as slow as f2() which is slightly slower (not 3 times slower as claimed by Robin) than f3(). So, I don't think that there is much to choose between the cleverer way and the most obvious way to compute integer powers. summary(f1time) Min. 1st Qu. MedianMean 3rd Qu.Max. 0.060 0.170 0.210 0.199 0.230 0.300 summary(f2time) Min. 1st Qu. MedianMean 3rd Qu.Max. 0.060 0.100 0.110 0.128 0.170 0.190 summary(f3time) Min. 1st Qu. MedianMean 3rd Qu.Max. 0. 0.0300 0.0950 0.0779 0.1100 0.1300 Ravi. -- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] -- -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Tuszynski, Jaroslaw W. Sent: Wednesday, June 29, 2005 8:25 AM To: 'Robin Hankin'; r-help Subject: Re: [R] x*x*x*... vs x^n I tried your code and got different results: system.time(ignore - f1(a)) [1] 0.83 0.09 1.08 NA NA system.time(ignore - f2(a)) [1] 0.38 0.01 0.41 NA NA system.time(ignore - f3(a)) [1] 0.32 0.04 0.43 NA NA So I tried it again but with a loop and got: for(i in 1:10) cat(system.time(ignore - f2(a)), \n) 0.36 0.04 0.44 NA NA 0.32 0.01 0.34 NA NA 0.28 0.03 0.32 NA NA 0.29 0.03 0.35 NA NA 0.3 0.02 0.32 NA NA 0.28 0.03 0.32 NA NA 0.3 0.02 0.32 NA NA 0.29 0.02 0.34 NA NA 0.23 0.03 0.32 NA NA 0.42 0 0.45 NA NA for(i in 1:10) cat(system.time(ignore - f3(a)), \n) 0.19 0.04 0.25 NA NA 0.17 0.04 0.25 NA NA 0.21 0.02 0.25 NA NA 0.21 0.02 0.23 NA NA 0.18 0.04 0.23 NA NA 0.18 0.05 0.23 NA NA 0.18 0.04 0.25 NA NA 0.17 0.06 0.23 NA NA 0.2 0.02 0.23 NA NA 0.14 0.06 0.25 NA NA It seems to me that f3 is 50% slower than f2 not 300%. My System is: - R version: R 2.1.1 - Operating System: Win XP - Compiler: mingw32-gcc-3.4.2 Jarek \=== Jarek Tuszynski, PhD. o / \ Science Applications International Corporation \__,| (703) 676-4192 \ [EMAIL PROTECTED] `\ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Robin Hankin Sent: Wednesday, June 29, 2005 7:32 AM To: r-help Subject: [R] x*x*x*... vs x^n Hi I have been wondering if there one can speed up calculating small powers of numbers such as x^8 using multiplication. In addition, one can be a bit clever and calculate x^8 using only 3 multiplies. look at this: f1 - function(x){x*x*x*x*x*x*x*x} f2 - function(x){x^8} f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)} [so f1() and f2() and f3() are algebraically identical] a - rnorm(100) system.time(ignore - f1(a)) [1] 0.50 0.17 2.88 0.00 0.00 system.time(ignore - f2(a)) [1] 0.31 0.03 1.40 0.00 0.00 system.time(ignore - f3(a)) [1] 0.10 0.07 0.18 0.00 0.00 [these figures show little variance from trial to trial] I was expecting f2() and f3() to be about the same. I was not expecting a factor of 3 there! anyone got any comments? -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-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] x*x*x*... vs x^n
In general, the Russian peasant algorithm, which requires only O(log n) multiplications, is very good. Section 4.6.3 of Knuth's The Art of Computer Programming. Volume 2: Seminumerical Algorithms has an in depth discussion. I have had to use this in the past, when computers were slower and compilers were not so clever. It is also better when x is not just a real number, say complex or matrix (as has been mentioned.) In many cases though, one needs many powers sequentially, and then it may be more efficient to just multiply the previous power by x and use the power, etc. (unless you have a parallel computer.) So pows - x^(1:1000) # use pows in calculations could be sped up by employing a faster algorithm, but probably a loop will be faster: pows - 1 for(i in 1:1000) { pows - pows * x # use this power } David L. Reiner, Ph.D. Rho Trading -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Robin Hankin Sent: Wednesday, June 29, 2005 9:13 AM To: Duncan Murdoch Cc: r-help; Robin Hankin Subject: Re: [R] x*x*x*... vs x^n On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote: On 6/29/2005 9:31 AM, Robin Hankin wrote: Hi Duncan library(gsl) system.time(ignore - pow_int(a,8)) [1] 1.07 1.11 3.08 0.00 0.00 why the slow execution time? Shouldn't you ask the gsl maintainer that? :-) well I did ask myself, but in this case the gsl maintainer told me to ask the gsl co-author, who is no doubt much better informed in these matters ;-) (Of course, I'm not suggesting that other programming tasks be suspended! All I'm pointing out is that there may exist a user to whom fast integer powers are very very important) Then that user should submit the patch, not you. But whoever does it should include an argument to convince an R core member that the change is worth looking at, and do it well enough that the patch is accepted. Changes to the way R does arithmetic affect everyone, so they had better be done right, and checking them takes time. yes, that's a fair point. But including a native R command pow.int(), say, wouldn't affect anyone, would it? One could even use the (tested) GSL code, as it is GPL'ed. This would just be a new function that users could use at their discretion, and no existing code would break. I assume that such a function would not suffer whatever performance disadvantage that the GSL package approach had, so it may well be quite a significant improvement over the method used by R_pow_di() in arithmetic.c especially for large n. best wishes rksh Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Selecting rows regarding the frequency of a factor variable.
Hi and sorry to disturb, I'll try to be as clear as possible: I want to select rows of a data frame called Data2.Iso regarding the frequency of a factor variable called Variete that I want =4. I used function table to have the frequency: FRAMEVARIETE-as.data.frame(table(Data2.Iso$Variete)) Then I selected the modalities with a frequency =4: FRAMEVARIETE2-FRAMEVARIETE[FRAMEVARIETE$Freq=4,] as.character(FRAMEVARIETE2[,Variete]) But then, how to select the rows with those modalities? Does anyone can help me? Thanks! Ghislain. __ 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] moving correlation coef ?
or ?rapply in package zoo. spencer graves Sundar Dorai-Raj wrote: vincent wrote: Hello, R gives us the correlation functions cor(). (Many thanks ;-)) Does it also exist a moving correlation coefficient ? (like the moving average). If not, could someone give me some infos or link on how to practically implement such a function in R. (I did search for moving correlation on the R homepage but didn't find anything.) Thank you. Vincent See ?running in the gtools package: library(gtools) X - rnorm(100); Y - rnorm(100) running(X, Y, cor) HTH, --sundar __ 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 -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ 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] Selecting rows regarding the frequency of a factor variab le.
See if this does what you want: dat - data.frame(f=factor(sample(letters[1:10], 100, replace=TRUE)), x=runif(100)) str(dat) `data.frame': 100 obs. of 2 variables: $ f: Factor w/ 10 levels a,b,c,d,..: 2 5 10 9 10 3 9 8 3 1 ... $ x: num 0.9162 0.0481 0.3048 0.0938 0.8599 ... g - names(which(table(dat$f) 11)) g [1] c j dat[dat$f %in% g,] f x 3 j 0.30477413 5 j 0.85992597 6 c 0.86881528 9 c 0.87317095 16 c 0.84252048 18 j 0.24039606 19 j 0.58927414 21 j 0.10077745 32 j 0.72275870 35 c 0.26001549 37 j 0.09608521 40 c 0.15481625 44 c 0.70203309 47 c 0.95907223 50 j 0.35258966 54 c 0.93422614 58 c 0.36546841 61 c 0.55123183 64 j 0.82995122 65 c 0.89104229 66 j 0.81661377 77 j 0.21134708 87 c 0.16602335 92 c 0.02175573 96 j 0.97864088 Andy From: Ghislain Vieilledent Hi and sorry to disturb, I'll try to be as clear as possible: I want to select rows of a data frame called Data2.Iso regarding the frequency of a factor variable called Variete that I want =4. I used function table to have the frequency: FRAMEVARIETE-as.data.frame(table(Data2.Iso$Variete)) Then I selected the modalities with a frequency =4: FRAMEVARIETE2-FRAMEVARIETE[FRAMEVARIETE$Freq=4,] as.character(FRAMEVARIETE2[,Variete]) But then, how to select the rows with those modalities? Does anyone can help me? Thanks! Ghislain. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] MLE with optim
Hello, I tried to fit a lognormal distribution by using optim. But sadly the output seems to be incorrect. Who can tell me where the bug is? test = rlnorm(100,5,3) logL= function(parm, x,...) -sum(log(dlnorm(x,parm,...))) start= list(meanlog=5, sdlog=3) optim(start,logL,x=test)$par Carsten. [[alternative HTML version deleted]] __ 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] x*x*x*... vs x^n
Looking at the code for gsl_pow_int, I see they do use that method. David L. Reiner -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of David Reiner [EMAIL PROTECTED] Sent: Wednesday, June 29, 2005 9:50 AM To: r-help Subject: Re: [R] x*x*x*... vs x^n In general, the Russian peasant algorithm, which requires only O(log n) multiplications, is very good. Section 4.6.3 of Knuth's The Art of Computer Programming. Volume 2: Seminumerical Algorithms has an in depth discussion. I have had to use this in the past, when computers were slower and compilers were not so clever. It is also better when x is not just a real number, say complex or matrix (as has been mentioned.) In many cases though, one needs many powers sequentially, and then it may be more efficient to just multiply the previous power by x and use the power, etc. (unless you have a parallel computer.) So pows - x^(1:1000) # use pows in calculations could be sped up by employing a faster algorithm, but probably a loop will be faster: pows - 1 for(i in 1:1000) { pows - pows * x # use this power } David L. Reiner, Ph.D. Rho Trading -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Robin Hankin Sent: Wednesday, June 29, 2005 9:13 AM To: Duncan Murdoch Cc: r-help; Robin Hankin Subject: Re: [R] x*x*x*... vs x^n On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote: On 6/29/2005 9:31 AM, Robin Hankin wrote: Hi Duncan library(gsl) system.time(ignore - pow_int(a,8)) [1] 1.07 1.11 3.08 0.00 0.00 why the slow execution time? Shouldn't you ask the gsl maintainer that? :-) well I did ask myself, but in this case the gsl maintainer told me to ask the gsl co-author, who is no doubt much better informed in these matters ;-) (Of course, I'm not suggesting that other programming tasks be suspended! All I'm pointing out is that there may exist a user to whom fast integer powers are very very important) Then that user should submit the patch, not you. But whoever does it should include an argument to convince an R core member that the change is worth looking at, and do it well enough that the patch is accepted. Changes to the way R does arithmetic affect everyone, so they had better be done right, and checking them takes time. yes, that's a fair point. But including a native R command pow.int(), say, wouldn't affect anyone, would it? One could even use the (tested) GSL code, as it is GPL'ed. This would just be a new function that users could use at their discretion, and no existing code would break. I assume that such a function would not suffer whatever performance disadvantage that the GSL package approach had, so it may well be quite a significant improvement over the method used by R_pow_di() in arithmetic.c especially for large n. best wishes rksh Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-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] MLE with optim
Carsten Steinhoff wrote: Hello, I tried to fit a lognormal distribution by using optim. But sadly the output seems to be incorrect. Who can tell me where the bug is? test = rlnorm(100,5,3) logL= function(parm, x,...) -sum(log(dlnorm(x,parm,...))) start= list(meanlog=5, sdlog=3) optim(start,logL,x=test)$par Carsten. You are only supplying the meanlog argument to dlnorm. Also you can set log = TRUE in dlnorm to avoid the log computation after calling dlnorm. set.seed(1) x - rlnorm(100, 5, 3) logL - function(par, x) -sum(dlnorm(x, par[1], par[2], TRUE)) start - list(meanlog = 5, sdlog = 2) optim(start, logL, x = x) Why not just use MASS::fitdistr instead? library(MASS) fitdistr(x, dlnorm, start) HTH, --sundar __ 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] MLE with optim
the following work for me: x - rlnorm(1000, 5, 3) fn - function(parms, dat) -sum(dlnorm(dat, parms[1], parms[2], log = TRUE)) optim(c(5, 3), fn, dat = x) library(MASS) fitdistr(x, log-normal, list(meanlog = 5, sdlog = 3)) I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/16/336899 Fax: +32/16/337015 Web: http://www.med.kuleuven.be/biostat/ http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm - Original Message - From: Carsten Steinhoff [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Wednesday, June 29, 2005 5:19 PM Subject: [R] MLE with optim Hello, I tried to fit a lognormal distribution by using optim. But sadly the output seems to be incorrect. Who can tell me where the bug is? test = rlnorm(100,5,3) logL= function(parm, x,...) -sum(log(dlnorm(x,parm,...))) start= list(meanlog=5, sdlog=3) optim(start,logL,x=test)$par Carsten. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] moving correlation coef ?
Thank you for your answers. In fact, i believe my question wasn't precise enough. I don't want to have a moving/sliding windows over the data to correlate (i am already doing that). If I have 2 vectors X = (x1, x2, x3, ..., xt) Y = (y1, y2, x3, ..., yt) I want the most recent elements (t) to have a heavier weight in the correlation calculus than the older ones (1). I think that one simple way to do that is for example to compute cor() over XP, YP where XP = (x1, x2, x2, x3, x3, x3, ..., xt, xt, xt) YP = (y1, y2, y2, y3, y3, y3, ..., yt, yt, yt) ie where each element is repeated several times according to its freshness. It's quite naive ! so if there is a cleaver idea, many thanks. Thanks Vincent __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] poly() in lm() leads to wrong coefficients (but correct residuals)
Dear all, I am using poly() in lm() in the following form. 1 DelsDPWOS.lm3 - lm(DelsPDWOS[,1] ~ poly(DelsPDWOS[,4],3)) 2 DelsDPWOS.I.lm3 - lm(DelsPDWOS[,1] ~ poly(I(DelsPDWOS[,4]),3)) 3 DelsDPWOS.2.lm3 - lm(DelsPDWOS[,1]~DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3)) 1 and 2 lead to identical but wrong results. 3 is correct. Surprisingly (to me) the residuals are the same (correct) in all cases although the coefficients of 1 (and 2) are wrong and do not in any way match the stated regression polynomial. (summaries below) QUESTION: Is there a correct way to use poly() in lm() since version 3 becomes quite tedious if used more often especially with higher order polynomials? I am using R Version 1.9.0 (2004-04-12). Thank you for your time, Andreas Neumann SUMMARIES: summary(DelsDPWOS.lm3) Call: lm(formula = DelsPDWOS[, 1] ~ poly(DelsPDWOS[, 4], 3)) Residuals: 12345678 -0.11414 0.12756 -0.21060 0.04636 -0.03244 0.16030 0.04290 -0.03944 9 0.01949 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 1.30.04861 27.430 1.21e-06 *** poly(DelsPDWOS[, 4], 3)1 -1.274640.14583 -8.741 0.000325 *** poly(DelsPDWOS[, 4], 3)2 0.274830.14583 1.885 0.118175 poly(DelsPDWOS[, 4], 3)3 0.115900.14583 0.795 0.462774 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 0.1458 on 5 degrees of freedom Multiple R-Squared: 0.9416, Adjusted R-squared: 0.9065 F-statistic: 26.86 on 3 and 5 DF, p-value: 0.001645 summary(DelsDPWOS.2.lm3) Call: lm(formula = DelsPDWOS[, 1] ~ DelsPDWOS[, 4] + I(DelsPDWOS[, 4]^2) + I(DelsPDWOS[, 4]^3)) Residuals: 12345678 -0.11414 0.12756 -0.21060 0.04636 -0.03244 0.16030 0.04290 -0.03944 9 0.01949 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) -1.486 7.242 -0.2050.846 DelsPDWOS[, 4] 8.640 14.161 0.6100.568 I(DelsPDWOS[, 4]^2) -6.526 8.799 -0.7420.492 I(DelsPDWOS[, 4]^3)1.375 1.730 0.7950.463 Residual standard error: 0.1458 on 5 degrees of freedom Multiple R-Squared: 0.9416, Adjusted R-squared: 0.9065 F-statistic: 26.86 on 3 and 5 DF, p-value: 0.001645 __ 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] poly() in lm() leads to wrong coefficients (but correct residuals)
On Wed, 2005-06-29 at 18:19 +0200, Andreas Neumann wrote: Dear all, I am using poly() in lm() in the following form. 1 DelsDPWOS.lm3 - lm(DelsPDWOS[,1] ~ poly(DelsPDWOS[,4],3)) 2 DelsDPWOS.I.lm3 - lm(DelsPDWOS[,1] ~ poly(I(DelsPDWOS[,4]),3)) 3 DelsDPWOS.2.lm3 - lm(DelsPDWOS[,1]~DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3)) 1 and 2 lead to identical but wrong results. 3 is correct. Surprisingly (to me) the residuals are the same (correct) in all cases although the coefficients of 1 (and 2) are wrong and do not in any way match the stated regression polynomial. (summaries below) QUESTION: Is there a correct way to use poly() in lm() since version 3 becomes quite tedious if used more often especially with higher order polynomials? The coefficients using 1 and 2 are not incorrect. poly() defines orthogonal polynomials, whereas your DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3 contruct defines an ordinary polynomial. You should be able to recover version 3 coefficients from 1 and 2. See help(poly) x - runif(10) x [1] 0.1878 0.2415 0.5834 0.6556 0.4112 0.3399 0.8144 0.1134 0.7360 0.0463 model.matrix(~ poly(x, 2)) (Intercept) poly(x, 2)1 poly(x, 2)2 11-0.27648 -0.0452 21-0.21052 -0.1899 31 0.20937 -0.2708 41 0.29799 -0.1021 51-0.00212 -0.4117 61-0.08970 -0.3621 71 0.49297 0.4968 81-0.36790 0.2148 91 0.39672 0.1620 10 1-0.45033 0.5082 attr(,assign) [1] 0 1 1 model.matrix(~ x + I(x^2)) (Intercept) x I(x^2) 11 0.1878 0.03528 21 0.2415 0.05834 31 0.5834 0.34040 41 0.6556 0.42982 51 0.4112 0.16911 61 0.3399 0.11554 71 0.8144 0.66320 81 0.1134 0.01286 91 0.7360 0.54169 10 1 0.0463 0.00214 attr(,assign) [1] 0 1 2 Regards, -- Markus Jantti Abo Akademi University [EMAIL PROTECTED] http://www.iki.fi/~mjantti __ 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] moving correlation coef ?
vincent wrote: Thank you for your answers. In fact, i believe my question wasn't precise enough. I don't want to have a moving/sliding windows over the data to correlate (i am already doing that). If I have 2 vectors X = (x1, x2, x3, ..., xt) Y = (y1, y2, x3, ..., yt) I want the most recent elements (t) to have a heavier weight in the correlation calculus than the older ones (1). I think that one simple way to do that is for example to compute cor() over XP, YP where XP = (x1, x2, x2, x3, x3, x3, ..., xt, xt, xt) YP = (y1, y2, y2, y3, y3, y3, ..., yt, yt, yt) ie where each element is repeated several times according to its freshness. It's quite naive ! so if there is a cleaver idea, many thanks. Thanks Vincent Perhaps ?cov.wt will work for you? Your example would be identical to: set.seed(1) X - rnorm(100); Y - rnorm(100) # using cov.wt rho1 - cov.wt(cbind(X, Y), 1:100, cor = TRUE)$cor[1, 2] # your weighting scheme rho2 - cor(X[rep(1:100, 1:100)], Y[rep(1:100, 1:100)]) all.equal(rho1, rho2) # [1] TRUE HTH, --sundar __ 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] moving correlation coef ?
One common weighting scheme is exponentially weighted, i.e., wt = L^(0:m) , where 0 L = 1 . David L. Reiner p.s. If your question is coming from a financial application, you might be interested in the R-sig-finance list, as well as reading the RiskMetrics (r) document Return to RiskMetrics: The Evolution of a Standard, by Jorge Mina and Jerry Yi Xiao (available at their site, after a free registration) where exponentially weighted moving statistics are discussed at length. -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of vincent Sent: Wednesday, June 29, 2005 11:02 AM Cc: [EMAIL PROTECTED] Subject: Re: [R] moving correlation coef ? Thank you for your answers. In fact, i believe my question wasn't precise enough. I don't want to have a moving/sliding windows over the data to correlate (i am already doing that). If I have 2 vectors X = (x1, x2, x3, ..., xt) Y = (y1, y2, x3, ..., yt) I want the most recent elements (t) to have a heavier weight in the correlation calculus than the older ones (1). I think that one simple way to do that is for example to compute cor() over XP, YP where XP = (x1, x2, x2, x3, x3, x3, ..., xt, xt, xt) YP = (y1, y2, y2, y3, y3, y3, ..., yt, yt, yt) ie where each element is repeated several times according to its freshness. It's quite naive ! so if there is a cleaver idea, many thanks. Thanks Vincent __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with regression modeling
Hello all, I'm using R version 2.0.1. I have been having trouble with my linear modeling. I have a table that looks something like this: T RSS DSS LSPFCOLS PS R RTT Actual Max COMM char 5 MSS 2 1251 1 4096 450 0.00164.300.13 64.17 char 5 MSS 50 1251 25 4096 450 0.00174.970.14 74.83 char 5 MSS 200 1251 100 4096 450 0.001111.12 0.13 110.99 char 5 MSS 3 1251 1 4096 450 0.00164.310.14 64.17 char 5 MSS 75 1251 25 4096 450 0.00175.990.14 75.85 char 5 MSS 300 1251 100 4096 450 0.001116.62 0.14 116.48 that I was modeling using the R command model.f - lm(comm ~ I(rtt * ceiling(rss/pf)) + I((s * rss) / ls) + ds*(I(((s+s^2)*rss)/r) + I(((s+s^2)*rss)/(r*ps R gives me a table like this for coefficients I(rtt * ceiling(rss/pf))I((s * rss)/ls) 8.174598e-01 1.353734e+00 dsASA dsDB2 3.101663e+00 2.587372e+00 dsMSS dsOracle 1.575822e+01 9.092041e+00 I(((s + s^2) * rss)/r) I(((s + s^2) * rss)/(r * ps)) 2.151170e-07 3.087502e-05 dsDB2:I(((s + s^2) * rss)/r) dsMSS:I(((s + s^2) * rss)/r) -6.972409e-08 1.812149e-07 dsOracle:I(((s + s^2) * rss)/r)dsDB2:I(((s + s^2) * rss)/(r * ps)) 1.251699e-08 6.687324e-05 dsMSS:I(((s + s^2) * rss)/(r * ps)) dsOracle:I(((s + s^2) * rss)/(r * ps)) 4.731049e-05 -9.243794e-05 What I want at the end of the day is a formula that says given this set of parameters, here is the predicted value. When I plug the above coefficients into Excel, I get a value that is WAY off the fitted value that R gives me. My questions are: 1) How do I use these coefficients in a formula? 2) How can I construct a formula from these coefficients that I can publish in my thesis? Thanks for your help! ty [[alternative HTML version deleted]] __ 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] x*x*x*... vs x^n
I was surprised to find that I was wrong about powers of complexes: seq.pow1 - function(x,n) { + y - rep(x,n) + for(i in 2:n) y[i] - y[i-1] * x + y + } seq.pow2 - function(x,n) x^(1:n) x - 1.001 + 1i * 0.999 # several reps of the following system.time(ignore - seq.pow1(x,10),gcFirst=TRUE) [1] 0.73 0.00 0.74 NA NA # several reps of the following system.time(ignore - seq.pow2(x,10),gcFirst=TRUE) [1] 0.35 0.00 0.35 NA NA I apologize for using probably below when not was correct (modulo grammar). David L. Reiner -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of David Reiner [EMAIL PROTECTED] Sent: Wednesday, June 29, 2005 10:24 AM To: r-help Subject: Re: [R] x*x*x*... vs x^n Looking at the code for gsl_pow_int, I see they do use that method. David L. Reiner -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of David Reiner [EMAIL PROTECTED] Sent: Wednesday, June 29, 2005 9:50 AM To: r-help Subject: Re: [R] x*x*x*... vs x^n In general, the Russian peasant algorithm, which requires only O(log n) multiplications, is very good. Section 4.6.3 of Knuth's The Art of Computer Programming. Volume 2: Seminumerical Algorithms has an in depth discussion. I have had to use this in the past, when computers were slower and compilers were not so clever. It is also better when x is not just a real number, say complex or matrix (as has been mentioned.) In many cases though, one needs many powers sequentially, and then it may be more efficient to just multiply the previous power by x and use the power, etc. (unless you have a parallel computer.) So pows - x^(1:1000) # use pows in calculations could be sped up by employing a faster algorithm, but probably a loop will be faster: pows - 1 for(i in 1:1000) { pows - pows * x # use this power } David L. Reiner, Ph.D. Rho Trading -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Robin Hankin Sent: Wednesday, June 29, 2005 9:13 AM To: Duncan Murdoch Cc: r-help; Robin Hankin Subject: Re: [R] x*x*x*... vs x^n On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote: On 6/29/2005 9:31 AM, Robin Hankin wrote: Hi Duncan library(gsl) system.time(ignore - pow_int(a,8)) [1] 1.07 1.11 3.08 0.00 0.00 why the slow execution time? Shouldn't you ask the gsl maintainer that? :-) well I did ask myself, but in this case the gsl maintainer told me to ask the gsl co-author, who is no doubt much better informed in these matters ;-) (Of course, I'm not suggesting that other programming tasks be suspended! All I'm pointing out is that there may exist a user to whom fast integer powers are very very important) Then that user should submit the patch, not you. But whoever does it should include an argument to convince an R core member that the change is worth looking at, and do it well enough that the patch is accepted. Changes to the way R does arithmetic affect everyone, so they had better be done right, and checking them takes time. yes, that's a fair point. But including a native R command pow.int(), say, wouldn't affect anyone, would it? One could even use the (tested) GSL code, as it is GPL'ed. This would just be a new function that users could use at their discretion, and no existing code would break. I assume that such a function would not suffer whatever performance disadvantage that the GSL package approach had, so it may well be quite a significant improvement over the method used by R_pow_di() in arithmetic.c especially for large n. best wishes rksh Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the
[R] Generalized Linear Mixed Models
Hello! I am trying to fit a Generalized Linear Mixed Model, ordinally I use GLIMMIX macros in SAS System, but I would like to fit this kind of models in R. Could anyone help me on what package I should use to? Thanks in advance Francisco [[alternative HTML version deleted]] __ 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] Generalized Linear Mixed Models
submit the following in R: RSiteSearch('glmm',restr='functions') -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA The business of the statistician is to catalyze the scientific learning process. - George E. P. Box -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: Wednesday, June 29, 2005 12:05 PM To: r-help@stat.math.ethz.ch Subject: [R] Generalized Linear Mixed Models Hello! I am trying to fit a Generalized Linear Mixed Model, ordinally I use GLIMMIX macros in SAS System, but I would like to fit this kind of models in R. Could anyone help me on what package I should use to? Thanks in advance Francisco [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] Generalized Linear Mixed Models
-Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of [EMAIL PROTECTED] Sent: 29 June 2005 18:05 To: r-help@stat.math.ethz.ch Subject: [R] Generalized Linear Mixed Models Hello! I am trying to fit a Generalized Linear Mixed Model, ordinally I use GLIMMIX macros in SAS System, but I would like to fit this kind of models in R. Could anyone help me on what package I should use to? Thanks in advance Francisco Si el modelo es de cualquier clase, revisa el package lme4. Si el modelo es espacial en las familias Poisson o binomial, revisa package geoRglm. R. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] return NA
A-c(1,2,NA,7,5) B-c(3,4,1,4,1) C-c(6,5,6,NA,9) D-c(8,7,4,6,2) df1-cbind(A,B,C,D) for(i in seq(1,ncol(df1)-1, by=2)) { ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] ) } Tried several variations but none worked. I wish to find any NA's in column's 1 or 3 and change the numerical value to the right of the NA's . In this case I wish to replace the 1 and 6 respectively with NA. Any help would be appreciated. __ 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] return NA
You need to use is.na(df1[,i]) to test for NA. Andy From: [EMAIL PROTECTED] A-c(1,2,NA,7,5) B-c(3,4,1,4,1) C-c(6,5,6,NA,9) D-c(8,7,4,6,2) df1-cbind(A,B,C,D) for(i in seq(1,ncol(df1)-1, by=2)) { ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] ) } Tried several variations but none worked. I wish to find any NA's in column's 1 or 3 and change the numerical value to the right of the NA's . In this case I wish to replace the 1 and 6 respectively with NA. Any help would be appreciated. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-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] return NA
[EMAIL PROTECTED] wrote: A-c(1,2,NA,7,5) B-c(3,4,1,4,1) C-c(6,5,6,NA,9) D-c(8,7,4,6,2) df1-cbind(A,B,C,D) for(i in seq(1,ncol(df1)-1, by=2)) { ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] ) } Tried several variations but none worked. I wish to find any NA's in column's 1 or 3 and change the numerical value to the right of the NA's . In this case I wish to replace the 1 and 6 respectively with NA. Any help would be appreciated. I think you want ?is.na. for(i in seq(1, ncol(df1) - 1, by = 2)) df1[, i + 1] - ifelse(is.na(df1[, i]), df1[, i], df1[, i + 1]) --sundar __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] predicted survival curve from a cox model
Hi there, I have a predictor varible class which is a categorical variable and a ' coxph' is used to find the coeffients. How can I plot the predicted survival proportion based on this model? Thanks Lisa Wang Princess Margaret Hospital Toronto tel 416 946 4501 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] plot (log scale on y-axis)
I am planning to plot my data on log scale (y-axis). There is a parameter in plot function, which is plot( ..., log=y, ...) While, the problem is that it is with base of e. Is there a way to let me change it to 10 instead of e? Thanks __ 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] return NA
Here is a way to do it without a loop that could save some time for a big dataset. df1 A B C D [1,] 1 3 6 8 [2,] 2 4 5 7 [3,] NA 1 6 4 [4,] 7 4 NA 6 [5,] 5 1 9 2 df2-cbind(0,ifelse(is.na(df1),NA,0))[,-ncol(df1)-1] df2 A B C [1,] 0 0 0 0 [2,] 0 0 0 0 [3,] 0 NA 0 0 [4,] 0 0 0 NA [5,] 0 0 0 0 df3-df1+df2 df3 A B C D [1,] 1 3 6 8 [2,] 2 4 5 7 [3,] NA NA 6 4 [4,] 7 4 NA NA [5,] 5 1 9 2 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: June 29, 2005 3:50 PM To: r-help@stat.math.ethz.ch Subject: [R] return NA A-c(1,2,NA,7,5) B-c(3,4,1,4,1) C-c(6,5,6,NA,9) D-c(8,7,4,6,2) df1-cbind(A,B,C,D) for(i in seq(1,ncol(df1)-1, by=2)) { ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] ) } Tried several variations but none worked. I wish to find any NA's in column's 1 or 3 and change the numerical value to the right of the NA's . In this case I wish to replace the 1 and 6 respectively with NA. Any help would be appreciated. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] sbrier (Brier score) and coxph
Hello I've decided to try and distill an earlier rather ill focused question to try and elicit a response. Any help is greatly appreciated. Why does mod.cox not work with sbrier whilst mod.km does? Can I make it work? data(DLBCL) DLBCL.surv-Surv(DLBCL$time,DLBCL$cens) mod.km-survfit(DLBCL.surv) mod.cox-survfit(coxph(DLBCL.surv~IPI, data=DLBCL)) sbrier(DLBCL.surv, mod.km) integrated Brier score 0.2076454 attr(,time) [1] 1.3 129.9 sbrier(DLBCL.surv, mod.cox) Error in switch(ptype, survfit = { : switch: EXPR must return a length 1 vector Thanks in advance Stephen Henderson ** This email and any files transmitted with it are confidentia...{{dropped}} __ 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] sbrier (Brier score) and coxph
Is this sbrier from package ipred? The short answer is that it contains ptype - class(pred) and assumes that is of length one. For a survfit.coxph fit it is of class c(survfit.cox, survfit). I suspect from the help page that this is not supported, but you need to contact the authors (as the posting guide suggests). On Wed, 29 Jun 2005, Stephen Henderson wrote: Hello I've decided to try and distill an earlier rather ill focused question to try and elicit a response. Any help is greatly appreciated. Why does mod.cox not work with sbrier whilst mod.km does? Can I make it work? data(DLBCL) DLBCL.surv-Surv(DLBCL$time,DLBCL$cens) mod.km-survfit(DLBCL.surv) mod.cox-survfit(coxph(DLBCL.surv~IPI, data=DLBCL)) sbrier(DLBCL.surv, mod.km) integrated Brier score 0.2076454 attr(,time) [1] 1.3 129.9 sbrier(DLBCL.surv, mod.cox) Error in switch(ptype, survfit = { : switch: EXPR must return a length 1 vector Thanks in advance Stephen Henderson ** This email and any files transmitted with it are confidentia...{{dropped}} __ 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 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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 c:\a\b to c:/a/b
Thank You, Prof. Ripley! Both test1.R and test2.R worked for me just now, as did the following minor modification: (x - readLines(stdin(), n=1)) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt Thanks again. spencer graves Prof Brian Ripley wrote: On Wed, 29 Jun 2005, David Duffy wrote: I couldn't resist adding a more literal answer This can only work for escapes which are preserved. The parser maps \n to a character (LF) and the deparser maps it back to \n. This happens to be true of \a \b \f \n \r \t \v \\ but no others. For example, \s is mapped to s, and there is no difference between \s and s in the parsed input. unback - function(x) { chars - unlist(strsplit(deparse(x),)) chars - chars[-c(1,length(chars))] paste(gsub(,/,chars),collapse=) } unback(\n) unback(\s) [1] s Spencer Graves keeps on insisting there is a better way, but all the solutions are to avoid sending the string to the parser, and hence avoiding having the string directly in an R script. This is common in shell scripts, which use 'here' documents to avoid 'quoting hell'. We can do that in R too. Here are two variants I have not seen in the thread test1.R: scan(, , allowEscapes=FALSE, n=1, quiet=TRUE) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt catIx, \n, sep=) R --slave --vanilla test1.R D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt (This one does not allow quoted strings.) test2.R: x - readLines(stdin(), n=1) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt x - gsub('^(.*)$', \\1, x) cat(x, \n) R --slave --vanilla test2.R D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt (This one allows surrounding quotes or not.) -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ 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 c:\a\b to c:/a/b?
One other comment. Ninotech Path Copy, which can be found at: http://home.worldonline.dk/ninotech/ is a free Windows utility that appears in the Windows Explorer context menu (i.e. it appears as the Copy Path menu entry when you right click any file in Windows Explorer). I had forgotten about it since even though I had installed it a long time ago I hardly ever use it. When one right clicks and chooses Copy Path one is given a number of different ways in which to copy the path. Using the Setup submenu (i.e. right click any file and choose Setup) one can define additional ways and its quite easy to define a new method which replaces backslashes with forward slashes as it copies -- in fact, there is a check box specifically for this. On 6/27/05, Spencer Graves [EMAIL PROTECTED] wrote: Hi, Gabor, James, et al.: Thanks for the help. I'm unable to get scan('clipboard', what='', allowEscapes=FALSE)) to work reliably on my system using Rgui under Windows XP, R 2.1.1 patched. However, file.choose() works like a charm. Thanks again, everyone. spencer graves Gabor Grothendieck wrote: Try using file.choose() to locate the file instead of Windows Explorer. That will return the name in a form useable within R. On 6/27/05, Spencer Graves [EMAIL PROTECTED] wrote: Thanks, Dirk, Gabor, Eric: You all provided appropriate solutions for the stated problem. Sadly, I oversimplified the problem I was trying to solve: I copy a character string giving a DOS path from MS Windows Explorer into an R script file, and I get something like the following: D:\spencerg\statmtds\R\Rnews I want to be able to use this in R with its non-R meaning, e.g., in readLine, count.fields, read.table, etc., after appending a file name. Your three solutions all work for my oversimplified toy example but are inadequate for the problem I really want to solve. Thanks, spencer graves Gabor Grothendieck wrote: On 6/27/05, Dirk Eddelbuettel [EMAIL PROTECTED] wrote: On 26 June 2005 at 20:30, Spencer Graves wrote: | How can one convert back slashes to forward slashes, e.g, changing | c:\a\b to c:/a/b? I tried the following: | | gsub(, /, c:\a\b) | [1] c:\a\b This does work, provided you remember that single backslashed don't exist as e.g. \a is a character in itself. So use doubles are you should be fine: gsub(, /, c:\\a\\b) [1] c:/a/b Also, if one finds four backslashes confusing one can avoid the use of four via any of these: gsub([\\], /, c:\\a\\b) gsub(\\, /, c:\\a\\b, fixed = TRUE) chartr(\\, /, c:\\a\\b) __ 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 -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ 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 c:\a\b to c:/a/b
Note that if you want to source it rather than than run it as a batch job from the command line you will something like this. The way it works is that one puts the file name into comments that are marked with tags and the script rereads itself as data picking out the tagged lines and removing the tags from those lines. script.name() returns the name of the script its running in and my.stdin sets up a connection to the script as data after stripping out all lines that do not match the tag and removing the tag from those lines that do. (This is a variation of some code that I previously posted.) # source the following from R script.name - function() showConnections()[as.character(eval.parent(quote(file), n=3)), description] my.stdin - function( tag = ^#, script.name. = script.name()) textConnection( sub(tag, , grep(tag,readLines(script.name.),value=TRUE)) ) x - readLines( my.stdin() ) #c:\a\b cat(x, \n) On 6/29/05, Spencer Graves [EMAIL PROTECTED] wrote: Thank You, Prof. Ripley! Both test1.R and test2.R worked for me just now, as did the following minor modification: (x - readLines(stdin(), n=1)) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt Thanks again. spencer graves Prof Brian Ripley wrote: On Wed, 29 Jun 2005, David Duffy wrote: I couldn't resist adding a more literal answer This can only work for escapes which are preserved. The parser maps \n to a character (LF) and the deparser maps it back to \n. This happens to be true of \a \b \f \n \r \t \v \\ but no others. For example, \s is mapped to s, and there is no difference between \s and s in the parsed input. unback - function(x) { chars - unlist(strsplit(deparse(x),)) chars - chars[-c(1,length(chars))] paste(gsub(,/,chars),collapse=) } unback(\n) unback(\s) [1] s Spencer Graves keeps on insisting there is a better way, but all the solutions are to avoid sending the string to the parser, and hence avoiding having the string directly in an R script. This is common in shell scripts, which use 'here' documents to avoid 'quoting hell'. We can do that in R too. Here are two variants I have not seen in the thread test1.R: scan(, , allowEscapes=FALSE, n=1, quiet=TRUE) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt catIx, \n, sep=) R --slave --vanilla test1.R D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt (This one does not allow quoted strings.) test2.R: x - readLines(stdin(), n=1) D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt x - gsub('^(.*)$', \\1, x) cat(x, \n) R --slave --vanilla test2.R D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt (This one allows surrounding quotes or not.) -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Extract fixed effects SE from lmer
Hi, Does anyone know how to extract fixed effects SE values from generalized linear mixed models estimated using the lmer function in the lme4 library? I searched attributes and structure with no luck. Thanks Frank A. La Sorte, Ph.D. Department of Fisheries and Wildlife Sciences University of Missouri Columbia, MO 65211 USA __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Fwd: Extract fixed effects SE from lmer
I forgot to cc: the list on this reply. -- Forwarded message -- From: Douglas Bates [EMAIL PROTECTED] Date: Jun 29, 2005 6:28 PM Subject: Re: [R] Extract fixed effects SE from lmer To: La Sorte, Frank A. [EMAIL PROTECTED] On 6/29/05, La Sorte, Frank A. [EMAIL PROTECTED] wrote: Hi, Does anyone know how to extract fixed effects SE values from generalized linear mixed models estimated using the lmer function in the lme4 library? I searched attributes and structure with no luck. Try vcov(myLmerModel). That returns a variance-covariance matrix for the fixed effects. The usual manipulations are used to get the standard errors. For the binomial and Poisson families you should use the optional argument useScale=FALSE as these families do not have a separately estimated scale parameter in the variance function. In the show method for the summary.lmer class the relevant piece of code is corF - as(as(vcov(object, useScale = useScale), pdmatrix), corrmatrix) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] deal package
Hi, I am wondering if anyone here used deal package in R to do the bayesian network. I am curious about its scalability: how many variables and how many observations can it handle in a reasonable time. If you have some good experience, please share your data configurations. thanks, -- Weiwei Shi, Ph.D Did you always know? No, I did not. But I believed... ---Matrix III __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Finding out collinearity in regression
Hi, I am trying to find out a collinearity in explanatory variables with alias(). I creat a dataframe: dat - ds[,sapply(ds,nlevels)=2] dat$Y - Response Explanatory variables are factor and response is continuous random variable. When I run a regression, I have the following error: fit - aov( Y ~ . , data = dat) Error in contrasts-(`*tmp*`, value = contr.treatment) : contrasts can be applied only to factors with 2 or more levels I think there is a dependency in explanatory variables. So, I wanted to use alias to find out a dependency in design matrix but I can't because I cannot create fit in the first place. One of examples I found is: carprice1.lm - lm(gpm100 ~ Type+Min.Price+Price+Max.Price+Range.Price,data=carprice) alias(carprice1.lm) But, what if I can create lm object ? Then is there a way to find out a dependency in design matrix? Thanks a lot for help in advance! -Young. __ 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] Finding out collinearity in regression
Hi, Young Cho wrote: fit - aov( Y ~ . , data = dat) Error in contrasts-(`*tmp*`, value = contr.treatment) : contrasts can be applied only to factors with 2 or more levels I think there is a dependency in explanatory variables. So, I wanted to use alias to find out a dependency in design matrix but I can't because I cannot create fit in the first place. The error message actually looks like you have got (at least) a variable that only has 1 level, e.g. a factor with only one level. Cheers, Kev -- Ko-Kang Kevin Wang PhD Student Centre for Bioinformation Science Building 27, Room 1004 Mathematical Sciences Institute (MSI) Australian National University Canberra, ACT 2601 Australia Homepage: http://wwwmaths.anu.edu.au/~wangk/ Ph (W): +61-2-6125-2431 Ph (H): +61-2-6125-7407 Ph (M): +61-40-451-8301 __ 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] Finding out collinearity in regression
At 01:45 PM 30/06/2005, Young Cho wrote: Hi, I am trying to find out a collinearity in explanatory variables with alias(). I creat a dataframe: dat - ds[,sapply(ds,nlevels)=2] dat$Y - Response Explanatory variables are factor and response is continuous random variable. When I run a regression, I have the following error: fit - aov( Y ~ . , data = dat) Error in contrasts-(`*tmp*`, value = contr.treatment) : contrasts can be applied only to factors with 2 or more levels 1. Sounds like at least one of your factors has only one level. This should be easy to spot. 2. Have you considered package perturb? HTH, Simon. I think there is a dependency in explanatory variables. So, I wanted to use alias to find out a dependency in design matrix but I can't because I cannot create fit in the first place. One of examples I found is: carprice1.lm - lm(gpm100 ~ Type+Min.Price+Price+Max.Price+Range.Price,data=carprice) alias(carprice1.lm) But, what if I can create lm object ? Then is there a way to find out a dependency in design matrix? Thanks a lot for help in advance! -Young. __ 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 Simon Blomberg, B.Sc.(Hons.), Ph.D, M.App.Stat. Centre for Resource and Environmental Studies The Australian National University Canberra ACT 0200 Australia T: +61 2 6125 7800 email: Simon.Blomberg_at_anu.edu.au F: +61 2 6125 0757 CRICOS Provider # 00120C __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] how to call egarch of sas in R
I use R to generate data and I need to estimate the data by egarch (that doesn't have in R). So how I can call egarch from SAS in R. Regards, luck __ 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] moving correlation coef ?
Sundar Dorai-Raj a écrit : Perhaps ?cov.wt will work for you? Your example would be identical to: set.seed(1) X - rnorm(100); Y - rnorm(100) # using cov.wt rho1 - cov.wt(cbind(X, Y), 1:100, cor = TRUE)$cor[1, 2] # your weighting scheme rho2 - cor(X[rep(1:100, 1:100)], Y[rep(1:100, 1:100)]) all.equal(rho1, rho2) # [1] TRUE HTH, --sundar Thank you very much Sundar for your advices. I will try this way. Have a good day. Vincent __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Dispersion parameter in Neg Bin GLM
Hi, Can someone tell me if it is possible to set the dispersion parameter constant when fitting a negative binomial glm in R? I've looked at the documentation and can't find the appropriate argument to pass. In STATA I can type: nbreg depvar [indepvar...], offset(offset) dispersion(constant). Thank you [[alternative HTML version deleted]] __ 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