Re: [R] Inverse of many small matrices

2015-04-02 Thread Feng Li

On 04/02/2015 10:40 PM, Duncan Murdoch wrote:

On 02/04/2015 9:31 AM, Feng Li wrote:

Dear all,

I am working with a likelihood function that requires the inverse of
many small covariance matrices for multivariate normal densities. When
the sample size is large, this calculation is really heavy. Those
matrices are independent but unfortunately I can hardly find a way to
vectorize them.

Can anyone give me a hint to speed this up? Thanks in advance!



Are you sure you need the inverses of those matrices?  For example, if
you are trying to compute x^t Ainv x,
where Ainv is A inverse, the naive calculation is t(x) %*% solve(A) %*%
x, but that's likely slower and less accurate than other equivalent
ones, such as x %*% solve(A, x), and I wouldn't be surprised if there
are better ones.

Duncan Murdoch


I agree what you suggested is good practice. But my likelihood function 
needs to calculate x %*% solve(A_i, x) for i=1,...,n, which is the 
bottom neck.



Feng

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Inverse of many small matrices

2015-04-02 Thread Feng Li

Dear all,

I am working with a likelihood function that requires the inverse of 
many small covariance matrices for multivariate normal densities. When 
the sample size is large, this calculation is really heavy. Those 
matrices are independent but unfortunately I can hardly find a way to 
vectorize them.


Can anyone give me a hint to speed this up? Thanks in advance!



Feng

--
Feng Li, Ph.D.
School of Statistics and Mathematics
Central University of Finance and Economics
100081 Beijing, China
http://feng.li/

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Use BatchJobs package to submit Rmpi job

2015-03-19 Thread Feng Li
Hi,
BatchJobs package is popular to submit R batch jobs on various job scheduler. 
Now we are facing the question is to submit a Rmpi jobs. So we want to know how 
we can leverage BatchJobs to do so. After a bit investigation/enhancement, we 
tested below approach on OpenLava job scheduler. we want to get some inputs 
from the group to see how this makes sense to you. thanks,

1. R script that submit a Rmpi job using BatchJobs:
library(BatchJobs)
conf = BatchJobs:::getBatchJobsConf()
conf$cluster.functions = 
makeClusterFunctionsOpenLava(/home/clusteradmin/rmpi.tmpl)  -- specify 
job template file that includes mpirun to start R job.
reg = makeRegistry(id = RMPIJobsExample)
f = function(x) {-- RMPI 
algorithm function. This is the normal Rmpi R script.
library(Rmpi)
...
mpi.spawn.Rslaves(nslaves=slaveno)
...
mpi.remote.exec(paste(I 
am,mpi.comm.rank(),of,mpi.comm.size()))
...
}
batchMap(reg, f, 1)   -- specify 1 Rmpi job
submitJobs(reg, np=2)  -- specify number of 
processors to use
showStatus(reg)
waitForJobs(reg)
showStatus(reg)

2. the job template: rmpi.tmpl
#BSUB -q mpi_job_queue
/opt/openlava/bin/openmpi-mpirun 
-np 1 R CMD BATCH %= rscript %


Please let us know your comments,

Thanks,
Feng Li
Teraproc Inc.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] The three-dot question

2013-01-14 Thread Feng Li
Dear all,

Why does not the three-dot accept arguments from the parent environment?
I am just confused with this error, can someone give me a hint? 

 rm(list=ls())
 testFun - function(a, ...)
+   {
+ if(a){
+ print(a)
+   }else
+   {
+ print(b)
+   }
+   }
 
 myTask - function(a)
+   {
+ b - 3
+ testFun(a, b = b)
+   }
 myTask(FALSE)
Error in print(b) : object 'b' not found


Thanks in advance!

Feng

-- 
Feng Li
Department of Statistics
Stockholm University
SE-106 91 Stockholm, Sweden
http://feng.li/

__
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] The three-dot question

2013-01-14 Thread Feng Li
Hi Michael,

Thanks for the reply.

On Mon, 2013-01-14 at 10:33 +, R. Michael Weylandt wrote: 
 Hi Feng,
 
 I'm afraid I don't entirely understansd your question -- the `...`
 construct only allows you to pass variable numbers of arguments, not
 to have arbitrary access to the parent frames. You need to manually
 extract b from the dots inside of testFun.
 
 Also, it's quite frowned upon to put ##rm(list = ls())## in your
 examples: it's the mailing list equivalent of asking a buddy to come
 over to help you move and then punching him in the face when he tries
 to lift your sofa.

The rm line just wants to declare there is no b in the global
environment. I did not mean people ought to try it. Also notice that
there is  in front of the rm. You won't get hurt if you just copy and
paste it. If it bites you, the mouse and keyboard did it:) 


Cheers,

Feng

 
 Cheers,
 Michael
 
 
 On Mon, Jan 14, 2013 at 10:21 AM, Feng Li m...@feng.li wrote:
  Dear all,
 
  Why does not the three-dot accept arguments from the parent environment?
  I am just confused with this error, can someone give me a hint?
 
  rm(list=ls())
  testFun - function(a, ...)
  +   {
  + if(a){
  + print(a)
  +   }else
  +   {
  + print(b)
  +   }
  +   }
 
  myTask - function(a)
  +   {
  + b - 3
  + testFun(a, b = b)
  +   }
  myTask(FALSE)
  Error in print(b) : object 'b' not found
 
 
  Thanks in advance!
 
  Feng
 
  --
  Feng Li
  Department of Statistics
  Stockholm University
  SE-106 91 Stockholm, Sweden
  http://feng.li/
 
  __
  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.

-- 
Feng Li
Department of Statistics
Stockholm University
SE-106 91 Stockholm, Sweden
http://feng.li/

__
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] The three-dot question

2013-01-14 Thread Feng Li
That makes sense. Thanks!


Feng


On Mon, 2013-01-14 at 09:09 -0500, Mark Leeds wrote:
 Hi: If you want testFun to know about b, then you would have to do
 b-list(...)$b inside
 TestFun itself.  But the dot dot dot argument is not really for that
 purpose. 
 
 The use of dotdotdot is for the case where a function INSIDE testFun
 has a formal argument named say b. Then you can pass the ... at the
 top level and the function inside will receive the ... but
 automatically slurp the b out of the dot dot dot and know about it.
 Below is an example of the use I'm talking about. It was created by a
 mentoR when I was trying to understand the dotdotdot concept. 
 hope it helps you. 
 
 #==
 
 Note that f does not directly know about x and g does.  That is g
 knows about x even though f does
 not know and g got it from f.
 
 g - function(x, ...) { cat(g: exists('x') = , exists(x), \n);
 list(...) }
 f - function(...) { cat(f: exists('x') = , exists(x), \n);
 g(...) }
 f(x = 3, y = 4)
 f: exists('x') =  FALSE
 g: exists('x') =  TRUE
 $y
 
 
 
 
 On Mon, Jan 14, 2013 at 5:21 AM, Feng Li m...@feng.li wrote:
 Dear all,
 
 Why does not the three-dot accept arguments from the parent
 environment?
 I am just confused with this error, can someone give me a
 hint?
 
  rm(list=ls())
  testFun - function(a, ...)
 +   {
 + if(a){
 + print(a)
 +   }else
 +   {
 + print(b)
 +   }
 +   }
 
  myTask - function(a)
 +   {
 + b - 3
 + testFun(a, b = b)
 +   }
  myTask(FALSE)
 Error in print(b) : object 'b' not found
 
 
 Thanks in advance!
 
 Feng
 
 --
 Feng Li
 Department of Statistics
 Stockholm University
 SE-106 91 Stockholm, Sweden
 http://feng.li/
 
 __
 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. 
 

-- 
Feng Li
Department of Statistics
Stockholm University
SE-106 91 Stockholm, Sweden
http://feng.li/

__
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] #!/usr/bin/env Rscript --vanilla ??

2012-04-13 Thread Feng Li

I guess you may try this

http://modules.sourceforge.net/

I have seen a lot of clusters using this to manage software of different 
versions.



Feng

On 04/13/2012 10:32 AM, Martin Maechler wrote:

I think that's my first true question (rather than answer)
to R-help.

As R has, for a long time, become my primary scripting and
programming language, I'm prefering at times to write  Rscript
files instead of shell scripts, notably when R has nice ways to
do some of the things.
On a standard standalone platform with standard R,
I would start such a script with
---
#! /usr/bin/Rscript --vanilla
---
(yes, the --vanilla is important to me, in this case)

However; as, at work, my scripts have to work correctly on quite a
few different (unixy : several flavors of Linux, Solaris, MacOS X) platforms,
*and* as an R developer, I have many different versions of R
installed simultaneously, using /usr/bin/Rscript  is not an
option.
Rather, I'd use the /usr/bin/env trick :

---
#! /usr/bin/env Rscript
---

which finds Rscript in the correct place, according to the
current PATH.  All fine till now.

PROBLEM:  It does not work with '--vanilla' or any other argument:
  If I start my script with 
#! /usr/bin/env Rscript --vanilla
  the error message simply is
/usr/bin/env: Rscript --vanilla: No such file or directory

I have tried a few variations on the theme, using quotes in
different places, but have not succeeded till now.
Any suggestions?

Martin Maechler, ETH Zurich



--
Feng Li
Department of Statistics
Stockholm University
SE-106 91 Stockholm, Sweden
http://feng.li/

__
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] #!/usr/bin/env Rscript --vanilla ??

2012-04-13 Thread Feng Li

I guess you may try this

http://modules.sourceforge.net/

I have seen a lot of clusters using this to manage software of different 
versions.



Feng

On 04/13/2012 10:32 AM, Martin Maechler wrote:

I think that's my first true question (rather than answer)
to R-help.

As R has, for a long time, become my primary scripting and
programming language, I'm prefering at times to write  Rscript
files instead of shell scripts, notably when R has nice ways to
do some of the things.
On a standard standalone platform with standard R,
I would start such a script with
---
#! /usr/bin/Rscript --vanilla
---
(yes, the --vanilla is important to me, in this case)

However; as, at work, my scripts have to work correctly on quite a
few different (unixy : several flavors of Linux, Solaris, MacOS X) platforms,
*and* as an R developer, I have many different versions of R
installed simultaneously, using /usr/bin/Rscript  is not an
option.
Rather, I'd use the /usr/bin/env trick :

---
#! /usr/bin/env Rscript
---

which finds Rscript in the correct place, according to the
current PATH.  All fine till now.

PROBLEM:  It does not work with '--vanilla' or any other argument:
  If I start my script with 
#! /usr/bin/env Rscript --vanilla
  the error message simply is
/usr/bin/env: Rscript --vanilla: No such file or directory

I have tried a few variations on the theme, using quotes in
different places, but have not succeeded till now.
Any suggestions?

Martin Maechler, ETH Zurich



--
Feng Li
Department of Statistics
Stockholm University
SE-106 91 Stockholm, Sweden
http://feng.li/

__
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] ess-tracebug to open a file

2012-03-13 Thread Feng Li

Dear all,

First I would like to thank the ESS people's all the hard work. I am 
watching the project closely and witnessing the improvements day by day.


Besides I found a strange situation using `ess-tracebug'. Please tell me 
if I am wrong or this is a bug.


Start Emacs with emacs -Q and load ESS and enable ess-tracebug.

Create a file named `testFun1.R' with the following contents

testFun1-function(x)
  {
y - x+2
browser()
return(y)
  }

testFun1(2)


Then start R under the same directory and run the command

source(testFun1.R)


which will show the following information

source(testFun1.R)

Called from: testFun1(2)
Browse[1] debug at testFun1.R#5: return(y)


Then if I left click `testFun1.R#5' I get this error

Wrong type argument: listp, [cl-struct-compilation--message (nil 5 
((testFun1.R nil) nil (5 #1)) nil nil) 2 nil]


Here are some system information might be useful

Emacs version

GNU Emacs 24.0.94.1 (x86_64-unknown-linux-gnu, GTK+ Version 3.2.3) of 
2012-02-27 on nova


ESS version

SVN trunk@4690


Linux version

3.2.0-2-amd64 (Debian 3.2.9-1) (debian-ker...@lists.debian.org) (gcc version 
4.6.3 (Debian 4.6.3-1) ) #1 SMP Sun Mar 4 22:48:17 UTC 2012



Best regards,

Feng

--
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

__
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] While loop working with TRUE/FALSE?

2012-02-01 Thread Feng Li

I guess you want something like

w=F
z - c(0,1,2,3,4,5,6,7,8,9)
r - 7

   while(w == F) {
for ( i in 1:10 ){
w - r == z[i]
print(w)
}
}

But this loop will run forever. The condition for while is checked 
when i jumps from 10 to 1. At that moment w is always FALSE.


Feng

On 02/01/2012 05:12 PM, Berend Hasselman wrote:


On 01-02-2012, at 16:55, Chris82 wrote:


Hi R users,

is there any possibilty that a while loop is working like that:

z- c(0,1,2,3,4,5,6,7,8,9)
r- 7

   while(w == T) {
for ( i in 1:10 ){
w- r == z[i]
print(w)
}
}


The loop should stop if w == TRUE


1. This won't work since the variable w doesn't exist at the start of the while 
loop.
2. The while loop condition is incorrect: it should be while( w == FALSE)
3. Don't use T and F. Use TRUE and FALSE.
4. If there is no element of z equal to r the loop will run forever (with the 
correct condition of course).

A much better way of doing what you seem to be wanting is

which(z == r)

or possibly

any(z == r)

Berend



--
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

 * English - detected
 * Chinese
 * English
 * Swedish

 * Chinese
 * English
 * Swedish

 javascript:void(0); #

__
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 help for linux

2011-04-19 Thread Feng Li
Dear R,

Usually when I want to read pdf help for R functions under Linux, I do
the following steps

 help(par,help_type=pdf)
 system(paste(getOption(pdfviewer),par.pdf))

which works well. First question, can I let help to save the pdf
file to a different place other than current working directory?

Now I want to write a small function to do the above together.

 viewpdfhelp - function(topic, reader = getOption(pdfviewer))
+   {
+ genpdf - help(topic = topic, help_type = pdf) # save pdf file
+ pdffile - paste(topic, .pdf, sep = ) # the pdf file name
+ system(paste(reader, pdffile)) #call the viewer
+   }

 viewpdfhelp(options)

To my surprise help does not generate the pdf file at all this time.
Can someone tell me what's going wrong here? Many thanks!


 sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.utf8   LC_NUMERIC=C
 [3] LC_TIME=en_US.utf8LC_COLLATE=en_US.utf8
 [5] LC_MONETARY=C LC_MESSAGES=en_US.utf8
 [7] LC_PAPER=en_US.utf8   LC_NAME=C
 [9] LC_ADDRESS=C  LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

loaded via a namespace (and not attached):
[1] tools_2.13.0


-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

__
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 gui problem for windows

2011-03-23 Thread Feng Li
Dear R,

I rare use the standard R-gui on Windows. Yesterday I tried the latest
stable release on Windows 7 and XP and found one thing interesting. Assume
currently I am running R code, say

 example(plot)

Then I click the close window button on the R main window. R asks me to
save workspace image or not. Then I click cancel. I suspect that my
program example(plot) should continue but it is interrupted eventually.

I would like to know the mechanism behind this and wonder whether this can
be improved or not. It is a minor problem but I believe it can be very
convenient for the window users especially when they are running time
consuming simulations. What do you think?


Best,

Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] R gui problem for windows

2011-03-23 Thread Feng Li
Thank you for spending time to point me to the posting guide. I
apology for the unclear statement for the R version. It is 32bit
2.12.2.

I love computer programs because I can always undo  what I did
before. For example I did not intend to close the program but I hit
the mouse by mistake. Then I thought the cancel command can help me.
And a lot programs have this option, so does R gui for windows.

What I found in R gui is that the click even terminates the running
task first  and then asks for saving image. But if the user regrets
what he did, it it too late to go back to the running task. I though
maybe it is possible to put a top layer to check if the user clicked
cancel, then do nothing, otherwise, do as usual did. I know how to
run R in a batched mode and all I said is based on user experience
point of view. A good user experience will let more people love R.



That is all my concern. Have a good day!

Feng

On Wed, Mar 23, 2011 at 11:30 AM, Prof Brian Ripley
rip...@stats.ox.ac.uk wrote:

 On Wed, 23 Mar 2011, Feng Li wrote:

 Dear R,

 I rare use the standard R-gui on Windows. Yesterday I tried the latest
 stable release on Windows 7 and XP and found one thing interesting. Assume

 You were asked in the posting guide (which clearly you have ignored as you 
 sent HTML) to be accurate about your version information.  We don't describe 
 any version as 'latest stable'.

 currently I am running R code, say

 example(plot)

 Then I click the close window button on the R main window. R asks me to
 save workspace image or not. Then I click cancel. I suspect that my
 program example(plot) should continue but it is interrupted eventually.

 I would like to know the mechanism behind this and wonder whether this can
 be improved or not. It is a minor problem but I believe it can be very
 convenient for the window users especially when they are running time
 consuming simulations. What do you think?

 How to run batch jobs on Windows is documented in the rw-FAQ (which the 
 posting guide pointed you to).  Rgui (sic) is not intended for batch jobs, 
 and like all other Windows GUIs I have ever seen, closing the main windows 
 terminates the process.



 Best,

 Feng

 --
 Feng Li
 Department of Statistics
 Stockholm University
 106 91 Stockholm, Sweden
 http://feng.li/

        [[alternative HTML version deleted]]

 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html

 That does mean you 

 --
 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



--
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

__
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] Singularity problem

2011-03-16 Thread Feng Li
Dear R,

If I have remembered correctly, a square matrix is singular if and only if
its determinant is zero. I am a bit confused by the following code error.
Can someone give me a hint?

 a - matrix(c(1e20,1e2,1e3,1e3),2)
 det(a)
[1] 1e+23
 solve(a)
Error in solve.default(a) :
  system is computationally singular: reciprocal condition number = 1e-17


Thanks in advance!


Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] Singularity problem

2011-03-16 Thread Feng Li
Thanks all of you for the very interesting discussion.  I think I get it 
now.


Feng

On 03/16/2011 09:25 PM, rex.dw...@syngenta.com wrote:

Feng,
Your matrix is *not* (practically) singular; its inverse is.
The message said that the *system* was singular, not the matrix.
Remember Cramer's Rule:  xi = |Ai| / |A|
The really, really large determinant of your matrix is going to appear in the 
denominator of your solutions, so, essentially, you get underflow. Try working 
out the entire solution with Cramer's Rule if you still don't see the problem.  
solve doesn't really use Cramer's Rule, but it will give you a feel for the 
issue.
HTH
Rex

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Berend Hasselman
Sent: Wednesday, March 16, 2011 1:33 PM
To: r-help@r-project.org
Subject: Re: [R] Singularity problem


Peter Langfelder wrote:

On Wed, Mar 16, 2011 at 8:28 AM, Feng Lilt;m...@feng.ligt; wrote:

Dear R,

If I have remembered correctly, a square matrix is singular if and only
if
its determinant is zero. I am a bit confused by the following code error.
Can someone give me a hint?


a- matrix(c(1e20,1e2,1e3,1e3),2)
det(a)

[1] 1e+23

solve(a)

Error in solve.default(a) :
  system is computationally singular: reciprocal condition number = 1e-17


You are right, a matrix is mathematically singular iff its determinant
is zero. However, this condition is useless in practice since in
practice one cares about the matrix being computationally singular,
i.e. so close to singular that it cannot be inverted using the
standard precision of real numbers. And that's what your matrix is
(and the error message you got says so).

You can write your matrix as

a = 1e20 * matrix (c(1, 1e-18, 1e-17, 1e-17), 2, 2)

Compared to the first element, all of the other elements are nearly
zero, so the matrix is numerically nearly singular even though the
determinant is 1e23. A better measure of how numerically unstable the
inversion of a matrix is is the condition number which IIRC is
something like the largest eigenvalue divided by the smallest
eigenvalue.


svd(a) indicates the problem.

largest singular value / smallest singular value=1e17  (condition number)
--  reciprocal condition number=1e-17
and the standard solve can't handle that.

(pivoted) QR decomposition does help. And so does SVD.

Berend

--
View this message in context: 
http://r.789695.n4.nabble.com/Singularity-problem-tp3382093p3382465.html
Sent from the R help mailing list archive at Nabble.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.



message may contain confidential information. If you are not the designated 
recipient, please notify the sender immediately, and delete the original and 
any copies. Any use of the message by you is prohibited.


--
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

__
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] Wired behavior of a 2-by-2 matrix indicies

2011-02-26 Thread Feng Li
Dear R,

I found a very wired behavior for a 2-by-2 matrix, see this example

 A - matrix(1:4, 2)
 idx4A - matrix(1:4, 2)
 A[idx4A]
Error in A[idx4A] : subscript out of bounds

But other matrices are fine,

 B - matrix(1:9, 3)
 idx4B - matrix(1:9, 3)
 B[idx4B]
[1] 1 2 3 4 5 6 7 8 9


I can reproduce this for both 32bit windows and 64bit linux with R 2.12.1
and some earlier versions.

Can someone tell me whether this a bug or a feature or something I don't
know?

Thanks!

Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] Weird behavior of a 2-by-2 matrix indicies

2011-02-26 Thread Feng Li
The real weird thing is I wrote weird and it comes out wired...

The background is that I did a lot of sparse matrices calculations.
Therefore need some indices tweaks, e.g subtract a block from a big matrix.
Then I just create an indicies matrix and use it directly as a vector since
a matrix is just a vector with dimensions attributed in R. I know I can
convert them to a vector but I did not do it anyway.

Thanks for this information!

Feng

On Sat, Feb 26, 2011 at 6:29 PM, Sarah Goslee sarah.gos...@gmail.comwrote:

 Hello,

 On Sat, Feb 26, 2011 at 12:11 PM, Feng Li m...@feng.li wrote:
  Dear R,
 
  I found a very wired behavior for a 2-by-2 matrix, see this example
 
  A - matrix(1:4, 2)
  idx4A - matrix(1:4, 2)
  A[idx4A]
  Error in A[idx4A] : subscript out of bounds

 This shouldn't work. From the help for [,

  When indexing arrays by ‘[’ a single argument ‘i’ can be a
  matrix with as many columns as there are dimensions of ‘x’;
  the result is then a vector with elements corresponding to
  the sets of indices in each row of ‘i’.

 Since there are two dimensions to A, and two rows to idx4A, then this
 form is used.
 The first row of idx4A is c(1, 3) so the first value extracted is
 A[1, 3]
 which doesn't exist, thus the error message.

  But other matrices are fine,
 
  B - matrix(1:9, 3)
  idx4B - matrix(1:9, 3)
  B[idx4B]
  [1] 1 2 3 4 5 6 7 8 9

 The part that surprises me is that this works.
 Apparently there is an implicit conversion to vector here because
 dim(B) != ncol(idx4B)

 What you seem to want is
 A[as.vector(idx4A)]
 and
 B[as.vector(idx4B)]

 but you neglected to tell us what you expected the result to be.

 Sarah
 --
 Sarah Goslee
 http://www.functionaldiversity.org




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] kmeans: number of cluster centres must lie between 1 and nrow(x)

2011-02-02 Thread Feng Li
Thank you for the suggestion and it is exactly as you said only one
observation in each cluster. I know I can avoid this anyway and I am just
out of curiosity of the error.

I am writing a special algorithm to cluster some datasets with different
numbers of observations. For some particular datasets, there is only one
observation(e.g. people died of a rare disease). kmeans() thus will not work
at this situation.


Feng

2011/2/2 Rafael Björk rafael.bj...@gmail.com

 If you change the algorithm, the function allow you to do this:

  kmeans(a, 20,algorithm=Lloyd)

 Maybe i'm missing something here, but why would you want to create as many
 clusters as there are observations? Won't the outcome just be one
 observation in each cluster?






-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] kmeans: number of cluster centres must lie between 1 and nrow(x)

2011-02-01 Thread Feng Li
Dear R,

Can't I cluster a dataset into k clusters where k is exactly the number of
observations? I have version 12.2 installed. See this example

 a - matrix(1:100, 20)
 kmeans(a, 20)
Error: number of cluster centres must lie between 1 and nrow(x)

This is a bit ad-hoc but I known R from version 2.12 allows number of
clusters to be one. So I guess allowing number of clusters to be nrow(x)
should be also possible in the future release?


Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] Sparse matrix multiplication

2011-01-30 Thread Feng Li
Thanks Martin. That is exactly what I want.

Feng

On Sat, Jan 29, 2011 at 11:59 PM, Martin Maechler 
maech...@stat.math.ethz.ch wrote:

  FL == Feng Li m...@feng.li
  on Sat, 29 Jan 2011 19:46:48 +0100 writes:

FL I meant sparse matrix, sorry for the typo.
 aha.. :-)


FL On Sat, Jan 29, 2011 at 7:02 PM, Feng Li m...@feng.li wrote:

 Dear R,

 I have a simple question concerning with a special case of
 sparse matrix multiplications. Say A is a 200-by-1 dense
 matrix. B is a 1-by-1 block- diagonal matrix, and each
 diagonal block B_i is 100-by-100. The usual way I did A%*%B
 will take about 30 seconds which is to time consuming because I
 have to do this thousands of times. I also tried to partition A
 into 100 small blocks and use mapply function to multiply by
 each B_i, but that is even slower.

 I am wondering if there is an efficient way to perform this
 type of multiplication with R?

 yes: e.g., via the (recommended, i.e. already installed) Matrix
 package's  bdiag() :

 require(Matrix)
 set.seed(1)
 A - matrix(rnorm(2e6), 200, 1)
 Blis - lapply(1:100, function(i)matrix(rnorm(1e4), 100,100))

 system.time(B. - .bdiag(Blis)) # 1.28 sec
 system.time(cc - A %*% B.) # 1.7  sec
 class(cc)# dgeMatrix .. i.e. dense
 ## and depending on the context you may revert to traditional unclassed
 matrices
 ## via
 c2 - as(cc, matrix)


 Thanks in advance!

 you are welcome.
 Martin Maechler, ETH Zurich


 Feng

 --
 Feng Li Department of Statistics Stockholm University 106 91
 Stockholm, Sweden http://feng.li/




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] Spare matrix multiplication

2011-01-29 Thread Feng Li
Dear R,

I have a simple question concerning with a special case of spare matrix
multiplications. Say A is a 200-by-1 dense matrix. B is a 1-by-1
block- diagonal matrix, and each diagonal block B_i is 100-by-100. The usual
way I did A%*%B will take about 30 seconds which is to time consuming
because I have to do this thousands of times. I also tried to partition A
into 100 small blocks and use mapply function to multiply by each B_i, but
that is even slower.

I am wondering if there is an efficient way to perform this type of
multiplication with R?

Thanks in advance!

Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] Spare matrix multiplication

2011-01-29 Thread Feng Li
I meant sparse matrix, sorry for the typo.

Feng

On Sat, Jan 29, 2011 at 7:02 PM, Feng Li m...@feng.li wrote:

 Dear R,

 I have a simple question concerning with a special case of sparse matrix
 multiplications. Say A is a 200-by-1 dense matrix. B is a 1-by-1
 block- diagonal matrix, and each diagonal block B_i is 100-by-100. The usual
 way I did A%*%B will take about 30 seconds which is to time consuming
 because I have to do this thousands of times. I also tried to partition A
 into 100 small blocks and use mapply function to multiply by each B_i, but
 that is even slower.

 I am wondering if there is an efficient way to perform this type of
 multiplication with R?

 Thanks in advance!

 Feng

 --
 Feng Li
 Department of Statistics
 Stockholm University
 106 91 Stockholm, Sweden
 http://feng.li/




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] Convert a matrix's columns to list

2011-01-18 Thread Feng Li
Dear R,

Is there an efficient way to make a list that each element is from the
corresponding column of a matrix. For example, if I have a matrix a

 a - matrix(1:10, 5, 2)
 a
 [,1] [,2]
[1,]16
[2,]27
[3,]38
[4,]49
[5,]5   10

I would like to have a list b like this

 b - list(a[, 1], a[, 2])
 b
[[1]]
[1] 1 2 3 4 5

[[2]]
[1]  6  7  8  9 10


Thanks in advance!


Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] How to disable using enter key to exit the browser in debugging mode

2011-01-13 Thread Feng Li
Thanks for this. I am an Emacs user but still having this problem with
ESS. Curious
to know how other editors manage this.

Feng



On Thu, Jan 13, 2011 at 1:29 AM, Gene Leynes
gleyne...@gmail.comgleynes%...@gmail.com
 wrote:

 That also drives me crazy!

 I don't have that problem when I use the StatEt plug-in for Eclipse.

 Of course, using a new IDE is a big undertaking, but I can assure you: it's
 worth it!  This is just one small benefit.

 On Wed, Jan 12, 2011 at 4:00 PM, Feng Li m...@feng.li wrote:

 Dear R,

 How can I disable using enter key to exit the browser() in debug mode? I
 would love to have this option because it is so annoying to jump out of
 the
 debugging mode unexpectedly when I don't want to. I guess some of us have
 encouraged at least one of these situations,

 1, Accidentally pressed the enter key within the browser.

 2, Copy and paste a piece of debugging code containing empty lines to the
 prompt within the debugging mode.

 3, If I paste a piece of code to the prompt to debug as follows, it will
 eventually jump out before I can do anything.

 ### copy starting from this line ##
 test - function()
 {
 x- 5
 browser()
 y-4
 }

 test()

  end of copy at this line 


 Any suggestions are most welcome!

 Feng

 --
 Feng Li
 Department of Statistics
 Stockholm University
 106 91 Stockholm, Sweden
 http://feng.li/

[[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.





-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] How to disable using enter key to exit the browser in debugging mode

2011-01-13 Thread Feng Li
Thanks for Prof. Ripley explaining this to me. I totally agree with you that
I can avoid this problem if I am debugging carefully enough. But it would be
very convenient to have this feature for the newbies like me.  I will submit
this as a feature request.

Feng


On Thu, Jan 13, 2011 at 12:20 PM, Prof Brian Ripley
rip...@stats.ox.ac.ukwrote:

 You do not mean the enter key: it is not that which exits the browser but
 rather a newline (which can be entered via 'return', and in other ways).

 This is part of the parser, and there is no way to turn it off.

 Somehow other experienced R users have never encountered this.  But if you
 wish you could submit this as a 'wishlist' request to R-bugs, and if enough
 users support it (and especially if someone submits a high-quality patch to
 do this), it might be added in a future release.



 On Wed, 12 Jan 2011, Feng Li wrote:

  Dear R,

 How can I disable using enter key to exit the browser() in debug mode? I
 would love to have this option because it is so annoying to jump out of
 the
 debugging mode unexpectedly when I don't want to. I guess some of us have
 encouraged at least one of these situations,

 1, Accidentally pressed the enter key within the browser.

 2, Copy and paste a piece of debugging code containing empty lines to the
 prompt within the debugging mode.

 3, If I paste a piece of code to the prompt to debug as follows, it will
 eventually jump out before I can do anything.

 ### copy starting from this line ##
 test - function()
 {
x- 5
browser()
y-4
 }

 test()

  end of copy at this line 


 Any suggestions are most welcome!

 Feng

 --
 Feng Li
 Department of Statistics
 Stockholm University
 106 91 Stockholm, Sweden
 http://feng.li/

[[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.


 --
 Brian D. Ripley,  rip...@stats.ox.ac.uk
 Professor of Applied Statistics,  
 http://www.stats.ox.ac.uk/~ripley/http://www.stats.ox.ac.uk/%7Eripley/
 University of Oxford, Tel:  +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UKFax:  +44 1865 272595




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] How to disable using enter key to exit the browser in debugging mode

2011-01-12 Thread Feng Li
Dear R,

How can I disable using enter key to exit the browser() in debug mode? I
would love to have this option because it is so annoying to jump out of the
debugging mode unexpectedly when I don't want to. I guess some of us have
encouraged at least one of these situations,

1, Accidentally pressed the enter key within the browser.

2, Copy and paste a piece of debugging code containing empty lines to the
prompt within the debugging mode.

3, If I paste a piece of code to the prompt to debug as follows, it will
eventually jump out before I can do anything.

### copy starting from this line ##
test - function()
{
 x- 5
 browser()
 y-4
}

test()

 end of copy at this line 


Any suggestions are most welcome!

Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] lapply to subsets

2010-10-13 Thread Feng Li
Yes, that is what I what...
Thanks.

Feng

On Wed, Oct 13, 2010 at 6:38 AM, Michael Bedward
michael.bedw...@gmail.comwrote:

 Hello Feng,

 I think you just want this...

 lapply(A, function(x) apply(x[,,-c(1,2)], c(1,2), mean))

 Michael


 On 13 October 2010 04:00, Feng Li feng...@stat.su.se wrote:
  Dear R,
 
  I have a silly question concerns with *apply. Say I have a list called A,
 
  A - list(a  =  array(1:20, c(2, 2, 5)), b  = array(1:30, c(2, 3, 5)))
 
  I wish to calculate the mean of A$a, and A$b w.r.t. their third dimension
 so
  I did
 
  lapply(A,apply,c(1,2),mean)
 
  Now if I still wish to do the above task but take away some burn-in, e.g.
 do
  not take A$a[,,1:2],and A$b[,,1:2] into account. How can I do then?
 
 
  Thanks!
 
 
 
  Feng
 
 
  --
  Feng Li
  Department of Statistics
  Stockholm University
  106 91 Stockholm, Sweden
  http://feng.li/
 
 [[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.
 




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] lapply to subsets

2010-10-12 Thread Feng Li
Dear R,

I have a silly question concerns with *apply. Say I have a list called A,

A - list(a  =  array(1:20, c(2, 2, 5)), b  = array(1:30, c(2, 3, 5)))

I wish to calculate the mean of A$a, and A$b w.r.t. their third dimension so
I did

lapply(A,apply,c(1,2),mean)

Now if I still wish to do the above task but take away some burn-in, e.g. do
not take A$a[,,1:2],and A$b[,,1:2] into account. How can I do then?


Thanks!



Feng


-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] Intel i7 utilization

2010-10-05 Thread Feng Li
 I guess you want a multi-threaded R. The cheapest way is to recompile 
R from the source and link it with an multi-threaded BLAS. See R-admin 
for details. See also Dirk's recent survey on this topic 
http://dirk.eddelbuettel.com/papers/gcbd.pdf I managed to do this under 
my linux box with Intel core 2 CPU which works well. But I haven't tried 
with i7 CPU.



Feng

On 10/05/2010 12:42 PM, Costas wrote:

 Hello,

Is there a way to force through R the amount of the cpu's cores (or 
the cpu's utilization level)  used under Windows 7 or Linux?


Thanks,
Costas



--
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

__
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] boundary check

2010-09-24 Thread Feng Li
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

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] boundary check

2010-09-24 Thread Feng Li
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 rk...@cam.ac.uk 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$soln1))
   }

 }
  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.

__
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] How to generate a particular sequence ?

2010-09-13 Thread Feng Li
Dear R,

I have a vector, say a = c(1,2,4,5,6,8). Can I generate a vector or array
(2-by-3-by-3) of this form c(1,2,1,2,1,2,4,5,4,5,4,5,6,8,6,8,6,8), in which
every two elements in a have been repeated twice?

I am to stupid today and could not figure this simple question out...  Many
many thanks!

Feng


-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] How to generate a particular sequence ?

2010-09-13 Thread Feng Li
Sorry. It was my typo. Should be three times as it in the example.


Feng

On Mon, Sep 13, 2010 at 11:32 AM, Jim Lemon j...@bitwrit.com.au wrote:

 On 09/13/2010 07:19 PM, Feng Li wrote:

 Dear R,

 I have a vector, say a = c(1,2,4,5,6,8). Can I generate a vector or array
 (2-by-3-by-3) of this form c(1,2,1,2,1,2,4,5,4,5,4,5,6,8,6,8,6,8), in
 which
 every two elements in a have been repeated twice?

 I am to stupid today and could not figure this simple question out...
  Many
 many thanks!

  Hi Feng,
 I would take a quick look at the help for rep and c, but I would first
 suggest that you count the number of times that the elements are to be
 repeated. While you have made the definition of the problem reasonably clear
 with your example, there are an awful lot of combinations of every two
 elements of a, whether you want to repeat them two or three times.

 Jim




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] How to generate a particular sequence ?

2010-09-13 Thread Feng Li
Thanks. That's what I want. Sorry for the typo:)

Feng

On Mon, Sep 13, 2010 at 11:42 AM, Ted Harding
ted.hard...@manchester.ac.ukwrote:

 On 13-Sep-10 09:19:21, Feng Li wrote:
  Dear R,
  I have a vector, say a = c(1,2,4,5,6,8). Can I generate a vector
  or array (2-by-3-by-3) of this form
  c(1,2,1,2,1,2,4,5,4,5,4,5,6,8,6,8,6,8)
  in which every two elements in a have been repeated twice?
 
  I am to stupid today and could not figure this simple question out...
  Many many thanks!
 
  Feng

 A possible solution (somewhat generalisable):

  a - c(1,2,4,5,6,8)
   Reps - 3

  pairs - matrix(a,nrow=2)
  as.vector(pairs[,rep(c(1,2,3),each=Reps)])
  # [1] 1 2 1 2 1 2 4 5 4 5 4 5 6 8 6 8 6 8

 (By the way, you have 3 repetitions but wrote twice -- I assume
 you meant thrice but the above generalises to 2 repetitions ... :)

 Ted.

 
 E-Mail: (Ted Harding) ted.hard...@manchester.ac.uk
 Fax-to-email: +44 (0)870 094 0861
 Date: 13-Sep-10   Time: 10:42:46
 -- XFMail --




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] a^c(1:3)

2010-09-07 Thread Feng Li
Dear R,

I have two small questions confused me recently. Now assume I have a matrix
a, like this,

 a - matrix(1:6, 2, 3)
 a
 [,1] [,2] [,3]
[1,]135
[2,]246

I sometimes need each row of a raised to a different exponent. So I do a
trick like this,

 a^c(2, 3)
 [,1] [,2] [,3]
[1,]19   25
[2,]8   64  216

My first question is that if it is possible to do this trick column wise?
Just out of curiosity, of course I know there are other ways of doing this.

And the second question is why I get such result when I put another element
in the exponent part like this,

 a^c(2, 3, 4)
 [,1] [,2] [,3]
[1,]1   81  125
[2,]8   16 1296



BTW, I have a 64bit R version (2.11) for Linux. Any advice would be
appreciated.



Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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.


Re: [R] a^c(1:3)

2010-09-07 Thread Feng Li
Very fruitful, thanks all of you:)

Feng

[[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] Speed up sparse matrices

2010-03-08 Thread Feng Li
Dear R,

I have three matrices like this

K:  pp-by-pp   commutation matrix,
I:  p-by-p diagonal matrix,
X:  p-by-q dense matrix,

and I wish to calculate

K(IoX)

where `o' denotes Kronecker product.

Can you give me any suggestion to speed it up when `p' and `q' are large?


Thanks in advance.


Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] R-package related to the topic of INARMA models

2010-01-05 Thread Feng Li
Dear R,

I am looking for R-package related to INARMA (Integer-valued ARMA). Can
anyone give me some information? I did not get information from task view.

Many thanks.


Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] rm(list-ls()) error

2009-11-05 Thread Feng Li
Dear R,

Why rm(list-ls()) gives an error but rm(list=ls()) not? I remember the
operator ‘-’ can be used anywhere...


Thanks!


Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] rm(list-ls()) error

2009-11-05 Thread Feng Li
Dear R,

Why rm(list-ls()) gives an error but rm(list=ls()) not?. I remember the
operator ‘-’ can be used anywhere...


Thanks!


Feng

-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden
http://feng.li/

[[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] How to handle large numbers?

2009-02-11 Thread Feng Li
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



-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden

[[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.


Re: [R] How to handle large numbers?

2009-02-11 Thread Feng Li
Hi Robin,

This works! Also your introduction to using S4 methods is very good!

Thanks!

Feng

On Wed, Feb 11, 2009 at 12:20 PM, Robin Hankin rk...@cam.ac.uk wrote:

 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




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden

[[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.


Re: [R] How to handle large numbers?

2009-02-11 Thread Feng Li
Thanks for your explanation. Now I am quite clear about that.

On Wed, Feb 11, 2009 at 12:05 PM, Thomas Lumley tlum...@u.washington.eduwrote:

 On Wed, 11 Feb 2009, Feng Li wrote:

  1, Why both R and Matlab give 0*Inf==NaN? To my knowledge, it should be
 zero
 mathematically. Am I right?


 No.  0*Inf is NaN according to the floating point arithmetic standards that
 R depends on.

  It's true that 0*x==0 for all finite x, but it's also true that x*Inf==Inf
 for all non-zero x, and you can't preserve both of these with 0*Inf.

  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.


 No. It should be close to 1. Try  1000/(1000+5) for a simpler example.

  My question is whether there
 is a general way to solve this kind of question or should I do some
 settings
 before computing?


 exp(1000)/(exp(1007)+5) is 1/(exp(7)+5*exp(-1000)), which is the same as
 1/exp(7) to more than 400 digits accuracy.

   -thomas

 Thomas Lumley   Assoc. Professor, Biostatistics
 tlum...@u.washington.eduUniversity of Washington, Seattle





-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden

[[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] Passing data among multiple instances

2009-02-04 Thread Feng Li
Dear R,

I have two R instances running at the same time, say instance A and instance
B. Is there a simpler way to pass the data in A to B?

More precise, I have a stupid example:

In instance A, I am running a function test1

test1 - function (x1)
{
 x2 - sin(x1)
 return(x2)
}

In instance B, another function test2

test2 - function (x2)
{
 x3 - cos(x2)
 return(x3)
}

where  test2 receives the input from test1's rueslt.  test1 and
test2 could be much more complex. They may take one minute each.

Now the whole procedure is instance A is running, while instance B is
waiting for the result of instance A. Once instance A is done (instance A
goes to run with new data), instance B should detect A is done, and instance
B receives the parameter from instance A. B begins to work. While B is done,
waiting for A's new results.

I want to repeat the loop many times and get x3 in the end.


Is it possible to do this job? Thanks !


Feng

[[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.


Re: [R] Passing data among multiple instances

2009-02-04 Thread Feng Li
Very very fruitful.

Now I only do the experiment on my single computer with a Quad CPU and more
than 2 G ram.

Let me have a try first.

Feng

On Wed, Feb 4, 2009 at 4:02 PM, Warren Young war...@etr-usa.com wrote:

 Feng Li wrote:


 I have two R instances running at the same time,


 On the same computer, or on different computers?

 Is the number of Rs likely to change, or will it always be just the two?

 Is this a simple one-off problem, or are you breaking the problem up into
 pieces so you can throw lots of hardware at it?

  Is there a simpler way to pass the data in A to B?


 Perhaps the simplest option is to write the data structure to a file, using
 any of the several R ways to do that.  When instance 2 sees that a file is
 available, it slurps its contents in and works on it.  The hard part is
 making the second instance wait until the whole file is written out by the
 first.  You wouldn't want it to read in half the file then hit the end
 because the first process hasn't finished writing out the file.  I don't see
 any good mechanism in R to fix this.

 A more robust option is to use sockets.  This is suitable even within a
 single machine.  See ?make.socket.  This solves the how do I know when I've
 got the full data structure problem because the second process can just
 keep reading until it gets an error indicating that the remote peer closed
 the connection.  Once you have the data structure in string form, you can
 eval() it to get an R object suitable for munching on.  Figuring out how to
 pass the data might be the hardest part.  deparse() might be the easiest
 way.

 If you're hoping to scale this up to lots of processes, look into Rmpi.
  This provides a very clean way for an R program on one computer to start
 slaves on other computers and then pass data to them in native R structures.
  Setting up MPI itself is not trivial, however.  It's best when you already
 have a cluster of computers linked with MPI.

 __
 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.




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden

[[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.


Re: [R] Passing data among multiple instances

2009-02-04 Thread Feng Li
Thanks your information any how.

I am using a Quad CPU. That's why I really want two or more than two
instances.


On Wed, Feb 4, 2009 at 3:46 PM, markle...@verizon.net wrote:

 There's no need to wait as long as you have the econd function fter the
 first ( the second won't
 start until the first one's finished )  but I think it's better to use
 lexical scoping by putting the second function right inside the first.
 Something like this:

 test1 - function (x1)
  {
  x2 - sin(x1)

  test2 - function(x2) {
   x3- cos(x2)
  }

   return(x3)
 }


 Then just call test1.

 But wait till someone else responds because I'm not an expeRt and someone
 else might
 say something more useful or different.








  In instance B, another function test2

 test2 - function (x2)
 {
 x3 - cos(x2)
 return(x3)
 }




 On Wed, Feb 4, 2009 at  9:08 AM, Feng Li wrote:

  Dear R,

 I have two R instances running at the same time, say instance A and
 instance
 B. Is there a simpler way to pass the data in A to B?

 More precise, I have a stupid example:

 In instance A, I am running a function test1

 test1 - function (x1)
 {
 x2 - sin(x1)
 return(x2)
 }

 In instance B, another function test2

 test2 - function (x2)
 {
 x3 - cos(x2)
 return(x3)
 }

 where  test2 receives the input from test1's rueslt.  test1 and
 test2 could be much more complex. They may take one minute each.

 Now the whole procedure is instance A is running, while instance B is
 waiting for the result of instance A. Once instance A is done (instance A
 goes to run with new data), instance B should detect A is done, and
 instance
 B receives the parameter from instance A. B begins to work. While B is
 done,
 waiting for A's new results.

 I want to repeat the loop many times and get x3 in the end.


 Is it possible to do this job? Thanks !


 Feng

[[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.




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden

[[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.


Re: [R] Passing data among multiple instances

2009-02-04 Thread Feng Li
On Wed, Feb 4, 2009 at 4:02 PM, Warren Young war...@etr-usa.com wrote:

 Feng Li wrote:


 I have two R instances running at the same time,


 On the same computer, or on different computers?

The first trial is only on my single computer with Quad CPU and more than 2G
mem.




 Is the number of Rs likely to change, or will it always be just the two?


I am planning to do three tasks at same time. That will be instance A,
instance B and instance C. There tasks are more or less the same. But one
always depends on others results.



 Is this a simple one-off problem, or are you breaking the problem up into
 pieces so you can throw lots of hardware at it?

This is just for one project. But if this is available, later I will try
more on this!




  Is there a simpler way to pass the data in A to B?


 Perhaps the simplest option is to write the data structure to a file, using
 any of the several R ways to do that.  When instance 2 sees that a file is
 available, it slurps its contents in and works on it.  The hard part is
 making the second instance wait until the whole file is written out by the
 first.  You wouldn't want it to read in half the file then hit the end
 because the first process hasn't finished writing out the file.  I don't see
 any good mechanism in R to fix this.

 A more robust option is to use sockets.  This is suitable even within a
 single machine.  See ?make.socket.  This solves the how do I know when I've
 got the full data structure problem because the second process can just
 keep reading until it gets an error indicating that the remote peer closed
 the connection.  Once you have the data structure in string form, you can
 eval() it to get an R object suitable for munching on.  Figuring out how to
 pass the data might be the hardest part.  deparse() might be the easiest
 way.

 If you're hoping to scale this up to lots of processes, look into Rmpi.
  This provides a very clean way for an R program on one computer to start
 slaves on other computers and then pass data to them in native R structures.
  Setting up MPI itself is not trivial, however.  It's best when you already
 have a cluster of computers linked with MPI.

 __
 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.




-- 
Feng Li
Department of Statistics
Stockholm University
106 91 Stockholm, Sweden

[[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.