Re: [R] iraq statistics - OT
It seems that as time goes by, people including statisticians, forget the past and must re-invent it. Anyone interested should read Richardson's The Statistics of Deadly Quarrels. Volume 2 of The World of Mathematics. The book used to be given out as sort of a cracker-jack prize by book clubs everywhere. Gabor Grothendieck wrote: I came across this one: http://www.nysun.com/article/32787 which says that the violent death rate in Iraq (which presumably includes violent deaths from the war) is lower than the violent death rate in major American cities. Does anyone have any insights from statistics on how to interpret this? __ 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 -- -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] OT: DOE - experiments for teaching
Since ECHIP no longer exists, I guess it is OK to reveal some of the stuff we did. We taught a course to engineers for more than 20 years and our statistics indicated that 9 out of every 10 attendees ran at least one experiment after the course. One of the things that was done was intended to break the mind set of the students, and we did this by posing a problem for them to solve before we attempted to teach them anything. We usually used the funnel, although in the early days it was a computer simulation. Many students had some idea of design to begin with, and tried to put what they knew to work with mixed success, while others rolled their own. Teams were used, and someone from each team reported the results. Criticism was always constructive. We followed this up by solving the problem using a design, and we explained the methodology as we went along with all students collecting and analyzing the data on their own computers. Almost all students (even those reluctantly required to attend by management) bought into the problem; becoming so interested, that we often had to chase them from the classroom at the end of the day. Berton Gunter wrote: I've had fun and luck with the apparatus described in my little paper: THROUGH A FUNNEL SLOWLY WITH BALL BEARING AND INSIGHT TO TEACH EXPERIMENTAL DESIGN The American Statistician, 47, 4 p. 265-269 (1993) We continue to use this in our industrial training. I also would strongly second Spencer's remarks re the difficulty of helping students see the big picture. For some reason, viewing experimentation as part of an overall learning process/strategy does not seem to be part of most scientist's or engineer's formal education. I suppose if you look at typical science or engineering labs where the goal is to come to a predetermined conclusion, it's not hard to see why. But we don't need to get into that imbroglio here. -- 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 Spencer Graves Sent: Friday, May 05, 2006 3:59 PM To: Thomas Kaliwe Cc: r-help@stat.math.ethz.ch Subject: Re: [R] OT: DOE - experiments for teaching I fully endorse Richard Heiberger's recommendation of the Bill Hunter articles on teaching experimental design. For a college-level semester D0E class, I had students do experiments in groups. I found it wise to have them do a preliminary presentation with a discussion of the experimental design plus their protocol for managing all the details of test materials, data collection, etc., then a final report with the results. Many students did fine, but some were clearly clueless about the whole process, which indicated a need for some adjustment in what I taught or in some individual assistance. If this is just a few hours or a 1-day thing, you might consider http://www.prodsyse.com/exped2b.pdf;. hope this helps. Spencer Graves Thomas Kaliwe wrote: Hi, I'm sorry for this not being related to R but I think this is a good place to ask. I'm looking for DOE examples(experiments) that can be done at home or in class, such as Paper Helicopter, Paper Towel etc.. I'm thankful for any comment. Thomas [[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 __ 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 -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] [R-pkgs] Zigurrat updated
I've updated the Ziggurat normal and exponential generator in SuppDists in accordance with Leong, et.al (2005). A comment on the implementation of the Ziggurat Method. Jour. Stat. Sci. 12-7, 1-4. The fix takes care of problems for very large sets of random values. I tested it a bit on small samples and it seems ok. It will be on CRAN shortly. -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R: fractional factorial design in R
I think you need to add factors=all to gen.factorial(), otherwise the model df will be less than what you expect. gen.orthogonal.design(c(2,2,3,3,3,3,2,2),numCards=16) [EMAIL PROTECTED] wrote: sorry, some small mistakes in the previuos syntax. This works! design.test - gen.orthogonal.design(c(2,4,3),numCards=16) design.test gen.orthogonal.design - function(listFactors,numCards){ library(AlgDesign) FactorsNames-c(A,B,C,D,E,F,G,H,J,K,L) numFactors-length(listFactors) dat-gen.factorial(listFactors,center=FALSE,varNames=FactorsNames[1:numFacto rs]) desPB-optFederov(~.,dat,nRepeats=20,approximate=FALSE,nTrials=numCards) design-desPB$design#[,2:(numFactors+1)] cat(Number of trials: , fill=T, length(design[,1]), append=T) print(cor(design)) return(design) } However, it is necessary to run the function and guess numCards until the correlation matrix is diagonal and all levels are selected for the final design. Any idea how to solve this problem without an iterative function? Roberto Furlan University of Turin, Italy La mia Cartella di Posta in Arrivo è protetta con SPAMfighter 188 messaggi contenenti spam sono stati bloccati con successo. Scarica gratuitamente SPAMfighter! -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R: fractional factorial design in R
If an orthogonal main effect plan exists for the number of trials you specify, optFederov() in AlgDesign will more than likely find it for you, since such a design should be an optimal design. Ulrike Grömping wrote: I think that there is an understandable wish to have the simple orthogonal plans (and be it only for non-experts to be able to analyse the results themselves). For mixed levels, there is e.g. the L36 that should be able to accomodate plans like 2x2x2x3x3x3. Unfortunately, R is not very strong in this arena. If I had more time, I would think about writing a package on comfortably designing experiments supported e.g. by the catalogues of Chen, J., Sun, D.X., and Wu, C.F.J. (1993). (A catalogue of two-level and three-level fractional factorial designs with small runs. International Statistical Review 61, 131-145.) Such a package should also provide the analysis facilities for any design generated with it, once it has been enriched with observed data. (This is a bit different from the typical R spirit, where users are often required to be experts themselves.) If anyone is planning a project like this or wants to make a diploma student work on it I would be interested in contributing. For the moment, if you want to implement main effects plans of the orthogonal sort (e.g. a Taguchi-plan like the L36) you have to use books or tables published on the internet, if you don't want to use expensive software like SPSS - not very comfortable, but possible. For example, you can find the L36 - which would be able to accomodate your 2x2x2x3x3x3 - in http://www.itl.nist.gov/div898/handbook/pri/section3/pri33a.htm. With kind regards, Ulrike In general, a main effects design need not be orthogonal -- the main effects merely need to be estimable. The trick is to estimate them with good efficiency, etc. I think you need to consult a local statistician for help to understand what these statistical concepts mean. In your example you could cross the 2^(3-1) with the 3^(3-1) to produce an orthogonal design to estimate main effects. But of course that's 72 runs, which I don't think you would consider small. As a previous poster commented, there are orthogonal mixed level arrays (Addleman, Kempthorne Youden -designs are a couple of phrases to try googling on) which stem from the 1960's. I doubt that, in general, they would satisfy your needs. I have not used the AlgDesign package myself. I suggest you direct questions about it to the author/maintainer, Bob Wheeler. -- 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: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of statistical.model at googlemail.com Sent: Monday, January 23, 2006 12:20 PM To: Berton Gunter; statistical.model at googlemail.com; r-help at stat.math.ethz.ch Subject: [R] R: fractional factorial design in R Yes, you're right. For, say, a 3 x 5 design, one can do this in as few as 7 runs -- but only in general by some version of one-factor-at-a-time (OFAT) designs, which are inefficient. It is easy, via, say model.matrix() to write a general function to produce these. But I think it's a bad idea; more efiicient algorithmic designs are better, IMO, which is why I suggested AlgDesign. You and others are free to disagree, of course. Hi Bert, thanks for your suggestion. However, let us say that i need a 2x2x2x3x3x3 design, which should not be too hard. I've loaded AlgDesign, and i am aware now that gen.factorial allows me to create a full desing. But how to create a main-effects-only factorial design (orthogonal)? I am still not able to produce what i need. The function model.matrix.formula is not very clear... :( Could you please indicate which syntax should i use? I'd really appreciate your help. Thanks in advance, Roberto Furlan University of Turin, Italy La mia Cartella di Posta in Arrivo è protetta con SPAMfighter 188 messaggi contenenti spam sono stati bloccati con successo. Scarica gratuitamente SPAMfighter! __ 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 -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] expected values of order statistics
normOrder() in SuppDists Anna Oganyan wrote: Hello, Could somebody point me, is there any function in R which returns expected values of order statistics for normal distribution? I have been looking and couldn't find it. Thanks! Anna __ 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 -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] Random numbers
You can use Marsaglia's multiply with carry. I haven't looked at the C code in R recently, but doubt if it has changed. The C code is very neat, using 6 #defines: static const double RANDCONST=2.32830643654e-10; unsigned long zSeed=362436069, wSeed=521288629; #define zNew ((zSeed=36969*(zSeed65535)+(zSeed16))16) #define wNew ((wSeed=18000*(wSeed65535)+(wSeed16))65535) #define IUNIFORM (zNew+wNew) #define UNIFORM ((zNew+wNew)*RANDCONST) #define setseed(A,B) zSeed=(A);wSeed=(B); #define getseed(A,B) A=zSeed;B=wSeed; See Marsaglia's DIEHARD page for more details: http://www.stat.fsu.edu/pub/diehard/ Carl wrote: Hi All. I have R code whose functionality is being replicated within a C+ program. The outputs are to be compared to validate the conversion somewhat - however (as is always the case) I have stuffed my code with random number calls. Random uniform numbers in C+ are being produced using the (Boost) mersenne-twister generators (mt11213b mt19937) - which is the default type of generator in R (if I read things correctly). If it was all within R I would just set the seed for reproducibility. Basically - how do I specify in C+ for a set of random uniform numbers such that they are the same as from R? I have considered the possibility of storing/using the R generated random numbers in the C+ version for validation purposes - but there are a lot of them, and that strikes me as a generally ugly way of doing things. thanks in advance C -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] Converting coordinates to actual distances
Paul Brewin wrote: Hello, I've been searching for a method of converting Lat/Lon decimal coordinates into actual distances between points, and taking into account the curvature of the earth. Is there such a package in R? I've looked at the GeoR package, but this does not seem to contain what I am looking for. Ideally the output would be a triangular matrix of distances. Thanks in advance, Paul Brewin Paul E Brewin (PhD) Center for Research in Biological Systems University of California San Diego 9500 Gilman Drive MC 0505 La Jolla CA, 92093-0505 USA Ph: 858-822-0871 Fax: 858-822-3631 __ 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 See http://www.meridianworlddata.com/Distance-Calculation.asp. The calculations are trivial. -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] insignificant factors in regression model
Weijie Cai wrote: Hi list, I am building a regression model with categorical predictor variable coded by treatment contrasts. The summary of the regression model shows that some levels are significant while others are not. The significant ones show that they are statistically significant from the basis factor (at level 0). How do we generally deal with those insignificant levels? If they are used for prediction, sometimes they will produce strange results because their coefficient estimates have large variances. Do we just simply ignore them assuming they are not different from level 0? Or do we exclude them in factors (by treating them as zero's) and refit the model? any suggestions? __ 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 The fact that an estimated coefficient is not significant does not mean that it is negligible. The lower tail of the F-distribution can be used to help you decide whether or not to omit a term. -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] power.anova.test for interaction effects
Your F value is so low as to make me suspect your model. Where did the 144 denominator degrees of freedom come from? Andrew Kniss wrote: This question will probably get me in trouble on theoretical grounds, but I will pose it anyway. The situation: I recently ran a field study looking for differences in sugarbeet cultivar tolerance to a specific herbicide. The study was set up so that 37 cultivars were treated with 4 different applications of the herbicide (37*4 factorial). In doing so, we found that the interaction effect was highly insignificant (ndf=108, ddf=144, F=0.28, p=1.). Now my problem is this... the study takes up an enormous amount of time, energy, and money (as you could guess with 37 cultivars in a field study). We need to determine weather it is worth the effort to repeat the study this summer (practically, it is not, but our funding source would like a more concrete demonstration). I decided to try using power.anova.test just as a starting point to see what our power was. My question is: is this valid to do on an interaction term? If I use power.anova.test with on the interaction term, this is what I get: ~ power.anova.test(groups=(37*4), n=3, between.var=12.06, ~+ within.var=21.23, sig.level=0.05) ~ ~ Balanced one-way analysis of variance power calculation ~ ~ groups = 148 ~ n = 3 ~between.var = 12.06 ~ within.var = 21.23 ~ sig.level = 0.05 ~ power = 1 ~ ~ NOTE: n is number in each group This would imply that given the variability we observed with 3 replications, we almost certainly would have found differences if they existed. But given what I have read on power analysis, a high p-value and wide confidence intervals nearly always suggest inadequate sample size. (Our 90% confidence intervals differed from the estimates by as much as 28%, when a 10% difference would be significant from a practical perspective.) So is this a valid approach? Or does the power.anova.test fall apart if using an interaction effect? Thank you in advance for any help or references you are willing to point me to. Best regards, Andrew Kniss Assistant Research Scientist University of Wyoming Department of Plant Sciences 1000 E. University Ave. Laramie, WY 82071 USA [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 -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] Johnson transformation
[EMAIL PROTECTED] wrote: Hello, I'm Carla, an italian student, I'm looking for a package to transform non normal data to normality. I tried to use Box Cox, but it's not ok. There is a package to use Johnson families' transormation? Can you give me any suggestions to find free software as R that use this trasform? Thank yuo very much Carla 6X velocizzare la tua navigazione a 56k? 6X Web Accelerator di Libero! Scaricalo su INTERNET GRATIS 6X http://www.libero.it __ 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 The Johnson system is in the SuppDists package. -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ 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] 2k-factorial design with 10 parameters
On Tue, 30 Nov 2004, Sven wrote: I'd like to apply a 2^k factorial design with k=10 parameters. Obviously this results in a quite long term for the model equation due to the high number of combinations of parameters. How can I specify the equation for the linear model (lm) without writing all combinations explicitly down by hand? Does a R command exist for this problematic? Use expand.formula() in the AlgDesign package. -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] inst directory
R CMD check now balks at my inst directory. It contains a single folder doc, but apparently any folder causes the problem. If inst is empty, the project checks OK. This was not a problem before 1.9. I've checked the documentation, but don't see a change. What am I missing. -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Density Estimation
Try fitting it with a Johnson function -- see SuppDists. If you can fit it you will then be able to use the functions in SuppDists just as you can for any other distribution supported by R. Brian Mac Namee wrote: Hi there, Sorry if this is a rather loing post. I have a simple list of single feature data points from which I would like to generate a probability that an unseen point comes from the same distribution. To do this I am trying to estimate the probability density of the list of points and use this to generate a probability for the new unseen points. I have managed to use the R density function to generate the density estimate but have not been able to do anything with this - i.e. generate a rpobability that a new point comes from the same distribution. Is there a function to do this, or am I way off the mark using the density function at all? Thanks in advance, Brian. __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Blom's approximation to rankits?
You need not use Blom's approximations. Exact values to several decimal places are given by normOrder() in SuppDists. Lisa Wang wrote: Hello, My name is Lisa and I'm a statistician at Princess Margare Hospital. I wonder if there is any function in R that calculate the Normal rankits based on Blom's approximation? Thank you very much Lisa Wang Msc. Princess Margaret hospital Toronto, Ca __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ghyper package
It is in SuppDists. Román Padilla Lizbeth wrote: Hello I am searching ghyper package (generalized hypergeometric distributions). Does anyone can send it to me? Regards from Mexico Lizbeth Román [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Spearman probabilities and SuppDists
cor.test is giving Pr(x0.6), SuppDists is giving Pr(x=0.6). My recollection is that it is the custom in R for discrete distributions to include the point in the left tail. Drew Hoysak wrote: cor.test and SuppDists give me different P-values for the same Spearman's rho. Which is correct, or am I doing something wrong? x - c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1) y - c( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8) cor.test(x,y,method=spearman) Spearman's rank correlation rho data: x and y S = 48, p-value = 0.0968 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 0.6 2*(1-pSpearman(.6,9)) [1] 0.08572531 Drew Hoysak __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] citation() doesn´t work
The entry type manual does not have a url field. In fact, I don't think any of the types do, but I haven't looked at this in a while. I usually expand the note field to include the url. Dr. Mahmoud K. Okasha wrote: No, I have just used it. the result is: citation() To cite R in publications, use R Development Core Team (2003). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-00-3, URL http://www.R-project.org. We have invested a lot of effort in creating R, please cite it when using it for data analysis. A BibTeX entry for LaTeX users is @Manual{, title= {R: A language and environment for statistical computing}, author = {{R Development Core Team}}, organization = {R Foundation for Statistical Computing}, address = {Vienna, Austria}, year = 2003, note = {ISBN 3-900051-00-3}, url = {http://www.R-project.org} } Best regards... -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] citing a package?
I faced this problem recently when documenting the AlgDesign package. It contains some stuff the isn't in the literature, so I added a citation statement in the AUTHOR section of each function. Even after the material is published, I think a citation to a working model is quite appropriate. Thomas Lumley wrote: On Mon, 9 Feb 2004, Ramon Diaz-Uriarte wrote: Dear Martin, I'd suggest you check the DESCRIPTION file and ask the author(s) of the package (e.g., a package might be related to a tech report which might, now, be in press, or whatever). The posted suggestions seem to be that you don't cite the package, you cite something else vaguely related to it instead. This violates both the purpose of a citation (a link to the original source) and the principle (which I hope R users support) that software is publishable in itself, not just as an appendage to text. Most citation styles give rules for citing software and rules for citing URIs. Even when the package author has been completely unhelpful in constructing a package title you can still put together a perfectly reasonable citation, eg, Lumley T (2003) Rmeta version 2.10. R package. http://cran.r-project.org Some publishers might want a download date, or an explicit statement that it is software (eg to make searching easier). -thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] number point under-flow
It's not the compiler. pghyper() in SuppDists does the same thing. Its just rounding error. Set the result to 0 if it bothers you. pghyper(24,514,53, 5961, lower.tail=F) [1] -3.325965e-12 Roger D. Peng wrote: Did you compile with gcc-2.96? I think there were some problems with the floating point arithmetic with that compiler (at least for the earlier versions released by Red Hat). -roger [EMAIL PROTECTED] wrote: Hello, I've come across the following situation in R-1.8.1 (compile + running under RedHat 7.1): phyper(24, 514, 5961-514, 53, lower.tail=T) [1] 1 phyper(24, 514, 5961-514, 53, lower.tail=F) [1] -1.037310e-11 I'd expect the later to be 0 or some very small positive number. Is this a number under-flow of the calculation? Do you think I'm safe if I just set the result to 0 in these cases? kind regards, Arne __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Call and memory
I use a large real matrix, X, in C code that is passed from R and transposed in place in the C code. I would like to conserve memory and, if possible, allocate space for only one copy of X -- hence I would like to pass a pointer to the data in the X object to the C code. The Writing R Extensions manual says that neither .Call nor .External copy their arguments. They also say that these arguments should be treated as read only. Fine, but in testing I seem to be able to transpose very large X's in place, in C code without an error. This leads me to assume that the manual was just giving good advice about treating arguments as read only. However, I find that I have done nothing to the X in R. It seems that a copy has been made after all. I may as well call the code with .C and avoid the use of macros. Could someone please point out the error in my thinking or suggest a way to accomplish my goal? My code follows: mListTest - function(X,N,k) { .Call(mList,as.double(X),as.integer(N),as.integer(k)); } SEXP mList( SEXP Xi, SEXP Ni, SEXP ki ) { double *pX=NUMERIC_POINTER(Xi); int N=INTEGER_POINTER(Ni)[0]; int k=INTEGER_POINTER(ki)[0]; SEXP alist; SEXP avector; SEXP nvector; SEXP rvector; SEXP kvector; int n; int i; transposeMatrix(pX,N,k); n=4; PROTECT(alist=NEW_LIST(n)); PROTECT(avector=NEW_NUMERIC(200)); for (i=0;i200;i++) { NUMERIC_POINTER(avector)[i]=pX[i]; } SET_ELEMENT(alist,0,avector); UNPROTECT(1); PROTECT(nvector=NEW_INTEGER(1)); INTEGER_POINTER(nvector)[0]=N; SET_ELEMENT(alist,1,nvector); UNPROTECT(1); PROTECT(rvector=NEW_NUMERIC(1)); NUMERIC_POINTER(rvector)[0]=0.5; SET_ELEMENT(alist,2,rvector); UNPROTECT(1); PROTECT(kvector=NEW_INTEGER(1)); INTEGER_POINTER(kvector)[0]=k; SET_ELEMENT(alist,3,kvector); UNPROTECT(1); UNPROTECT(1); return alist; } -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] random number generator?
Well if the base generator is to be changed, it might be a good idea to consider implementing one of Marsaglia's really long period generators. His latest, using titanic primes, has a period of 2^131086. It is extremely fast. Even more could be done, by replacing the two phase normal generator by faster generator using a single phase. Liaw, Andy wrote: Might I suggest taking a poll (even though unscientific) of how many people will be affected by a change in default RNG? My totally arbitrary guess is very few, if any. If I'm not mistaken, Python had only recently changed the default RNG to Mersenne-Twister. If Python can do it, I should think R can, too, without too much pain... Just my $0.02... Andy -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]] Sent: Tuesday, January 28, 2003 3:53 PM To: Charles Annis, P.E. Cc: [EMAIL PROTECTED] Subject: RE: [R] random number generator? Can I suggest RNGkind(Mersenne-Twister, Inversion) and especially the use of Inversion where tail behaviour of the normal is important. Were it not for concerns about reproducibility we would have switched to Inversion a while back. On Tue, 28 Jan 2003, Charles Annis, P.E. wrote: Earlier today I reported finding an unbalanced number of observations in the p=0.0001 tails of rnorm. Many thanks to Peter Dalgaard who suggested changing the normal.kind generator. Using RNGkind(kind = NULL, normal.kind =Box-Muller) seems to have provided the remedy. For example: observed.fraction.below [1] 0.000103 observed.fraction.above [1] 0.000101 Thank you, Peter! Charles Annis, P.E. [EMAIL PROTECTED] phone: 561-352-9699 eFAX: 503-217-5849 http://www.StatisticalEngineering.com -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]] On Behalf Of Peter Dalgaard BSA Sent: Tuesday, January 28, 2003 2:36 PM To: Charles Annis, P.E. Cc: [EMAIL PROTECTED] Subject: Re: [R] random number generator? Charles Annis, P.E. [EMAIL PROTECTED] writes: Dear R-Aficionados: I realize that no random number generator is perfect, so what I report below may be a result of that simple fact. However, if I have made an error in my thinking I would greatly appreciate being corrected. I wish to illustrate the behavior of small samples (n=10) and so generate 100,000 of them. n.samples - 100 sample.size = 10 p - 0.0001 z.normal - qnorm(p) # generate n.samples of sample.size each from a normal(mean=0, sd=1) density # small.sample - matrix(rnorm(n=sample.size*n.samples, mean=0, sd=1), nrow=n.samples, ncol=sample.size) # Verify that from the entire small.sample matrix, p sampled values are below, p above. # observed.fraction.below - sum(small.sample z.normal)/length(small.sample) observed.fraction.above - sum(small.sample -z.normal)/length(small.sample) observed.fraction.below [1] 6.3e-05 observed.fraction.above [1] 0.000142 I've checked the behavior of the entire sample's mean and median and they seem fine. The total fraction in both tails is 0.0002, as it should be. However in every instance about 1/3 are in the lower tail, 2/3 in the upper. I also observe the same 1/3:2/3 ratio for one million samples of ten. Is this simply because random number generators aren't perfect? Or have I stepped in something? Thank you for your kind counsel. You stepped in something, I think, but I probably shouldn't elaborate on the metaphor ... There's an unfortunate interaction between the two methods that are used for generating uniform and normal variables (the latter uses the former). This has been reported a couple of times before and typically gives anomalous tail behaviour. Changing one of the generators (see help(RNGkind)) usually helps. -- 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 __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help -- __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- (302) 239-6620, voice FAX 724 Yorklyn Rd., Hockessin, DE 19707 Randomness comes in bunches