[R] [R-pkgs] new package on CRAN: multivator
Dear list, The new package 'multivator' is now available on CRAN. This presents a multivariate generalization of the emulator package. The corresponding JSS article is: Robin K. S. Hankin (2012), "Introducing multivator: A Multivariate Emulator", Journal of Statistical Software, 46(8), 1-20. URL http://www.jstatsoft.org/v46/i08/ best wishes Robin -- Robin Hankin Uncertainty Analyst hankin.ro...@gmail.com ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org 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] Trouble installing gsl wrapper
An interesting exchange. The package has some installation tips in Misc.Rd, but perhaps these are not sufficiently prominent. I think it'd be a good idea to include a READ.ME file in inst/doc with installation information, and include a pointer to this in the DESCRIPTION file; I'll adapt parts of this thread for the next release. Robin On Sat, Oct 30, 2010 at 8:14 AM, Gang Chen wrote: > You nailed it, Prof. Ripley! Thanks a lot... > > Gang > > On Sat, Oct 30, 2010 at 2:58 PM, Prof Brian Ripley > wrote: >> On Sat, 30 Oct 2010, Gang Chen wrote: >> >>> Hi, >>> >>> I'm trying to install the gsl wrapper source code >>> (http://cran.r-project.org/src/contrib/gsl_1.9-8.tar.gz) on a Linux >>> system (OpenSuse 11.1), but encountering the following problem. I've >>> already installed 'gsl' version 1.14 >>> (ftp://ftp.gnu.org/gnu/gsl/gsl-1.14.tar.gz) on the system. What's >>> missing? Thanks a lot... >> >> Installing the gsl library correctly? I need to guess quite a bit here (and >> a clean install attempt would have given a few more clues). >> >> I suspect you installed into /usr/local/lib whereas your OS probably uses >> /usr/local/lib64 (most x86_64 Linuxen do, and you seem to be using lib64 for >> R). In that case ld.so most likely will not find the dynamic library in >> /usr/local/lib. >> >> You can avoid such problems by installing auxiliary software such as gsl >> from RPMs -- I would be very surprised if OpenSuse did not have gsl and >> gsl-devel RPMs. Otherwise you need to install from the sources by something >> like >> >> make install libdir=/usr/local/lib64 >> >>> >>>> R CMD INSTALL gsl >>> >>> * installing to library ‘/usr/lib64/R/library’ >>> * installing *source* package ‘gsl’ ... >>> checking for gsl-config... /usr/local/bin/gsl-config >>> checking if GSL version >= 1.12... checking for gcc... gcc >>> checking for C compiler default output file name... a.out >>> checking whether the C compiler works... yes >>> checking whether we are cross compiling... no >>> checking for suffix of executables... >>> checking for suffix of object files... o >>> checking whether we are using the GNU C compiler... yes >>> checking whether gcc accepts -g... yes >>> checking for gcc option to accept ISO C89... none needed >>> yes >>> configure: creating ./config.status >>> config.status: creating src/Makevars >>> ** libs >>> make: Nothing to be done for `all'. >> >> This was not a clean install >> >>> installing to /usr/lib64/R/library/gsl/libs >>> ** R >>> ** inst >>> ** preparing package for lazy loading >>> ** help >>> *** installing help indices >>> ** building package indices ... >>> ** testing if installed package can be loaded >>> Error in dyn.load(file, DLLpath = DLLpath, ...) : >>> unable to load shared object '/usr/lib64/R/library/gsl/libs/gsl.so': >>> libgsl.so.0: cannot open shared object file: No such file or directory >>> ERROR: loading failed >>> * removing ‘/usr/lib64/R/library/gsl’ >>> >>>> sessionInfo() >>> >>> R version 2.12.0 (2010-10-15) >>> Platform: x86_64-unknown-linux-gnu (64-bit) >> >> >> -- >> Brian D. Ripley, rip...@stats.ox.ac.uk >> 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, UK Fax: +44 1865 272595 > -- Robin Hankin Uncertainty Analyst hankin.ro...@gmail.com __ R-help@r-project.org 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] boundary check
Feng
thanks for this. The problem you report is
reproducible; it originates from simplex()
of the boot packge.
It is ultimately due to
the fact that x_i is precisely *on* the convex hull,
which is evidently causing problems. I'll
investigate it.
In the short term, you can break the degeneracy:
> isin.chull(X[1,]+1e-10,X)
[1] FALSE
>
but I don't know if that is acceptable to you.
best
rksh
On 24/09/10 13:17, Feng Li wrote:
> Thanks. I agree with you that the speed and memory issues might be
> (actually is) a big problem for big dimensions. It is interesting to
> know to solve this by using linear programming. Buy the way, it seems
> a potential bug in your function if you try this
>
>
>
>> X <- matrix(rnorm(50), 10, 5)
>> x_i<-X[1,,drop=FALSE]
>> isin.chull(x_i,X)
>>
> Error in A.out[, basic] <- iden(M) : subscript out of bounds
>
> Or I must be wrong somewhere.
>
>
> Feng
>
>
> On Sep 24, 12:39 pm, Robin Hankin wrote:
>
>> Hello
>>
>> convex hulls in large numbers of dimensions are hard.
>>
>> For your problem, though, one can tell whether a given
>> point is inside or outside by using linear programming:
>>
>>
>>> X <- matrix(rnorm(50), 10, 5)
>>> x_i <- matrix(rnorm(5), 1, 5)
>>> isin.chull
>>>
>> function(candidate,p,plot=FALSE,give.answers=FALSE,
>> ...){
>> if(plot){
>> plot(p,...)
>> p(candidate[1],candidate[2], pch=16)
>> }
>> n <- nrow(p) # number of points
>> d <- ncol(p) # number of dimensions
>>
>> p <- t(sweep(p,2,candidate))
>> jj <- simplex(a=rep(1,n),A3=rbind(p,1),b3=c(0*candidate,1))
>> if(give.answers){
>> return(jj)
>> } else {
>> return((jj$solved >= 0) & all(jj$soln<1))
>> }
>>
>> }
>>
>>> isin.chull(x_i,X)
>>>
>> [1] FALSE
>>
>> (we can discuss offline; I'll summarize)
>>
>> HTH
>>
>> rksh
>>
>> On 24/09/10 10:44, Feng Li wrote:
>>
>>
>>
>>
>>> Dear R,
>>>
>>
>>> I have a covariates matrix with 10 observations, e.g.
>>>
>>
>>>> X <- matrix(rnorm(50), 10, 5)
>>>> X
>>>>
>>
>>> [,1][,2][,3][,4] [,5]
>>> [1,] 0.24857135 0.30880745 -1.44118657 1.10229027 1.0526010
>>> [2,] 1.24316806 0.36275370 -0.40096866 -0.24387888 -1.5324384
>>> [3,] -0.33504014 0.42996246 0.03902479 -0.84778875 -2.4754644
>>> [4,] 0.06710229 1.01950917 -0.09325091 -0.03222811 0.4127816
>>> [5,] -0.13619141 1.33143821 -0.79958805 2.08274102 0.6901768
>>> [6,] -0.45060357 0.19348831 -1.23793647 -0.72440163 0.5057326
>>> [7,] -1.20740516 0.20231086 1.15584485 0.8170 -1.2719855
>>> [8,] -1.81166284 -0.07913113 -0.91080581 -0.34774436 0.9552182
>>> [9,] 0.19131383 0.14980569 -0.37458224 -0.09371273 -1.7667203
>>> [10,] -0.85159276 -0.66679528 1.63019340 0.56920196 -2.4049600
>>>
>>
>>> And I define a boundary of X: The smallest "ball" that nests all the
>>> observations of X. I wish to check if a particular point x_i
>>>
>>
>>>> x_i <- matrix(rnorm(5), 1, 5)
>>>> x_i
>>>>
>>
>>>[,1] [,2] [,3] [,4] [,5]
>>> [1,] -0.1525543 0.4606419 -0.1011011 -1.557225 -1.035694
>>>
>>
>>> is inside the boundary of X or not. I know it's easy to do it with 1-D or
>>> 2-D, but I don't knot how to manage it when the dimension is large.
>>>
>>
>>> Can someone give a hint? Thanks in advance!
>>>
>>
>>> Feng
>>>
>> --
>> Robin K. S. Hankin
>> Uncertainty Analyst
>> University of Cambridge
>> 19 Silver Street
>> Cambridge CB3 9EP
>> 01223-764877
>>
>> __
>> r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] boundary check
Hello
convex hulls in large numbers of dimensions are hard.
For your problem, though, one can tell whether a given
point is inside or outside by using linear programming:
> X <- matrix(rnorm(50), 10, 5)
> x_i <- matrix(rnorm(5), 1, 5)
> isin.chull
function(candidate,p,plot=FALSE,give.answers=FALSE,
...){
if(plot){
plot(p,...)
p(candidate[1],candidate[2], pch=16)
}
n <- nrow(p) # number of points
d <- ncol(p) # number of dimensions
p <- t(sweep(p,2,candidate))
jj <- simplex(a=rep(1,n),A3=rbind(p,1),b3=c(0*candidate,1))
if(give.answers){
return(jj)
} else {
return((jj$solved >= 0) & all(jj$soln<1))
}
}
> isin.chull(x_i,X)
[1] FALSE
>
(we can discuss offline; I'll summarize)
HTH
rksh
On 24/09/10 10:44, Feng Li wrote:
> Dear R,
>
> I have a covariates matrix with 10 observations, e.g.
>
>
>> X <- matrix(rnorm(50), 10, 5)
>> X
>>
> [,1][,2][,3][,4] [,5]
> [1,] 0.24857135 0.30880745 -1.44118657 1.10229027 1.0526010
> [2,] 1.24316806 0.36275370 -0.40096866 -0.24387888 -1.5324384
> [3,] -0.33504014 0.42996246 0.03902479 -0.84778875 -2.4754644
> [4,] 0.06710229 1.01950917 -0.09325091 -0.03222811 0.4127816
> [5,] -0.13619141 1.33143821 -0.79958805 2.08274102 0.6901768
> [6,] -0.45060357 0.19348831 -1.23793647 -0.72440163 0.5057326
> [7,] -1.20740516 0.20231086 1.15584485 0.8170 -1.2719855
> [8,] -1.81166284 -0.07913113 -0.91080581 -0.34774436 0.9552182
> [9,] 0.19131383 0.14980569 -0.37458224 -0.09371273 -1.7667203
> [10,] -0.85159276 -0.66679528 1.63019340 0.56920196 -2.4049600
>
> And I define a boundary of X: The smallest "ball" that nests all the
> observations of X. I wish to check if a particular point x_i
>
>
>> x_i <- matrix(rnorm(5), 1, 5)
>> x_i
>>
>[,1] [,2] [,3] [,4] [,5]
> [1,] -0.1525543 0.4606419 -0.1011011 -1.557225 -1.035694
>
> is inside the boundary of X or not. I know it's easy to do it with 1-D or
> 2-D, but I don't knot how to manage it when the dimension is large.
>
> Can someone give a hint? Thanks in advance!
>
>
> Feng
>
>
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] simple table/matrix problem
Hi
Given three vectors
x <- c(fish=3, dogs=5, bats=2)
y <- c(dogs=1, hogs=3)
z <- c(bats=3, dogs=5)
How do I create a multi-way table like the following?
> out
x y z
bats 2 0 3
dogs 5 1 5
fish 3 0 0
hogs 0 3 0
('out' is a matrix).
See how the first line shows 'x' has 2 bats, 'y' has zero bats, and 'z'
has 3 bats
and so on for each line.
The real application would have a matrix of size ~10 by ~1.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] Function to compute the multinomial beta function?
It's usually better to build vectorization in to functions:
> beta3<- function (n1, n2, n3)
exp(lgamma(n1)+lgamma(n2)+lgamma(n3)-lgamma(n1+n2+n3))
> f <- function(x){exp(sum(lgamma(x))-lgamma(sum(x)))}
> beta3(5,3,8)
[1] 1.850002e-07
> f(c(5,3,8))
[1] 1.850002e-07
>
rksh
On 07/06/2010 01:54 AM, Robert A LaBudde wrote:
At 05:10 PM 7/5/2010, Gregory Gentlemen wrote:
Dear R-users,
Is there an R function to compute the multinomial beta function? That
is, the normalizing constant that arises in a Dirichlet distribution.
For example, with three parameters the beta function is
Beta(n1,n2,n2) = Gamma(n1)*Gamma(n2)*Gamma(n3)/Gamma(n1+n2+n3)
> beta3<- function (n1, n2, n3)
exp(lgamma(n1)+lgamma(n2)+lgamma(n3)-lgamma(n1+n2+n3))
> beta3(5,3,8)
[1] 1.850002e-07
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] Apply a shift which is a function of the array.
Hello again Jim
It seems that ashift() from the same package *doesn't* do what you want.
But you can use shift() as follows:
> myshift <- function(x){shift(x,1-which.max(x))}
> a <- matrix(runif(30),5,6)
> a
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.268955362 0.18446939 0.6489957 0.51000680 0.05877378 0.6671306
[2,] 0.007827401 0.03128131 0.7508847 0.27802680 0.93626913 0.2460778
[3,] 0.936496465 0.96065059 0.5452295 0.04976117 0.91768925 0.3285821
[4,] 0.244059311 0.31785472 0.7398967 0.31343310 0.67126086 0.2168343
[5,] 0.083498454 0.24745674 0.3611509 0.74341994 0.45215530 0.7454390
> t(apply(a,1,myshift))
[,1] [,2][,3] [,4] [,5] [,6]
[1,] 0.6671306 0.26895536 0.184469391 0.64899567 0.5100068 0.05877378
[2,] 0.9362691 0.24607779 0.007827401 0.03128131 0.7508847 0.27802680
[3,] 0.9606506 0.54522946 0.049761166 0.91768925 0.3285821 0.93649646
[4,] 0.7398967 0.31343310 0.671260864 0.21683431 0.2440593 0.31785472
[5,] 0.7454390 0.08349845 0.247456736 0.36115095 0.7434199 0.45215530
> apply(a,2,myshift)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.936496465 0.96065059 0.7508847 0.74341994 0.93626913 0.7454390
[2,] 0.244059311 0.31785472 0.5452295 0.51000680 0.91768925 0.6671306
[3,] 0.083498454 0.24745674 0.7398967 0.27802680 0.67126086 0.2460778
[4,] 0.268955362 0.18446939 0.3611509 0.04976117 0.45215530 0.3285821
[5,] 0.007827401 0.03128131 0.6489957 0.31343310 0.05877378 0.2168343
>
rksh
On 07/02/2010 12:05 PM, Jim Hargreaves wrote:
Dear List,
I have a 2,000x10,000 array of time domain data which when plotted
draws a distinct pulse. The matrix is 10,000 pulses of length 2000. I
would like the pulse to be shifted so that the peak (which.max of the
pulse data) is consistently at point 400.
I don't want to use loops as it'll take a long time. The problem with
using apply is the pulse peak varies for each pulse, so I want to use
something like:
apply(really_mean_pulse, 2, shift, which.max(really_mean_pulse))
(Using the shift function from the "magic" library)
This doesn't work, so I need to modify this so that the argument to
shift is dependent on the column the shift is being applied to.
Any suggestions would be greatly appreciated!
Kind Regards,
Jim Hargreaves
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] Apply a shift which is a function of the array.
Hello Jim you can use ashift() from the same library which does (I think) what you want. HTH, Robin On 07/02/2010 12:05 PM, Jim Hargreaves wrote: Dear List, I have a 2,000x10,000 array of time domain data which when plotted draws a distinct pulse. The matrix is 10,000 pulses of length 2000. I would like the pulse to be shifted so that the peak (which.max of the pulse data) is consistently at point 400. I don't want to use loops as it'll take a long time. The problem with using apply is the pulse peak varies for each pulse, so I want to use something like: apply(really_mean_pulse, 2, shift, which.max(really_mean_pulse)) (Using the shift function from the "magic" library) This doesn't work, so I need to modify this so that the argument to shift is dependent on the column the shift is being applied to. Any suggestions would be greatly appreciated! Kind Regards, Jim Hargreaves __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] table() of a factor
thanks everyone.
I think the motto should be "always specify the levels of a factor when
you create it
if you possibly can".
best wishes
Robin
On 06/29/2010 12:39 PM, Felix Andrews wrote:
Just use factor(), not levels(); you can pass a factor to factor() too.
x<- factor(c(rep("a",3),"b","d"), levels = letters[1:5])
table(x)
x
a b c d e
3 1 0 1 0
Cheers,
-Felix
On 29 June 2010 20:59, Robin Hankin wrote:
Hi
suppose I have a factor 'x':
x<- as.factor(c(rep("a",3),"b","d"))
table(x)
x
a b d
3 1 1
But this is not what I want because
I need to include the fact that the count of "c" is zero.
I can't just change the levels of x:
levels(x)<- c("a","b","c","d")
table(x)
x
a b c d
3 1 1 0
because this records the single "d" in the original 'x' as a "c".
What I want is:
a b c d
3 1 0 1
How to get this from 'x'?
(my real application has dozens of levels with complicated names).
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] table() of a factor
Hi
suppose I have a factor 'x':
> x <- as.factor(c(rep("a",3),"b","d"))
> table(x)
x
a b d
3 1 1
>
>
But this is not what I want because
I need to include the fact that the count of "c" is zero.
I can't just change the levels of x:
> levels(x) <- c("a","b","c","d")
> table(x)
x
a b c d
3 1 1 0
>
because this records the single "d" in the original 'x' as a "c".
What I want is:
a b c d
3 1 0 1
How to get this from 'x'?
(my real application has dozens of levels with complicated names).
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] add points to 3D plot using p3d {onion}
Hello Bradley I don't think there's an easy way to do what you want because the viewing angles are internal to p3d(). Frankly p3d() tries to be all things to all men (the arguments are a mess) and inevitably isn't as flexible as one might wish. I take it you want to do this: data(bunny) p3d(head(bunny,100),d0=2,theta=3) points(tail(bunny), col='blue') You'd want the call to points() to "remember" theta=3, and possibly d0=2 as well. Although I can see a hack I'd be very happy to help you offline. best wishes Robin Bradley Christoffersen wrote: Hi, Can anyone guide me as to how I can add points to a p3d() plot from the onion package? I want to plot points with different colors on the same 3D plot. Perhaps I can do this without adding points but somehow directing the 'h' parameter to give different color to points based on a factor I assign to them? FYI, I can do this using using scatterplot3d() and points3d(), but these plots lack perspective and hence it is hard to sense depth without the use of color. Thanks, Brad __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] Generating data from Null Distribution
Jim the 2x2 case is reasonably straightforward because the support is quite a small set. With the aylmer package you could do this: > a <- matrix(c(1,5,7,8),2,2) > sample(seq_along(allprobs(a)),100,replace=TRUE,prob=allprobs(a)) [1] 3 2 4 1 3 4 4 4 4 3 3 4 4 2 3 4 5 4 3 3 4 4 2 1 4 3 3 3 4 2 3 4 3 4 4 2 3 [38] 3 4 3 4 4 3 4 3 2 3 4 2 3 3 4 2 3 2 3 1 3 5 1 3 2 3 3 4 5 2 5 1 3 3 4 3 5 [75] 4 2 3 4 5 2 3 4 4 3 3 3 4 3 4 6 3 3 4 2 4 4 2 3 2 3 > For bigger tables the matter becomes quite complicated but you can use randomboards() of the aylmer package, which generates random samples using the Metropolis-Hastings algorithm. HTH Robin Jim Silverton wrote: Hello everyone, Can someone tell me exactly how to generate data from a null distribution for the fisher exact test? I know I have to use the hypergrometric but exactly what commands do I use? Jim [[alternative HTML version deleted]] __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] expand.grid game
Hi Ted.
you've found a bug in the documentation for blockparts().
It should read "0 <= ai <= yi". I'll fix it before the next
major release (which will include sampling without replacement from
a multiset, Insha'Allah).
Best wishes
rksh
(Ted Harding) wrote:
I wonder whether this answers Baptiste's question as asked.
1: An 8-digit number can have some digits equal to 0;
see Baptiste's comment "maxi <- 9 # digits from 0 to 9"
2: According to the man-page fror blockparts in partitions,
"all sets of a=(a1,...,an) satisfying Sum[ai] = n subject
to 0 < ai <= yi are given in lexicographical order."
So it would seem that blockparts would not count 8-digit numbers
which have some zero digits.
One could presumably "fake" it by looping over the number of
non-zero digits, from 2 to 8 -- something like:
all <- 0
for(i in (2:8)){
jj <- blockparts(rep(9,8),17)
all <- all + dim(jj)
}
Or am I missing something?!
Ted.
On 21-Dec-09 07:57:32, Robin Hankin wrote:
Hi
library(partitions)
jj <- blockparts(rep(9,8),17)
dim(jj)
gives 318648
HTH
rksh
baptiste auguie wrote:
Dear list,
In a little numbers game, I've hit a performance snag and I'm not sure
how to code this in C.
The game is the following: how many 8-digit numbers have the sum of
their digits equal to 17?
The brute-force answer could be:
maxi <- 9 # digits from 0 to 9
N <- 5 # 8 is too large
test <- 17 # for example's sake
sum(rowSums(do.call(expand.grid, c(list(1:maxi), rep(list(0:maxi),
N-1 == test)
## 3675
Now, if I make N = 8, R stalls for some time and finally gives up
with:
Error: cannot allocate vector of size 343.3 Mb
I thought I could get around this using Reduce() to recursively apply
rowSum to intermediate results, but it doesn't seem to help,
a=list(1:maxi)
b=rep(list(0:maxi), N-1)
foo <- function(a, b, fun="sum", ...){
switch(fun,
'sum' = rowSums(do.call(expand.grid, c(a, b))),
'mean' = rowMeans(do.call(expand.grid, c(a, b))),
apply(do.call(expand.grid, c(a, b)), 1, fun, ...)) # generic
case
}
sum(Reduce(foo, list(b), init=a) == test)
## 3675 # OK
Same problem with N=8.
Now, if N was fixed I could write a little C code to do this
calculation term-by-term, something along those lines,
test = 0;
for (i1=1, i1=9, i1++) {
for (i2=0, i2=9, i2++) {
[... other nested loops ]
test = test + (i1 + i2 + [...] == 17);
} [...]
}
but here the number of for loops, N, should remain a variable.
In despair I coded this in R as a wicked eval(parse()) construct, and
it does produce the expected result after an awfully long time.
makeNestedLoops <- function(N=3){
startLoop <- function(ii, start=1, end=9){
paste("for (i", ii, " in seq(",start,", ",end,")) {\n", sep="")
}
first <- startLoop(1)
inner.start <- lapply(seq(2, N), startLoop, start=0)
calculation <- paste("test <- test + (", paste("i", seq(1, N),
sep="", collapse="+"), "==17 )\n")
end <- replicate(N, "}\n")
code.to.run <- do.call(paste, c(list(first), inner.start,
calculation, end))
cat(code.to.run)
invisible(code.to.run)
}
test <- 0
eval(parse(text = makeNestedLoops(8)) )
## 229713
I hope I have missed a better way to do this in R. Otherwise, I
believe what I'm after is some kind of C or C++ macro expansion,
because the number of loops should not be hard coded.
Thanks for any tips you may have!
Best regards,
baptiste
E-Mail: (Ted Harding)
Fax-to-email: +44 (0)870 094 0861
Date: 21-Dec-09 Time: 08:45:09
-- XFMail --
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] expand.grid game
Hello again everybody.
I fired off my reply before reading the correspondence about
the leading zeros.
You can also assume that there is at least one block
at the leading position, [so that position can take 0,1,2,...,8 additional
blocks] and distribute the remaining 16 blocks amongst all 8
places:
system.time(dim(blockparts(c(8,rep(9,7)),16)))
user system elapsed
0.056 0.016 0.069
>
this is faster! :-)
best wishes
rksh
baptiste auguie wrote:
Wow!
system.time({
all = blockparts(rep(9,8),17)
print( dim(all[,all[1,]!=0])[2] ) # remove leading 0s
})
## 229713
user system elapsed
0.160 0.068 0.228
In some ways I think this is close to Hadley's suggestion, though I
didn't know how to implement it.
Thanks a lot to everybody who participated, I have learned interesting
things from a seemingly innocuous question.
Best regards,
baptiste
2009/12/21 Robin Hankin :
Hi
library(partitions)
jj <- blockparts(rep(9,8),17)
dim(jj)
gives 318648
HTH
rksh
baptiste auguie wrote:
Dear list,
In a little numbers game, I've hit a performance snag and I'm not sure
how to code this in C.
The game is the following: how many 8-digit numbers have the sum of
their digits equal to 17?
The brute-force answer could be:
maxi <- 9 # digits from 0 to 9
N <- 5 # 8 is too large
test <- 17 # for example's sake
sum(rowSums(do.call(expand.grid, c(list(1:maxi), rep(list(0:maxi),
N-1 == test)
## 3675
Now, if I make N = 8, R stalls for some time and finally gives up with:
Error: cannot allocate vector of size 343.3 Mb
I thought I could get around this using Reduce() to recursively apply
rowSum to intermediate results, but it doesn't seem to help,
a=list(1:maxi)
b=rep(list(0:maxi), N-1)
foo <- function(a, b, fun="sum", ...){
switch(fun,
'sum' = rowSums(do.call(expand.grid, c(a, b))),
'mean' = rowMeans(do.call(expand.grid, c(a, b))),
apply(do.call(expand.grid, c(a, b)), 1, fun, ...)) # generic case
}
sum(Reduce(foo, list(b), init=a) == test)
## 3675 # OK
Same problem with N=8.
Now, if N was fixed I could write a little C code to do this
calculation term-by-term, something along those lines,
test = 0;
for (i1=1, i1=9, i1++) {
for (i2=0, i2=9, i2++) {
[... other nested loops ]
test = test + (i1 + i2 + [...] == 17);
} [...]
}
but here the number of for loops, N, should remain a variable.
In despair I coded this in R as a wicked eval(parse()) construct, and
it does produce the expected result after an awfully long time.
makeNestedLoops <- function(N=3){
startLoop <- function(ii, start=1, end=9){
paste("for (i", ii, " in seq(",start,", ",end,")) {\n", sep="")
}
first <- startLoop(1)
inner.start <- lapply(seq(2, N), startLoop, start=0)
calculation <- paste("test <- test + (", paste("i", seq(1, N),
sep="", collapse="+"), "==17 )\n")
end <- replicate(N, "}\n")
code.to.run <- do.call(paste, c(list(first), inner.start, calculation,
end))
cat(code.to.run)
invisible(code.to.run)
}
test <- 0
eval(parse(text = makeNestedLoops(8)) )
## 229713
I hope I have missed a better way to do this in R. Otherwise, I
believe what I'm after is some kind of C or C++ macro expansion,
because the number of loops should not be hard coded.
Thanks for any tips you may have!
Best regards,
baptiste
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] expand.grid game
Hi
library(partitions)
jj <- blockparts(rep(9,8),17)
dim(jj)
gives 318648
HTH
rksh
baptiste auguie wrote:
Dear list,
In a little numbers game, I've hit a performance snag and I'm not sure
how to code this in C.
The game is the following: how many 8-digit numbers have the sum of
their digits equal to 17?
The brute-force answer could be:
maxi <- 9 # digits from 0 to 9
N <- 5 # 8 is too large
test <- 17 # for example's sake
sum(rowSums(do.call(expand.grid, c(list(1:maxi), rep(list(0:maxi),
N-1 == test)
## 3675
Now, if I make N = 8, R stalls for some time and finally gives up with:
Error: cannot allocate vector of size 343.3 Mb
I thought I could get around this using Reduce() to recursively apply
rowSum to intermediate results, but it doesn't seem to help,
a=list(1:maxi)
b=rep(list(0:maxi), N-1)
foo <- function(a, b, fun="sum", ...){
switch(fun,
'sum' = rowSums(do.call(expand.grid, c(a, b))),
'mean' = rowMeans(do.call(expand.grid, c(a, b))),
apply(do.call(expand.grid, c(a, b)), 1, fun, ...)) # generic case
}
sum(Reduce(foo, list(b), init=a) == test)
## 3675 # OK
Same problem with N=8.
Now, if N was fixed I could write a little C code to do this
calculation term-by-term, something along those lines,
test = 0;
for (i1=1, i1=9, i1++) {
for (i2=0, i2=9, i2++) {
[... other nested loops ]
test = test + (i1 + i2 + [...] == 17);
} [...]
}
but here the number of for loops, N, should remain a variable.
In despair I coded this in R as a wicked eval(parse()) construct, and
it does produce the expected result after an awfully long time.
makeNestedLoops <- function(N=3){
startLoop <- function(ii, start=1, end=9){
paste("for (i", ii, " in seq(",start,", ",end,")) {\n", sep="")
}
first <- startLoop(1)
inner.start <- lapply(seq(2, N), startLoop, start=0)
calculation <- paste("test <- test + (", paste("i", seq(1, N),
sep="", collapse="+"), "==17 )\n")
end <- replicate(N, "}\n")
code.to.run <- do.call(paste, c(list(first), inner.start, calculation, end))
cat(code.to.run)
invisible(code.to.run)
}
test <- 0
eval(parse(text = makeNestedLoops(8)) )
## 229713
I hope I have missed a better way to do this in R. Otherwise, I
believe what I'm after is some kind of C or C++ macro expansion,
because the number of loops should not be hard coded.
Thanks for any tips you may have!
Best regards,
baptiste
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] Fishers exact test at < 2.2e-16
The aylmer package has some functionality in this regard which you may find useful. In particular, you can use good() to get a feel for the number of tableaux that are consistent with the specified marginal totals: > good(dat2) [1] 42285210 > good(dat3) [1] 2.756286e+12 > HTH rksh Søren Faurby wrote: In an effort to select the most appropriate number of clusters in a mixture analysis I am comparing the expected and actual membership of individuals in various clusters using the Fisher?s exact test. I aim for the model with the lowest possible p-value, but I frequently get p-values below 2.2e-16 and therefore does not get exact p-values with standard Fisher?s exact tests in R. Does anybody know if there is a version of Fisher?s exact test in any package which can handle lower probabilities, or have other suggestions as to how I can compare the probabilities? I am for instance comparing the following two: dat2<-matrix(c(29,0,29,0,12,0,18,0,0,29,0,16,0,19), nrow=2) fisher.test(dat2, workspace=3000) dat3<-matrix(c(29,0,0,29,0,0,12,0,0,17,0,1,0,29,0,0,15,1,0,0,19), nrow=3) fisher.test(dat3, workspace=3000) Which both result in p-value < 2.2e-16 Kind regards, Søren __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to creat a matrix
Hi try R> library(magic) R> ashift(diag(5),1) HTH rksh enrique Dallazuanna wrote: Try this: N <- 5 diag(1, N)[c(N, 1:(N - 1)),] On Fri, Dec 11, 2009 at 1:47 PM, Moohwan Kim wrote: Dear R family I am attempting to create a matrix. e.g., 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 How could I write a R program? Later I want to extend it to a N by N case. Thanks in advance best Moohwan __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] access elements of a named list using a factor
Hello Dimitris
thanks for this. It works! I guess I was fixated on the dollar sign.
I must confess that I don't really understand any of the error
messages below. Can anyone help me interpret them?
rksh
Dimitris Rizopoulos wrote:
do you mean:
f <- factor(c("pigs", "pigs", "slugs"))
jj <- list(pigs = 1:10, slugs = 1:3)
jj[levels(f)[1]]
jj[[levels(f)[1]]]
Best,
Dimitris
Robin Hankin wrote:
Hi
I have a factor 'f' and a named list 'jj'.
I want names(jj) to match up with levels(f).
How do I use levels(f) to access elements of jj?
> f <- factor(c("pigs","pigs","slugs"))
> f
[1] pigs pigs slugs
Levels: pigs slugs
>
> jj <- list(pigs=1:10,slugs=1:3)
My attempts to produce jj$pigs all give errors:
> jj$levels(f)[1]
Error: attempt to apply non-function
> do.call("$",jj,levels(f)[1])
Error in if (quote) { : argument is not interpretable as logical
> "$"(jj,levels(f)[1])
Error in jj$levels(f)[1] : invalid subscript type 'language'
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] access elements of a named list using a factor
Hi
I have a factor 'f' and a named list 'jj'.
I want names(jj) to match up with levels(f).
How do I use levels(f) to access elements of jj?
> f <- factor(c("pigs","pigs","slugs"))
> f
[1] pigs pigs slugs
Levels: pigs slugs
>
> jj <- list(pigs=1:10,slugs=1:3)
My attempts to produce jj$pigs all give errors:
> jj$levels(f)[1]
Error: attempt to apply non-function
> do.call("$",jj,levels(f)[1])
Error in if (quote) { : argument is not interpretable as logical
> "$"(jj,levels(f)[1])
Error in jj$levels(f)[1] : invalid subscript type 'language'
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] reference on fisher.test()
Hi
fexact.c points you to the original ACM paper:
/*
ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM.
THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488.
-
You may find the discussion in the vignette("fishervig")
in the aylmer package helpful.
HTH
Robin
Peter Dalgaard wrote:
Peng Yu wrote:
On Thu, Oct 15, 2009 at 4:19 PM, RICHARD M. HEIBERGER
wrote:
On Thu, Oct 15, 2009 at 4:56 PM, Peng Yu wrote:
Can somebody point me a book on Fisher's exact test? I looked a few
webpages. But the descriptions on the webpages are not very complete.
Is there a book on that covers all the aspect of Fisher's exact test
that is implemented in R?
Section 15.2 of my book (Statistical Analysis and Data Display, with
Burt Holland and published by Springer)
shows a detailed example.
It doesn't mention odd ratio.
The general idea of basing the inference on the noncentral
hypergeometric distribution is something I have first seen in
Breslow&Day's famous 1980 book on case-control studies, including the
fact that the conditional MLE differs from the ordinary OR. (I'm sure
there's an earlier reference, but I happened to be a grad student when
that book came out...)
The rest of what R does is "carbon copied" from similar procedures for
the binomial distribution. I wouldn't know what kind of book to look
for for that sort of minutiae. Alan Agresti is a possible source.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] generalization of tabulate()
Hi
I want a generalization of tabulate() which works on rows of a matrix.
Suppose I have an integer matrix 'observation':
> observation
y1 y2 y3
1 4 0
1 4 0
2 0 3
4 1 0
0 5 0
0 1 4
2 0 3
Each row corresponds to a (multivariate) observation. Note that the
first two rows are identical: this means that data "c(1,4,0)" was
observed twice.
Now suppose I can list the sample space:
> S
[1,] 5 0 0
[2,] 4 1 0
[3,] 3 2 0
[4,] 2 3 0
[5,] 1 4 0
[6,] 0 5 0
[7,] 4 0 1
[8,] 3 1 1
[9,] 2 2 1
[10,] 1 3 1
[11,] 0 4 1
[12,] 3 0 2
[13,] 2 1 2
[14,] 1 2 2
[15,] 0 3 2
[16,] 2 0 3
[17,] 1 1 3
[18,] 0 2 3
[19,] 1 0 4
[20,] 0 1 4
[21,] 0 0 5
(thus each row corresponds to a point in my sample space).
Now what I need to do is to construct a new matrix, which uses the
'observation' matrix above, which is a sort of table:
> desired
y1 y2 y3 d
[1,] 5 0 0 0
[2,] 4 1 0 1
[3,] 3 2 0 0
[4,] 2 3 0 0
[5,] 1 4 0 2
[6,] 0 5 0 1
[7,] 4 0 1 0
[8,] 3 1 1 0
[9,] 2 2 1 0
[10,] 1 3 1 0
[11,] 0 4 1 0
[12,] 3 0 2 0
[13,] 2 1 2 0
[14,] 1 2 2 0
[15,] 0 3 2 0
[16,] 2 0 3 2
[17,] 1 1 3 0
[18,] 0 2 3 0
[19,] 1 0 4 0
[20,] 0 1 4 1
[21,] 0 0 5 0
Thus the 'd' column counts the number of times that each row occurs in
variable 'observation'. So desired[5,4]=2 because the observation
corresponding to desired[5,1:3] (viz c(1,4,0)) occurred twice. And
desired[1,4]=0 because the observation corresponding to desired[1,1:3]
(viz c(5,0,0)) did not occur once (it was not observed).
In my application I have dim(S) ~= c(5,4e6).
I've tried merge(), stack(), reshape(), but the best I can do
is the (derisory):
require(partitions)
obs <- matrix(as.integer(c(
1, 4, 0,
1, 4, 0,
2, 0, 3,
4, 1, 0,
0, 5, 0,
0, 1, 4,
2, 0, 3
)),ncol=3,byrow=TRUE)
S <- t(compositions(5,3))
d <- rep(0,nrow(S))
for(i in seq_len(nrow(obs))){
for(j in seq_len(nrow(S))){
if(all(obs[i,,drop=TRUE] == S[j,,drop=TRUE])){
d[j] <- d[j]+1
}
}
}
S <- cbind(S,d)
Anyone got anything better before I try C?
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] S4 tutorial
Peng the Brobdingnag package includes a vignette that gives a step-by-step guide to creating a simple package that uses S4. best wishes Robin Peng Yu wrote: I'm looking for some tutorial on S4. I only find the following one, which is not in English. Can somebody let me know if there is any introductory material? I'm very familiar with OO and C++. If there is some material that suits my background, it will be great. https://stat.ethz.ch/pipermail/r-help/2009-January/184108.html __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] sparse vectors
Hi guys thanks for this, it works fine, but I'm not sure the Matrix package does what I want: > a = sparseMatrix(i=c(20, 30, 10), j=rep(1, 3), x=c(2.2, 3.3, 4.4)) Error in asMethod(object) : Cholmod error 'out of memory' at file:../Core/cholmod_memory.c, line 148 Surely an efficient storage mechanism would need only six pieces of information? I've been pondering the solution that Henrique suggested, that uses merge(). This seems to be fine, although it might be possible to squeeze some efficiency gains by using the fact that the index vector is always sorted, which migh save some searching time. Any thoughts anyone? best wishes Robin Benilton Carvalho wrote: library(Matrix) a = sparseMatrix(i=c(20, 30, 1), j=rep(1, 3), x=c(2.2, 3.3, 4.4)) b = sparseMatrix(i=c(3, 30), j=rep(1, 2), x=c(0.1, 0.1), dims=dim(a)) theSum = a+b summary(theSum) hth, b On Sep 8, 2009, at 10:19 AM, Henrique Dallazuanna wrote: Try this: abMerge <- merge(a, b, by = 'index', all = TRUE) list(index = abMerge$index, val = rowSums(abMerge[,2:3], na.rm = TRUE)) On Tue, Sep 8, 2009 at 10:06 AM, Robin Hankin wrote: Hi I deal with long vectors almost all of whose elements are zero. Typically, the length will be ~5e7 with ~100 nonzero elements. I want to deal with these objects using a sort of sparse vector. The problem is that I want to be able to 'add' two such vectors. Toy problem follows. Suppose I have two such objects, 'a' and 'b': a $index [1]20 30 1 $val [1] 2.2 3.3 4.4 b $index [1] 3 30 $val [1] 0.1 0.1 What I want is the "sum" of these: AplusB $index [1]3 20 30 1 $val [1] 0.1 2.2 3.4 4.4 See how the value for index=30 (being common to both) is 3.4 (=3.3+0.1). What's the best R idiom to achieve this? -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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. -- Henrique Dallazuanna Curitiba-Paraná-Brasil 25° 25' 40" S 49° 16' 22" O [[alternative HTML version deleted]] -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] sparse vectors
Hi I deal with long vectors almost all of whose elements are zero. Typically, the length will be ~5e7 with ~100 nonzero elements. I want to deal with these objects using a sort of sparse vector. The problem is that I want to be able to 'add' two such vectors. Toy problem follows. Suppose I have two such objects, 'a' and 'b': > a $index [1]20 30 1 $val [1] 2.2 3.3 4.4 > b $index [1] 3 30 $val [1] 0.1 0.1 > What I want is the "sum" of these: > AplusB $index [1]3 20 30 1 $val [1] 0.1 2.2 3.4 4.4 > See how the value for index=30 (being common to both) is 3.4 (=3.3+0.1). What's the best R idiom to achieve this? -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] Goldbach partitions code
Hi interesting blog! not strictly relevant, but there are various number-theoretic functions implemented in the elliptic package which you might find useful. best wishes Robin murali.me...@fortisinvestments.com wrote: Folks, I put up a brief note describing my naive attempts to compute Goldbach partitions, starting with a brute-force approach and refining progressively. http://jostamon.blogspot.com/2009/02/goldbachs-comet.html I'd welcome your suggestions on improvements, alternatives, other optimisations, esp. to do with space vs time tradeoffs. Is this an example interesting enough for pedagogical purposes, do you think? Please advise. Cheers, MM [[alternative HTML version deleted]] __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R parser for If-else
I too have had many problems with if-else.
My solution is to always always always
use the line
"} else {"
in any if-else construction. This guarantees that
there won't be problems of the sort discussed here.
HTH
rksh
Martin Maechler wrote:
"D" == Dani
on Tue, 24 Feb 2009 14:09:36 -0800 writes:
D> Hi list,
D> I don't know if somebody has spent a lot of time debugging strange
D> problems with if else positioning - the parser seems to recognize only
D> the syntax bellow - this is the only way of making these pieces of
D> code to work.
D> As far as i'm concerned, no examples were available (it would be so
D> awesome to have them in the introductory manual!)
D> #Try to change the placement of the keywords and you are dead!
["dead"?]
Oh dear...
Note this has nothing to do with if( ) .. else ..
but indeed with how things are parsed.
I think this is FAQ (or should become one):
?if [the help page you really should read before spending too
much time or even post to R-help]
has the following section
> Note that it is a common mistake to forget to put braces ('{ .. }')
> around your statements, e.g., after 'if(..)' or 'for()'.
> In particular, you should not have a newline between '}' and
> 'else' to avoid a syntax error in entering a 'if ... else'
> construct at the keyboard or via 'source'. For that reason, one
> (somewhat extreme) attitude of defensive programming is to always
> use braces, e.g., for 'if' clauses.
Regards,
Martin Maechler, ETH Zurich
D> Ex1:
D> if (1==1){
D> print('if')
D> print('if again')
D> }else
D> print('else')
D> Ex2:
D> if (2==2) print('if') else print('else')
D> Ex3:
D> if (2==2){
D> print('if')
D> print('if again')
D> }else
D> {
D> print('else')
D> print('else2')
D> }
D> Ex4:
D> if (2==2){
D> print('if')
D> print('if again')
D> }else print('else')
D> cheers,
D> -
D> Daniela
D> __
D> R-help@r-project.org mailing list
D> https://stat.ethz.ch/mailman/listinfo/r-help
D> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
D> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to handle large numbers?
Feng checkout the Brobdingnag package: > library(Brobdingnag) > exp(1000)/(exp(1007)+5) [1] NaN > as.numeric(exp(as.brob(1000))/(exp(as.brob(1007))+5)) [1] 0.000911882 > Feng Li wrote: Dear R, I have two questions: 1, Why both R and Matlab give 0*Inf==NaN? To my knowledge, it should be zero mathematically. Am I right? 2, I need to calculate e.g. exp(a)/(exp(b)+c), where both a and b are very large numbers (>>1000, e.g a=1000, b=1007, and c=5). R gives me NaN when I use the following command: exp(1000)/(exp(1007)+5) [1] NaN I am pretty sure this should be close to zero. My question is whether there is a general way to solve this kind of question or should I do some settings before computing? Thanks in advance! Feng -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] pdf() and pch problems
Hi R-2.8.1, Suse 11.1 I'm having problems with pdf(). In the following transcript, file 'f.pdf' does not use the expected symbols for the plot. It uses a 'q' letter instead of the open circle I get when viewing the graphics window. I also get the same under r47678. Does anyone else get this? le112:~/scratch/R-2.8.1% R --vanilla --quiet > pdf(file='~/f.pdf') > plot(1:10 , pch=1) > dev.off() null device 1 > sessionInfo() R version 2.8.1 (2008-12-22) i686-pc-linux-gnu locale: LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base > q() le112:~/scratch/R-2.8.1% -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] reshape() problems
Hi I have a data frame with timeseries information like this: year cell Q1Q2 Q3 Q4 1940 1 1.2 1.4 1.41.9 1941 1 2.9 2.1 3.4 2.4 1942 1 2.7 3.2 1.52.6 1940 2 1.4 2.1 2.62.4 1941 2 2.4 1.4 1.43.4 1942 2 1.4 2.4 2.54.4 where the Qs mean 'quarter'. I want to extract from this a dataframe with a timeseries for each cell: year quarter cell1 cell2 1940 1 1.2 1.4 1940 2 1.4 2.1 1940 3 1.4 2.6 1940 4 1.9 2.4 1941 1 2.9 2.4 1941 2 2.1 1.4 1941 3 3.4 1.4 1942 4 2.4 3.4 1942 1 2.7 1.4 1942 2 3.2 2.4 1942 3 1.5 2.5 1942 4 2.6 4.4 Thus the third and fourth columns are the timeserieses for cell 1 and cell 2. Is there a nice vectorized way to do this? I can't quite make reshape() do what I want. [the real dataset is months, not quarters, has ~2000 cells and ~60 years] -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R package tests
I think the OP was asking about test suites that test the software.
The R package structure includes a test/ directory which you can use
to put tests.
For example, in the onion package I check that I have got my
signs and multiplication table correctly implemented:
stopifnot(Hi*Hj == Hk)
stopifnot(Hj*Hi == -Hk)
stopifnot(Hj*Hk == Hi)
stopifnot(Hk*Hj == -Hi)
stopifnot(Hk*Hi == Hj)
stopifnot(Hi*Hk == -Hj)
[and a whole lot of others]
and the elliptic package includes a whole bunch of code
that verifies identities that appear in AMS-55. It also
includes numerical verification that the functions,
using randomish arguments, match the output of
mathematica or maple.
Philippe Grosjean wrote:
There is a mechanism for testing code in R packages (R CMD check), see
the Writing R extensions manual. If you need more flexibility for your
tests, you could look at RUnit on CRAN, or svUnit on R-Forge
(http://r-forge.r-project.org, on CRAN soon). For the later one, you
install it by:
install.packages("svUnit", repos="http://R-Forge.R-project.org";)
These is a vignette associated with svUnit:
vignette("svUnit")
Note that RUnit and svUnit are "test suite code" compatible, but they
use very different mechanisms internally.
Best,
Philippe Grosjean
..<°}))><
) ) ) ) )
( ( ( ( (Prof. Philippe Grosjean
) ) ) ) )
( ( ( ( (Numerical Ecology of Aquatic Systems
) ) ) ) ) Mons-Hainaut University, Belgium
( ( ( ( (
..
Nathan S. Watson-Haigh wrote:
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1
I was wondering if anyone could point me in the right direction for
reading up on writing tests in
R. I'm writing some functions for inclusion into a package and would
like to test them to ensure
they're doing what I expect them to do.
Are these approaches used for testing packages in CRAN?
Cheers,
Nathan
- --
-
Dr. Nathan S. Watson-Haigh
OCE Post Doctoral Fellow
CSIRO Livestock Industries
Queensland Bioscience Precinct
St Lucia, QLD 4067
Australia
Tel: +61 (0)7 3214 2922
Fax: +61 (0)7 3214 2900
Web: http://www.csiro.au/people/Nathan.Watson-Haigh.html
-
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.9 (MingW32)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org
iEYEARECAAYFAkluio8ACgkQ9gTv6QYzVL5X9QCgwvg5xjwZW2A2Z5G41iADu1Kz
hIkAoI5ISuAtHyQ+JwJSRBAc9q/oyeEt
=cqm4
-END PGP SIGNATURE-
__
R-help@r-project.org 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@r-project.org 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@r-project.org 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] Construct All Possible Strings from 4 Bases (ATCG)
Gundala
f <- function(n){expand.grid(rep(list(seq_len(4)),n))}
HTH
Robin
Gundala Viswanath wrote:
Dear all,
Is there an efficient way in R to construct all strings from 4 bases (ATCG).
If we want a length L string, there are 4 ^ L possible strings of such.
e . g with L = 2 we have AA, AT, AC, AG, .. GC, GA, GT, GG as many as
4 ^ 2 = 16 strings,
with L = 3 we have as many as 4 ^ 3 = 64 strings
- Gundala Viswanath
Jakarta - Indonesia
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] Complex integration in R
Hi Borja
library(elliptic)
?myintegrate
HTH
rksh
Borja Soto Varela wrote:
Dear R-user
I need a function to approximate a complex integration. My function is:
aprox2=function(s,x,rate){
dexp(x,rate)*exp(-s*x)
}
where argument s is a complex number. I can't use the integrate function
because it's only used with "numeric" arguments
Does anyone know some function to approximate complex integrals?
Thanks
Borja
[[alternative HTML version deleted]]
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] for loop query
Hi
start simple!
Work out *each* row combined with *each* row,
to give (in your case) a 26-by-26 matrix.
Only after you have got this working, start thinking about
making it run faster [eg by
only evaluating the upper triangular entries]
To do a nested loop, do
M <- matrix(0,n,n)
for(i in seq_len(n)){
for(j in seq_len(n)){
M[i,j] <- f(i,j)
}
}
which will fill in matrix M for you.
HTH
rksh
Gerard M. Keogh wrote:
Hi all,
apologies if this is obvious - but I can't see it and would appreciate some
quick help!
the matrix mhouse is 26x3 and I'm computing odds ratios. The simple code
below "should" compute the odds vector for every pair (325) i.e. 26C2 in
cols 1 and 2.
On the first i=1 outer loop the inner j loop runs from 2 to 26 ok and then
I get the error (Error: subscript out of bounds)
Why isn't my loop incrementing i - the outer loop to 2 and then resetting
j=3?
Am I missing something obvious?
thanks
Gerard
> mhouse
[,1] [,2] [,3]
[1,] 275 81949
[2,] 593 1323 192
[3,] 813 1181 292
[4,] 2177 5189 1320
[5,] 1651 2243 270
[6,] 1061 5629 11035
[7,] 1690 2302 589
[8,] 1130 1203 345
[9,] 565 1898 655
[10,] 580 730 234
[11,] 343 176173
[12,] 372 53667
[13,] 666 1713 397
[14,] 382 918 279
[15,] 486 921 247
[16,] 1141 988 313
[17,] 626 1135 666
[18,] 438 436 168
[19,] 425 691 101
[20,] 609 71699
[21,] 467 661 141
[22,] 879 137379
[23,] 444 1101 130
[24,] 459 898 351
[25,] 995 1801 398
[26,] 396 1107 201
> # set up the odds vector by declaring it to be null
> odds=NULL
>
> # compute the odds ratios for Individual House vs Scheme House
> for (i in 1:25)
+ {
+ for (j in i+1:26)
+ {
+ todds = (mhouse[i,1]*mhouse[j,2])/(mhouse[j,1]*mhouse[i,2])
+ # compute the todds for row i with row j:
j>i
+ odds = c(odds,todds)
+ # append todds to the odds vector
+ }
+ }
Error: subscript out of bounds
> odds
[1] 0.7491244 0.4877622 0.8003391 0.4561745 1.7814132 0.4573697
0.3574670 1.1279674 0.4226138 1.7239078 0.4838053
[12] 0.8636384 0.8069156 0.6363150 0.2907502 0.6087939 0.3342421
0.5459312 0.3947703 0.4752623 0.5244818 0.8326321
[23] 0.6569199 0.6077702 0.9386447
>
**
The information transmitted is intended only for the person or entity to which
it is addressed and may contain confidential and/or privileged material. Any
review, retransmission, dissemination or other use of, or taking of any action
in reliance upon, this information by persons or entities other than the
intended recipient is prohibited. If you received this in error, please contact
the sender and delete the material from any computer. It is the policy of the
Department of Justice, Equality and Law Reform and the Agencies and Offices
using its IT services to disallow the sending of offensive material.
Should you consider that the material contained in this message is offensive
you should contact the sender immediately and also mailminder[at]justice.ie.
Is le haghaidh an duine nó an eintitis ar a bhfuil sí dírithe, agus le haghaidh
an duine nó an eintitis sin amháin, a bheartaítear an fhaisnéis a tarchuireadh
agus féadfaidh sé go bhfuil ábhar faoi rún agus/nó faoi phribhléid inti.
Toirmisctear aon athbhreithniú, atarchur nó leathadh a dhéanamh ar an
bhfaisnéis seo, aon úsáid eile a bhaint aisti nó aon ghníomh a dhéanamh ar a
hiontaoibh, ag daoine nó ag eintitis seachas an faighteoir beartaithe. Má fuair
tú é seo trí dhearmad, téigh i dteagmháil leis an seoltóir, le do thoil, agus
scrios an t-ábhar as aon ríomhaire. Is é beartas na Roinne Dlí agus Cirt,
Comhionannais agus Athchóirithe Dlí, agus na nOifígí agus na nGníomhaireachtaí
a úsáideann seirbhísí TF na Roinne, seoladh ábhair cholúil a dhícheadú.
Más rud é go measann tú gur ábhar colúil atá san ábhar atá sa teachtaireacht seo is ceart duit dul i dteagmháil leis an seoltóir láithreach agus le mailminder[ag]justice.ie chomh maith.
***
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Ca
Re: [R] Lexical Permutation Algorithm in R
Rory
there are several packages that perform this.
I would use permn() of the combinat library, then, if
lexicographical order is important, sort it explicitly.
HTH
rksh
[EMAIL PROTECTED] wrote:
Hi all
Here is a rather naive implementation of the SEPA algorithm for generating
lexical permutations:
lexperm3 <- function(x, n=length(x)) {
perms <- list()
k <- 1
perms[[k]] <- x
k <- k + 1
for (y in 1:(factorial(n)-1)) {
i <- n-1
while (x[i] > x[i+1] && i > 0) {
i <- i - 1
}
# i is largest index st x[i] > x[i+1]
j <- n
# find min{ x[j], st n>=j>=i+1 and x[j]>x[i] }
while (x[j] <= x[i] && j > i) {
j <- j - 1
}
# swap x[i] and x[j]
tmp <- x[i];x[i] <- x[j]; x[j] <- tmp
# now sort everything from x[i+1]..x[n]
# by reversing the elements within
p <- i + 1
q <- n
while (p < q) {
tmp <- x[p]; x[p] <- x[q]; x[q] <- tmp
p <- p + 1
q <- q - 1
}
perms[[k]] <- x
k <- k + 1
}
perms
}
This, as you can imagine, is severely slow. I would like to speed up this
function if possible, which I guess would involve vectorizing the inner
loop..does anyone have any ideas about how to improve this code's runtime?
One small potential optimization I tried was to shorten the "sort by reverse
ordering" near the end of the inner loop : I tried replacing it with rev() and
sort(decreasing=TRUE) over a partial subset of the x vector: however rev() was much
slower, and I dont think sort() supports lexicographic ordering, so that didnt work.
Thanks
rory
Rory Winston
RBS Global Banking & Markets
280 Bishopsgate, London, EC2M 4RB
Office: +44 20 7085 4476
***
The Royal Bank of Scotland plc. Registered in Scotland No 90312. Registered Office: 36 St Andrew Square, Edinburgh EH2 2YB.
Authorised and regulated by the Financial Services Authority
This e-mail message is confidential and for use by the=2...{{dropped:25}}
__
R-help@r-project.org 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.
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org 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] small numbers
Hi use the logarithmic representation for your problem. The Brobdingnag package uses such a form in a (more-or-less) user-transparent manner. HTH rksh Marc Jekel wrote: Dear R Fans, I have a simple probem but cannot find any reference to the soultion. I want to do calculations with small numbers (for max likelihood estimations). The lowest value R is storing by default is 1*10^-323, a smaller numer like 1*10^-324 is stored as a 0. How can I circumvent this problem? Is there a way to define how small a number can be in R. Thanks for a reply in advance, Marc -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] hypergeometric
The hypergeo package should be able to deal with this, although the function you specify below looks like a degenerate case (if I understand it correctly) so the convergence rate is likely to be slow. Let me know how you get on best wishes Robin (author of hypergeo) Jarle Brinchmann wrote: The dhyper etc deal with the hypergeometric _distribution_ while what you appear to want have is the hypergeometric special function (the connection is that the regular hypergeometric function is the generating function for the hypergeometric distribution if I recall correctly). Anyway, what you need, I believe, is the hypergeo library from CRAN. Cheers, Jarle. On Wed, Dec 3, 2008 at 9:17 AM, Zakaria, Roslinazairimah - zakry001 <[EMAIL PROTECTED]> wrote: Hi, I hope somebody can help me on how to use the hypergeometric function. I did read through the R documentation on hypergeometric but not really sure what it means. I would like to evaluate the hypergeometric function as follows: F((2*alpha+1)/2, (2*alpha+2)/2 , alpha+1/2, betasq/etasq). I'm not sure which function should be used- either phyper or qhyper or dhyper Where alpha <- .75; beta1 <- 7 ; beta2 <- 5.5; etasq <- ((beta1+beta2)/(2*beta1*beta2*(1-rho))) ^2 betasq <- ((beta1-beta2)^2+4*beta1*beta2*rho)/(4*beta1^2*beta2^2*(1-rho)^2) Thank you so much for your help. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] Multidimensional array with R
Hello a good place to start is R-and-octave.txt, in the contributed docs section of CRAN. This translates between common matlab and R commands HTH rksh Michael Zak wrote: Hi there I know, I'm sure you discussed this stuff 100 times, but I really have a basic understanding problem, if and how do I create a multidimensional array in R. I'm coming from MATLAB and there it's as easy as you ever could imagine. Ok, so, I want to have an array, where I can fill in data from a Excel spreadsheet. The array should be addressed like this: data["Cacao"]["Open"] or data["Cotton"]["Close"], you get me? Ok, that's fine. Otherwise feel free to ask questions. Because I tried to build such an array but without success (even with googling)... Thank you, Michael __ R-help@r-project.org 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. -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ R-help@r-project.org 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] Generating unique permutations of a vector
Annette
I understand your problem.
I think you may find 'blockparts(rep(5,5),5)' helpful.
I'm working on permutations of multisets right now
and expect to have functionality in partitions package
as soon as I finish the Other Ten Thousand Things
On My Things To Do List.arts
Perhaps we could talk offline.
best wishes
Robin
Annette Heisswolf wrote:
Hi all,
I try to generate sets of strategies that contain probability
distributions for a defined number of elements, e.g. imagine an
animal that can produce 5 different types of offspring and I want to
figure out which percentage of each type it should produce in order to
maximize its fitness. In order to do so, I need to calculate the fitness
for all potential strategies. As an example, if I assume the probability
levels to be 0, 0.2, 0.4, 0.6, 0.8, 1, I can generate all possible
probability distributions for the 5 types of offspring by using the
following R-code:
library(partitions)
library(crank)
n=5
P=parts(n)
L=length(P[1,])
for (i in 1:L){
ma=unique(permute(P[,i]))
if (i==1) m=ma else m=rbind(m,ma)
}
m=m/n
m
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0 0.0 0.0 0.0 0.0
[2,] 0.0 1.0 0.0 0.0 0.0
[3,] 0.0 0.0 1.0 0.0 0.0
[4,] 0.0 0.0 0.0 1.0 0.0
[5,] 0.0 0.0 0.0 0.0 1.0
[6,] 0.8 0.2 0.0 0.0 0.0
[7,] 0.8 0.0 0.2 0.0 0.0
[8,] 0.8 0.0 0.0 0.2 0.0
[9,] 0.8 0.0 0.0 0.0 0.2
[10,] 0.2 0.8 0.0 0.0 0.0
[11,] 0.2 0.0 0.8 0.0 0.0
[12,] 0.2 0.0 0.0 0.8 0.0
[13,] 0.2 0.0 0.0 0.0 0.8
[14,] 0.0 0.8 0.2 0.0 0.0
[15,] 0.0 0.8 0.0 0.2 0.0
[16,] 0.0 0.8 0.0 0.0 0.2
[17,] 0.0 0.2 0.8 0.0 0.0
[18,] 0.0 0.2 0.0 0.8 0.0
[19,] 0.0 0.2 0.0 0.0 0.8
[20,] 0.0 0.0 0.8 0.2 0.0
[21,] 0.0 0.0 0.8 0.0 0.2
[22,] 0.0 0.0 0.2 0.8 0.0
[23,] 0.0 0.0 0.2 0.0 0.8
[24,] 0.0 0.0 0.0 0.8 0.2
[25,] 0.0 0.0 0.0 0.2 0.8
[26,] 0.6 0.4 0.0 0.0 0.0
[27,] 0.6 0.0 0.4 0.0 0.0
.
.
.
[106,] 0.4 0.2 0.2 0.2 0.0
[107,] 0.4 0.2 0.2 0.0 0.2
[108,] 0.4 0.2 0.0 0.2 0.2
[109,] 0.4 0.0 0.2 0.2 0.2
[110,] 0.2 0.4 0.2 0.2 0.0
[111,] 0.2 0.4 0.2 0.0 0.2
[112,] 0.2 0.4 0.0 0.2 0.2
[113,] 0.2 0.2 0.4 0.2 0.0
[114,] 0.2 0.2 0.4 0.0 0.2
[115,] 0.2 0.2 0.2 0.4 0.0
[116,] 0.2 0.2 0.2 0.0 0.4
[117,] 0.2 0.2 0.0 0.4 0.2
[118,] 0.2 0.2 0.0 0.2 0.4
[119,] 0.2 0.0 0.4 0.2 0.2
[120,] 0.2 0.0 0.2 0.4 0.2
[121,] 0.2 0.0 0.2 0.2 0.4
[122,] 0.0 0.4 0.2 0.2 0.2
[123,] 0.0 0.2 0.4 0.2 0.2
[124,] 0.0 0.2 0.2 0.4 0.2
[125,] 0.0 0.2 0.2 0.2 0.4
[126,] 0.2 0.2 0.2 0.2 0.2
Using the "unique()" function is necessary, as the "permute()" function
doesn't check whether some elements within the vector to be permuted are
identical, e.g. permute(c(1,0,0)) gives:
[,1] [,2] [,3]
[1,]100
[2,]100
[3,]010
[4,]001
[5,]010
[6,]001
In my above example, this means that permuting the first column of P,
P[,1]
[1] 5 0 0 0 0
gives 5! = 120 permutations, although there are in fact only 5
unique ones, namely:
5 0 0 0 0
0 5 0 0 0
0 0 5 0 0
0 0 0 5 0
0 0 0 0 5
On my computer (Windows XP SP3, Pentium(R) 4 CPU 2.40GHz, 504 MB RAM -
but I've also tried it on another computer with 1 GB RAM) this leads to
the problem that as soon as the vector to be permuted has more than 9
elements, the resulting matrix get's too big and R won't execute the
permutation.
Thus, even when I aim to reduce the total number of unique permutations,
e.g. by using less probability levels than there are elements, the above
code doesn't work. For instance, I've used 10 elements but only the
probability levels 0, 0.2, 0.4, 0.6, 0.8, 1 as above, thus, P looks like
this:
P
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]5433221
[2,]0121211
[3,]0001111
[4,]0000011
[5,]0000001
[6,]0000000
[7,]0000000
[8,]0000000
[9,]0000000
[10,]0000000
In this case, permute(P[,1]) would give 10! = 3 628 800 permutations,
although there are only 10 unique ones. There are more then 10 unique
permutations for the other columns, but still their numbers should be
small enough.
Thus: Has anyone been working on a similar problem and found a better
way to generate such probability distributions? I'd appreciate any kind
of hint how to solve my problem. Thanks!
Annette
--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/p
Re: [R] Computational problems in R
Hello. The Brobdingnag package uses that identity for a logarithmic representation and also has a hack for negative numbers. HTH rksh A.Noufaily wrote: Many thanks for your suggestions... I am still checking which one is the most useful for my simulations. Concerning using logs, this might be very helpful, but I am not sure if I can use the following: A+B = A*(1 + B/A) = exp(log(A) + log(1 + B/A)) because unfortunately B can be negative. However, I might still use logs in case (1 + B/A)>0. Regards, Amy -Original Message- From: Duncan Murdoch [mailto:[EMAIL PROTECTED] Sent: Saturday, October 25, 2008 11:36 AM To: Steven McKinney Cc: A.Noufaily; r-help@r-project.org Subject: Re: [R] Computational problems in R On 24/10/2008 9:50 PM, Steven McKinney wrote: I suspect there's a deeper issue here. sum(exp(yi)) when large yi occur is problematic. exp(yi) for yi>710 is just a huge number, and summing additional values only makes the overall sum larger as all components of the summation are positive. There's no way around that. Sure there is, and you quoted it below. Work on a log scale. The log of exp(yi) is yi, and it sounds as though the yi values are manageable. You might end up knowing that the log of the final answer is 2 and not be able to evaluate exp(2) in R, but you still know that the answer is exp(2). Duncan Murdoch You could try this with Robin Hankins' package "brobdingnag" which can handle bunches of bizarrely large numbers. What kind of process are you studying? What kind of process generates values such as exp(8/0.01) when other values are much smaller? Are these outliers in an otherwise well-behaved data set? Perhaps then they need to be set aside and investigated separately, etc. Steven McKinney Statistician Molecular Oncology and Breast Cancer Program British Columbia Cancer Research Centre email: smckinney +at+ bccrc +dot+ ca tel: 604-675-8000 x7561 BCCRC Molecular Oncology 675 West 10th Ave, Floor 4 Vancouver B.C. V5Z 1L3 Canada -Original Message- From: [EMAIL PROTECTED] on behalf of Duncan Murdoch Sent: Fri 10/24/2008 4:04 PM To: A.Noufaily Cc: r-help@r-project.org Subject: Re: [R] Computational problems in R On 24/10/2008 12:42 PM, A.Noufaily wrote: Dear all, I would be grateful if anyone can help me with the following: My aim is to compute explicitely the sum S=A+B where A=sum(exp(c_i/d)), i=1,...,n; B, c_i, and d are real numbers with -Inf0. The problem is that when c_i/d >710 (for some i) R is setting exp(c_i/d) to be equal to +Inf and hence the whole summation S. So in simple cases where for example c_i=8 (for some i), and d=0.01 the whole summation is turning out to be infinite. Is there a way to get round that in R? Can anyone suggest any computational trick to calculate S when c_i/d>710 (for some i)? Work on a log scale. Use the identity that A+B = A*(1 + B/A) = exp(log(A) + log(1 + B/A)) (where you chose A to be the biggest term in the sum). Duncan Murdoch Any suggestions would be much appreciated. Regards, Amy - The Open University is incorporated by Royal Charter (RC 000391), an exempt charity in England & Wales and a charity registered in Scotland (SC 038302). [[alternative HTML version deleted]] __ R-help@r-project.org 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@r-project.org 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@r-project.org 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. - The Open University is incorporated by Royal Charter (RC 000391), an exempt charity in England & Wales and a charity registered in Scotland (SC 038302). __ R-help@r-project.org 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@r-project.org 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] weird behavior with the 3rd root....
This comes up from time to time. The problem is that one needs complex numbers to address taking the third root: there are three cube roots for any nonzero number (real or complex). To wit: > (-0.084121928394+0i)^(1/3) [1] 0.2190818+0.3794609i > (-0.084121928394-0i)^(1/3) [1] 0.2190818+0.3794609i > (-0.084121928394+1e-100i)^(1/3) [1] 0.2190818+0.3794609i > (-0.084121928394-1e-100i)^(1/3) [1] 0.2190818-0.3794609i > Note the first two are identical but the second two differ. Anyone care to start discussing signed zero again? [you probably want the *real* cube root, in which case it is best to take minus the unique real cube root of the absolute value: > -(0.084121928394)^(1/3) [1] -0.4381637 > (which is what you did, of course!)] HTH rksh Juan Manuel Barreneche wrote: Well, this is what i got... -0.084121928394^(1/3) [1] -0.438163696867656 (-0.084121928394)^(1/3) [1] NaN and i don't have a clue of why this happens or how to avoid it, any suggestions? thank you, Juan __ R-help@r-project.org 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. -- Robin K. S. Hankin Senior Research Associate Cambridge Centre for Climate Change Mitigation Research (4CMR) Department of Land Economy University of Cambridge [EMAIL PROTECTED] 01223-764877 __ R-help@r-project.org 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] Generating a valid covariance matrix
Megh corr.matrix() in the 'emulator' package can calculate P-D variance matrices using any of a very broad class of methods. HTH rksh Megh Dal wrote: I want to generate a valid variance-covariance matrix. One way could be to generate some random sample from multivariate normal distribution and then calculate cov. matrix. Another way could be to sample from wishart distribution itself. However both cases need a valid i.e. PD covariance matrix. As I need to generate that covariance matrix only, I am not interested those two methods. Can anyone suggest me some other way out? Regards, __ R-help@r-project.org 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. -- Robin K. S. Hankin Senior Research Associate Cambridge Centre for Climate Change Mitigation Research (4CMR) Department of Land Economy University of Cambridge [EMAIL PROTECTED] 01223-764877 __ R-help@r-project.org 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] Generalising to n-dimensions
Laura Bonnett wrote:
Sorry to hassle you, but I really need to get my code up and running.
Please can you therefore explain what a and v are?
Hi Laura. I've been away (in Norwich). Sorry not to give an example.
Variable 'a' is an array and variable 'd' is the same as in your
original email.
> a <- array(1:4,c(3,2,2,4))
> d <- c(1,2)
> f(a,d)
[,1] [,2]
[1,]14
[2,]21
[3,]32
>
Thus in this case f(a,d) gives a[,,1,2]. Function f() is necessary because
you specified that the length of 'd' is not known in advance.
You can get 'd' from a single row of expand.grid() [but you will have to
coerce it to a matrix]
HTH
rksh
Thank you,
Laura
On Wed, Sep 24, 2008 at 8:27 PM, Laura Bonnett
<[EMAIL PROTECTED] <mailto:[EMAIL PROTECTED]>> wrote:
Can I ask what a and v are?
Thanks,
Laura
On Sat, Aug 23, 2008 at 11:41 AM, Robin Hankin <[EMAIL PROTECTED]
<mailto:[EMAIL PROTECTED]>> wrote:
Laura Bonnett wrote:
crosstable[,,expand[d,1],expand[d,2],expand[d,3],...expand[d,n]]
crosstable is just a crosstabulation of an
n+2-dimensional dataset and I am trying to pick out those
that are in combination 'd' of expand. So for example, for
5-dimensional data using your example:
Var1 Var2 Var3
1 111
2 211
3 311
4 121
5 221
6 321
7 112
8 212
9 312
10122
11222
12322
d refers to the row of the matrix above - d=2 is 2,1,1 so
crosstable[,,2,1,1] would retrieve all the data where Var1
=2, Var2=1, Var3=1 and the two remaining variables are
given in the crosstabulations for all values.
Is that any better?
OK I think I understand. The magic package uses this type of
construction extensively, but not this particular one.
It's trickier than I'd have expected.
Try this:
f <- function(a,v){
jj <-
sapply(dim(a)[seq_len(length(dim(a))-length(v))],seq_len,simplify=FALSE)
jj <- c(jj , as.list(v))
do.call("[" , c(list(a) , jj, drop=TRUE))
}
[you will have to coerce the output from expand.grid() to a
matrix in order to extract a row from it]
HTH
rksh
--
Robin K. S. Hankin
Senior Research Associate
Cambridge Centre for Climate Change Mitigation Research (4CMR)
Department of Land Economy
University of Cambridge
[EMAIL PROTECTED] <mailto:[EMAIL PROTECTED]>
01223-764877
--
Robin K. S. Hankin
Senior Research Associate
Cambridge Centre for Climate Change Mitigation Research (4CMR)
Department of Land Economy
University of Cambridge
[EMAIL PROTECTED]
01223-764877
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to draw the graph of f(x,y) = x * y ?
Paul you might find the view() function in the 'elliptic' package useful. This function implements various methods to visualize functions over the complex plane. HTH rksh Paul Smith wrote: Dear All, The function curve() draws the graph of functions from R to R. Is there some homologous function to curve() to draw functions from R^2 to R? Thanks in advance, Paul __ R-help@r-project.org 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. -- Robin K. S. Hankin Senior Research Associate Cambridge Centre for Climate Change Mitigation Research (4CMR) Department of Land Economy University of Cambridge [EMAIL PROTECTED] 01223-764877 __ R-help@r-project.org 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] Generalising to n-dimensions
Laura Bonnett wrote:
crosstable[,,expand[d,1],expand[d,2],expand[d,3],...expand[d,n]]
crosstable is just a crosstabulation of an n+2-dimensional dataset and
I am trying to pick out those that are in combination 'd' of expand.
So for example, for 5-dimensional data using your example:
Var1 Var2 Var3
1 111
2 211
3 311
4 121
5 221
6 321
7 112
8 212
9 312
10122
11222
12322
d refers to the row of the matrix above - d=2 is 2,1,1 so
crosstable[,,2,1,1] would retrieve all the data where Var1 =2, Var2=1,
Var3=1 and the two remaining variables are given in the
crosstabulations for all values.
Is that any better?
OK I think I understand. The magic package uses this type of
construction extensively, but not this particular one.
It's trickier than I'd have expected.
Try this:
f <- function(a,v){
jj <-
sapply(dim(a)[seq_len(length(dim(a))-length(v))],seq_len,simplify=FALSE)
jj <- c(jj , as.list(v))
do.call("[" , c(list(a) , jj, drop=TRUE))
}
[you will have to coerce the output from expand.grid() to a matrix in
order to extract a row from it]
HTH
rksh
--
Robin K. S. Hankin
Senior Research Associate
Cambridge Centre for Climate Change Mitigation Research (4CMR)
Department of Land Economy
University of Cambridge
[EMAIL PROTECTED]
01223-764877
__
R-help@r-project.org 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] Generalising to n-dimensions
First bit: > x <- c(3,2,2) > expand.grid(sapply(x,seq_len)) Var1 Var2 Var3 1 111 2 211 3 311 4 121 5 221 6 321 7 112 8 212 9 312 10122 11222 12322 > second bit I'm not sure about. I didn't quite get why d=2 implied the order is 2,1. Could you post a small self-contained example? HTH rksh Laura Bonnett wrote: Hi R-helpers, I have two queries relating to generalising to n dimensions: What I want to do in the first one is generalise the following statement: expand<-expand.grid(1:x[1],1:x[2],...1:x[n]) where x is a vector of integers and expand.grid gives every combination of the set of numbers, so for example, expand.grid(1:2, 1:3) takes 1,2 and 1,2,3 and gives 1,1 2,1 1,2 2,2 1,3 2,3 My x vector has varying lengths and I can't find a way of giving it every set without stating each set individually. Secondly and similarly, I want to get the table within crosstable that has the elements defined by the combinations given in expand above crosstable[,,expand[d,1],expand[d,2],expand[d,3],...expand[d,n]] where crosstable is just a crosstabulation of an n+2-dimensional dataset and I am trying to pick out those that are in combination 'd' of expand. So for example, using x[1]=2 and x[2]=3 as above example, if d =2 then the order is 2,1 so I take crosstable[,,2,1]. Can anyone suggest a way to give the code every set without stating each set individually? Thank you [[alternative HTML version deleted]] __ R-help@r-project.org 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. -- Robin K. S. Hankin Senior Research Associate Cambridge Centre for Climate Change Mitigation Research (4CMR) Department of Land Economy University of Cambridge [EMAIL PROTECTED] 01223-764877 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to keep up with R?
Adaikalavan Ramasamy wrote: I agree! The best way to learn (and remember for longer) is to teach someone else about it. And there is not reason not to repeat some of the anlysis done on SAS with R. That way you can verify your outputs or compare the presentations. If you consistently find differences in the outputs, then trying to figure out the reason may lead you to better understand the methods (e.g. different optimization or estimation procedures). My take on this: I have repeatedly found that it is surprisingly easy to improve on existing (non-R) implementations of statistical and non-statistical computation, when working in R. Something about the structure of the language, something about the package mechanism, something about R-help, something about R-core, something about open-source, something about JSS or R-news, whatever it is, there is SOMETHING ABOUT R which lends itself to straightforward production of quality software. And that something is missing from other programming languages, IMO. rksh Regards, Adai Barry Rowlingson wrote: 2008/9/19 Wensui Liu <[EMAIL PROTECTED]>: Dear Listers, I've been a big fan of R since graduate school. After working in the industry for years, I haven't had many opportunities to use R and am mainly using SAS. However, I am still forcing myself really hard to stay close to R by reading R-help and books and writing R code by myself for fun. But by and by, I start realizing I have hard time to keep up with R and am afraid that I would totally forget how to program in R. I really like it and am very unwilling to give it up. Is there any idea how I might keep touch with R without using it in work on daily basis? I really appreciate it. -- Robin K. S. Hankin Senior Research Associate Cambridge Centre for Climate Change Mitigation Research (4CMR) Faculty of Economics The University of Cambridge [EMAIL PROTECTED] 01223-764877 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to keep up with R?
Hi Wensei. Why not do as I do? Find an interesting area of numerical computation (perhaps not statistical) that has not been implemented in open-source. Then write an R package for it, under GPL-2, then write an article about the new package in R-news or JSS. works for me. Best wishes Robin Wensui Liu wrote: Dear Listers, I've been a big fan of R since graduate school. After working in the industry for years, I haven't had many opportunities to use R and am mainly using SAS. However, I am still forcing myself really hard to stay close to R by reading R-help and books and writing R code by myself for fun. But by and by, I start realizing I have hard time to keep up with R and am afraid that I would totally forget how to program in R. I really like it and am very unwilling to give it up. Is there any idea how I might keep touch with R without using it in work on daily basis? I really appreciate it. [[alternative HTML version deleted]] __ R-help@r-project.org 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. -- Robin K. S. Hankin Senior Research Associate Cambridge Centre for Climate Change Mitigation Research (4CMR) Faculty of Economics The University of Cambridge [EMAIL PROTECTED] 01223-764877 __ R-help@r-project.org 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] give all combinations
Hi Yuan, Lucien, List.
try this:
f <-
function (...)
{
args <- list(...)
if(length(args)==0){
return(NULL)
}
if (length(args) == 1) {
return(args[[1]])
}
if (length(args) > 2) {
jj <- do.call("Recall", c(args[-1]))
return(do.call("Recall", c(list(args[[1]]), list(jj) )))
}
a <- args[[1]]
b <- args[[2]]
if (is.null(b)) {
return(a)
}
jj <- outer(a,b,paste)
return(jj[!lower.tri(jj)])
}
[the difficult bit (IMO) is to make f() work with any number of
arguments; thus
f(a,a) and f(a,a,a,a,a,b,b,a,b) and whatever should also work and this
is why
the Recall bit is needed].
Comments anyone?
HTH
rksh
Lucien Lemmens wrote:
Another solution requiring also a bit of programming is:
l<-letters[1:3]
c2<-c()
for(i in 1:3){c2<-c(c2,paste(letters[i],letters[i:3],sep=""))}
c2
[1] "aa" "ab" "ac" "bb" "bc" "cc"
n<-length(c2)
c3<-c();for(i in 1:n){c3<-c(c3,paste(c2[i],letters[ceiling(i/2):3],sep=""))}
c3
[1] "aaa" "aab" "aac" "aba" "abb" "abc" "acb" "acc" "bbb" "bbc" "bcc" "ccc"
__
R-help@r-project.org 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@r-project.org 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] [R-pkgs] new package multipol
Hello List please find a new package, multipol, recently uploaded to CRAN. This package generalizes the polynom package (which handles univariate polynomials) to the multivariate case. A short article discussing the package will appear in the next issue of Rnews, Insha'Allah enjoy -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org 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] S4 : package creation
One other reference is the vignette in the Brobgingnag
package; this is a cookbook approach that gives step-by-step
instructions on how to build a simple package using S4 methods.
rksh
On 21 Mar 2008, at 17:29, Martin Morgan wrote:
> Hi Christophe --
>
> In terms of documentation, see ?promptClass, ?promptMethods.
>
> I don't think the description of package creation in 'S4 Classes in 15
> pages, more or less' is the way things are generally done these
> days. A package might normally look like
>
> DESCRIPTION
> NAMESPACE
> R/AllClasses.R
> R/methods-SomeClass.R
>
> The content and organization of the files in the R directory are up to
> the package author, but consist of the usual
>
> setClass("A", ...
> setMethod("foo", ...
>
> as well as other code. The content of DESCRIPTION and NAMESPACE are
> described in the 'Writing R Extensions' manual; important points are
> to Depends: methods and LazyLoad: yes in the DESCRIPTION file. It is
> also good to arrange that the files in the R directory are sourced in
> such a way that classes are defined before the methods that use them
> (using Collate: in the DESCRIPTION).
>
> There are a number of example packages available, including those
> written specifically for illustrating these issues. Several of these
> were referenced in recent threads on this news group, or perhaps it
> was the R-devel news group.
>
> Martin
>
> Christophe Genolini <[EMAIL PROTECTED]> writes:
>
>> Hi the list,
>>
>> Using S4, how can we create a package? In "S4 Classes in 15 pages,
>> more
>> or less", they put all the classes definition in a function that
>> will be
>> called at the opening of the library and they add "by hand" a Rd
>> file.
>> Is it the only way ? Is there something like "S4.package.skeleton"?
>>
>> Thanks
>>
>> Christophe
>>
>> __
>> R-help@r-project.org 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.
>
> --
> Martin Morgan
> Computational Biology / Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N.
> PO Box 19024 Seattle, WA 98109
>
> Location: Arnold Building M2 B169
> Phone: (206) 667-2793
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst and Neutral Theorist,
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] empty array
Hello everyone
I know other, more knowledgeable, people
have replied to Christophe's question, but
perhaps the List would be interested to know
that zero-extent arrays are useful (to me at least)
because although such an array has no content, the
dimname are nevertheless retained:
> a <- array(0,dim=c(0,3))
> dimnames(a) <- list(size=c() , fish=c("cod","skate","crab"))
> b <- array(0,dim=c(2,0))
> dimnames(b) <- list(size=c("huge","small"),depth=c())
>
We can attach these arrays ---both of which are of length 0---using
adiag():
> library(magic)
> adiag(a,b)
fish
sizecod skate crab
huge0 00
small 0 00
>
Note how the dimnames of "a" and "b" are retained in the output.
The contents of this array are the default "pad" value of adiag().
This is terribly useful in the humble workaday world
of high-dimensional magic hypercubes.
rksh
On 15 Mar 2008, at 15:33, Christophe Genolini wrote:
> Hi the list
>
> Is it possible to create an empty matrix ? I do not mean an matrix
> with
> a single value that is NA (which is not empty) but a real empty one,
> with length=0.
>
> I do not understand why we have length(numeric()), length(factor())
> and
> length(character()) to zero, and length(array()) to one... Any rason
> for
> that ?
>
> Thanks
>
> Christophe
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst and Neutral Theorist,
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] Bessel functions of complex argument
Hello Baptiste Bessel functions with complex arguments are not supported in R. Neither matlab nor the Gnu Scientific Library support them either. . . . but . . . the pari/gp system (released on the GPL) does: ? besselj(1+I,3) %3 = 0.6919067491368555819808728680 + 0.4484268613977010268818252591*I ? You can access some pari/gp functionality from within R by using the elliptic package, although unfortunately its wrapper function, P.pari(), is not quite flexible enough to deal with besselj(). I'd be happy to discuss this offline; P.pari() will need only minor changes to accommodate besselj(). HTH Robin On 9 Mar 2008, at 13:44, baptiste Auguié wrote: > Dear R users, > > > I'm porting a piece of Matlab code to R, but I'm now stuck with the > following: I need an equivalent of besselJ(x, nu) that can handle a > complex argument x. I couldn't find any R implementation. I did find > a possible fortran solution in SLATEC (< http://www.netlib.org/slatec/ >> , CBESJ-C), however I've never tried to use external C or Fortran > code together with my R code, so I'm not sure where to go for a > simple solution. > > Any advice welcome, > > Best regards > > baptiste > > _ > > Baptiste Auguié > > Physics Department > University of Exeter > Stocker Road, > Exeter, Devon, > EX4 4QL, UK > > Phone: +44 1392 264187 > > http://newton.ex.ac.uk/research/emag > http://projects.ex.ac.uk/atto > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] [OT] "normal" (as in "Guassian")
On 4 Mar 2008, at 08:20, Ingmar Visser wrote: > Wikipedia says: > > Stigler attributes the discovery of Stigler's Law to Robert K. Merton > (which makes the law self-referencing). > Stigler's law certainly applies in mathematics, where standard procedure is to name a concept in honour of the first person after Euler to have (re)discovered it. rksh > (Working as a historian of science he should have proceeded to name > the law Merton's law only to find out > later that actually someone had discovered it even earlier.) > > Ingmar > > [edit] > On 3 Mar 2008, at 19:17, Douglas Bates wrote: > >>>> >>>> -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] PDF with computationally expensive normalizing constant
Hi I am writing some functionality for a multivariate PDF. One problem is that evaluating the normalizing constant (NC) is massively computationally intensive [one recent example took 4 hours and bigger examples would take much much longer] and it would be good allow for this in the design of the package somehow. For example, the likelihood function doesn't need the NC but (eg) the moment generating function does. So a user wanting a maximum-likelihood estimate shouldn't have to evaluate the NC but a user wanting a mean has to. Some simple forms of the PDF have an easily-evaluated analytical expression for the NC. And once the NC is evaluated, it would be good to store it somehow. I thought perhaps I could define an S4 class with a slot for the parameters and a slot for the NC; and if the NC is unknown this would have an "NA" entry. Then a user could execute something like a <- CalculateNormalizingConstant(a) and after this, object "a" would then have the numerically computed NC in place. Is this a Good Idea? Are there any PDFs implemented in R in which this is an issue? -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] learning S4
Christophe
you might find the Brobdingnag package on CRAN helpful here.
I wrote the package partly to teach myself S4; it includes a
vignette that builds the various S4 components from scratch,
in a step-by-step annotated "cookbook".
HTH
rksh
On 8 Feb 2008, at 15:30, [EMAIL PROTECTED] wrote:
> Hi the list.
>
> I try to learn the S4 programming. I find the wiki and several doc.
> But
> I still have few questions...
>
> 1. To define 'representation', we can use two syntax :
>- representation=list(temps = 'numeric',traj = 'matrix')
>- representation(temps = 'numeric',traj = 'matrix')
> Is there any difference ?
> 2. 'validityMethod' check the intialisation of a new object, but not
> the latter
> modifications. Is it possible to set up a validation that check
> every
> modifications ?
> 3. When we use setMethod('initialize',...) does the validityMethod
> become un-used ?
> 4. Is it possible to set up several initialization processes ? One
> that build an objet from a data.frame, one from a matrix...
>
> Thanks
>
> Christophe
>
>
> Ce message a ete envoye par IMP, grace a l'Universite Paris 10
> Nanterre
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst and Neutral Theorist,
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] Appell Hypergeometric function
hits=-2.6 tests=BAYES_00 X-USF-Spam-Flag: NO Hi the closest I can offer is genhypergeo() from the Davies package HTH rksh (actually, I have a few items on my Things To Do List for the Davies package, and adding Appell hypergeometric functions is one of them. But don't hold your breath ;-) On 7 Feb 2008, at 10:10, Giovanni Parrinello wrote: > Dear All, > I am looking for an implementation in R of the Appell Hypergeometric > function. > > Any suggestions will be more than appreciated! > > GP > > > -- > dr. Giovanni Parrinello > External Lecturer > Medical Statistics Unit > Department of Biomedical Sciences > Viale Europa, 11 - 25123 Brescia Italy > Tel: +390303717528 > Fax: +390303717488 > email: [EMAIL PROTECTED] > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] "=" in functions and matlab (was: Multiplying each row of a big matrix with a vector)
hits=-2.6 tests=BAYES_00
X-USF-Spam-Flag: NO
[this after a mistakenly private email to Megh]
As others have said, "%*%" is sufficiently vectorized to do
what you want.
Speaking as a reformed Matlab user,
it might be as well to add that the definition of function fu()
that you give looks like matlab code to me.
Check this out:
> f <- function(x){f = x^2}
> f(4)
>
> dput(f(4))
16
>
Thus fu() returns 16 but invisibly. This can be terribly
confusing to those from a matlab background.
In R, one can just write
> f<- function(x){x^2}
because functions return the last evaluated expression.
Compare matlab, where(IIRC) one has to state in the function
the name of the variable whose value will be returned.
HTH
rksh
On 30 Jan 2008, at 07:20, Megh Dal wrote:
> I have a big matrix 'ret'. I want to multiply each row of it with a
> 2nd vector 'pos', resulting result, I want to save in a vector named
> 'port'. I wrote following code:
>
>> pos
> [1] 2593419 2130220 6198197 1673888 198 1784732 2052120
> -7490228 -5275000
>
>
>> dim(ret)
> [1] 500 9
>
>> fu# user defined function
> function(x)
> {
>fu = x %*% t(pos)
> }
> port = apply(ret, 1, fu)
>
>> dim(port)
> [1] 81 500
>
> My desire is to get port as a vector with length 500. However I am
> not getting that?
>
> Can anyone tell me how to correct that?
>
> Regards,
>
>
>
> -
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst and Neutral Theorist,
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] An "R is slow"-article
Hello Gustaf, List. Thanks Gustaf for your post! well I am working pretty intensively with fisher.test() right now, as some of you will know. The comparison is not fair: R's fisher.test() does a whole bunch of error checking and testing for the size of the input matrix and assessing of other arguments, and puts together a nice little list of class "htest". The C routine does none of this. The clincher is that fisher.test() as called gives an estimate for the odds ratio using uniroot() to numerically solve an equation in terms of the hypergeometric probability distribution. This takes a lo time, but one doesn't notice it in a standard R session. Sorry, but the time comparison is simply not worth reporting. On 9 Jan 2008, at 15:25, Gustaf Rydevik wrote: > Hi all, > > Reading the wikipedia page on R, I stumbled across the following: > http://fluff.info/blog/arch/0172.htm > > It does seem interesting that the C execution is that much slower from > R than from a native C program. Could any of the more technically > knowledgeable people explain why this is so? > > The author also have some thought-provoking opinions on R being > no-good and that you should write everything in C instead (mainly > because R is slow and too good at graphics, encouraging data > snooping). See http://fluff.info/blog/arch/0041.htm > While I don't agree (granted, I can't really write C), it was > interesting to read something from a very different perspective than > I'm used to. > > Best regards, > > Gustaf > > _ > Department of Epidemiology, > Swedish Institute for Infectious Disease Control > work email: gustaf.rydevik at smi dot ki dot se > skype:gustaf_rydevik > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] array addition
[snip snip snip]
> suppose I have two arrays x1,x2 of dimensions a1,b1,c1 and
> a2,b2,c2 respectively.
>
> I want x = x1 "+" x2 with dimensions c(max(a1,a2), max(b1,b2),max
> (c1,c2))
[snip snip snip]
perhaps it wouldn't be too much to ask for you to
check the most recent version of the "magic" package?
[and we *really* don't want any whingeing about
magic_1.3-31 not being available. If I were you, I'd
email the package maintainer and tell him to release
updates in a more timely manner. . . ]
> library(magic)
> aplus
function (...)
{
args <- list(...)
if (length(args) == 1) {
return(args[[1]])
}
if (length(args) > 2) {
jj <- do.call("Recall", c(args[-1]))
return(do.call("Recall", c(list(args[[1]]), list(jj
}
a <- args[[1]]
b <- args[[2]]
dima <- dim(a)
dimb <- dim(b)
stopifnot(length(dima) == length(dimb))
out <- array(0, pmax(dima, dimb))
return(do.call("[<-", c(list(out), lapply(dima, seq_len),
list(a))) + do.call("[<-", c(list(out), lapply(dimb,
seq_len), list(b
}
>
--
Robin Hankin
Uncertainty Analyst and Neutral Theorist,
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] array addition
Hi suppose I have two arrays x1,x2 of dimensions a1,b1,c1 and a2,b2,c2 respectively. I want x = x1 "+" x2 with dimensions c(max(a1,a2), max(b1,b2),max (c1,c2)) with x[a,b,c] = x1[a1,b1,c1] + x2[a2,b2,c2] ifa <=min(a1,a2) , b<=min (b1,b2), c<=min(c1,c2) and the other bits either x1 or x2 or zero according to whether the coordinates are "in range" for x1 or x2 or neither. The answer has to work for arbitrary-dimensioned arrays. toy example follows (matrices): > x1 [,1] [,2] [,3] [,4] [,5] [1,]13579 [2,]2468 10 > x2 [,1] [,2] [,3] [1,]147 [2,]258 [3,]369 > x [,1] [,2] [,3] [,4] [,5] [1,]27 1279 [2,]49 148 10 [3,]36900 > Note the zeros at lower-right. Is there a ready-made solution to this? -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] can R solve these paired equations
Hello [An answer was posted just now using numerical ideas; here is an answer from a symbolic perspective] These equations involve x^y in more than one unknown, so inverse functions cannot be used. I do not think you will be able to characterize even the number of solutions, let alone their nature. To see the difficulty, look at the Lambert W function. My advice would be to simplify, simplify, simplify your problem as far as possible; remove terms successively until you are left with a trivial system, then add *one* term and see if this produces any insights. HTH rksh On 17 Dec 2007, at 05:43, Xin wrote: > Dear: > > > > I have a paired equation below. Can I solve (x,y) using R. > > > > Thanks! > > > > Xin > > > > A=327.727 > > B=9517.336 > > p=0.114^10 > > > > (1-p)*y*(1-x)/x/(1-x^y)=A > > A(1+(1-x)*(1+y)/x-A))=B > > > > > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst and Neutral Theorist, National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] Matrix Inversion
Hello Wang
matrix bb is symmetric positive semidefinite, so
algebraically the eigenvalues are nonnegative.
I would use
bb <- crossprod(b)
to calculate bb (faster and possibly more accurate)
Look at eigen(bb,TRUE,TRUE)$values
(see ?eigen for the meaning of the arguments) to see how
many very small eigenvalues you have. The number of zero
eigenvalues is equal to the number of linear relations
in the columns of b.
HTH
rksh
On 12 Dec 2007, at 10:59, Wang Chengbin wrote:
> I got the following error:
>
> a = read.csv("mat.csv")
> b = as.matrix(a)
> tb = t(b)
> bb = tb %*% b
> dim(bb)
> ibb = solve(bb)
> bb %*% ibb
>
>> ibb = solve(bb)
> Error in solve.default(bb) :
> system is computationally singular: reciprocal condition number =
> 1.77573e-19
>>
> Are there any ways to find more information about why it is singular?
>
> Thanks.
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst and Neutral Theorist,
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] rowSums() and is.integer()
On 21 Nov 2007, at 08:30, Prof Brian Ripley wrote: > On Tue, 20 Nov 2007, Tim Hesterberg wrote: > >> I wrote the original rowSums (in S-PLUS). >> There, rowSums() does not coerce integer to double. > > Actaully, neither does R. It computes a double answer but does no > coercion per se. > >> However, one advantage of coercion is to avoid integer overflow. > > Indeed, as I told Robin Hankin privately, that was the design reason. > Brian Ripley also reminded me that the sum() of integers is an integer, behaviour that I find desirable. The reason for my starting this thread is that sometimes I actually *want* sums of integers to overflow: my interest is in exact computations where I must be absolutely certain that there can be no rounding error. If the sum cannot be represented in integers, I want this fact to be flagged with extreme vigour as it signals what might be catastrophic loss of precision. At least, that's my current thinking. best wishes rksh >> >> Tim Hesterberg >> >>> ... So, why does rowSums() coerce to double (behaviour >>> that is undesirable for me)? >> >> __ >> R-help@r-project.org 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. >> > > -- > 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 -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] All nonnegative integer solution
Hello Amin The partitions library does this. If N=4 and k=3: > library(partitions) > blockparts(rep(4,3),4) [1,] 4 3 2 1 0 3 2 1 0 2 1 0 1 0 0 [2,] 0 1 2 3 4 0 1 2 3 0 1 2 0 1 0 [3,] 0 0 0 0 0 1 1 1 1 2 2 2 3 3 4 > The solutions are enumerated in the columns. HTH rksh On 19 Nov 2007, at 09:57, [EMAIL PROTECTED] wrote: > Dear all, > Is there any method in R to find all possible nonnegative integer > solutions to the linear equation with unit coefficients as follow: > X1+X2+...+Xk=N > Thank you, > Amin Zollanvari > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R det
Hello Kalin
det() does not take a complex matrix as an argument.
To get the determinant of a complex matrix, use eigen():
mydet <- function(a){prod(eigen(a,only.values=TRUE)$values)}
a <- matrix(1:9,3,3)
a[1,1] <- 1i
mydet(a)
[List: can we not add the above, or something like it,
to the definition of det() so that it can deal with complex matrices?]
HTH
Robin
On 16 Nov 2007, at 19:27, kalin lagno wrote:
> Hi,
> Which R function I should use to obtain determinant of a matrix
> with real(and complex) numbers?
>
> Kalin
>
> -
> Never miss a thing. Make Yahoo your homepage.
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] creating discretized data
Hi
The "trick" is to define a function f() that does
what you want elementwise, then use lapply():
> f <- function(i){c(rep(0,i-1),1)}
> x <- c(2,1,3,5)
> c(lapply(x,f),recursive=T)
[1] 0 1 1 0 0 1 0 0 0 0 1
>
HTH
rksh
> Hi, Ia m working in discretized data. Here my data:
>
> x <- c(2,1,3, 5), and I want to make (0,1) data based on the
> length of
> each component in x.
> So the new data should like: y = (0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1).
> I spent
> too much time with
> "seq", "rep". Still didn't get it. Any help? Thanks
>
> Ilham
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] rowSums() and is.integer()
On 10 Nov 2007, at 07:32, Prof Brian Ripley wrote: > On Fri, 9 Nov 2007, Robin Hankin wrote: > >> Hi >> >> [R-2.6.0, macOSX 10.4.10]. >> >> The helppage says that rowSums() and colSums() >> are equivalent to 'apply' with 'FUN = sum'. >> >> But I came across this: >> >> > a <- matrix(1:30,5,6) >> > is.integer(apply(a,1,sum)) >> [1] TRUE >> > is.integer(rowSums(a)) >> [1] FALSE >> > > > 'equivalent' does not mean 'identical': the wording was deliberate. > >> so rowSums() returns a float. > > And that is what the help page says it does (albeit more > accurately: there is no 'float' type, but there is numeric aka > double and the result could be complex). > >> Why is this? > > You seem to be asking why R works as documented! > Yes, that's exactly what I was asking [perhaps this should have been R-devel?]. What is the thinking behind converting to double? I expect that part of the answer is speed: # First define an integer matrix: a <- matrix(as.integer(rpois(1e6,3)),1000,1000) > system.time(rowSums(a)) user system elapsed 0.049 0.000 0.050 > system.time(rowSums(a)) user system elapsed 0.050 0.000 0.051 > system.time(rowSums(a)) user system elapsed 0.050 0.001 0.052 > system.time(colSums(a)) user system elapsed 0.043 0.001 0.046 > system.time(colSums(a)) user system elapsed 0.043 0.000 0.044 About the same speed. Now use apply() to see whether integer summation is faster than double summation for this kind of problem: > system.time(ignore <- apply(a,1,sum)) user system elapsed 0.085 0.009 0.094 > system.time(ignore <- apply(a,1,sum)) user system elapsed 0.086 0.010 0.095 > system.time(ignore <- apply(a,1,sum)) user system elapsed 0.089 0.010 0.104 > system.time(ignore <- apply(a,2,sum)) user system elapsed 0.071 0.008 0.078 > system.time(ignore <- apply(a,2,sum)) user system elapsed 0.069 0.007 0.076 > system.time(ignore <- apply(a,2,sum)) user system elapsed 0.070 0.008 0.081 # Now convert to double: > a <- a+0 > system.time(ignore <- apply(a,1,sum)) user system elapsed 0.127 0.019 0.151 > system.time(ignore <- apply(a,1,sum)) user system elapsed 0.121 0.017 0.139 > system.time(ignore <- apply(a,1,sum)) user system elapsed 0.130 0.022 0.175 > system.time(ignore <- apply(a,2,sum)) user system elapsed 0.084 0.015 0.098 > system.time(ignore <- apply(a,2,sum)) user system elapsed 0.085 0.015 0.105 > system.time(ignore <- apply(a,2,sum)) user system elapsed 0.087 0.016 0.107 [can anyone comment on the difference between the first three and the last three double precision summations?] perhaps a little bit faster for the integers, but there's not much in it. So, why does rowSums() coerce to double (behaviour that is undesirable for me)? -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] rowSums() and is.integer()
Hi [R-2.6.0, macOSX 10.4.10]. The helppage says that rowSums() and colSums() are equivalent to 'apply' with 'FUN = sum'. But I came across this: > a <- matrix(1:30,5,6) > is.integer(apply(a,1,sum)) [1] TRUE > is.integer(rowSums(a)) [1] FALSE > so rowSums() returns a float. Why is this? -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] "flipping" vector and matrix
Sorry misread the question. Function arev() takes a second argument that specifies which dimensions to flip: > b <- matrix(1:9,3,3,byrow=T) > b [,1] [,2] [,3] [1,]123 [2,]456 [3,]789 > arev(b,1) [,1] [,2] [,3] [1,]789 [2,]456 [3,]123 > On 23 Oct 2007, at 11:45, Rainer M Krug wrote: > Hi > > I have a vector > > x <- c(1, 2, 3, 4, 5) > > and I want to "flip" it around, i.e. I need > > 5, 4, 3, 2, 1 > > Is there a ssolution apart from > > y <- x[length(x):1] > > > I am also looking for the same for a matrix M, i.e. > > 1 2 3 > 4 5 6 > 7 8 9 > > should become > > 7 8 9 > 4 5 6 > 1 2 3 > > again, I am using > > M[1:dim(M)[1], dim(M)[2]:1] > > > Thanks > > Rainer > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] "flipping" vector and matrix
Hello for vectors, use rev() For matrices and arbitrary-dimensioned arrays, use arev() of the magic package: > library(magic) > b <- matrix(1:9,3,3) > b [,1] [,2] [,3] [1,]147 [2,]258 [3,]369 > arev(b) [,1] [,2] [,3] [1,]963 [2,]852 [3,]741 > On 23 Oct 2007, at 11:45, Rainer M Krug wrote: > Hi > > I have a vector > > x <- c(1, 2, 3, 4, 5) > > and I want to "flip" it around, i.e. I need > > 5, 4, 3, 2, 1 > > Is there a ssolution apart from > > y <- x[length(x):1] > > > I am also looking for the same for a matrix M, i.e. > > 1 2 3 > 4 5 6 > 7 8 9 > > should become > > 7 8 9 > 4 5 6 > 1 2 3 > > again, I am using > > M[1:dim(M)[1], dim(M)[2]:1] > > > Thanks > > Rainer > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] create diagonal matrix within a for cycle
Hello
I suspect that adiag() of the magic package is what you need:
> library(magic)
Loading required package: abind
Attaching package: 'magic'
The following object(s) are masked from package:mgcv :
magic
> a <- matrix(1:4,2,2)
> b <- matrix(1:9,3,3)
> adiag(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,]13000
[2,]24000
[3,]00147
[4,]00258
[5,]00369
>
HTH
rksh
On 23 Oct 2007, at 14:14, Niccolò Bassani wrote:
> Dear R users, I'm trying to build a diagonal matrix from a group of
> matrices. These matrices have been built in a for loop. That is, I've
> subsetted a preliminar matrix to obtain a certain number of square
> sub-matrices, and now I need to create a diagonal out of these.
> Suppose my matrix is a square matrix 21x21, and that I want to
> diagonalize
> 616 submatrices of this one. To do this I do:
>
> diag <- rep(1,21)
> work.cov <- matrix(0,21,21)
> diag(work.cov) <- diag
> samp <- sample(seq(1:21),616,replace=T)
> for (i in 1:length(samp)){
> A <- work.cov[1:samp[i],1:samp[i]]
> }
>
> Now, I want to put these 616 square matrices on the diagonal of a new
> matrix. The question is: how I can I do this? I know there's a
> function
> bdiag that creates diagonal matrix with various elements (vectors and
> matrices), but the problem here's that my matrices exists only in
> the loop,
> not outside of it, and they all correspond to the same matrix V
> computed
> under different values of i.
> Thansk in advance
> niccolò
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] Multiply a 3D-array by a vector (weighted combination of matrices)
Hi you need the tensor library: > library(tensor) > z <- array(runif(27),rep(3,3)) > w <- runif(3) > w[1]* z[,,1] + w[2]*z[,,2] + w[3]*z[,,3] [,1] [,2] [,3] [1,] 1.2700333 1.1920113 0.8015904 [2,] 0.5175217 0.7808569 0.6306053 [3,] 0.8386015 0.6143882 0.6382314 > tensor(z,w,3,1) [,1] [,2] [,3] [1,] 1.2700333 1.1920113 0.8015904 [2,] 0.5175217 0.7808569 0.6306053 [3,] 0.8386015 0.6143882 0.6382314 > > w[1]* z[,,1] + w[2]*z[,,2] + w[3]*z[,,3] - tensor(z,w,3,1) [,1] [,2] [,3] [1,] -2.220446e-1600 [2,] 0.00e+0000 [3,] 0.00e+0000 > HTH rksh On 17 Oct 2007, at 08:58, Yvonnick NOEL wrote: > Hello, > > I would like to compute a weighted combination of matrices. > > I have a number of matrices, arranged in a 3D-array, say: > > z = array(rep(1:3,c(9,9,9)),c(3,3,3)) > > so that z[,,1] is my first matrix, and z[,,2] and z[,,3] the second > and > third one, and a vector of coefficients: > > w = rep(1/3,3) > > I would like to compute: > > w[1]* z[,,1] + w[2]*z[,,2] + w[3]*z[,,3] > > I could of course do this using a for() loop, but would like to > know if > there is a way to do it in a "vectorized" manner, or any other way > that > is likely to result in faster computation. > > Any hint ? > > Thank you very much in advance, > > YNOEL > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] reference for logistic regression
> [snip] > > Whenever you use wikipedia you should be cautious of the quality of > the information in the articles. Generally the articles are good as a > brief introduction but they can and do contain errors so you should > check important facts and not take them at face value. A person in > one of my classes asked about the standard deviation and I suggested > that they look at the wikipedia article on the topic. Then I looked > at it myself and saw that one of the things mentioned is that the > standard deviation of the Cauchy distribution is undefined, which is > true, but the reason given is because E[X] is undefined, which is not > true. > The variance page has Many distributions, such as the Cauchy distribution, do not have a variance because the relevant integral diverges. In particular, if a distribution does not have an expected value, it does not have a variance either. The converse is not true: there are distributions for which the expected value exists, but the variance does not. which seems to be right. I'd say that A => B but B =/> A together with A being true would be expressed as "B is true because of A" which is pretty much what the standard deviation page says. Is this what you meant? -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] rearrange data columns
On 11 Oct 2007, at 12:55, Martin Ivanov wrote: > Dear R users, > I need to to the the following. Let a= 1 2 3 > 4 5 6 > and b= -1 -2 -3 be (2x3) matrices. > -4 -5 -6 > I need to combine the two matrices into a new (2x6) matrix like this: > > ab = ( 1 -1 2 -2 3 -3 ) > 4 -4 5 -5 6 -6 > > How can this be done in R? > > > a [,1] [,2] [,3] [1,]123 [2,]456 > b [,1] [,2] [,3] [1,] -1 -2 -3 [2,] -4 -5 -6 > x <- cbind(a,b)+NA > x [,1] [,2] [,3] [,4] [,5] [,6] [1,] NA NA NA NA NA NA [2,] NA NA NA NA NA NA > x[,seq(from=1,by=2,len=3)] <- a > x[,seq(from=2,by=2,len=3)] <- b > x [,1] [,2] [,3] [,4] [,5] [,6] [1,]1 -12 -23 -3 [2,]4 -45 -56 -6 > HTH rksh > > - > Крайна цел - Да оцелееш! www.survivor.btv.bg > > __ > R-help@r-project.org 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. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org 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] Multivariate chi-square distribution function
Hello
If you want the multivariate t-distribution,
use rmvt() of the mvtnorm package.
If you want the Wishart distribution, one
can write a little nonce function:
library(mvtnorm)
rwis <- function(n,sigma){
crossprod(rmvnorm(n,sigma=sigma))
}
I'm not sure what you mean by the multivariate gamma
distribution, but a good candidate might be the
Dirichlet. This is implemented in several packages;
see ?rdiric of the VGAM package for one example.
HTH
rksh
On 9 Oct 2007, at 18:39, [EMAIL PROTECTED] wrote:
> Dear All,
>
> Is there any function in R for computing "multivariate chi-square
> distribution"?
> How about "multivariate gamma distribution"?
> I appreciate any comment on this subject.
>
> Thank you,
>
> Amin Zollanvari
> PhD student
> Department of Electrical and Computer Engineering,
> Texas A&M University,
> College Station, TX
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Question about Distinguishable Permutation
On 6 Oct 2007, at 20:05, house-ball wrote: > Hi, > > How can I list all possilities of Distinguishable Permutation by > R? For example, N=6, n1=n2=n3=2, the total possible answers are 6!/ > (2!2!2!)=90. > > Please help me. > Hello packages "gregmisc", "combinat", and "partitions" include related functionality. However, none of these (AFAICS) solves your problem, which is now on my Things To Do List [don't hold your breath]. What is your application? > > Thank you So much. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] error installing gsl pkg
Hello Randy
I get emails like this quite a lot.
The most likely problem is in the installation of
the gsl library. To verify that it is in fact installed,
try to compile and run the little Bessel function
example given in gsl-ref, section 2.1.
Get this working first. If it works, this means
that everything is where it should be, and the
configuration script of the R package
should detect this fact.
While I'm writing, version 1-10 of
GSL came out the other day, and
the current configure script fails with this version.
I'll upload a fixed version of the gsl
library to CRAN when I get a minute (patch
from Matt Clegg gratefully acknowledged).
HTH
rksh
On 3 Oct 2007, at 14:24, Randy Heiland wrote:
> Newbie here (to R) and running Linux...
>
>> install.packages("gsl","~/R")
> ...
> trying URL 'http://cran.wustl.edu/src/contrib/gsl_1.8-4.tar.gz'
> Content type 'application/x-tar' length 57051 bytes
> opened URL
> ==
> downloaded 55Kb
>
> * Installing *source* package 'gsl' ...
> checking for gcc... gcc
> checking for C compiler default output... a.out
> checking whether the C compiler works... yes
> checking whether we are cross compiling... no
> checking for suffix of executables...
> checking for suffix of object files... o
> checking whether we are using the GNU C compiler... yes
> checking whether gcc accepts -g... yes
> checking for gcc option to accept ANSI C... none needed
> checking for gsl_sf_airy_Ai_e in -lgsl... no
> configure: error: Cannot find Gnu Scientific Library.
> ERROR: configuration failed for package 'gsl'
>
>
> and I have gsl:
> LD_LIBRARY_PATH=/N/soft/linux-sles9-ppc64/gsl-1.8-xlc/lib
>
> and, fwiw:
> /N/soft/linux-sles9-ppc64/R-2.5.0-ibm-64/lib/R/bin/exec> file R
> R: ELF 64-bit MSB executable, cisco 7500, version 1 (SYSV), for GNU/
> Linux 2.4.21, dynamically linked (uses shared libs), not stripped
>
> tia, Randy
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Sweave: tables vs matrices
Hi Gavin
thanks for that. . . it does 99% of what I wanted.
I'd forgotten about the na.print argument.
It's considerably nicer than my other solution
which converted to character, then jj[is.na(jj)] <- "-"
then noquote(jj).
But sometimes I just need nice LaTeX tables
and I can't think of a way to arrange things
so that: (i) I have only one set of numbers to maintain,
and (ii) an NA appears as a "-" in the LaTeX table.
best wishes
rksh
On 14 Sep 2007, at 09:52, Gavin Simpson wrote:
> On Fri, 2007-09-14 at 09:34 +0100, Robin Hankin wrote:
>> Hello everyone
>>
>>
>> I am preparing a document using Sweave in which I want my matrices
>> to appear as tables. I am running into problems because as my
>> Rnw files stand, I have to change table entries twice, once for
>> the matrix and once for the typeset table.
>>
>> I have lots of material like the following. How can I arrange
>> my Rnw file so that I only have to change one set of figures
>> when my numbers change?
>>
>> One reason I prefer tables here is that the NA entries
>> appear as "-" in the table, but as "NA" in the Schunk.
>> Is there a way to make the Schunk typeset NAs
>> as minuses?
>
> See ?print.default and its argument na.print:
>
>> print.default(jj, na.print = "-")
> [,1] [,2] [,3] [,4] [,5]
> [1,]2341 10
> [2,]057- 12
> [3,]37-4 14
> [4,]2--24
> [5,]7 15 117 40
>
> Is that what you meant? It still prints the [1,] bits...
>
> HTH
>
> G
>
>>
>>
>>
>> \begin{table}
>> \centering
>> \begin{tabular}{||c|}\hline
>> \multicolumn{4}{|c|}{brand}&\\ \hline
>> A&B&C&D&total\\ \hline
>> 2 & 3 & 4 & 1& 10 \\
>> 0 & 5 & 7 & -& 12 \\
>> 3 & 7 & - & 4& 14 \\
>> 2 & - & - & 2& 4\\ \hline
>> 7&15&11&7&40\\ \hline
>> \end{tabular}
>> \caption{snipped caption}
>> \end{table}
>>
>>
>> <<>>=
>> jj <- matrix(c(2, 3, 4, 1,
>> 0, 5, 7, NA,
>> 3, 7, NA, 4,
>> 2, NA, NA, 2
>> ),byrow=TRUE,nrow=4)
>> jj <- rbind(jj,apply(jj,2,sum,na.rm=TRUE))
>> jj <- cbind(jj,apply(jj,1,sum,na.rm=TRUE))
>> jj
>> @
>>
>>
>>
>
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] Sweave: tables vs matrices
Hello everyone
I am preparing a document using Sweave in which I want my matrices
to appear as tables. I am running into problems because as my
Rnw files stand, I have to change table entries twice, once for
the matrix and once for the typeset table.
I have lots of material like the following. How can I arrange
my Rnw file so that I only have to change one set of figures
when my numbers change?
One reason I prefer tables here is that the NA entries
appear as "-" in the table, but as "NA" in the Schunk.
Is there a way to make the Schunk typeset NAs
as minuses?
\begin{table}
\centering
\begin{tabular}{||c|}\hline
\multicolumn{4}{|c|}{brand}&\\ \hline
A&B&C&D&total\\ \hline
2 & 3 & 4 & 1& 10 \\
0 & 5 & 7 & -& 12 \\
3 & 7 & - & 4& 14 \\
2 & - & - & 2& 4\\ \hline
7&15&11&7&40\\ \hline
\end{tabular}
\caption{snipped caption}
\end{table}
<<>>=
jj <- matrix(c(2, 3, 4, 1,
0, 5, 7, NA,
3, 7, NA, 4,
2, NA, NA, 2
),byrow=TRUE,nrow=4)
jj <- rbind(jj,apply(jj,2,sum,na.rm=TRUE))
jj <- cbind(jj,apply(jj,1,sum,na.rm=TRUE))
jj
@
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] vectorize a matrix conversion
Hello
I have X, an n-by-n matrix and want to convert it to Y, an
n(n-1)/2 -by- n matrix such that each row of Y
corresponds to an element of the upper diagonal
of X. Say row k of Y corresponds to [i,j] with i\neq j.
Then Y[i,k] = X[i,j] and Y[j,k] = X[j,i].
and Y[-c(i,j),k] = NA.
How to do this vectorizedly?
Example follows:
> X
[,1] [,2] [,3] [,4]
[1,] NA 1087
[2,] 10 NA7 12
[3,] 12 13 NA8
[4,] 138 12 NA
> Y
[,1] [,2] [,3] [,4]
[1,] 10 10 NA NA
[2,] 12 NA8 NA
[3,] 13 NA NA7
[4,] NA 137 NA
[5,] NA8 NA 12
[6,] NA NA 128
>
[matrix X corresponds to an all-play-all competition amongst 4
individuals,
entry [i,j] corresponding to the number of times individual "i" won
when competing against individual "j". Thus individual 2 beat
individual
3 seven times and individual 3 beat individual 2 thirteen times.
Note X[i,j] + X[j,i]=20 as there were 20 trials for each pair]
Pitiful nonvectorized code follows.
n <- nrow(X)
Y <- matrix(NA,n*(n-1)/2,n)
k <- 1
for(i in 1:(n-1)){
for(j in (i+1):n){
if( !(i==j)){
print(c(i,j,k))
Y[k,i] <- X[i,j]
Y[k,j] <- X[j,i]
}
k <- k+1
}
}
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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] Passing parameters to 'optim' fn function
Hi
Use the dots argument at the end to pass in named parameters.
This is documented as the last in the argument list under ?optim.
Here's a toy example using x and p:
> f <- function(x,p){sum(x^2)+p*sum(x)}
> gr <- function(x,p){2*x+p*rep(1,length(x))}
> optim(par=c(10,10),fn=f,gr=gr,p=13)
$par
[1] -6.498188 -6.499652
$value
[1] -84.5
$counts
function gradient
63 NA
$convergence
[1] 0
$message
NULL
HTH
rksh
On 12 Sep 2007, at 09:32, José Luis Aznarte M. wrote:
> Hi again! I'm using the 'optim' method to fix the parameters of a
> model. I have written the function to be minimised and another
> function
> which returns the gradient of the error. Now my problem is that, in
> order to compute that gradient, it would be extremely convenient to be
> able to pass some parameters to the gradient function. I don't see how
> to do it given the fixed syntax of 'optim', which does not allow for
> parameters being passed to fn and gr:
>
>> optim(par, fn, gr = NULL, ...)
>
> Of course the first idea would be to "pack" the extra
> parameters in
> the vector 'par', but in that case the extra parameters would be
> optimized too.
> Does anyone have an idea on how to pass parameters to gr in optim?
> Thanks for your time!
>
> __
> R-help@r-project.org 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.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@r-project.org 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.
