Re: [R] filling a list faster
another solution: l - vector(mode = list, length = 10^5) system.time(for(i in (1:10^5)) l[[i]] - c(i,i+1,i)) On my system this version is even slightly faster than the matrix version ... Best, Matthias - original message Subject: Re: [R] filling a list faster Sent: Fri, 13 Jul 2007 From: Philippe Grosjean[EMAIL PROTECTED] If all the data coming from your iterations are numeric (as in your toy example), why not to use a matrix with one row per iteration? Also, do preallocate the matrix and do not add row or column names before the end of the calculation. Something like: m - matrix(rep(NA, 3*10^5), ncol = 3) system.time(for(i in (1:10^5)) m[i, ] - c(i,i+1,i)) user system elapsed 1.362 0.033 1.424 That is, about 1.5sec on my Intel Duo Core 2.33Mhz MacBook Pro, compared to: l - list(1-c(1,2,3)) system.time(for(i in (1:10^5)) l[[length(l)+1]] - c(i,i+1,i)) user system elapsed 191.629 49.110 248.454 ... more than 4 minutes for your code. By the way, what is your very fast machine, that is actually four times faster than mine (gr!)? Best, Philippe Grosjean ..?})) ) ) ) ) ) ( ( ( ( (Prof. Philippe Grosjean ) ) ) ) ) ( ( ( ( (Numerical Ecology of Aquatic Systems ) ) ) ) ) Mons-Hainaut University, Belgium ( ( ( ( ( .. Balazs Torma wrote: hello, first I create a list: l - list(1-c(1,2,3)) then I run the following cycle, it takes over a minute(!) to complete on a very fast mashine: for(i in (1:10^5)) l[[length(l)+1]] - c(i,i+1,i) How can I fill a list faster? (This is just a demo test, the elements of the list are calculated iteratively in an algorithm) Are there any packages and documents on how to use more advanced and fast data structures like linked-lists, hash-tables or trees for example? Thank you, Balazs Torma __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. --- original message end __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Generic distributions
Hello Alberto, hello Greg, in distr you can do: library(distr) N - Norm(mean = 1, sd = 2) p(N)(0.5) r(N)(100) !!! not: p(N, 0.5) or r(N, 100) !!! A detailed description of package distr is given in package distrDoc. library(distrDoc) vignette(distr) hth Matthias - original message Subject: Re: [R] Generic distributions Sent: Tue, 06 Mar 2007 From: Greg Snow[EMAIL PROTECTED] I think the distr package does this. There are also packages that link to winbugs if that is what you really want to do. -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare [EMAIL PROTECTED] (801) 408-8111 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Alberto Monteiro Sent: Tuesday, March 06, 2007 1:38 PM To: r-help@stat.math.ethz.ch Subject: [R] Generic distributions Is there any class that generalizes distributions? For example, I could say x - generic_distribution(normal, list(mean=1, sigma=0.5)) and then use it like rgeneric_distribution(100, x) to get a sample of 100, or pgeneric_distribution(0.5, x) to get the pdf at (x = 0.5). In the openbugs/winbugs package, that uses a language that looks like R/S, we can do things like x ~ dnorm(mu, tau), forget that x is a normal with mean mu and variance 1/tau, and then treat it generically. Alberto Monteiro PS: this is noise... but due to spam invasion, anything that increases the nonspam/spam ratio should be welcome :-) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. --- original message end __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] RE gap statistic in cluster analysis
there is an implementation in package SLmisc and also in the bioconductor package SAGx. hth Matthias GRAHAM LEASK schrieb: Has anyone implemented Tibrishani's gap statistic in R or S plus? If so I would greatly appreciate the relevant script file. Any help will be much appreciated Kind regards Dr Graham Leask Economics and Strategy Group Aston Business School Aston University Aston Triangle Birmingham B4 7ET Tel: Direct line 0121 204 3150 email [EMAIL PROTECTED] [[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 and provide commented, minimal, self-contained, reproducible code. -- Dr. rer. nat. Matthias Kohl E-Mail: [EMAIL PROTECTED] Home: www.stamats.de __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Integration + Normal Distribution + Directory Browsing Processing Questions
Hi Nils,
I would say, pnorm is faster and has a higher precision.
Best,
Matthias
- original message
Subject: Re: [R] Integration + Normal Distribution + Directory Browsing
Processing Questions
Sent: Mon, 22 Jan 2007
From: Nils Hoeller[EMAIL PROTECTED]
Thank you,
both work fine. Why is pnorm to prefer?
Nils
Matthias Kohl schrieb:
Hi,
why don't you use pnorm?
E.g.,
pnorm(1, mean = 0.1, sd = 1.2) - pnorm(0, mean = 0.1, sd = 1.2)
Matthias
- original message
Subject: Re: [R] Integration + Normal Distribution + Directory Browsing
Processing Questions
Sent: Sun, 21 Jan 2007
From: Dimitris Rizopoulos[EMAIL PROTECTED]
you can use the `...' argument of integrate, e.g.,
integrate(dnorm, 0, 1)
integrate(dnorm, 0, 1, mean = 0.1)
integrate(dnorm, 0, 1, mean = 0.1, sd = 1.2)
look at ?integrate for more info.
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/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
http://www.student.kuleuven.be/~m0390867/dimitris.htm
Quoting Nils Hoeller [EMAIL PROTECTED]:
Hi everyone,
I am new to R, but it's really great and helped me a lot!
But now I have 2 questions. It would be great, if someone can help me:
1. I want to integrate a normal distribution, given a median and sd.
The integrate function works great BUT the first argument has to be a
function
so I do integrate(dnorm,0,1) and it works with standard m. and sd.
But I have the m and sd given.
So for fixed m and sd I work around with a new function mynorm
mynorm - function(n) {
ret - dnorm(n,0.6,0.15)
ret
}
for example.
BUT what can I do for dynamic m and sd?
I want something like integrate(dnorm(,0.6,0.15),0,1), with the first
dnorm parameter open for the
integration but fixed m and sd.
I hope you can help me.
2. I am working with textfiles with rows of measure data.
read.table(file) works fine.
Now I want R to read.table all files within a given directory and
process them one by the other.
for(all files in dir xy) {
x - read.table(nextfile)
process x
}
Is that possible with R? I hope so. Can anyone give me a link to
examples.
Thanks for your help
Nils
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--- original message end
--
Dr. rer. nat. Matthias Kohl
[EMAIL PROTECTED]
www.stamats.de
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--- original message end
--
Dr. rer. nat. Matthias Kohl
[EMAIL PROTECTED]
www.stamats.de
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Integration + Normal Distribution + Directory Browsing Processing Questions
Hi,
why don't you use pnorm?
E.g.,
pnorm(1, mean = 0.1, sd = 1.2) - pnorm(0, mean = 0.1, sd = 1.2)
Matthias
- original message
Subject: Re: [R] Integration + Normal Distribution + Directory Browsing
Processing Questions
Sent: Sun, 21 Jan 2007
From: Dimitris Rizopoulos[EMAIL PROTECTED]
you can use the `...' argument of integrate, e.g.,
integrate(dnorm, 0, 1)
integrate(dnorm, 0, 1, mean = 0.1)
integrate(dnorm, 0, 1, mean = 0.1, sd = 1.2)
look at ?integrate for more info.
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/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
http://www.student.kuleuven.be/~m0390867/dimitris.htm
Quoting Nils Hoeller [EMAIL PROTECTED]:
Hi everyone,
I am new to R, but it's really great and helped me a lot!
But now I have 2 questions. It would be great, if someone can help me:
1. I want to integrate a normal distribution, given a median and sd.
The integrate function works great BUT the first argument has to be a
function
so I do integrate(dnorm,0,1) and it works with standard m. and sd.
But I have the m and sd given.
So for fixed m and sd I work around with a new function mynorm
mynorm - function(n) {
ret - dnorm(n,0.6,0.15)
ret
}
for example.
BUT what can I do for dynamic m and sd?
I want something like integrate(dnorm(,0.6,0.15),0,1), with the first
dnorm parameter open for the
integration but fixed m and sd.
I hope you can help me.
2. I am working with textfiles with rows of measure data.
read.table(file) works fine.
Now I want R to read.table all files within a given directory and
process them one by the other.
for(all files in dir xy) {
x - read.table(nextfile)
process x
}
Is that possible with R? I hope so. Can anyone give me a link to
examples.
Thanks for your help
Nils
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--- original message end
--
Dr. rer. nat. Matthias Kohl
[EMAIL PROTECTED]
www.stamats.de
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] package dependency tree
Hello, http://bioconductor.org/packages/1.9/bioc/html/pkgDepTools.html resp. http://bioconductor.org/packages/2.0/bioc/html/pkgDepTools.html may help you. Best regards, Matthias Gabor Grothendieck schrieb: Try this, noting that available.packages() returns a matrix whose columns include Depends and Suggests and whose rownames are the package names: AP - available.packages() rownames(AP)[grep(quantreg, AP[, Depends])] [1] cobsemplik lss rankreg rqmcmb2 rownames(AP)[grep(quantreg, AP[, Suggests])] [1] apTreeshape diveMovedyn ggplot gsubfn [6] np On 1/2/07, roger koenker [EMAIL PROTECTED] wrote: Is there a painless way to find the names of all packages on CRAN that Depend on a specified package? url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Dr. rer. nat. Matthias Kohl E-Mail: [EMAIL PROTECTED] Home: www.stamats.de __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Questions regarding integrate function
Hello,
your integrand needs to be a function which accepts a numeric vector as
first argument and returns a vector of the same length (see ?integrate).
Your function does not fulfill this requirement. Hence, you have to
rewrite your function or use sapply, apply or friends; something like
newintegrand - function(x) sapply(x, integrand)
By the way you use a very old R version and I would recommend to update
to R 2.4.0.
hth
Matthias
Le Wang schrieb:
Hi there. Thanks for your time in advance.
I am using R 2.2.0 and OS: Windows XP.
My final goal is to calculate 1/2*integral of
(f1(x)^1/2-f2(x)^(1/2))^2dx (Latex codes:
$\frac{1}{2}\int^{{\infty}}_{\infty}
(\sqrt{f_1(x)}-\sqrt{f_2(x)})^2dx $.) where f1(x) and f2(x) are two
marginal densities.
My problem:
I have the following R codes using adapt package. Although adapt
function is mainly designed for more than 2 dimensions, the manual
says it will also call up integrate if the number of dimension
equals one. I feed in the data x1 and x2 and bandwidths h1 and h2.
These codes worked well when my final goal was to take double
integrals.
integrand - function(x) {
# x input is evaluation point for x1 and x2, a 2x1 vector
x1.eval - x[1]
x2.eval - x[2]
# n is the number of observations
n - length(x1)
# x1 and x2 are the vectors read from data.dat
# Compute the marginal densities
f.x1 - sum(dnorm((x1.eval-x1)/h1))/(n*h1)
f.x2 - sum(dnorm((x2.eval-x2)/h2))/(n*h2)
# Return the integrand #
return((sqrt(f.x1)-sqrt(f.x2))**2)
}
estimate-0.5*adapt(1, lo=lo.default, up=up.default,
minpts=minpts.default, maxpts=maxpts.default,
functn=integrand, eps=eps.default, x1, x2,h1,h2)$value
But when I used it for one-dimension, it failed. Some of my
colleagues suggested getting rid of x2.eval in the integrand
because it is only one integral. But after I changed it, it still
didn't work. R gave the error msg: evaluation of function gave a
result of wrong length
I am not a frequent R user..although I looked up the mailing list
for a while and there were few postings asking similar questions, I
can't still figure out why my codes won't work. Any help will be
appreciated.
Le
-
~~
Le Wang, Ph.D
Population Center
University of Minnesota
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Dr. rer. nat. Matthias Kohl
E-Mail: [EMAIL PROTECTED]
Home: www.stamats.de
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Computing time for matrix addition or subtraction
not concerning your subject line, but function crossprod may be useful, too Matthias - original message Subject: Re: [R] Computing time for matrix addition or subtraction Sent: Mon, 13 Nov 2006 From: Prof Brian Ripley[EMAIL PROTECTED] On Sun, 12 Nov 2006, YONGWAN CHUN wrote: I wonder by chance if there is a way to reduce computing time for matrix addition or subtraction. With a lot of iterations, it would be helpful to reduce a little amount time. Yes, by making use of an optimized BLAS: see the R-admin manual. On my (dual CPU) system this reduced your example from 36 to 6 seconds. BTW, it is the calculation of PP that is taking the most of time, not as in your subject line. Simple example is as below n - 2000 P - matrix(rnorm(n*n),n,n) PP - P %*% P M - diag(n) - P R - M + t(M) - diag(n) + PP I would like to reduce time in calculating R. -- 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 and provide commented, minimal, self-contained, reproducible code. --- original message end __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] median of gamma distribution
Hi, to compute the median (or expectation, var, sd, IQR, mad, ...) you can also use package distrEx. library(distrEx) (G - Gammad()) median(G) Matthias - original Nachricht Betreff: Re: [R] median of gamma distribution Gesendet: Fri, 30. Jun 2006 Von: [EMAIL PROTECTED] On 30-Jun-06 Philip He wrote: Doese anyone know a R function to find the median of a gamma distribution? qgamma will do it. Test: -log(0.5) [1] 0.6931472 qgamma(0.5,1) [1] 0.6931472 Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 30-Jun-06 Time: 16:53:16 -- XFMail -- __ 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 --- original Nachricht Ende -- Dr. rer. nat. Matthias Kohl [EMAIL PROTECTED] www.stamats.de __ 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] Generic code for simulating from a distribution.
Hi,
have a look at the packages distr and distrSim which are on CRAN.
hth
Matthias
- original Nachricht
Betreff: [R] Generic code for simulating from a distribution.
Gesendet: Mo, 10. Apr 2006
Von: [EMAIL PROTECTED]
Hello all,
I have the code below to simulate samples of certain size from a
particular distribution (here,beta distribution) and compute some
statistics for the samples.
betasim2-function(nsim,n,alpha,beta)
{
sim-matrix(rbeta(nsim*n,alpha,beta),ncol=n)
xmean-apply(sim,1,mean)
xvar-apply(sim,1,var)
xmedian-apply(sim,1,median)
simset-data.frame(sampleno=seq1:nsim),means=xmean,vars=xvar,medians=xmedian
)
return(simset)
}
I can write a similar coding for any distribution individually.
Now, I would like to have a generic code, say if I specify the
distribution with the parameters and simulation and sample size I would
like to have the simulations done for the mentioned distribution and the
statistics performed.
I would appreciate any help in doing so?
Thanks for your time.
Mathangi
Mathangi Gopalakrishnan
Graduate student
Dept of Mathematics and Statistics
University of Maryland, Baltimore County
Baltimore, MD
__
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
--- original Nachricht Ende
__
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] discrete probability distribution
Hi, have a look at http://mathworld.wolfram.com/GeometricDistribution.html respectively http://mathworld.wolfram.com/NegativeBinomialDistribution.html with r = 1. In R have a look at ?rnbinom with n = 1 and in your case: prob = 1-p hth Matthias Hi I want to sample from a discrete random variable X, defined on the non-negative integers, with prob(X=0) = (1-p) prob(X=1) = (1-p)*p ... prob(X=i)=(1-p)*p^i ... Before reinventing the wheel, is there a ready-made R function to give me such a thing? -- 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 -- StaMatS - Statistik + Mathematik Service Dr. rer. nat. Matthias Kohl Gottlieb-Keim-Straße 60 95448 Bayreuth Phone: +49 921 50736457 E-Mail: [EMAIL PROTECTED] Home: www.stamats.de __ 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] New user: Custom probability distribution
Hi, you can use our package distr respectively distrEx. require(distrEx) D1 - DiscreteDistribution(supp=0:3, prob = 12/(25*(1:4))) plot(D1) r(D1)(50) hth Matthias Hello, Given a probability function: p(x) = 12 / (25(x+1)) , x=0, 1, 2, 3 we generate the following values: C1 C2 0 0.48 1 0.24 2 0.16 3 0.12 Now, I'm supposed to create 50 random values using this table. In MiniTab, I simply selected Calc - Random Data - Discrete, and selected the columns, and it created 50 random values in a new column.[1] How do I do the same thing in R? [1] You may be wondering why I'm telling you this. Well, it's because if I were in your shoes, I'd think Oh, someone wants me to solve his homework.. I have already solved it using MiniTab, but I want to be able to use R instead (since I prefer NetBSD). __ 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 -- StaMatS - Statistik + Mathematik Service Dr. rer. nat. Matthias Kohl Gottlieb-Keim-Straße 60 95448 Bayreuth Phone: +49 921 50736457 E-Mail: [EMAIL PROTECTED] Home: www.stamats.de __ 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] Extending a data frame with S4
the help page on setOldClass might help you. In particular the section
Register or Convert?.
Matthias
hadley wickham schrieb:
I'm trying to create an extension to data.frame with more complex row
and column names, and have run into some problems:
setClass(new-data.frame, representation(data.frame))
[1] new-data.frame
Warning message:
old-style ('S3') class data.frame supplied as a superclass of
new-data.frame, but no automatic conversion will be peformed for S3
classes in: .validDataPartClass(clDef, name)
Do I need to be worried about this?
new(new-data.frame, data.frame())
Error in initialize(value, ...) : initialize method returned an object
of class data.frame instead of the required class new-data.frame
I guess this is related to the warning above. I presume I can fix
this with an initialize function, but I'm not sure how to go about
referring to the data frame that is the object.
Is there a way to extend a data.frame, or do I need to create an
object that contains the data frame in a slot?
Thanks for your help,
Hadley
__
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
--
StaMatS - Statistik + Mathematik Service
Dipl.Math.(Univ.) Matthias Kohl
www.stamats.de
__
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] convolution of the double exponential distribution
Hi David,
you can even increase the precision of these computations by changing
the variables
DefaultNrFFTGridPointsExponent (default: 12 - 2^12 grid points) and
TruncQuantile
(default: 1e-5) in our package distr; see
distroptions()
You can for instance use
distroptions(DefaultNrFFTGridPointsExponent, 14)
and
distroptions(TruncQuantile, 1e-8)
We checked our algorithm in case of Binomial, Poisson, Normal and
Exponential distributions and found a very high precision; confer
Section 5 of
http://www.stamats.de/comp.pdf
Moreover, for n-fold convolutions you can also use the function
convpow which can be found under
http://www.uni-bayreuth.de/departments/math/org/mathe7/DISTR/nFoldConvolution.R
hth,
Matthias
Ravi, Duncan, and Matthias,
Thank you very much for the helpful replies. For convolutions with 2 or
3 copies, I found that the CDFs from the distr package closely match the
analytic results from this paper:
K. Singh, M. Xie, and W. E. Strawderman, Combining Information From
Independent Sources Through Confidence Distributions, Annals of
Statistics 33, no. 1 (2005): 159-183.
That gives me confidence that the package will also work well for higher
copy numbers. At least for me, using the package is much more convenient
than programming all the needed integrals into R.
David
___
David R. Bickel http://davidbickel.com
Research Scientist
Pioneer Hi-Bred International (DuPont)
Bioinformatics and Exploratory Research
7200 NW 62nd Ave.; PO Box 184
Johnston, IA 50131-0184
515-334-4739 Tel
515-334-4473 Fax
[EMAIL PROTECTED], [EMAIL PROTECTED]
| -Original Message-
| From: Matthias Kohl [mailto:[EMAIL PROTECTED]
| Sent: Friday, December 23, 2005 9:09 AM
| To: Bickel, David
| Cc: Duncan Murdoch; r-help@stat.math.ethz.ch
| Subject: Re: [R] convolution of the double exponential distribution
|
| Duncan Murdoch schrieb:
|
| On 12/22/2005 7:56 PM, Bickel, David wrote:
|
|
| Is there any R function that computes the convolution of the double
| exponential distribution?
|
| If not, is there a good way to integrate ((q+x)^n)*exp(-2x)
| over x from
| 0 to Inf for any value of q and for any positive integer n?
| I need to
| perform the integration within a function with q and n as
| arguments. The
| function integrate() is giving me this message:
|
| evaluation of function gave a result of wrong length
|
|
|
| Under the substitution of y = q+x, that looks like a gamma integral.
| The x = 0 to Inf range translates into y = q to Inf, so
| you'll need an
| incomplete gamma function, such as pgamma. Be careful to get the
| constant multiplier right.
|
| 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
|
|
|
| Hi,
|
| you can use our package distr.
|
| require(distr)
| ## define double exponential distribution
| loc - 0 # location parameter
| sca - 1 # scale parameter
|
| rfun - function(n){ loc + scale * ifelse(runif(n) 0.5, 1,
| -1) * rexp(n) }
| body(rfun) - substitute({ loc + scale * ifelse(runif(n)
| 0.5, 1, -1) *
| rexp(n) },
| list(loc = loc, scale = sca))
|
| dfun - function(x){ exp(-abs(x-loc)/scale)/(2*scale) }
| body(dfun) - substitute({ exp(-abs(x-loc)/scale)/(2*scale)
| }, list(loc
| = loc, scale = sca))
|
| pfun - function(x){ 0.5*(1 +
| sign(x-loc)*(1-exp(-abs(x-loc)/scale))) }
| body(pfun) - substitute({ 0.5*(1 +
| sign(x-loc)*(1-exp(-abs(x-loc)/scale))) },
| list(loc = loc, scale = sca))
|
| qfun - function(x){ loc - scale*sign(x-0.5)*log(1 - 2*abs(x-0.5)) }
| body(qfun) - substitute({ loc - scale*sign(x-0.5)*log(1 -
| 2*abs(x-0.5)) },
| list(loc = loc, scale = sca))
|
| D1 - new(AbscontDistribution, r = rfun, d = dfun, p =
| pfun, q = qfun)
| plot(D1)
|
| D2 - D1 + D1 # convolution based on FFT
| plot(D2)
|
| hth,
| Matthias
|
| --
| StaMatS - Statistik + Mathematik Service
| Dipl.Math.(Univ.) Matthias Kohl
| www.stamats.de
|
This communication is for use by the intended recipient and ...{{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
--
StaMatS - Statistik + Mathematik Service
Dipl.Math.(Univ.) Matthias Kohl
Gottlieb-Keim-Straße 60
95448 Bayreuth
Phone: +49 921 50736 457
E-Mail: [EMAIL PROTECTED]
www.stamats.de
__
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] convolution of the double exponential distribution
Duncan Murdoch schrieb:
On 12/22/2005 7:56 PM, Bickel, David wrote:
Is there any R function that computes the convolution of the double
exponential distribution?
If not, is there a good way to integrate ((q+x)^n)*exp(-2x) over x from
0 to Inf for any value of q and for any positive integer n? I need to
perform the integration within a function with q and n as arguments. The
function integrate() is giving me this message:
evaluation of function gave a result of wrong length
Under the substitution of y = q+x, that looks like a gamma integral.
The x = 0 to Inf range translates into y = q to Inf, so you'll need an
incomplete gamma function, such as pgamma. Be careful to get the
constant multiplier right.
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
Hi,
you can use our package distr.
require(distr)
## define double exponential distribution
loc - 0 # location parameter
sca - 1 # scale parameter
rfun - function(n){ loc + scale * ifelse(runif(n) 0.5, 1, -1) * rexp(n) }
body(rfun) - substitute({ loc + scale * ifelse(runif(n) 0.5, 1, -1) *
rexp(n) },
list(loc = loc, scale = sca))
dfun - function(x){ exp(-abs(x-loc)/scale)/(2*scale) }
body(dfun) - substitute({ exp(-abs(x-loc)/scale)/(2*scale) }, list(loc
= loc, scale = sca))
pfun - function(x){ 0.5*(1 + sign(x-loc)*(1-exp(-abs(x-loc)/scale))) }
body(pfun) - substitute({ 0.5*(1 +
sign(x-loc)*(1-exp(-abs(x-loc)/scale))) },
list(loc = loc, scale = sca))
qfun - function(x){ loc - scale*sign(x-0.5)*log(1 - 2*abs(x-0.5)) }
body(qfun) - substitute({ loc - scale*sign(x-0.5)*log(1 - 2*abs(x-0.5)) },
list(loc = loc, scale = sca))
D1 - new(AbscontDistribution, r = rfun, d = dfun, p = pfun, q = qfun)
plot(D1)
D2 - D1 + D1 # convolution based on FFT
plot(D2)
hth,
Matthias
--
StaMatS - Statistik + Mathematik Service
Dipl.Math.(Univ.) Matthias Kohl
www.stamats.de
__
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] Loading namespaces
BBK schrieb: Just noticed the mssing ) at the end of the setClass statement, it is there in the orginal Phineas -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of BBK Sent: Thursday, December 08, 2005 8:18 PM To: 'R-Help Subject: [R] Loading namespaces I'm creating a package for my own use that uses some S4 classes but no methods. I have a file called NAMESPACE it contains the line: exportClasses(foo) and at the top of the R file I have setClass(foo, representation(x=numeric) and the line: .onLoad-function(libname,pkgname) Do you mean .onLoad - function(lib, pkg) require(methods) as given in Section 1.6.6 of Writing R Extensions When I run R CMD check I get Syntax error in the only R file. If I comment out the .onLoad function I get a package/namespace load failed error. Are libname and pkgname parameters for the function .onLoad that need to explicitly stated, or does R populate them when the package is loaded? Does .onLoad as defined above do enough to ensure that the namesapce is loaded? see Section 1.6.6 (ibid.): There needs to be an .onLoad action to ensure that the methods package is loaded and attached hth Matthias Phineas Campbell version _ platform sparc-sun-solaris2.9 arch sparc os solaris2.9 system sparc, solaris2.9 status major2 minor1.0 year 2005 month04 day 18 language 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-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 -- StaMatS - Statistik + Mathematik Service Dipl.Math.(Univ.) Matthias Kohl www.stamats.de __ 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] Vectors of S4 Classes
Phineas schrieb: I have a function from which I wish to return two vectors of equal length of a class eg setClass(ClassOne, representation(x=numeric)) [1] ClassOne first-new(ClassOne, x=1) second-new(ClassOne,x=2) Do you want list(first = first, second = second) or something like this setClass(ClassOneList, contains = list) new(ClassOneList, list(first = first, second = second)) hth Matthias first-rbind(first,second) first first second Is it possible to create vector or list of an S4 class? Phineas version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor1.0 year 2005 month04 day 18 language 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 -- StaMatS - Statistik + Mathematik Service Dipl.Math.(Univ.) Matthias Kohl Gottlieb-Keim-Straße 60 95448 Bayreuth Phone: +49 921 50736 457 E-Mail: [EMAIL PROTECTED] www.stamats.de __ 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] solution of convolution equation
Anna Oganyan wrote: Hello, May be somebody can help me... I am trying to find a solution of a convolution equation using fft (and unfortunately I do not have a good background for this). So I am just trying to figure out how it can be implemented in R. I have two multidimensional independent variables X and Z and I know their densities fx and fz, which are multidimensional arrays. So I have to find the density of Y, such that X+Y=Z. So, first I tried on a simple example, where the variables are one-dimensional, say X is N(0,1) and Z is N(7,3). So I want to find the density of Y (which should be N(7, sqrt(8)). I did: x - seq(-6, 20, length=500) fx - dnorm(x) z - seq(-6, 20, length=500) fz - dnorm(z, mean=7, sd=3) ffty - fft(fz)/fft(fx) fy - fft(ffty, inverse=T)/length(ffty) plot(Re(fy)) But the plot is far from being normal. May be the problem is with the lengths of fx and fz? should they be of different lengths and fx padded with zeros? May be somebody could point out that…? 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 Hello Anna, in our R package distr (on CRAN) we have implemented a convolution algorithm via fft. See also: http://www.uni-bayreuth.de/departments/math/org/mathe7/KOHL/pubs/comp.pdf respectively library(distr) getMethod(+, signature=c(AbscontDistribution,AbscontDistribution)) (or see the sources) Unfortunatelly, we haven't implemented our algorithms for multidimensional distributions so far. hope that helps, Matthias __ 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] two-tailed exact binomail test
Peter Dalgaard wrote:
katrina smith [EMAIL PROTECTED] writes:
I am trying to find a definition for the two-tailed exact binomial
test but have been unsuccessful. Can you help?
Just read binom.test. The relevant bit is this:
(m is the mean == n*p)
else if (x m) {
i - seq(from = ceiling(m), to = n)
y - sum(dbinom(i, n, p) = d * relErr)
pbinom(x, n, p) + pbinom(n - y, n, p, lower = FALSE)
}
i.e. we take the lower tail, including the value observed + the part
of the upper tail where the binomial density is less than or equal to
that of x (with a little fuzz added in). Symmetrically for observations
in the upper tail of course.
If you were looking for an official definition of the two sided
exact test, I don't think one exists. R's version is equivalent to the
likelihood ratio test, but there are alternatives (tail-balancing,
doubling the one-sided p, and maybe more).
there is a reference:
Section 2.4.2 (Zweiseitige Tests in einparametrigen
Exponentialfamilien - two sided tests in one-parameter exponential
families) in
H. Witting (1985): Mathematische Statistik I. Teubner. Stuttgart
confer Satz 2.70, Korollar 2.73 (in case of symmetry)
and Beispiel 2.74 (application of Korollar 2.73 to binomial model for p=
0.5).
Matthias
__
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] rlm/M/MM-estimator questions
Christian Hennig wrote: Hi list, 1) How can the MM-estimator method=MM in function rlm be tuned to 85% efficiency? It seems that there is a default tuning to 95%. I presume, but am not sure, that the MM-estimator uses phi=phi.bisquare as default and the tuning constant could be set by adding a parameter c=... Is this true? Which value to use for 85%? (In principle I should be able to figure that out theoretically, but it would be much easier if somebody already knew the constant or a straightforward way to compute it.) Hi Christian, I have not calculated the efficiency myself ... But the thesis of Matias Salibian-Barrera (SB 2000) might help you to find the answer (cf. Chapter 4). See: http://mathstat.math.carleton.ca:16080/~matias/thesis.pdf As far as I understand the choice k0=1.548 is to obtain a breakdown point 0.5 whereas k0=1.988 leads to a breakdown point of 0.4 - at least in the location case; confer p. 60 of SB 2000. In the article Optimal robust $M$-estimates of location by Fraiman, Yohai and Zamar (Ann. Stat. 29(1): 194 - 223) which is, of course, concerned with the location case, the authors recommend to use k0=1.988 instead of k0=1.548 (cf. p. 206). Hope that helps! Matthias 2) The M-estimator with bisquare uses rescaled MAD of the residuals as scale estimator according to the rlm help page. Does this mean that a scale estimator is used which is computed from least squares residuals? Are M-estimator and residual scale estimator iterated until convergence of them both? (Does this converge?) Or what else? What does rescaled mean? Thank you, Christian *** NEW ADDRESS! *** Christian Hennig University College London, Department of Statistical Science Gower St., London WC1E 6BT, phone +44 207 679 1698 [EMAIL PROTECTED], www.homepages.ucl.ac.uk/~ucakche __ 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] User-defined random variable
Achim Zeileis wrote:
On Sat, 30 Apr 2005, Peter Dalgaard wrote:
Paul Smith [EMAIL PROTECTED] writes:
Dear All
I would like to know whether it is possible with R to define a
discrete random variable different from the ones already defined
inside R and generate random numbers from that user-defined
distribution.
Yes. One generic way is to specify the quantile function (as in
qpois() etc.) and do qfun(runif(N)).
If the support discrete but also finite, you can also use sample(), e.g.
sample(myset, N, replace = TRUE, prob = myprob)
Z
--
O__ Peter Dalgaard Blegdamsvej 3
c/ /'_ --- Dept. of Biostatistics 2200 Cph. N
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
__
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
Hi,
one can also use our R package distr to generate discrete random
variables. The subsequent code provides a function which generates an
object of class DiscreteDistribution based on a finite support supp.
If prob is missing all elements in supp are equally weighted.
Matthias
# generating function
DiscreteDistribution - function(supp, prob){
if(!is.numeric(supp))
stop('supp' is no numeric vector)
if(any(!is.finite(supp))) # admit +/- Inf?
stop(inifinite or missing values in supp)
len - length(supp)
if(missing(prob)){
prob - rep(1/len, len)
}else{
if(len != length(prob))
stop('supp' and 'prob' must have equal lengths)
if(any(!is.finite(prob)))
stop(inifinite or missing values in prob)
if(!identical(all.equal(sum(prob), 1,
tolerance = 2*distr::TruncQuantile), TRUE))
stop(sum of 'prob' has to be (approximately) 1)
if(!all(prob = 0))
stop('prob' contains values 0)
}
if(length(usupp - unique(supp)) len){
warning(collapsing to unique support values)
prob - as.vector(tapply(prob, supp, sum))
supp - sort(usupp)
len - length(supp)
}else{
o - order(supp)
supp - supp[o]
prob - prob[o]
}
if(len 1){
if(min(supp[2:len] - supp[1:(len - 1)]) distr::DistrResolution)
stop(grid too narrow -- change DistrResolution)
}
rfun - function(n){
sample(x = supp, size = n, replace = TRUE, prob = prob)
}
intervall - distr::DistrResolution / 2
supp.grid - as.numeric(matrix(
rbind(supp - intervall, supp + intervall), nrow = 1))
prob.grid - c(as.numeric(matrix(rbind(0, prob), nrow = 1)), 0)
dfun - function(x){ stepfun(x = supp.grid, y = prob.grid)(x) }
cumprob - cumsum(prob)
pfun - function(x){ stepfun(x = supp, y = c(0, cumprob))(x) }
qfun - function(x){ supp[sum(cumprobx)+1] }
return(new(DiscreteDistribution, r = rfun, d = dfun, p = pfun,
q = qfun, support = supp))
}
# some examples
supp - rpois(20, lambda=7) # some support vector
D1 - DiscreteDistribution(supp = supp) # prob is missing
r(D1)(10) # 10 random numbers of Distribution D1
support(D1) # support
d(D1)(support(D1)) # pdf
p(D1)(5) # cdf
q(D1)(0.5) # quantile (here: median)
plot(D1) # plot of pdf, cdf and quantile
D2 - lgamma(D1) # apply member of group generic Math
plot(D2)
D3 - D1 + Binom(size=10) # convolution with object of class
DiscreteDistribution
plot(D3)
D4 - D1 + Norm() # convolution with object of class AbscontDistribution
plot(D4)
__
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] Using R to illustrate the Central Limit Theorem
(Ted Harding) wrote:
On 21-Apr-05 [EMAIL PROTECTED] wrote:
Here's a bit of a refinement on Ted's first suggestion.
[ corrected from runif(M*k), N, k) to runif(N*k), N, k) ]
N - 1
graphics.off()
par(mfrow = c(1,2), pty = s)
for(k in 1:20) {
m - (rowMeans(matrix(runif(N*k), N, k)) - 0.5)*sqrt(12*k)
hist(m, breaks = FD, xlim = c(-4,4), main = k,
prob = TRUE, ylim = c(0,0.5), col = lemonchiffon)
pu - par(usr)[1:2]
x - seq(pu[1], pu[2], len = 500)
lines(x, dnorm(x), col = red)
qqnorm(m, ylim = c(-4,4), xlim = c(-4,4), pch = ., col = blue)
abline(0, 1, col = red)
Sys.sleep(1)
}
Very nice! (I can better keep up with it mentally, though, with
Sys.sleep(2) or Sys.sleep(3), which moght be better for classroom
demo).
One thing occurred to me, watching it: people might say Yes,
we can see how the distribution - Normal, nice and smooth,
especially in the tails and side-arms; but the peaks always look
a bit rough.
Which could be the cue for introducing SD(ni) = sqrt(E[ni]),
and the following hack of the above code seems to show this OK
in the rootograms:
N - 1
graphics.off()
par(mfrow = c(1,2), pty = s)
for(k in 1:20) {
m - (rowMeans(matrix(runif(N*k), N, k)) - 0.5)*sqrt(12*k)
hm - hist(m, breaks = FD, xlim = c(-4,4), main = k, plot=FALSE,
prob = TRUE, ylim = c(0,0.5), col = lemonchiffon)
hm$counts-sqrt(hm$counts) ;
plot(hm,xlim = c(-4,4),main = k,ylab=sqrt(Frequency))
pu - par(usr)[1:2]
x - seq(pu[1], pu[2], len = 500)
lines(x, sqrt(N*dnorm(x)*(hm$breaks[2]-hm$breaks[1])), col = red)
qqnorm(m, ylim = c(-4,4), xlim = c(-4,4), pch = ., col = blue)
abline(0, 1, col = red)
Sys.sleep(2)
}
(and also shows clearly how the tails of the sample move outwards
into the tails of the Normal, as in fact you expect from the finite
range of mean(runif(k)), especially initially: very visible for
k up to about 5, and not really settled down for k10).
Next stop: Hanging rootograms!
Best wishes,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 22-Apr-05 Time: 13:10:19
-- XFMail --
__
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
Hello,
here is another idea to illustrate the Central Limit Theorem which is
based on our package distr.
# Illustration of the Central Limit Theorem
# using package distr
require(distr)
CLT - function(Distr, n, sleep = 1){
# Distr: object of class AbscontDistribution
# n: iterations
# sleep: time to sleep
graphics.off()
par(mfrow = c(1,2))
# expectation of Distr
fun1 - function(x, Distr){x*d(Distr)(x)}
E - try(integrate(fun1, lower = q(Distr)(0), upper = q(Distr)(1),
Distr = Distr)$value,
silent = TRUE)
if(!is.numeric(E))
E - try(integrate(fun1, lower = q(Distr)(distr::TruncQuantile),
upper = q(Distr)(1-distr::TruncQuantile),
Distr = Distr)$value,
silent = TRUE)
# standard deviation of Distr
fun2 - function(x, Distr){x^2*d(Distr)(x)}
E2 - try(integrate(fun2, lower = q(Distr)(0), upper = q(Distr)(1),
Distr = Distr)$value,
silent = TRUE)
if(!is.numeric(E2))
E2 - try(integrate(fun2, lower = q(Distr)(distr::TruncQuantile),
upper = q(Distr)(1-distr::TruncQuantile),
Distr = Distr)$value,
silent = TRUE)
std - sqrt(E2 - E^2)
Sn - 0
N - Norm()
for(k in 1:n) {
Sn - Sn + Distr
Tn - (Sn - k*E)/(std*sqrt(k))
x - seq(-5,5,0.01)
dTn - d(Tn)(x)
ymax - max(1/sqrt(2*pi), dTn)
plot(x, d(Tn)(x), ylim = c(0, ymax), type = l, ylab =
densities, main = k, lwd = 4)
lines(x, d(N)(x), col = orange, lwd = 2)
plot(x, p(Tn)(x), ylim = c(0, 1), type = l, ylab = cdfs,
main = k, lwd = 4)
lines(x, p(N)(x), col = orange, lwd = 2)
Sys.sleep(sleep)
}
}
# some examples
distroptions(DefaultNrFFTGridPointsExponent, 13)
CLT(Distr = Unif(), n = 20, sleep = 0)
CLT(Distr = Exp(), n = 20, sleep = 0)
CLT(Distr = Chisq(), n = 20, sleep = 0)
CLT(Distr = Td(df = 5), n = 20, sleep = 0)
CLT(Distr = Beta(), n = 20, sleep = 0)
distroptions(DefaultNrFFTGridPointsExponent, 14)
CLT(Distr = Lnorm(), n = 20, sleep = 0)
__
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] usage and behavior of 'setIs'
Hello,
am I using 'setIs' in the correct way in the subsequent (artifical) example?
Do I have to specify explicit 'setAs' for 'list' and 'vector' or
should this work automatically, since getClass(List1) states
an explicit coerce also for these classes.
I'm working with R 2.0.0 Patched (2004-10-06) on windows 2000.
Thanks for your advice,
Matthias
# example
setClass(Class = List1, representation(List = list))
setClass(Class = List2, contains = list)
setIs(class1 = List1, class2 = List2,
coerce = function(obj){ new(List2, [EMAIL PROTECTED]) },
replace = function(obj, value){
[EMAIL PROTECTED] - value
})
getClass(List1)
# states explicit coerce for 'list' and 'vector'
getClass(List2)
L1 - new(List1, List = list(a))
# all TRUE
is(L1, List2)
is(L1, list)
is(L1, vector)
as(L1, List2) # works
# both return 'list()'
# why not a 'list' with entry a?
# Is there an additional 'setAs' needed?
as(L1, list)
as(L1, vector)
L2 - as(L1, List2)
as(L2, list) # works
as(L2, vector) # works
__
[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] usage and behavior of 'setIs'
Thank you,
Matthias
Hi Matthias,
A similar problem to yours (with one level of inheritance less) was
disccussed this month on the r-devel list.
You find an answer from JChambers here:
https://stat.ethz.ch/pipermail/r-devel/2004-October/030980.html
And yes specifying _setAs_ to each _setIs_ with the coerce and replace
is a _hack_ which is with this version of methods necessary when
inherting from Old Classes.
/E
[EMAIL PROTECTED] wrote:
Hello,
am I using 'setIs' in the correct way in the subsequent (artifical)
example?
Do I have to specify explicit 'setAs' for 'list' and 'vector' or
should this work automatically, since getClass(List1) states
an explicit coerce also for these classes.
I'm working with R 2.0.0 Patched (2004-10-06) on windows 2000.
Thanks for your advice,
Matthias
# example
setClass(Class = List1, representation(List = list))
setClass(Class = List2, contains = list)
setIs(class1 = List1, class2 = List2,
coerce = function(obj){ new(List2, [EMAIL PROTECTED]) },
replace = function(obj, value){
[EMAIL PROTECTED] - value
})
getClass(List1)
# states explicit coerce for 'list' and 'vector'
getClass(List2)
L1 - new(List1, List = list(a))
# all TRUE
is(L1, List2)
is(L1, list)
is(L1, vector)
as(L1, List2) # works
# both return 'list()'
# why not a 'list' with entry a?
# Is there an additional 'setAs' needed?
as(L1, list)
as(L1, vector)
L2 - as(L1, List2)
as(L2, list) # works
as(L2, vector) # works
__
[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
--
Dipl. bio-chem. Witold Eryk Wolski
MPI-Moleculare Genetic
Ihnestrasse 63-73 14195 Berlin
tel: 0049-30-83875219 __(_
http://www.molgen.mpg.de/~wolski \__/'v'
http://r4proteomics.sourceforge.net||/ \
mail: [EMAIL PROTECTED]^^ m m
[EMAIL PROTECTED]
__
[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
__
[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] Gumbel distribution
please, take a look at package evd Matthias Does R have built in Gumbel distribution (pdf, ecdf, hazard, parameters estimation) for the minimum case? Thanks Anne [[alternative HTML version deleted]] __ [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 __ [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] setClassUnion
Hello,
I have a question concerning setClassUnion.
I'm working with R 2.0.0 Patched (2004-10-06) on windows 2000.
I tried to use setClassUnion in a package I am currently working on. The
situation is similar to the following example:
The DESCRIPTION file has entries:
Depends: R (= 2.0.0), methods
Imports: methods
LazyLoad: yes
The NAMESPACE file has entries:
importClassesFrom(methods, NULL, numeric)
exportClass(OptionalNumeric, class1, class2)
The example R code is:
.onLoad - function(lib, pkg){
require(methods, character = TRUE, quietly = TRUE)
}
setClassUnion(OptionalNumeric, c(numeric, NULL))
setClass(class1,
representation(test1 = OptionalNumeric),
prototype(test1 = numeric(1)))
# why does this not work?
# The error I get is:
# Error in makePrototypeFromClassDef(properties, ClassDef, immediate,
# where) :In making the prototype for class class1 elements of the
# prototype failed to match the corresponding slot class: test1
# (class OptionalNumeric )
# Sourcing this into R gives no error for me
# but instead using
prototype(test1 = NULL)
# works
# Moreover, using the second version (with test1 = NULL)
# the following works, too
setClass(class2,
representation(test2 = class1),
prototype(test2 = new(class1, test1 = numeric(1
What am I doing wrong?
Can someone please explain this to me?
Thanks for your help,
Matthias
__
[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] setClassUnion
sorry, please replace exportClass by exportClasses
Hello,
I have a question concerning setClassUnion.
I'm working with R 2.0.0 Patched (2004-10-06) on windows 2000.
I tried to use setClassUnion in a package I am currently working on.
The situation is similar to the following example:
The DESCRIPTION file has entries:
Depends: R (= 2.0.0), methods
Imports: methods
LazyLoad: yes
The NAMESPACE file has entries:
importClassesFrom(methods, NULL, numeric)
exportClass(OptionalNumeric, class1, class2)
exportClasses(OptionalNumeric, class1, class2)
The example R code is:
.onLoad - function(lib, pkg){
require(methods, character = TRUE, quietly = TRUE)
}
setClassUnion(OptionalNumeric, c(numeric, NULL))
setClass(class1,
representation(test1 = OptionalNumeric),
prototype(test1 = numeric(1)))
# why does this not work?
# The error I get is:
# Error in makePrototypeFromClassDef(properties, ClassDef, immediate, #
where) :In making the prototype for class class1 elements of the #
prototype failed to match the corresponding slot class: test1
# (class OptionalNumeric )
# Sourcing this into R gives no error for me
# but instead using
prototype(test1 = NULL)
# works
# Moreover, using the second version (with test1 = NULL)
# the following works, too
setClass(class2,
representation(test2 = class1),
prototype(test2 = new(class1, test1 = numeric(1
What am I doing wrong?
Can someone please explain this to me?
Thanks for your help,
Matthias
__
[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
__
[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] How to pass strings to functions? [once once more, now With content I hope...]
Dear expeRts,
I fail to succesfully pass strings to functions. It comes down to the
observation that
plot(someVariable,anotherVariable)
works fine, but
x - someVariable
y - anotherVariable
plot(x,y)
does not.
Does this have something to do with the returned value of x being
/someVariable/ and not /someVariable/, i.e. without the quotation
marks? Is there any way to work around this?
Ultimately I'd like to make multiple graphs by looping throught the
values in vectors. Something like:
var-c(var1,var2...n)
for (v in var)
{
plot(var, x))
}
what about plot(get(v), x)?
I've looked around for help on this but am stuck.
Hope you can help,
Gijs Plomp
__
[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
[R] questions about setMethod(Arith, ...)
Hi,
we have some questions concerning the definition of new arithmetic methods.
In our package distr (on CRAN) we define some new arithmetic methods
for +, -, *, /.
After loading distr the corresponding arithmetic methods work. Now, if
we define a new
class and also a new method for one of the arithmetic methods +, -,
*, /, still everything works fine. But, if we then define new
methods for the whole group Arith, the arithmetic methods of distr
no longer work. (for more details see example code below)
What are we missing?
Moreover, there is a certain precendence; i.e., if one defines a single
arithmetic method (e.g., /) and alterwards defines a method for the
whole group Arith, the old method / remains valid.
However, if we first define a method for the whole group Arith and
afterwards define a new single arithmetic method (e.g., +) the new one
is valid. (for more details see example code below).
Is this intended?
Thanks for your help,
Matthias, Thomas
###
## Example code
###
require(distr)
getMethods(/) # shows the corresponding methods of distr
## now define a new class track (see Chambers (1998))
## and define /
setClass(track, representation(x = numeric, y = numeric))
setMethod(/, signature(track, numeric),
function(e1, e2){ [EMAIL PROTECTED] = [EMAIL PROTECTED]/e2; e1 })
getMethods(/) # shows the corresponding methods
# of distr and class track
(N - Norm()) # creates an object of standard normal distribution
(N1 - N/3) # works
(tr - new(track, x = 1:3, y = 4:6))
(tr1 - tr/3) # works
## now define new methods for Arith
setMethod(Arith, signature(track, numeric),
function(e1, e2){
[EMAIL PROTECTED] = callGeneric([EMAIL PROTECTED], e2)
e1
})
getMethods(/) # / for distr is lost
N2 - N/3 # fails
(tr2 - tr/3) # works, but still the old method
tr + 2 # works
## now a new method +
setMethod(+, signature(track, numeric),
function(e1, e2){ [EMAIL PROTECTED] = [EMAIL PROTECTED]; e1 })
tr + 2 # works, but with the new method
__
[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] Fw: Evaluating strings as variables
Robin Gruna [EMAIL PROTECTED] writes: I have the following problem: I have a list as follows, values - list(red = 1, yellow = 2, blue = 3) values $red [1] 1 $yellow [1] 2 $blue [1] 3 There is also a vector containing the diffrent colors as character strings: colors - c(red, red, blue, yellow, red, blue) colors [1] redredblue yellow redblue Now i can attach the list values to R: attach(values) red [1] 1 etc... Now to my problem: How can I make R in a simple way to evaluate the strings red, blue etc. as variables, returning their numeric values ? As result I want to get a vector containing the values of the colors like this one: values.colors [1] 1 1 3 2 1 3 Forget about the attach() business, and do values - unlist(values) # or values - c(red = 1, yellow = 2, blue = 3) values.colors - values[colors] If you insist on going via variables, try sapply(colors, get) or? unlist(mget(colors, envir = as.environment(-1), inherits = TRUE)) Matthias -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ [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
[R] setValidity() changes Extends?
Hi,
I'm using Version 1.9.0 (2004-04-12) on Windows NT/98/2000 and found
the following difference between using setClass(..., valdity = ),
respectively using setValidity() afterwards:
setValidity() changes the Extends-part of a derived class, is this
intended or a bug or am I missing something?
##
## Example code
##
setClass(Class1, representation(name = character))
if(!isGeneric(name)) setGeneric(name, function(object)
standardGeneric(name))
setMethod(name, Class1, function(object) [EMAIL PROTECTED])
setClass(Class2, representation(Class1))
setClass(Class3, representation(Class2))
getClass(Class3) # as I expected
#Slots:
#
#Name: name
#Class: character
#
#Extends:
#Class Class2, directly
#Class Class1, by class Class2
validClass3 - function(object){TRUE}
setValidity(Class3, validClass3)
#Slots:
#
#Name: name
#Class: character
#
#Extends: Class2 # has been changed???
getClass(Class3) # why does setValidity change Extends?
# am I missing something?
# This doesn't happen if I use
# setClass(..., validity = )
# It, of course, also works if I
explicitly use
# setClass(Class3, contains =
c(Class2, Class1)
C32 - new(Class3)
name(C32) # generates an error?!
#Error in name(C32) : No direct or inherited method for function name
for this call
is(C32, Class1) # however
#TRUE
__
[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] privileged slots
Martin Maechler schrieb:
Matthias == Matthias Kohl [EMAIL PROTECTED]
on Thu, 27 May 2004 14:01:51 +0100 writes:
Matthias Hi all, in the help for RClassUtils I found the
Matthias expression privileged slots in function
Matthias checkSlotAssignment with the explanation:
Matthias /privileged slots (those that can only be set by
Matthias accesor functions defined along with the class
Matthias itself)/
RClassUtils ???
help.search(RClassUtils)
your right, sorry
but, at least a R Site search in Functions
gives me one match: Utilities for Managing Class Definitions
which hast the title: RClassUtils{methods} R Documentation
No help files found with alias or concept or title matching
'RClassUtils' using fuzzy matching.
-
So I guess that's not something in a standard R document.
You should rather keep to the 'official documentation' ...
I thought this is a official documentation ...
Matthias I thought all slots of a (not private) class can
Matthias be a accessed by a user via the @ Operator.
I tend to agree with your thoughts...
Matthias Is there a way to make a single slot of a class (not
Matthias the whole class) private, so that you can access
Matthias this slot only via an accessor function (not via @)?
I'd rather guess not.
Matthias Thanks, for your help Matthias
Martin
__
[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
[R] privileged slots
Hi all, in the help for RClassUtils I found the expression privileged slots in function checkSlotAssignment with the explanation: /privileged slots (those that can only be set by accesor functions defined along with the class itself)/ I thought all slots of a (not private) class can be a accessed by a user via the @ Operator. Is there a way to make a single slot of a class (not the whole class) private, so that you can access this slot only via an accessor function (not via @)? Thanks, for your help Matthias __ [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] (offtopic) I need two sets of 5 different color scales
maybe, the RColorBrewer package does what you want? see also: ColorBrewer.org Hi, I am plotting a policy function (result from a dynamic stochastic optimization problem, discretized approximation). The policy function maps from an 2 x 2 x 2 x 3 x B x F state space to a B x F state space (B and F are usually between 4-6, and represent domestic and foreign savings. The other variables are income (Y), inflation (Pi), domestic and foreign interest rates (R and Z)). I actually wrote a plotting function to represent all this, the result is attached -- please have a look at it and help me... I need advice in the following: I need two sets of colors for B and F which are easy to distinguish (when printed on a color laser printer), represent cardinality (ie have an intuitive mapping to an interval) or at least ordinality. I have experimented with the following: Bcolors - hsv(.6, seq(0.2, 1, length=5), 1) Fcolors - hsv(seq(.1,0, length=5), seq(0.2, 1, length=5) this is what you see in the plot. What colors would you use? Do you think that varying both brightness and hue helps to distinguish colors? Should I change saturation, too? Thanks, Tamas PS.: The plot is simply gzipped. If you need a zipped version, or the source code, contact me. -- Tamás K. Papp E-mail: [EMAIL PROTECTED] Please try to send only (latin-2) plain text, not HTML or other garbage. __ [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] Exactness of ppois
Hello,
by checking the precision of a convolution algorithm, we found the
following inexactness:
We work with R Version 1.8.1 (2003-11-21) on Windows systems (NT, 2000,
XP).
Try the code:
## Kolmogorov distance between two methods to
## determine P(Poisson(lambda)=x)
Kolm.dist - function(lam, eps){
x - seq(0,qpois(1-eps, lambda=lam), by=1)
max(abs(ppois(x, lambda=lam)-cumsum(dpois(x, lambda=lam
}
erg-optimize(Kolm.dist, lower=900, upper=1000, maximum=TRUE, eps=1e-15)
erg
Kolm1.dist - function(lam, eps){
x - seq(0,qpois(1-eps, lambda=lam), by=1)
which.max(abs(ppois(x, lambda=lam)-cumsum(dpois(x, lambda=lam
}
Kolm1.dist(lam=erg$max, eps=1e-15)
So for lambda=977.8 and x=1001 we get a distance of about 5.2e-06.
(This inexactness seems to hold for all lambda values greater than about
900.)
BUT, summing about 1000 terms of exactness around 1e-16,
we would expect an error of order 1e-13.
We suspect algorithm AS 239 to cause that flaw.
Do you think this could cause other problems apart from
that admittedly extreme example?
Thanks for your attention!
Matthias
__
[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] var on a vector
[EMAIL PROTECTED] schrieb: Hello I am an R newbie. I have a problem with computing a variance on a vector. data(cars) variance - function (x) mean(x^2)-mean(x)^2; variance(cars[,1]) [1] 27.4 var(cars[,1]) [1] 27.95918 What did I assume/understand wrong ? TIA Hello, help(var) says: The denominator n - 1 is used which gives an unbiased estimator of the (co)variance for i.i.d. observations. Try: n - length(cars[,1]) var(cars[,1])*(n-1)/n Matthias Kohl __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
