[R] factors and ops

2006-05-26 Thread Erin Hodgess
Dear R People

I have a variable which is an ID number that is a factor.

I would like to look for ID numberbs that are greeater than
a particular value.

However, factors are a big ugly for these operations.

I was messing with the HR data set from the SASmixed package.

HR$Patient is a factor like that.

Now if I used:
subset(HR,as.integer(as.character(Patient))  214)

that will work.

If seems to me that there may be a better way.
Is that true?

R Version 2.3.0 Windows

Thanks in advance!


sincerely,
Erin Hodgess
Associate Professor
Department of Computer and Mathematical Sciences
University of Houston - Downtown
mailto: [EMAIL PROTECTED]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Is there a way to draw 3d plot?

2006-05-26 Thread Michael
Hi all,

I have a 2D matrix, which has 100 rows, and 100 columns,

I have a 2D matrix, with 100 rows and 100 columns,

I want to display it using 3D plot, much like plot3d and mesh/surf functions
in matlab.

Specifically, in matlab, I just need to do the following:


[X, Y]=meshgrid([0:0.01:0.99, 0:0.01:0.99]);
% Z is my 2D matrix,
surf(X, Y, Z);


Note that X and Y are created so that I can associate physical meaning onto
the x and y axis of the 3D plot.

For example, my 100 rows represent 0, 0.01, 0.02, ... 0.99 here.

In Matlab I can also drag in the graphic window and see from different
visual angle and perspective of the 3D plot...

Are there similar functions in R that (1) show 3D plot; (2) let me
manipulate view angles easily?

Thanks a lot!

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] factors and ops

2006-05-26 Thread Petr Pikal
Hi


On 26 May 2006 at 1:33, Erin Hodgess wrote:

Date sent:  Fri, 26 May 2006 01:33:34 -0500
From:   Erin Hodgess [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Subject:[R]  factors and ops

 Dear R People
 
 I have a variable which is an ID number that is a factor.
 
 I would like to look for ID numberbs that are greeater than
 a particular value.
 
 However, factors are a big ugly for these operations.
 
 I was messing with the HR data set from the SASmixed package.
 
 HR$Patient is a factor like that.
 
 Now if I used:
 subset(HR,as.integer(as.character(Patient))  214)
 
 that will work.

AFAIK no other way. If you know you do not want HR$Patient as factor 
you can use 
HR$Patient - as.integer(as.character(HR$Patient))
or
 HR$Patient - as.numeric(as.character(HR$Patient)) 
and then
subset(HR,Patient  214)

but beware of attaching data.frame as if it is attached you need to 
detach it first and than to change factor to numeric and only then to 
attach it again.

HTH
Petr


 
 If seems to me that there may be a better way.
 Is that true?
 
 R Version 2.3.0 Windows
 
 Thanks in advance!
 
 
 sincerely,
 Erin Hodgess
 Associate Professor
 Department of Computer and Mathematical Sciences
 University of Houston - Downtown
 mailto: [EMAIL PROTECTED]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html

Petr Pikal
[EMAIL PROTECTED]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Distribution Fitting

2006-05-26 Thread Lorenzo Isella
Hi,
I know this is a bit off-topic, but I am quite puzzled. I am going
through several papers about aerosol physics and in this field you
often have determine the parameters of a distribution to match your
experimental data (one typically uses a Gaussian mixture).
However, in many cases people plot a normalized empirical distribution
function and then perform some least-square fitting rather than using
likelihood functions.
As an undergrad, I was told that the former approach is correct only
if you have a model for the dynamics (e.g. Ohm law and you perform a
least-square fitting), but not if you deal with a distribution and you
pick random draws from it (in that case, one should maximize the
probability of drawing the data which were actually observed and this
leads to likelihood functions).
The two approaches do not seem equivalent to me, but I cannot believe
that this distinction is ignored in practice...
Many thanks

Lorenzo

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Is there a way to draw 3d plot?

2006-05-26 Thread Uwe Ligges
Michael wrote:

 Hi all,
 
 I have a 2D matrix, which has 100 rows, and 100 columns,
 
 I have a 2D matrix, with 100 rows and 100 columns,
 
 I want to display it using 3D plot, much like plot3d and mesh/surf functions
 in matlab.
 
 Specifically, in matlab, I just need to do the following:
 
 
 [X, Y]=meshgrid([0:0.01:0.99, 0:0.01:0.99]);
 % Z is my 2D matrix,
 surf(X, Y, Z);
 
 
 Note that X and Y are created so that I can associate physical meaning onto
 the x and y axis of the 3D plot.
 
 For example, my 100 rows represent 0, 0.01, 0.02, ... 0.99 here.
 
 In Matlab I can also drag in the graphic window and see from different
 visual angle and perspective of the 3D plot...
 
 Are there similar functions in R that (1) show 3D plot; (2) let me
 manipulate view angles easily?

(1) See ?persp

(1) *and* (2): See package rgl.

Uwe Ligges



 Thanks a lot!
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Is there a way to draw 3d plot?

2006-05-26 Thread Michael
On 5/26/06, Uwe Ligges [EMAIL PROTECTED] wrote:

 Michael wrote:

  Hi all,
 
  I have a 2D matrix, which has 100 rows, and 100 columns,
 
  I have a 2D matrix, with 100 rows and 100 columns,
 
  I want to display it using 3D plot, much like plot3d and mesh/surf
 functions
  in matlab.
 
  Specifically, in matlab, I just need to do the following:
 
  
  [X, Y]=meshgrid([0:0.01:0.99, 0:0.01:0.99]);
  % Z is my 2D matrix,
  surf(X, Y, Z);
  
 
  Note that X and Y are created so that I can associate physical meaning
 onto
  the x and y axis of the 3D plot.
 
  For example, my 100 rows represent 0, 0.01, 0.02, ... 0.99 here.
 
  In Matlab I can also drag in the graphic window and see from different
  visual angle and perspective of the 3D plot...
 
  Are there similar functions in R that (1) show 3D plot; (2) let me
  manipulate view angles easily?

 (1) See ?persp

 (1) *and* (2): See package rgl.

 Uwe Ligges



 Thanks a lot,

But a glance at rgl seems requireing shape, etc... and very
complicated...

Any easier approaches?

persp does not allow me to use mouse to rotate

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] missed ylim from plot.default

2006-05-26 Thread Antonio, Fabio Di Narzo
2006/5/26, Peter Ehlers [EMAIL PROTECTED]:


 Antonio, Fabio Di Narzo wrote:

  Hi all.
  On my R-2.3.0,  calling:
 
 
 plot(1:5, rep(1,5), xlim=c(-0.5,7), ylim=c(-0.5,2.8), asp=1)
 
 
  gives to me a plot (tried pdf and X11) whose y axis goes from about -2
 to
  about 4.5. What I've missed? How can I show some pre-specified x and y
  ranges from a plot (keeping fixed the aspect ratio)?
 
  Tnx all,
  Antonio, Fabio Di Narzo.
 
[[alternative HTML version deleted]]

 Do you want a wide, but not very tall plot? If so, you
 could specify the width/height in your X11() call.

 Peter Ehlers


Tnx for your indication, now I see: x and y plotting ranges are influenced
by the actual window size (at least when fixing aspect ratio). I think this
is obvious... What I missed/forgot above, was that default plot width and
height are fixed indipendently from the 'x-ylim' and 'asp' arguments.

However now I have a problem: I can write a wrapper to set device width and
height depending on the (xlim,ylim,asp) triple. But there's often
extra-space in the device due to plot title, axis labels, etc. So, I fear I
can't simply set width=height*asp. Some suggestion?
(Anyway, by now I only have to produce few plots, so today I can go try by
try :-) ).

Bests,
Antonio.

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Indexing vector with repeated values

2006-05-26 Thread Andris Jankevics
Hi all,
I have a vector which contains many repeated values.

Z - c(1,2,3,4,5,1,2,3,4,5,1,2,3,5,5,2,3,4,5)

I need to know how many times each number in this vecetor is repeated.
I can use a command like this:

 table (cut(Z,seq(1,5,1)))

(1,2] (2,3] (3,4] (4,5] 
4 4 3 5 

But how can I find a break points for vector with random values and random 
number sequence intervals? Why result for interval (1,2) is 4, if there is 
only 3 values of 1?

Can I get from vector Z a smal vector Zs 1,2,3,4,5 ?


Thanks a lot!

Andris Jankevics

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Labels in cluster plots

2006-05-26 Thread Sasha Pustota
I'm trying to follow an example from ?clara, which plots two clusters:

  x - rbind(cbind(rnorm(10,0,8), rnorm(10,0,8)),
  cbind(rnorm(20,50,8), rnorm(20,50,8)))
 clarax - clara(x, 2)
 clusplot(clarax)

Suppose I'd like to label these 30 observations on the plot according to a
preconceived classification, e.g. c(rep(A,15),rep(B,7),rep(C,8))

Is there a way to plot observations with these pre-specified labels?
That is, I'd like the first 15 of 'x' to be labeled as 'A', the next
seven as 'B', the next eight as 'C'.

(P.S. For principal component plots it was helpfully suggested to use
the following as an argument to the plot:
pch=c('A','B','C')[factor(c(rep(A,15),rep(B,7),rep(C,8)))]. I
can't figure out if there is a similar way here, and what is a general
way to put repeated row labels on various graphs...)

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Is there a way to draw 3d plot?

2006-05-26 Thread Antonio, Fabio Di Narzo
2006/5/26, Michael [EMAIL PROTECTED]:

 On 5/26/06, Uwe Ligges [EMAIL PROTECTED] wrote:
 
  Michael wrote:
 
   Hi all,
  
   I have a 2D matrix, which has 100 rows, and 100 columns,
  
   I have a 2D matrix, with 100 rows and 100 columns,
  
   I want to display it using 3D plot, much like plot3d and mesh/surf
  functions
   in matlab.
  
   Specifically, in matlab, I just need to do the following:
  
   
   [X, Y]=meshgrid([0:0.01:0.99, 0:0.01:0.99]);
   % Z is my 2D matrix,
   surf(X, Y, Z);
   
  
   Note that X and Y are created so that I can associate physical meaning
  onto
   the x and y axis of the 3D plot.
  
   For example, my 100 rows represent 0, 0.01, 0.02, ... 0.99 here.
  
   In Matlab I can also drag in the graphic window and see from different
   visual angle and perspective of the 3D plot...
  
   Are there similar functions in R that (1) show 3D plot; (2) let me
   manipulate view angles easily?
 
  (1) See ?persp
 
  (1) *and* (2): See package rgl.
 
  Uwe Ligges
 
 
 
  Thanks a lot,

 But a glance at rgl seems requireing shape, etc... and very
 complicated...

 Any easier approaches?

 persp does not allow me to use mouse to rotate

 [[alternative HTML version deleted]]

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html



Have you seen the example in 'rgl' man page?
It seems sufficient something like:
rgl.surface(x,y,z)

In that example, is also showed how to colorize the surface.

Antonio.

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Indexing vector with repeated values

2006-05-26 Thread Uwe Ligges
Andris Jankevics wrote:

 Hi all,
 I have a vector which contains many repeated values.
 
 
Z - c(1,2,3,4,5,1,2,3,4,5,1,2,3,5,5,2,3,4,5)
 
 
 I need to know how many times each number in this vecetor is repeated.
 I can use a command like this:
 
 
table (cut(Z,seq(1,5,1)))


Why not simply
  table(Z)
?


 
 (1,2] (2,3] (3,4] (4,5] 
 4 4 3 5 
 
 But how can I find a break points for vector with random values and random 
 number sequence intervals? Why result for interval (1,2) is 4, if there is 
 only 3 values of 1?


There is no 1 counted, but 4 2s, you have not specified an interval 
with 1 in it.
You were looking for
  table(cut(Z,seq(0,5,1)))
but don't need cut() here at all.



 Can I get from vector Z a smal vector Zs 1,2,3,4,5 ?

seq(min(Z), max(Z))

Uwe Ligges


 
 Thanks a lot!
 
 Andris Jankevics
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] changing font size in plot(effect())

2006-05-26 Thread John Fox
Dear Peter,

Won't that wipe out the other components of axis.text, etc.?

Regards,
 John

On Thu, 25 May 2006 06:22:22 -0600
 P Ehlers [EMAIL PROTECTED] wrote:
 
 John Fox wrote:
 
  Dear Emilie,
  
  This is, I guess, the effect() function in the effects package. If
 so,
  note that the plot.effect() method uses trellis graphics (via the
  lattice package), not standard R graphics, so you have to control
  aspects of the plot in a manner consistent with trellis graphics.
  
  Unfortunately, you can't just specify the scales argument when you
 call
  plot(), since plot.effect() already includes a scales argument when
 it
  calls xyplot(). You can, however, set trellis parameters globally,
  e.g., via something like
  
  axis.text - trellis.par.get(axis.text)
  par.ylab.text - trellis.par.get(par.ylab.text)
  par.xlab.text - trellis.par.get(par.xlab.text)
  axis.text$cex - 1.5
  par.ylab.text$cex - 1.5
  par.xlab.text$cex - 1.5
  trellis.par.set(axis.text, axis.text)
  trellis.par.set(par.ylab.text, par.ylab.text)
  trellis.par.set(par.xlab.text, par.xlab.text)
  
  Perhaps this can be done more efficiently -- I'm far from a trellis
  whiz -- but the above should work.
  
  More generally, in designing the plot method for effect objects, it
  wasn't my goal to produce publication-quality plots; for that, I
  suggest that you build a custom graph using the information
 included in
  the effect object.
  
  I hope this helps,
   John
 
 Possibly a bit more efficient:
 
 trellis.par.set(list(axis.text = list(cex = 2),
   par.ylab.text = list(cex = 1.5),
   par.xlab.text = list(cex = 0.5)))
 
 Also good to know, Emilie:
 
trellis.par.get()
 
 to see all the things that can be (re)set. Most are self-explanatory.
 
   - Peter Ehlers
 
  
  On Wed, 24 May 2006 13:24:03 -0400
   Emilie Berthiaume [EMAIL PROTECTED] wrote:
  
 
 I can't seem to be able to change the font size in an effect
 display.
  I've
 tried the following:
 
 
 par(cex.lab=4)
 plot(effect (alti,reg8), ylab=detection probability)
 
 and
 
 
 plot(effect (alti,reg8), ylab=detection probability, cex=4)
 
 but nothing changes.  Can anyone help me?
 thanks.
 
 



John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] factors and ops

2006-05-26 Thread Gabor Grothendieck
Try an ordered factor:

  subset(HR, ordered(Patient)  214)

or

  HR$Patient - ordered(HR$Patient)
  subset(HR, Patient  214)

You might also check that it orders them in the way you want.
ordered(HR$Patient) in the first case or just HR$Patient in the
second case will display the data followed by the ordering.

On 5/26/06, Erin Hodgess [EMAIL PROTECTED] wrote:
 Dear R People

 I have a variable which is an ID number that is a factor.

 I would like to look for ID numberbs that are greeater than
 a particular value.

 However, factors are a big ugly for these operations.

 I was messing with the HR data set from the SASmixed package.

 HR$Patient is a factor like that.

 Now if I used:
 subset(HR,as.integer(as.character(Patient))  214)

 that will work.

 If seems to me that there may be a better way.
 Is that true?

 R Version 2.3.0 Windows

 Thanks in advance!


 sincerely,
 Erin Hodgess
 Associate Professor
 Department of Computer and Mathematical Sciences
 University of Houston - Downtown
 mailto: [EMAIL PROTECTED]

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Vector elements and ratios

2006-05-26 Thread Andrej Kastrin
Dear useRs,

I have two different length vectors: one column (1...m) and one row 
vector (1...n):

20
40
20
60

5 4 2

Now I have to calculate ratios between column vector elements and each 
row vector elements:

4 5 10
8 10 20
4 5 20
15 12 30

Thank's in advance for any suggestions,

Andrej

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Vector elements and ratios

2006-05-26 Thread Chuck Cleland
Andrej Kastrin wrote:
 Dear useRs,
 
 I have two different length vectors: one column (1...m) and one row 
 vector (1...n):
 
 20
 40
 20
 60
 
 5 4 2
 
 Now I have to calculate ratios between column vector elements and each 
 row vector elements:
 
 4 5 10
 8 10 20
 4 5 20
 15 12 30

A - c(20,40,20,60)
B - c(5,4,2)

A %*% t(1/B)
  [,1] [,2] [,3]
[1,]45   10
[2,]8   10   20
[3,]45   10
[4,]   12   15   30

 Thank's in advance for any suggestions,
 
 Andrej
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Chuck Cleland, Ph.D.
NDRI, Inc.
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Vector elements and ratios

2006-05-26 Thread jim holtman
?outer
 v1 - c(20,40,20,60)
 v2 - c(5,4,2)
 outer(v1,v2,'/')
 [,1] [,2] [,3]
[1,]45   10
[2,]8   10   20
[3,]45   10
[4,]   12   15   30




On 5/26/06, Andrej Kastrin [EMAIL PROTECTED] wrote:

 Dear useRs,

 I have two different length vectors: one column (1...m) and one row
 vector (1...n):

 20
 40
 20
 60

 5 4 2

 Now I have to calculate ratios between column vector elements and each
 row vector elements:

 4 5 10
 8 10 20
 4 5 20
 15 12 30

 Thank's in advance for any suggestions,

 Andrej

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html




-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390 (Cell)
+1 513 247 0281 (Home)

What is the problem you are trying to solve?

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Vector elements and ratios

2006-05-26 Thread David Hugh-Jones
 outer(c(20,40,20,60), c(5,4,2), /)
 [,1] [,2] [,3]
[1,]45   10
[2,]8   10   20
[3,]45   10
[4,]   12   15   30


cheers
D

On 26/05/06, Andrej Kastrin [EMAIL PROTECTED] wrote:
 Dear useRs,

 I have two different length vectors: one column (1...m) and one row
 vector (1...n):

 20
 40
 20
 60

 5 4 2

 Now I have to calculate ratios between column vector elements and each
 row vector elements:

 4 5 10
 8 10 20
 4 5 20
 15 12 30

 Thank's in advance for any suggestions,

 Andrej

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Simulate Strauss process in 3D

2006-05-26 Thread Adrian Baddeley
Apple Ho writes:

 I am searching some information about simulating Strauss process in 3D 

If you need this immediately, I suggest you follow Brian Ripley's advice.
This probably involves downloading a source tar file of the R package
and finding the Fortran code in library/spatial/src/. The Fortran should be
edited to generate 3D patterns (it's not hard: everywhere you see 'x' and 'y'
just include 'z'). Then also edit the R interface to the Fortran.

If you can wait about 3 months, we will be releasing an upgrade of 
the R package `spatstat' that will deal with 3D point patterns, 
including simulation of the Strauss process in 3D.

regards
Adrian Baddeley

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Is there a way to draw 3d plot?

2006-05-26 Thread Duncan Murdoch
On 5/26/2006 3:39 AM, Antonio, Fabio Di Narzo wrote:
 2006/5/26, Michael [EMAIL PROTECTED]:
 On 5/26/06, Uwe Ligges [EMAIL PROTECTED] wrote:
 Michael wrote:

 Hi all,

 I have a 2D matrix, which has 100 rows, and 100 columns,

 I have a 2D matrix, with 100 rows and 100 columns,

 I want to display it using 3D plot, much like plot3d and mesh/surf
 functions
 in matlab.

 Specifically, in matlab, I just need to do the following:

 
 [X, Y]=meshgrid([0:0.01:0.99, 0:0.01:0.99]);
 % Z is my 2D matrix,
 surf(X, Y, Z);
 

 Note that X and Y are created so that I can associate physical meaning
 onto
 the x and y axis of the 3D plot.

 For example, my 100 rows represent 0, 0.01, 0.02, ... 0.99 here.

 In Matlab I can also drag in the graphic window and see from different
 visual angle and perspective of the 3D plot...

 Are there similar functions in R that (1) show 3D plot; (2) let me
 manipulate view angles easily?
 (1) See ?persp

 (1) *and* (2): See package rgl.

 Uwe Ligges



 Thanks a lot,
 But a glance at rgl seems requireing shape, etc... and very
 complicated...

 Any easier approaches?

 persp does not allow me to use mouse to rotate

 [[alternative HTML version deleted]]

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html

 
 
 Have you seen the example in 'rgl' man page?
 It seems sufficient something like:
 rgl.surface(x,y,z)
 
 In that example, is also showed how to colorize the surface.

And an in-development update to rgl adds a number of new routines with 
interfaces similar to ones from the graphics package.  For example, 
persp3d is probably what Michael wants.

If you like playing with experimental versions and don't mind having 
defaults changing with new releases, etc., you can download a recent 
build from

http://www.stats.uwo.ca/faculty/murdoch/software

Bug reports and suggestions about the new stuff would be welcome.

Duncan

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Thanks for help

2006-05-26 Thread Thomas Wutzler
Hello David and Spencer,

thank you for your replies!
The systemfit-package was exactly what I was looking for.
The only little drawback is, that the nonlinear systemfit does not 
support predifined weights yet (in oder to ensure homegenity of 
variances when predictios are hetescedastic)

Thanks also to Katharina, who mailed me offline and drew my attention to 
the nls and the nlm function. This answerd my second question.

Thomas

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] lme, best model without convergence

2006-05-26 Thread Thomas Wutzler
Dear R-help list readers,

I am fitting mixed models with the lme function of the nlme package.
If I get convergence depends on how the method (ML/REM) and which (and 
how much) parameters will depend randomly on the cluster-variable.

How get the bist fit without convergence?




I set the parameters msVerbose and returnObject to TRUE:

lmeControl(maxIter=5, msMaxIter=200, tolerance=1e-4, niter=50, 
msTol=1e-5, nlmStepMax=500, 
,msVerbose=TRUE
,returnObject=TRUE
)

However, the lme-functions does not produce verbose output, nor does it 
return the best fit if lme is not converging.
It returns only an error:
Error in lme.formula(y ~ lndbh + I(lndbh^2) + lnh + I(lnh^2), random = 
~lndbh +  :
 iteration limit reached without convergence (9)



Best regards
Thomas

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Computing a reliability index of a statistic with missing data

2006-05-26 Thread Chaouch, Aziz
 Thanks Spencer, that is interesting but I must say I'm a bit lost with
the terminology. I'll try to catch up but I'm not sure I need a
complicated model (MC sounds complicated to me but it may not be...). I
plan to use this reliability index just as an indication and I need to
compute it in batch for several different charts so I try to keep the
statistic as simple as possible but yet efficient.

Aziz

-Original Message-
From: Spencer Graves [mailto:[EMAIL PROTECTED] 
Sent: May 25, 2006 8:12 PM
To: Chaouch, Aziz
Cc: R-help@stat.math.ethz.ch
Subject: Re: [R] Computing a reliability index of a statistic with
missing data

  Have you considered some kind of binary time series model? 
'RSiteSearch(binary time series)' produced 150 hits.  One of the first
20 mentioned continuous-time hidden Markov chains 
(http://finzi.psych.upenn.edu/R/library/repeated/html/chidden.html).  I
don't know if this will help you or not, but it might be worth
examining.

  hope this helps.
  Spencer Graves

Chaouch, Aziz wrote:
 Hi All,
 
 I'd like to compute a kind of reliability index (RI) that would in a 
 sense stand as a measure of reliability of a statistic (histogram etc)

 computed on a time serie with missing values. The final goal is that:
 
 RI=1 for a perfect reliability
 RI=0 for a total unreliability (no data at all as an extreme case...)
 
 The percentage of missing data is one indication: the more missing 
 data, the less confidence we can have in the statistic. But the 
 distribution of missing data throughout the data serie is important as
well:
 independently of the number of missing data, if available data are 
 regularily spaced in time the RI should be higher than if available 
 data are irregulary spaced. As a measure of sampling regularity, I 
 thought about computing the time to next record and then take its 
 variance over the time interval on which the statistic is computed. 
 The variance of the time to next record would be a measure of sampling

 regularity so that the final RI could be of the form:
 
 RI=1 when n=0
 RI~1/n*var(T)
 
 with
 n=% of missing data
 T=time to next record (in hours)
 
 However I need to normalize var(T) to use it to compute the RI. Does

 someone have an idea on how to do this (or another proposal to compute

 the RI)?
 
 Thanks,
 
 Aziz
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! 
 http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Try to debbug R script

2006-05-26 Thread LAURENS, Guillaume

Hi,

I have a bug and I don't find it. Can u help me please.

In attachement :
- My script file : CPU_study.txt
- My data sample file : a.txt

When I use it, the error message is :
 source(H:/R/_workspace/CPU_study_1.R)
Erreur dans [-(`*tmp*`, i, k, value = numeric(0)) :   rien à remplacer
 (= error in ... nothing to replace)

I try a lot of thing but it doesn't work.
The probleme is line 39 :
val - (val + t[(((i*step)-step)+j),k])

Thanks

Guillaume LAURENS - trainee student
Simulation Industry - EYYSIDV
Airbus France


This e-mail is intended only for the above addressee. It may contain
privileged information. If you are not the addressee you must not copy,
distribute, disclose or use any of the information in it. If you have
received it in error please delete it and immediately notify the sender.
Security Notice: all e-mail, sent to or from this address, may be
accessed by someone other than the recipient, for system management and
security reasons. This access is controlled under Regulation of
Investigatory Powers Act 2000, Lawful Business Practises.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] query: sample size calculation

2006-05-26 Thread ManuelPerera-Chang

Dear all,

I am doing something wrong.

I am trying to apply a formula for  sample size calculation as in the book
Design and Analysis of Clinical Trials, from Chow et all.

There a suggested sample size strategy uses the formula

0.30/0.45=F(0.80,2n,2n)/F(0.025,2n,2n)

which gives n=96, as in the book.

(here F(alfa,k,n) is the upper (alfa)th quantile of an F distribution with
k,n degrees of freedom)

I have been trying to get n=96 using the following code in R.

val-rep(NA,200)
 for (n in 1:200) {
  val[n]- qf(0.8, 2*n, 2*n,lower.tail = TRUE, log.p = FALSE)/qf(0.025,
2*n, 2*n,lower.tail = TRUE, log.p = FALSE);
 }

but val doesnot get any close to 0.30/0.45.

Thank you,

Manuel

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Try to debbug R script

2006-05-26 Thread Duncan Murdoch
On 5/26/2006 8:45 AM, LAURENS, Guillaume wrote:
 Hi,
 
 I have a bug and I don't find it. Can u help me please.
 
 In attachement :
 - My script file : CPU_study.txt
 - My data sample file : a.txt
 
 When I use it, the error message is :
 source(H:/R/_workspace/CPU_study_1.R)
 Erreur dans [-(`*tmp*`, i, k, value = numeric(0)) :   rien à remplacer
  (= error in ... nothing to replace)
 
 I try a lot of thing but it doesn't work.
 The probleme is line 39 :
 val - (val + t[(((i*step)-step)+j),k])

The error says that your indexing returned a zero length vector.  I'd 
suggest printing out the values of (((i*step)-step)+j) and k just 
before this line (or using the debugger to examine them).  They probably 
aren't in the range you were expecting.

Duncan Murdoch

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] changing font size in plot(effect())

2006-05-26 Thread P Ehlers
No, I don't think so, John. I think what I proposed is
Deepayan's newer way of setting parameters. I used to do it the
other way, too (S-PLUS may still need it), but I find this new
way very handy, since it avoids the get/save/alter/set process.

Try this:

trellis.device()
trellis.par.set(list(axis.text = list(col = blue))) # reset colour
trellis.par.get(axis.text)  # inspect
trellis.par.set(list(axis.text = list(cex = 2)))  # reset cex
trellis.par.get(axis.text)  # inspect; colour is still 'blue'

Best wishes,
Peter

John Fox wrote:

 Dear Peter,
 
 Won't that wipe out the other components of axis.text, etc.?
 
 Regards,
  John
 
 On Thu, 25 May 2006 06:22:22 -0600
  P Ehlers [EMAIL PROTECTED] wrote:
 
John Fox wrote:


Dear Emilie,

This is, I guess, the effect() function in the effects package. If

so,

note that the plot.effect() method uses trellis graphics (via the
lattice package), not standard R graphics, so you have to control
aspects of the plot in a manner consistent with trellis graphics.

Unfortunately, you can't just specify the scales argument when you

call

plot(), since plot.effect() already includes a scales argument when

it

calls xyplot(). You can, however, set trellis parameters globally,
e.g., via something like

axis.text - trellis.par.get(axis.text)
par.ylab.text - trellis.par.get(par.ylab.text)
par.xlab.text - trellis.par.get(par.xlab.text)
axis.text$cex - 1.5
par.ylab.text$cex - 1.5
par.xlab.text$cex - 1.5
trellis.par.set(axis.text, axis.text)
trellis.par.set(par.ylab.text, par.ylab.text)
trellis.par.set(par.xlab.text, par.xlab.text)

Perhaps this can be done more efficiently -- I'm far from a trellis
whiz -- but the above should work.

More generally, in designing the plot method for effect objects, it
wasn't my goal to produce publication-quality plots; for that, I
suggest that you build a custom graph using the information

included in

the effect object.

I hope this helps,
 John

Possibly a bit more efficient:

trellis.par.set(list(axis.text = list(cex = 2),
  par.ylab.text = list(cex = 1.5),
  par.xlab.text = list(cex = 0.5)))

Also good to know, Emilie:

   trellis.par.get()

to see all the things that can be (re)set. Most are self-explanatory.

  - Peter Ehlers


On Wed, 24 May 2006 13:24:03 -0400
 Emilie Berthiaume [EMAIL PROTECTED] wrote:


I can't seem to be able to change the font size in an effect

display.

I've
tried the following:



par(cex.lab=4)
plot(effect (alti,reg8), ylab=detection probability)

and



plot(effect (alti,reg8), ylab=detection probability, cex=4)

but nothing changes.  Can anyone help me?
thanks.


 
 
 
 
 John Fox
 Department of Sociology
 McMaster University
 Hamilton, Ontario, Canada
 http://socserv.mcmaster.ca/jfox/
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Try to debbug R script

2006-05-26 Thread Cuvelier Etienne

LAURENS, Guillaume a écrit :
 Hi,
 
 I have a bug and I don't find it. Can u help me please.
 
 In attachement :
 - My script file : CPU_study.txt
 - My data sample file : a.txt
 
 When I use it, the error message is :
 source(H:/R/_workspace/CPU_study_1.R)
 Erreur dans [-(`*tmp*`, i, k, value = numeric(0)) :   rien à remplacer
  (= error in ... nothing to replace)
 
 I try a lot of thing but it doesn't work.
 The probleme is line 39 :
 val - (val + t[(((i*step)-step)+j),k])
 

Have you try the debug library ?

With the function mtrace() you can follow the work step by step and in 
the same time with the function bp() you can add break points...

Etienne

 Thanks
 
 Guillaume LAURENS - trainee student
 Simulation Industry - EYYSIDV
 Airbus France
 
 
 This e-mail is intended only for the above addressee. It may contain
 privileged information. If you are not the addressee you must not copy,
 distribute, disclose or use any of the information in it. If you have
 received it in error please delete it and immediately notify the sender.
 Security Notice: all e-mail, sent to or from this address, may be
 accessed by someone other than the recipient, for system management and
 security reasons. This access is controlled under Regulation of
 Investigatory Powers Act 2000, Lawful Business Practises.
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
 
 
 
 


--

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] query: sample size calculation

2006-05-26 Thread P Ehlers
[EMAIL PROTECTED] wrote:

 Dear all,
 
 I am doing something wrong.
 
 I am trying to apply a formula for  sample size calculation as in the book
 Design and Analysis of Clinical Trials, from Chow et all.
 
 There a suggested sample size strategy uses the formula
 
 0.30/0.45=F(0.80,2n,2n)/F(0.025,2n,2n)
 
 which gives n=96, as in the book.
 
 (here F(alfa,k,n) is the upper (alfa)th quantile of an F distribution with
 k,n degrees of freedom)
 
 I have been trying to get n=96 using the following code in R.
 
 val-rep(NA,200)
  for (n in 1:200) {
   val[n]- qf(0.8, 2*n, 2*n,lower.tail = TRUE, log.p = FALSE)/qf(0.025,
 2*n, 2*n,lower.tail = TRUE, log.p = FALSE);
  }

Don't you want 'lower.tail = FALSE'?

Peter Ehlers

 
 but val doesnot get any close to 0.30/0.45.
 
 Thank you,
 
 Manuel
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] changing font size in plot(effect())

2006-05-26 Thread John Fox
Dear Peter,

OK -- I see that this works as advertized.

Thanks,
 John

On Fri, 26 May 2006 07:12:52 -0600
 P Ehlers [EMAIL PROTECTED] wrote:
 No, I don't think so, John. I think what I proposed is
 Deepayan's newer way of setting parameters. I used to do it the
 other way, too (S-PLUS may still need it), but I find this new
 way very handy, since it avoids the get/save/alter/set process.
 
 Try this:
 
 trellis.device()
 trellis.par.set(list(axis.text = list(col = blue))) # reset colour
 trellis.par.get(axis.text)  # inspect
 trellis.par.set(list(axis.text = list(cex = 2)))  # reset cex
 trellis.par.get(axis.text)  # inspect; colour is still 'blue'
 
 Best wishes,
 Peter
 
 John Fox wrote:
 
  Dear Peter,
  
  Won't that wipe out the other components of axis.text, etc.?
  
  Regards,
   John
  
  On Thu, 25 May 2006 06:22:22 -0600
   P Ehlers [EMAIL PROTECTED] wrote:
  
 John Fox wrote:
 
 
 Dear Emilie,
 
 This is, I guess, the effect() function in the effects package. If
 
 so,
 
 note that the plot.effect() method uses trellis graphics (via the
 lattice package), not standard R graphics, so you have to control
 aspects of the plot in a manner consistent with trellis graphics.
 
 Unfortunately, you can't just specify the scales argument when you
 
 call
 
 plot(), since plot.effect() already includes a scales argument
 when
 
 it
 
 calls xyplot(). You can, however, set trellis parameters globally,
 e.g., via something like
 
 axis.text - trellis.par.get(axis.text)
 par.ylab.text - trellis.par.get(par.ylab.text)
 par.xlab.text - trellis.par.get(par.xlab.text)
 axis.text$cex - 1.5
 par.ylab.text$cex - 1.5
 par.xlab.text$cex - 1.5
 trellis.par.set(axis.text, axis.text)
 trellis.par.set(par.ylab.text, par.ylab.text)
 trellis.par.set(par.xlab.text, par.xlab.text)
 
 Perhaps this can be done more efficiently -- I'm far from a
 trellis
 whiz -- but the above should work.
 
 More generally, in designing the plot method for effect objects,
 it
 wasn't my goal to produce publication-quality plots; for that, I
 suggest that you build a custom graph using the information
 
 included in
 
 the effect object.
 
 I hope this helps,
  John
 
 Possibly a bit more efficient:
 
 trellis.par.set(list(axis.text = list(cex = 2),
   par.ylab.text = list(cex = 1.5),
   par.xlab.text = list(cex = 0.5)))
 
 Also good to know, Emilie:
 
trellis.par.get()
 
 to see all the things that can be (re)set. Most are
 self-explanatory.
 
   - Peter Ehlers
 
 
 On Wed, 24 May 2006 13:24:03 -0400
  Emilie Berthiaume [EMAIL PROTECTED] wrote:
 
 
 I can't seem to be able to change the font size in an effect
 
 display.
 
 I've
 tried the following:
 
 
 
 par(cex.lab=4)
 plot(effect (alti,reg8), ylab=detection probability)
 
 and
 
 
 
 plot(effect (alti,reg8), ylab=detection probability, cex=4)
 
 but nothing changes.  Can anyone help me?
 thanks.
 
 
  
  
  
  
  John Fox
  Department of Sociology
  McMaster University
  Hamilton, Ontario, Canada
  http://socserv.mcmaster.ca/jfox/
  
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] MART(tm) vs. gbm

2006-05-26 Thread manuel martin
Hello,
I tried two different implementations of the stochastic gradient 
boosting (Friedman 2002) : the MART(tm) with R tool 
(http://www-stat.stanford.edu/~jhf/R-MART.html) and the gbm R package. 
To me, both seemed fairly comparable, except maybe regarding the 
different loss criterion proposed and the fact that the gbm tool is 
sightly more convenient to use. However, it seemed that the MART with R 
tool systematically outperforms the gbm tool in terms of goodness of fit 
(whatever the way of choosing the best iteration for the gbm package).
I tried to find out if there were specific options that could have 
explained it but nothing came out.  See below for an example of how I 
compare both implementations. Did any one had the same experience, and 
can anyone give me hints about such performance differences or tell me 
if I am missing something obvious?

Thank you in advance,Manuel

Here are the arguments and options I used for comparison purposes, 
working on a 1600 records * 15 variables dataset :

# the MART with R tool
lx -
mart( as.matrix(x), y, c(1,1,1,1,1,1,1,1,1,1,1,1,1,2,2)
 niter=1000, tree.size=6, learn.rate=0.01,
 loss.cri=2 # gaussian
)

# for gbm
gbm1 - gbm(y ~ v1 + v2 + v3 + v4 + v5+ v6+ v7+ v8 + v9 + v10 + v11 + 
v12 + v13 + v14 + v15,
data=data, var.monotone=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
distribution=gaussian, n.trees=1000, shrinkage=0.01,
interaction.depth=6, bag.fraction = 0.5, train.fraction = 0.5,
n.minobsinnode = 10, cv.folds = 1, keep.data=TRUE)

# I then do predictions on the same dataset, and further perform 
goodness of fit comparisons
#...

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] I cannot load the package Rcmdr

2006-05-26 Thread zhang jian
I can load the package in other R file, but I can not do it in the file.
How to solve the question?

 local({pkg - select.list(sort(.packages(all.available = TRUE)))
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
Loading required package: tcltk
Loading required package: car
Error in parse(file, n, text, prompt) : syntax error in *
Error: .onAttach failed in 'attachNamespace'
Error: package/namespace load failed for 'Rcmdr'

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] MART(tm) vs. gbm

2006-05-26 Thread manuel.martin
Hello,
I have been using two different implementations of the stochastic 
gradient boosting (Friedman 2002) : MART(tm) with R and the gbm package. 
Both are fairly comparable except that the MART with R systematically 
strongly (depending on the dataset though) outperforms the gbm tool in 
terms of goodness of fit.
For instance, a

# gbm package
   gbm1 - gbm(Y~X2+X3+X4+X5+X6,
  data=data,
  var.monotone=c(0,0,0,0,0),   #  0: no monotone restrictions
  distribution=gaussian, # bernoulli, adaboost, gaussian,
   # poisson, and coxph available
  n.trees=3000,# number of trees
  shrinkage=0.005, # shrinkage or learning rate,
   # 0.001 to 0.1 usually work
  interaction.depth=6, # 1: additive model, 2: two-way 
interactions, etc.
  bag.fraction = 0.5,  # subsampling fraction, 0.5 is 
probably best
  train.fraction = 0.5,# fraction of data for training,
   # first train.fraction*N used for 
training
  n.minobsinnode = 10, # minimum total weight needed in 
each node
  cv.folds = 5,# do 5-fold cross-validation
  keep.data=TRUE,  # keep a copy of the dataset with 
the object
  verbose=TRUE)# print out progress

# MART with R
X - as.matrix(cbind(data$X2,as.numeric(data$X3), 
as.numeric(data$X4),as.numeric(data$X5),data$X6))
Y - data$Y
mart(X, Y, c(1,2,2,2,1) , niter=3000, tree.size=6,learn.rate=0.005, 
loss.cri=2 #gaussian too
)

leads to very different goodnesses of fit (I can provide the dataset if 
needed).

Did anyone already encountered this, is there an explanation, am I 
missing something obvious in the argument settings?
Thank you in advance,

Manuel

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Function as.Date leading to error implying that strptime requires 3 arguments

2006-05-26 Thread Rob Balshaw
Thank you to Prof Ripley, Gabor Grothendieck, Dirk Eddelbuettel and
others who emailed with suggestions and comments.

Dr Ripley's suggestion that I reinstall appears to have fixed the
problem, though I confess I have upgraded only to 2.3.0 rather than the
beta 2.3.1.  I use Windows XP Pro, Service Pack 2.0 -- my apologies for
not mentioning this.

While this instance was not serious, in itself, is it a symptom of
something more troubling?  I worry that some feature of my usage has led
to this problem.  However, if it were something I did regularly, then it
should recur.  At that time, I will certainly have to try to correct
whatever it happens to be.  

In the meantime, however, if anyone has any suggestions about what I
might have done to create this problem I'd be pleased to hear them.  I
do (or rather, did) have multiple versions of R installed, and I do make
extensive use of several libraries, including the Bioconductor set.  My
suspicion is that this would not occur if I were to remove old versions
when I upgrade.

Cheers,

Rob


-Original Message-
From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] 
Sent: Friday, May 19, 2006 11:39 PM
To: Rob Balshaw
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Function as.Date leading to error implying that
strptime requires 3 arguments

Your system (unstated) is corrupted.  R runs its examples as part of it
test suite, so this is not an error in R 2.2.1 (which is not current:
see the posting guide which asked you to update *before* posting).  One
possibility is that you have strptime from a much earlier version of R
in use somehow.

I suggest you remove all versions of R from your system and reinstall
one latest version (2.3.1 beta).

On Fri, 19 May 2006, Rob Balshaw wrote:


 I'm using R V 2.2.1.  When I try an example from the as.Date help 
 page, I get an error.

 x - c(1jan1960, 2jan1960, 31mar1960, 30jul1960) z - 
 as.Date(x, %d%b%Y)
 Error in strptime(x, format) : 2 arguments passed to 'strptime' which 
 requires 3


 Any suggestions would be appreciated.

 Thanks,

 Rob

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! 
 http://www.R-project.org/posting-guide.html


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

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Is there a way to draw 3d plot?

2006-05-26 Thread Greg Snow
There is also the rotate.persp function in the TeachingDemos package.
It creates a set of slider bars that you can move with the mouse to
change the different options to persp including the angles.  It is not
quite as convenient as clicking in the plot and dragging like rgl
allows, but it does have the advantage that when you have the plot you
like, you can see what the current settings are and recreate the same
plot using persp and those settings rather than having to remember how
you rotated the graph interactively (you can query the current rotation
in rgl as well and set it by command in future plots). 


-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111
 

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Michael
Sent: Friday, May 26, 2006 1:25 AM
To: Uwe Ligges
Cc: R-help@stat.math.ethz.ch
Subject: Re: [R] Is there a way to draw 3d plot?

On 5/26/06, Uwe Ligges [EMAIL PROTECTED] wrote:

 Michael wrote:

  Hi all,
 
  I have a 2D matrix, which has 100 rows, and 100 columns,
 
  I have a 2D matrix, with 100 rows and 100 columns,
 
  I want to display it using 3D plot, much like plot3d and mesh/surf
 functions
  in matlab.
 
  Specifically, in matlab, I just need to do the following:
 
  
  [X, Y]=meshgrid([0:0.01:0.99, 0:0.01:0.99]); % Z is my 2D matrix, 
  surf(X, Y, Z);
  
 
  Note that X and Y are created so that I can associate physical 
  meaning
 onto
  the x and y axis of the 3D plot.
 
  For example, my 100 rows represent 0, 0.01, 0.02, ... 0.99 here.
 
  In Matlab I can also drag in the graphic window and see from 
  different visual angle and perspective of the 3D plot...
 
  Are there similar functions in R that (1) show 3D plot; (2) let me 
  manipulate view angles easily?

 (1) See ?persp

 (1) *and* (2): See package rgl.

 Uwe Ligges



 Thanks a lot,

But a glance at rgl seems requireing shape, etc... and very
complicated...

Any easier approaches?

persp does not allow me to use mouse to rotate

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] how to multiply a constant to a matrix?

2006-05-26 Thread Tony Plate
I still can't see why this is a problem.  If a 1x1 matrix should be 
treated as a scalar, then it can just be wrapped in drop(), and the 
arithmetic will be computed correctly by R.

Are there any cases where this cannot be done?  More specifically, are 
there any matrix algebra expressions where, depending on the particular 
dimensions of the variables used, drop() must be used in some cases, and 
not in other cases?

A related but different behavior is the default dropping dimensions with 
extent equal to one by indexing operations.  This can be problematic 
because if one is not careful, incorrect results can be obtained for 
particular values used in the expression.

For example, consider the following, in which we are trying to compute 
the cross product of some columns of x with some rows of y.  If x has n 
rows and y has n columns, then the result should always be an nxn 
matrix.  However, if we are not careful with using drop=F in the 
indexing expressions, we can inadvertently end up with a 1x1 inner 
product matrix result for the case where we just use one column of x and 
one row of y.  The solution to this is to always use drop=F in indexing 
in situations where this can occur.

  x - matrix(1:9, ncol=3)
  y - matrix(-(1:9), ncol=3)
  i - 1:2
  x[,i] %*% y[i,]
  [,1] [,2] [,3]
[1,]   -9  -24  -39
[2,]  -12  -33  -54
[3,]  -15  -42  -69
  i - 1:3
  x[,i] %*% y[i,]
  [,1] [,2] [,3]
[1,]  -30  -66 -102
[2,]  -36  -81 -126
[3,]  -42  -96 -150
  # i has just one element -- the expression without drop=F
  # no longer computes an outer product
  i - 2
  x[,i] %*% y[i,]
  [,1]
[1,]  -81
  x[,i,drop=F] %*% y[i,,drop=F]
  [,1] [,2] [,3]
[1,]   -8  -20  -32
[2,]  -10  -25  -40
[3,]  -12  -30  -48
 

Cannot all cases in the situations you mention be handled in an 
analogous manner, by always wrapping appropriate quadratic expressions 
in drop(), or are there some cases where the result of the quadratic 
expression must be treated as a matrix, and other cases where the result 
of the quadratic expression must be treated as a scalar?

-- Tony Plate

Michael wrote:
 imagine when you have complicated matrix algebra computation using R,
 
 you cannot prevent some middle-terms become quadratic and absorbs into one
 scalar, right?
 
 if R cannot intelligently determine this, and you  have to manually add
 drop everywhere,
 
 do you think it is reasonable?
 
 On 5/23/06, Patrick Burns [EMAIL PROTECTED] wrote:
 
I think

drop(B/D) * solve(A)

would be a more transparent approach.

It isn't that R can not do what you want, it is that
it is saving you from shooting yourself in the foot
in your attempt.  What you are doing is not really
a matrix computation.


Patrick Burns
[EMAIL PROTECTED]
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of S Poetry and A Guide for the Unwilling S User)

Michael wrote:


This is very strange:

I want compute the following in R:

g = B/D * solve(A)

where B and D are  quadratics so they are just a scalar number, e.g.

B=t(a)

%*% F %*% a;

I want to multiply B/D to A^(-1),

but R just does not allow me to do that and it keeps complaining that
nonconformable array, etc.


I tried the following two tricks and they worked:

as.numeric(B/D) * solve(A)

diag(as.numeric(B/D), 5, 5) %*% solve (A)



But if R cannot intelligently do scalar and matrix multiplication, it is
really problemetic.

It basically cannot be used to do computations, since in complicated

matrix

algebras, you have to distinguish where is scalar, and scalars obtained

from

quadratics cannot be directly used to multiply another matrix, etc. It is
going to a huge mess...

Any thoughts?

  [[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!

http://www.R-project.org/posting-guide.html




 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] upgrade or recompile

2006-05-26 Thread Mike Wolfgang
Dear list,

My current R 2.2.0 in my linux desktop was installed by compiling sources,
now I want to upgrade it to 2.3.0, is it possible to upgrade it without
re-compiling sources? Thanks,

mike

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] multiple comparisons of time series data

2006-05-26 Thread Kyle Hall
I am interested in a statistical comparison of multiple (5) time series' 
generated from modeling software (Hydrologic Simulation Program Fortran). The 
model output simulates daily bacteria concentration in a stream. The multiple 
time series' are a result of varying our representation of the stream within 
the model.

Our main question is: Do the different methods used to represent a stream 
produce different results at a statistically significant level?

We want to compare each otput time series to determine if there is a 
difference before looking into the cause within the model.  In a previous 
study, the Kolmogorov-Smirnov k-sample test was used to compare multiple time 
series'.

I am unsure about the strength of the Kolmogorov-Smirnov test and I have set 
out to determine if there are any other tests to compare multiple time 
series'.

I know htat R has the ks.test but I am unsure how this test handles multiple 
comparisons.  Is there something similar to a pairwise.t.test with a 
bonferroni corection, only with time series data?

Does R currently (v 2.3.0) have a comparison test that takes into account the 
strong serial correlation of time series data?


Kyle Hall

Graduate Research Assistant
Biological Systems Engineering
Virginia Tech

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Maximum likelihood estimate of bivariate vonmises-weibulldistribution

2006-05-26 Thread Chaouch, Aziz
Hi,

I'm still strugling with this copula model but this seems to be the way
to go. I'm now trying to model the marginal distributions and and for
wind direction I use a mixture of 2 von mises. I'd like to estimate all
the parameters (m1,m1,kappa1,kappa2,p) by maximizing the likelihood but
I don't know how to define the likelihood (or log-likelihood) of a
mixture of 2 Von Mises to use it with the function fitDistr in the MASS
package. Can you help me define this likelihood function and use it
through the fitDistr function?

Thanks,

Aziz

-Original Message-
From: Dimitrios Rizopoulos [mailto:[EMAIL PROTECTED] 
Sent: May 12, 2006 4:35 PM
To: Chaouch, Aziz
Subject: RE: [R] Maximum likelihood estimate of bivariate
vonmises-weibulldistribution

look at the following code:

library(copula)
par(mfrow = c(2, 2))
x - mvdc(normalCopula(sin(0.5 * pi /2)), c(norm, norm),
list(list(mean = 0, sd = 1), list(mean = 0, sd = 1))) contour(x, dmvdc,
xlim = c(-2.7, 2.7), ylim = c(-2.7, 2.7))

x - mvdc(frankCopula(5.736276), c(norm, norm), list(list(mean = 0,
sd = 1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7),
ylim = c(-2.7, 2.7))

x - mvdc(gumbelCopula(2), c(norm, norm), list(list(mean = 0, sd =
1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7), ylim
= c(-2.7, 2.7))

x - mvdc(claytonCopula(2), c(norm, norm), list(list(mean = 0, sd =
1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7), ylim
= c(-2.7, 2.7))


the values of the association parameter I've chosen in each copula 
correspond to Kendall's tau 0.5; assuming also standard normal 
marginal distributions look at the different shapes you get!

If possible try something similar for you case (i.e., using von Mises 
and Weibull marginals) and check if the association shape for a 
specific copula is more appropriate for your application. If this is 
not possible fit models assumig different copulas and check which one 
provides a better fit to your data.

I hope it helps.

Best,
Dimitris



Quoting Chaouch, Aziz [EMAIL PROTECTED]:

 Hi Dimitris,
 
 I'm not sure to understand your suggestion. How would you build that
 contour plot for a particular copula and on what is computed the
 Kendall's tau? 
 
 Thanks,
 
 Aziz
 
 -Original Message-
 From: Dimitris Rizopoulos
 [mailto:[EMAIL PROTECTED] 
 Sent: May 12, 2006 9:57 AM
 To: Chaouch, Aziz; [EMAIL PROTECTED]
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] Maximum likelihood estimate of bivariate
 vonmises-weibulldistribution
 
 the choice of the copula is, in fact, a model selection problem. 
 First, you could have a look at the contour plots of different
 copulas
 (preferably for the same value of Kendall's tau), and decide if some
 of
 them assume a more appropriate association structure for your
 application, compared to the others. Afterwards, you may fit various
 copula functions, check the fit on the data, compute AIC (since
 these
 are typically not nested models), etc.
 
 regarding the Von Mises distribution and if could be used in mvdc(),
 that I don't know. It'd be better to contact the copula package
 maintainer and ask.
 
 I hope it helps.
 
 Best,
 Dimitirs
 
 
 Dimitris Rizopoulos
 Ph.D. Student
 Biostatistical Centre
 School of Public Health
 Catholic University of Leuven
 
 Address: Kapucijnenvoer 35, Leuven, Belgium
 Tel: +32/(0)16/336899
 Fax: +32/(0)16/337015
 Web: http://www.med.kuleuven.be/biostat/
  http://www.student.kuleuven.be/~m0390867/dimitris.htm
 
 
 - Original Message -
 From: Chaouch, Aziz [EMAIL PROTECTED]
 To: Dimitris Rizopoulos [EMAIL PROTECTED];
 [EMAIL PROTECTED]
 Cc: r-help@stat.math.ethz.ch
 Sent: Friday, May 12, 2006 3:13 PM
 Subject: RE: [R] Maximum likelihood estimate of bivariate
 vonmises-weibulldistribution
 
 
 Thanks a lot! I wasn't aware of that copula package and it could well
 be
 appropriate to use it for my application. However if I read the
 copula
 help correctly, I still need to know what kind of copula to use to
 link
 the distribution of wind speeds and directions. Is there a way to
 determine this in R?
 
 Moreover can I use the Von Mises distribution from the circular or
 CircStats package into the mvdc function of the copula package or
 does
 the mvdc function only recognize distributions available natively
 within R?
 
 Thanks again to all, your help is highly appreciated for a newbie
 like
 me!
 
 Regards,
 
 Aziz
 
 -Original Message-
 From: Dimitris Rizopoulos
 [mailto:[EMAIL PROTECTED]
 Sent: May 12, 2006 3:01 AM
 To: Philip He; Chaouch, Aziz
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] Maximum likelihood estimate of bivariate
 vonmises-weibulldistribution
 
 
 - Original Message -
 From: Philip He [EMAIL PROTECTED]
 To: Chaouch, Aziz [EMAIL PROTECTED]
 Cc: r-help@stat.math.ethz.ch
 Sent: Thursday, May 11, 2006 11:21 PM
 Subject: Re: [R] Maximum likelihood estimate of bivariate
 vonmises-weibulldistribution
 
 
  On 5/11/06, Chaouch, Aziz [EMAIL PROTECTED] 

Re: [R] upgrade or recompile

2006-05-26 Thread Uwe Ligges
Mike Wolfgang wrote:
 Dear list,
 
 My current R 2.2.0 in my linux desktop was installed by compiling sources,
 now I want to upgrade it to 2.3.0, is it possible to upgrade it without
 re-compiling sources? Thanks,

No, except if installing a binary distribution.

BTW: You might want to try an R-2.3.1 beta release instead anyway.

Uwe Ligges



 
 mike
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Maximum likelihood estimate of bivariatevonmises-weibulldistribution

2006-05-26 Thread Ravi Varadhan
Hi Aziz,

I am attaching a file that contains the functions that I wrote to compute
the MLE for a binary vonMises mixture, using the EM algorithm.  It also
contains a simulation example.  Let me know if you have any trouble.

Hope this is helpful,
Ravi.

--
Ravi Varadhan, Ph.D.
Assistant Professor,  The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email:  [EMAIL PROTECTED]
Webpage:http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html 
--
 -Original Message-
 From: [EMAIL PROTECTED] [mailto:r-help-
 [EMAIL PROTECTED] On Behalf Of Chaouch, Aziz
 Sent: Friday, May 26, 2006 12:02 PM
 To: Dimitrios Rizopoulos
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] Maximum likelihood estimate of bivariatevonmises-
 weibulldistribution
 
 Hi,
 
 I'm still strugling with this copula model but this seems to be the way
 to go. I'm now trying to model the marginal distributions and and for
 wind direction I use a mixture of 2 von mises. I'd like to estimate all
 the parameters (m1,m1,kappa1,kappa2,p) by maximizing the likelihood but
 I don't know how to define the likelihood (or log-likelihood) of a
 mixture of 2 Von Mises to use it with the function fitDistr in the MASS
 package. Can you help me define this likelihood function and use it
 through the fitDistr function?
 
 Thanks,
 
 Aziz
 
 -Original Message-
 From: Dimitrios Rizopoulos [mailto:[EMAIL PROTECTED]
 Sent: May 12, 2006 4:35 PM
 To: Chaouch, Aziz
 Subject: RE: [R] Maximum likelihood estimate of bivariate
 vonmises-weibulldistribution
 
 look at the following code:
 
 library(copula)
 par(mfrow = c(2, 2))
 x - mvdc(normalCopula(sin(0.5 * pi /2)), c(norm, norm),
 list(list(mean = 0, sd = 1), list(mean = 0, sd = 1))) contour(x, dmvdc,
 xlim = c(-2.7, 2.7), ylim = c(-2.7, 2.7))
 
 x - mvdc(frankCopula(5.736276), c(norm, norm), list(list(mean = 0,
 sd = 1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7),
 ylim = c(-2.7, 2.7))
 
 x - mvdc(gumbelCopula(2), c(norm, norm), list(list(mean = 0, sd =
 1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7), ylim
 = c(-2.7, 2.7))
 
 x - mvdc(claytonCopula(2), c(norm, norm), list(list(mean = 0, sd =
 1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7), ylim
 = c(-2.7, 2.7))
 
 
 the values of the association parameter I've chosen in each copula
 correspond to Kendall's tau 0.5; assuming also standard normal
 marginal distributions look at the different shapes you get!
 
 If possible try something similar for you case (i.e., using von Mises
 and Weibull marginals) and check if the association shape for a
 specific copula is more appropriate for your application. If this is
 not possible fit models assumig different copulas and check which one
 provides a better fit to your data.
 
 I hope it helps.
 
 Best,
 Dimitris
 
 
 
 Quoting Chaouch, Aziz [EMAIL PROTECTED]:
 
  Hi Dimitris,
 
  I'm not sure to understand your suggestion. How would you build that
  contour plot for a particular copula and on what is computed the
  Kendall's tau?
 
  Thanks,
 
  Aziz
 
  -Original Message-
  From: Dimitris Rizopoulos
  [mailto:[EMAIL PROTECTED]
  Sent: May 12, 2006 9:57 AM
  To: Chaouch, Aziz; [EMAIL PROTECTED]
  Cc: r-help@stat.math.ethz.ch
  Subject: Re: [R] Maximum likelihood estimate of bivariate
  vonmises-weibulldistribution
 
  the choice of the copula is, in fact, a model selection problem.
  First, you could have a look at the contour plots of different
  copulas
  (preferably for the same value of Kendall's tau), and decide if some
  of
  them assume a more appropriate association structure for your
  application, compared to the others. Afterwards, you may fit various
  copula functions, check the fit on the data, compute AIC (since
  these
  are typically not nested models), etc.
 
  regarding the Von Mises distribution and if could be used in mvdc(),
  that I don't know. It'd be better to contact the copula package
  maintainer and ask.
 
  I hope it helps.
 
  Best,
  Dimitirs
 
  
  Dimitris Rizopoulos
  Ph.D. Student
  Biostatistical Centre
  School of Public Health
  Catholic University of Leuven
 
  Address: Kapucijnenvoer 35, Leuven, Belgium
  Tel: +32/(0)16/336899
  Fax: +32/(0)16/337015
  Web: http://www.med.kuleuven.be/biostat/
   http://www.student.kuleuven.be/~m0390867/dimitris.htm
 
 
  - Original Message -
  From: Chaouch, Aziz [EMAIL PROTECTED]
  To: Dimitris Rizopoulos [EMAIL PROTECTED];
  [EMAIL PROTECTED]
  Cc: r-help@stat.math.ethz.ch
  Sent: Friday, May 12, 2006 3:13 PM
  Subject: RE: [R] Maximum likelihood estimate of bivariate
  vonmises-weibulldistribution
 
 
  Thanks a lot! I wasn't aware of that copula package and it could well
  be
  appropriate to use it for my application. 

[R] how to pick a value from AR equation

2006-05-26 Thread Lorenzo Bencivelli
i need to compute several (hundreds) of regression AR and PP.test for 
untary roots.
is there an easy way to pick the values of interest from the output of 
these operations?
i would like to fill a matrix (in a for cycle) with all the values of 
coefficients, standard deviations, PP statistics and the relevant 
critical value for the series i have.
thanks in advance L

-- 



credo nella ragione umana, 
e nella liberta' e nella giustizia che dalla ragione scaturiscono. (sciascia)

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] I cannot load the package Rcmdr

2006-05-26 Thread John Fox
Dear Zhang Jian,

Is it possible that you're using a version of the Rcmdr package that's
incompatible with your version of R?

(I won't have access to email this weekend, so will be unable to reply
until Monday to your response.)

I hope this helps,
 John

On Fri, 26 May 2006 09:52:31 -0400
 zhang jian [EMAIL PROTECTED] wrote:
 I can load the package in other R file, but I can not do it in the
 file.
 How to solve the question?
 
  local({pkg - select.list(sort(.packages(all.available = TRUE)))
 + if(nchar(pkg)) library(pkg, character.only=TRUE)})
 Loading required package: tcltk
 Loading required package: car
 Error in parse(file, n, text, prompt) : syntax error in *
 Error: .onAttach failed in 'attachNamespace'
 Error: package/namespace load failed for 'Rcmdr'
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html


John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] combinatorial programming problem

2006-05-26 Thread Kjetil Brinchmann Halvorsen
Hola!

I am programming a class (S3) symarray for
storing the results of functions symmetric in its
k arguments. Intended use is for association indices
for more than two variables, for instance coresistivity
against antibiotics.

There is one programming problem I haven't solved, making an inverse
of the index function indx() --- se code below. It could for instance 
return the original k indexes in strictly increasing order, to make 
indx() formally invertible.

Any ideas?

Kjetil


Code:


# Implementing an S3 class for symarrays with array rank r for dimension
# [k, k, ..., k] with k=r repeated r times. We do not store the diagonal.

# Storage requirement is given by {r, k}=  choose(k, r)
# where r=array rank, k=maximum index

symarray - function(data=NA, dims=c(1,1)){
  r - dims[1]
  k - dims[2]
  if(r  k) stop(symarray needs dimension larger than array rank)
  len - choose(k, r)
  out - data[1:len]
  attr(out, dims) - dims
  class(out) - symarray
  out
}

# Index calculation:

indx - function(inds, k){
r - length(inds)
if(r==1) return(inds) else {
   if(inds[1]==1) {
  return( indx(inds[-1]-1, k-1 ) ) } else {
  return( indx(c(inds[1]-1, seq(from=k-r+2, by=1, to=k)), k) +
  indx( inds[-1]-inds[1], k-inds[1] ))
  }
}
   } # end indx

# Methods for assignment and indexing:

[.symarray - function(x, inds, drop=FALSE){
 dims - attr(x, dims)
 k - dims[2]
 inds - indx(inds, k)
 res - NextMethod([, x)
 res
}

[-.symarray - function(x, inds, value){
 dims - attr(x, dims)
 k - dims[2]
 inds - indx(inds, k)
 res - NextMethod([-, x)
 res
}

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Returned mail: Data format error

2006-05-26 Thread hci
Hello,

Your mail to [EMAIL PROTECTED] was caught by the
SpamAssassin filter running on the bcs.org.uk mail system.

To confirm that your mail is genuine, please click this
link, or paste it into your browser:
https://bcsnet.bcs.org.uk/approve.php?c=66e3c8769c1611445312f4

You will not have to do this again for any mail sent
to this recipient ([EMAIL PROTECTED]).

Thank you.

-- 
British Computer Society - www.bcs.org.uk
Email Services from gradwell dot com - www.gradwell.com

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] how to pick a value from AR equation

2006-05-26 Thread Rogerio Porto
Lorenzo,

in order to help you more, you should consider sending to the
list some relevant code with the functions and packages you
are using.

Rogerio.
- Original Message - 
From: Lorenzo Bencivelli [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Sent: Friday, May 26, 2006 1:21 PM
Subject: [R] how to pick a value from AR equation


i need to compute several (hundreds) of regression AR and PP.test for
 untary roots.
 is there an easy way to pick the values of interest from the output of
 these operations?
 i would like to fill a matrix (in a for cycle) with all the values of
 coefficients, standard deviations, PP statistics and the relevant
 critical value for the series i have.
 thanks in advance L

 -- 


 
 credo nella ragione umana,
 e nella liberta' e nella giustizia che dalla ragione scaturiscono. 
 (sciascia)

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! 
 http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] how to pick a value from AR equation

2006-05-26 Thread Lorenzo Bencivelli
ok, i'll try to be more precise.
both functions i need (ar() and PP.test()) are in the package stats.
what i would like to do:

omega-array(0, dim=c(N+1,7))
omega[1,]-(series,mean,std,max,min,ar(1)coeff,stdar(1))
for (i in (2:N+1)){
omega[i,1]-i-1
omega[i,2]-mean(y[(i-1),])
(...)
omega[i,6]- (???)
omega[i,7]- (???)}

for the PP.test, i would like to do something similar. N is the number 
of the series for which i would like to compute the AR(1) model and the 
unitary root test. thanks in advance. L

-- 



credo nella ragione umana, 
e nella liberta' e nella giustizia che dalla ragione scaturiscono. (sciascia)

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] I cannot load the package Rcmdr

2006-05-26 Thread Prof Brian Ripley
On Fri, 26 May 2006, John Fox wrote:

 Dear Zhang Jian,

 Is it possible that you're using a version of the Rcmdr package that's
 incompatible with your version of R?

That ought to have generated an error much earlier, but of course the 
dependencies in an old version of Rcmdr might not have been prescient.
I would suggest running update.packages() to be sure.

 (I won't have access to email this weekend, so will be unable to reply
 until Monday to your response.)

 I hope this helps,
 John

 On Fri, 26 May 2006 09:52:31 -0400
 zhang jian [EMAIL PROTECTED] wrote:
 I can load the package in other R file, but I can not do it in the
 file.
 How to solve the question?

 local({pkg - select.list(sort(.packages(all.available = TRUE)))
 + if(nchar(pkg)) library(pkg, character.only=TRUE)})
 Loading required package: tcltk
 Loading required package: car
 Error in parse(file, n, text, prompt) : syntax error in *
 Error: .onAttach failed in 'attachNamespace'
 Error: package/namespace load failed for 'Rcmdr'

Please look at the section on debugging in the `Writing R Extensions' 
manual.  (If you do not see it, you need to update your R to = 2.3.0.) At 
a minimum, running traceback() at that point would have helped both you 
and us.  (I am guessing this occurs in the call to Commander(), but would 
like not to be guessing.)

Also, please check out the posting guide and tell us the information 
requested, e.g. sessionInfo(), and Sys.getlocale() in case it is a charset 
issue.

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

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] how to get a series of results by loop

2006-05-26 Thread zhang jian
I donot know how to do it remainly.
Perhaps my question is not clear.
For example,

test=data.frame(sp=paste(c(sp),1:100,sep=),x=rnorm(100),y=rnorm(100))

# I want to get the data in the circle of radius = 1 in every point, and I
want to add a sign (e.g. sp1,sp2) in order that I can distinguish
different groups
# if I want to choose a point, I can do it like this

r = 1 # radius
attach(test)
point100=test[(x-x[100])^2+(y-y[100])^2=r^2,]
result100=data.frame
(note=paste(c(p),rep(100,length(point100$sp)),sep=),point100)
detach(test)

But if I want to do all points in test, I donot know how to do it. I try
to do it by using for, but it doesnot work.
Can you give me some advices?
thanks!
  jian zhang

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] how to pick a value from AR equation

2006-05-26 Thread Marco Geraci
Ciao,
you didn't provide lots of details so I'll try to
answer your questions with examples you can find
typing ?PP.test and ?ar. In general, when you run a
command in R you get an output. You just need to
understand what kind of object you get. For the
PP.test

 x - rnorm(1000)
 pp - PP.test(x)

The command 'str' is very useful

 str(pp)
List of 5
 $ statistic: Named num -33.2
  ..- attr(*, names)= chr Dickey-Fuller
 $ parameter: Named num 7
  ..- attr(*, names)= chr Truncation lag parameter
 $ p.value  : num 0.01
 $ method   : chr Phillips-Perron Unit Root Test
 $ data.name: chr x
 - attr(*, class)= chr htest

It's a list. So you can select the element you need,
e.g.

 pp$statistic
Dickey-Fuller 
-33.22576 

or

 as.numeric(pp$statistic)
[1] -33.22576

so you get rid of the attributes.
Now, same thing for 'ar'.

 fit.ar - ar(lh)
 str(fit.ar)
List of 14
 $ order   : int 3
 $ ar  : num [1:3]  0.6534 -0.0636 -0.2269
 $ var.pred: num 0.196
 $ x.mean  : num 2.4
 $ aic : Named num [1:17] 18.307  0.996  0.538
 0.000  1.490 ...
  ..- attr(*, names)= chr [1:17] 0 1 2 3 ...
 $ n.used  : int 48
 $ order.max   : num 16
 $ partialacf  : num [1:16, 1, 1]  0.576 -0.223 -0.227
 0.103 -0.076 ...
 $ resid   : Time-Series [1:48] from 1 to 48:
NA NA NA -0.200 -0.169 ...
 $ method  : chr Yule-Walker
 $ series  : chr lh
 $ frequency   : num 1
 $ call: language ar(x = lh)
 $ asy.var.coef: num [1:3, 1:3]  0.02156 -0.01518 
0.00482 -0.01518  0.03117 ...
 - attr(*, class)= chr ar
 

This is for what concerns your questions. I'd like to
give you a hint for a simpler programming code. You
use the first row of your matrix as header
 omega-array(0, dim=c(N+1,7))

omega[1,]-(series,mean,std,max,min,ar(1)coeff,stdar(1))

you can simply type

 omega-array(0, dim=c(N,7))
 colnames(omega) -
c(series,mean,std,max,min,ar(1)coeff,stdar(1))
so your matrix has Nx7 empty cells and avoid an
headache trying to adjust the index 'i' in the loop
'for (i in (2:N+1))'.
Also, if you really need to specify the 'id' of the
record
 omega[i,1]-i-1
you just type before the loop
 omega[,1] - 1:N
or
 rownames(omega) - 1:N

In summary,

 omega-array(0, dim=c(N,7))
 colnames(omega) -
c(series,mean,std,max,min,ar(1)coeff,stdar(1))
 omega[,1] - 1:N
 for (i in 1:N){
 omega[i,2]-mean(y[i,])
 (...)
 omega[i,6]- (???) # see before
 omega[i,7]- (???)} # see before

Finally, if 'y' is a prespecified matrix I suggest you
to take a look at ?rowMeans, ?rowSums, and ?apply in
order to execute 'mean' or 'sd' or any other function
applied to 'y[i,]' at once instead of N times

Hope this helps,
Marco Geraci


--- Lorenzo Bencivelli [EMAIL PROTECTED]
wrote:

 ok, i'll try to be more precise.
 both functions i need (ar() and PP.test()) are in
 the package stats.
 what i would like to do:
 
 omega-array(0, dim=c(N+1,7))

omega[1,]-(series,mean,std,max,min,ar(1)coeff,stdar(1))
 for (i in (2:N+1)){
 omega[i,1]-i-1
 omega[i,2]-mean(y[(i-1),])
 (...)
 omega[i,6]- (???)
 omega[i,7]- (???)}
 
 for the PP.test, i would like to do something
 similar. N is the number 
 of the series for which i would like to compute the
 AR(1) model and the 
 unitary root test. thanks in advance. L
 
 -- 
 
 


 credo nella ragione umana, 
 e nella liberta' e nella giustizia che dalla ragione
 scaturiscono. (sciascia)
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] find position to last number in a vector

2006-05-26 Thread Anders Bjørgesæter
Hello

Is there a function to extract the position (i.e. row number or more 
desirable the value from another column with the same rownumber ) of 
last number in a vector? This is done in a loop with 1000s of columns.

e.g.

vector/column.n: na,na,10,na,2,na,na,1,na,na
fixed column:10,20,30,40,50,60,70,80,85,100

result: 8 and/or 80

Thanks...

Best Regards
Anders

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] curve peeling

2006-05-26 Thread Greg Tarpinian
R2.3.0, WinXP.

I am trying to fit the model

   Y(t) = C1*exp(-exp(a)*t) + C2*exp(-exp(b)*t) - (C1+C2)*exp(-exp(c)*t)

to some pharmacokinetic data.  The data derives from a single oral dose
of a drug.  Does anyone know how to use peeling and/or other methods to 
obtain good starting estimates of the parameters in the triexponential?

The book Mixed Effects Models in S and SPLUS describes self-starting 
models for the biexponential and first-order one-compartment PK models,
but not for the triexponential.


Much thanks in advance,

   Greg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] find position to last number in a vector

2006-05-26 Thread Liaw, Andy
Is this what you want?

 x
 [1] NA NA 10 NA  2 NA NA  1 NA NA
 y
 [1]  10  20  30  40  50  60  70  80  85 100
 y[max(which(!is.na(x)))]
[1] 80

Andy 

From: Anders Bjørgesæter
 
 Hello
 
 Is there a function to extract the position (i.e. row number or more 
 desirable the value from another column with the same rownumber ) of 
 last number in a vector? This is done in a loop with 1000s of columns.
 
 e.g.
 
 vector/column.n: na,na,10,na,2,na,na,1,na,na
 fixed column:10,20,30,40,50,60,70,80,85,100
 
 result: 8 and/or 80
 
 Thanks...
 
 Best Regards
 Anders
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! 
 http://www.R-project.org/posting-guide.html
 


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] pdf output incompatible with latest ghostscript

2006-05-26 Thread ivo welch
I just upgraded to the latest ghostscript, AFPL 8.54.  When I run the
R output through gs, I get

AFPL Ghostscript 8.54: Set UseCUEColor for UseDeviceIndependentColor
to work properly.

followed by a segfault.  the resulting ghostscript file is incomplete
(no %%EOF).  probably a ghostscript bug, but caused by R output.

regards,

/iaw

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Maximum likelihood estimate of bivariatevonmises-weibulldistribution

2006-05-26 Thread Chaouch, Aziz
Thanks Ravi! That's probably more than what I need. I suppose I'll have
to get back to you for some more explanations but for now I'm trying to
apply the log-likelihood (likelihood) that you defined to the fitdistr
function or use the log-likelihood with optim().

Thanks again, that's really helpful!

Regards,

Aziz

-Original Message-
From: Ravi Varadhan [mailto:[EMAIL PROTECTED] 
Sent: May 26, 2006 12:18 PM
To: Chaouch, Aziz
Cc: r-help@stat.math.ethz.ch
Subject: RE: [R] Maximum likelihood estimate of
bivariatevonmises-weibulldistribution

Hi Aziz,

I am attaching a file that contains the functions that I wrote to
compute the MLE for a binary vonMises mixture, using the EM algorithm.
It also contains a simulation example.  Let me know if you have any
trouble.

Hope this is helpful,
Ravi.


--
Ravi Varadhan, Ph.D.
Assistant Professor,  The Center on Aging and Health Division of
Geriatric Medicine and Gerontology Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email:  [EMAIL PROTECTED]
Webpage:http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

--
 -Original Message-
 From: [EMAIL PROTECTED] [mailto:r-help- 
 [EMAIL PROTECTED] On Behalf Of Chaouch, Aziz
 Sent: Friday, May 26, 2006 12:02 PM
 To: Dimitrios Rizopoulos
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] Maximum likelihood estimate of bivariatevonmises- 
 weibulldistribution
 
 Hi,
 
 I'm still strugling with this copula model but this seems to be the 
 way to go. I'm now trying to model the marginal distributions and and 
 for wind direction I use a mixture of 2 von mises. I'd like to 
 estimate all the parameters (m1,m1,kappa1,kappa2,p) by maximizing the 
 likelihood but I don't know how to define the likelihood (or 
 log-likelihood) of a mixture of 2 Von Mises to use it with the 
 function fitDistr in the MASS package. Can you help me define this 
 likelihood function and use it through the fitDistr function?
 
 Thanks,
 
 Aziz
 
 -Original Message-
 From: Dimitrios Rizopoulos 
 [mailto:[EMAIL PROTECTED]
 Sent: May 12, 2006 4:35 PM
 To: Chaouch, Aziz
 Subject: RE: [R] Maximum likelihood estimate of bivariate 
 vonmises-weibulldistribution
 
 look at the following code:
 
 library(copula)
 par(mfrow = c(2, 2))
 x - mvdc(normalCopula(sin(0.5 * pi /2)), c(norm, norm), 
 list(list(mean = 0, sd = 1), list(mean = 0, sd = 1))) contour(x, 
 dmvdc, xlim = c(-2.7, 2.7), ylim = c(-2.7, 2.7))
 
 x - mvdc(frankCopula(5.736276), c(norm, norm), list(list(mean = 
 0, sd = 1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 
 2.7), ylim = c(-2.7, 2.7))
 
 x - mvdc(gumbelCopula(2), c(norm, norm), list(list(mean = 0, sd =

 1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7), 
 ylim = c(-2.7, 2.7))
 
 x - mvdc(claytonCopula(2), c(norm, norm), list(list(mean = 0, sd 
 = 1), list(mean = 0, sd = 1))) contour(x, dmvdc, xlim = c(-2.7, 2.7), 
 ylim = c(-2.7, 2.7))
 
 
 the values of the association parameter I've chosen in each copula 
 correspond to Kendall's tau 0.5; assuming also standard normal 
 marginal distributions look at the different shapes you get!
 
 If possible try something similar for you case (i.e., using von Mises 
 and Weibull marginals) and check if the association shape for a 
 specific copula is more appropriate for your application. If this is 
 not possible fit models assumig different copulas and check which one 
 provides a better fit to your data.
 
 I hope it helps.
 
 Best,
 Dimitris
 
 
 
 Quoting Chaouch, Aziz [EMAIL PROTECTED]:
 
  Hi Dimitris,
 
  I'm not sure to understand your suggestion. How would you build that

  contour plot for a particular copula and on what is computed the 
  Kendall's tau?
 
  Thanks,
 
  Aziz
 
  -Original Message-
  From: Dimitris Rizopoulos
  [mailto:[EMAIL PROTECTED]
  Sent: May 12, 2006 9:57 AM
  To: Chaouch, Aziz; [EMAIL PROTECTED]
  Cc: r-help@stat.math.ethz.ch
  Subject: Re: [R] Maximum likelihood estimate of bivariate 
  vonmises-weibulldistribution
 
  the choice of the copula is, in fact, a model selection problem.
  First, you could have a look at the contour plots of different 
  copulas (preferably for the same value of Kendall's tau), and decide

  if some of them assume a more appropriate association structure for 
  your application, compared to the others. Afterwards, you may fit 
  various copula functions, check the fit on the data, compute AIC 
  (since these are typically not nested models), etc.
 
  regarding the Von Mises distribution and if could be used in mvdc(),

  that I don't know. It'd be better to contact the copula package 
  maintainer and ask.
 
  I hope it helps.
 
  Best,
  Dimitirs
 
  
  Dimitris Rizopoulos
  Ph.D. Student
  Biostatistical Centre
  School of Public Health
  Catholic University of Leuven
 
  Address: Kapucijnenvoer 35, Leuven, 

Re: [R] find position to last number in a vector

2006-05-26 Thread apjaworski
I am sure there are several ways of doing this but how about something like
this:

y - c(NA,NA,10,NA,2,NA,NA,1,NA,NA)
x - c(10,20,30,40,50,60,70,80,85,100)

tail(x[!is.na(y)], 1)

or

rev(x[!is.na(y)])[1]

Rudimentary checks indicate that the second expression is much faster than
the first.  Also, if y is contains NAs only, the first expression returns
numeric(0) and the second one NA.

Hope this helps,

Andy

__
Andy Jaworski
518-1-01
Process Laboratory
3M Corporate Research Laboratory
-
E-mail: [EMAIL PROTECTED]
Tel:  (651) 733-6092
Fax:  (651) 736-3122


   
 Anders
 Bjørgesæter   
 anders.bjorgesat  To 
 [EMAIL PROTECTED]r-help@stat.math.ethz.ch
 Sent by:   cc 
 [EMAIL PROTECTED] 
 at.math.ethz.ch   Subject 
   [R] find position to last number in 
   a vector
 05/26/2006 01:30  
 PM
   
   
   
   




Hello

Is there a function to extract the position (i.e. row number or more
desirable the value from another column with the same rownumber ) of
last number in a vector? This is done in a loop with 1000s of columns.

e.g.

vector/column.n: na,na,10,na,2,na,na,1,na,na
fixed column:10,20,30,40,50,60,70,80,85,100

result: 8 and/or 80

Thanks...

Best Regards
Anders

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] how to get a series of results by loop

2006-05-26 Thread Gabor Grothendieck
Try this:

set.seed(1)
r - .1
n - 5
test - data.frame(sp = paste(sp, 1:n, sep = ), x= rnorm(n), y = rnorm(n))
r - 1 # radius
f - function(i) cbind(orig = test[i,1], subset(test, (x - x[i])^2 +
(y - y[i])^2  r^2))
do.call(rbind, lapply(1:n, f))


On 5/26/06, zhang jian [EMAIL PROTECTED] wrote:
 I donot know how to do it remainly.
 Perhaps my question is not clear.
 For example,

 test=data.frame(sp=paste(c(sp),1:100,sep=),x=rnorm(100),y=rnorm(100))

 # I want to get the data in the circle of radius = 1 in every point, and I
 want to add a sign (e.g. sp1,sp2) in order that I can distinguish
 different groups
 # if I want to choose a point, I can do it like this

 r = 1 # radius
 attach(test)
 point100=test[(x-x[100])^2+(y-y[100])^2=r^2,]
 result100=data.frame
 (note=paste(c(p),rep(100,length(point100$sp)),sep=),point100)
 detach(test)

 But if I want to do all points in test, I donot know how to do it. I try
 to do it by using for, but it doesnot work.
 Can you give me some advices?
 thanks!
  jian zhang

[[alternative HTML version deleted]]

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] PC rotation question

2006-05-26 Thread Sasha Pustota
I wrote:
 On p. 48 of Statistics Complements to the 3rd MASS edition,
 http://www.stats.ox.ac.uk/pub/MASS3/VR3stat.pdf
 I read that the orthogonal rotations of Z Lambda^-1 remain
 uncorrelated, where Z is the PC and Lambda is the diag matrix of
 singular values. However, the example below that text is
 A - loadings(ir.pca) %*% diag(ir.pca$sdev)
 If ir.pca$sdev are the singular values, should that be diag(1 /
 ir.pca$sdev), or is it some discrepancy between S+ and R that I'm
 missing?


To dwell some more on this and in hopes to get replies I did a small
experimantation (below) that makes me to suspect that the correct
syntax
is A - loadings(ir.pca) %*% diag(1/ir.pca$sdev). Any comments?

library(MASS)
a - 1; b - 0.2; c - -.3; d - .7
x - scale(mvrnorm(n=1000,c(0,0,0),matrix(c(a,b,c, b,a,d, c,d,a),3,3)))
e - eigen(cov(x))

cat(expect identity corr for orthogonal rotations:\n)
zs - e$vectors %*% diag(1/sqrt(e$values))
pdx - x %*% varimax(zs, normalize = FALSE)$loadings
print.table(cor(pdx), digits=2)

cat(expect zero: , prcomp(x)$sdev - sqrt(e$values), \n)

cat(Now zero corr between projections is not preserved:\n)
zs - e$vectors %*% diag(sqrt(e$values))
pdx - x %*% varimax(zs, normalize = FALSE)$loadings
print.table(cor(pdx), digits=2)

Output:
expect identity corr for orthogonal rotations:
 [,1] [,2] [,3]
[1,]  1.0e+00 -3.0e-16 -4.5e-16
[2,] -3.0e-16  1.0e+00  2.3e-16
[3,] -4.5e-16  2.3e-16  1.0e+00
expect zero:  4.440892e-16 2.220446e-16 2.775558e-16
Now zero corr between projections is not preserved:
 [,1]  [,2]  [,3]
[1,]  1.00 -0.35 -0.88
[2,] -0.35  1.00 -0.11
[3,] -0.88 -0.11  1.00

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] missed ylim from plot.default

2006-05-26 Thread Greg Snow
Try the following function to see if it does what you want.  The basic
syntax for your example would be:

 tmp - squishplot( xlim=c(-0.5,7), ylim=c(-0.5,2.8), asp=1 )
 plot(1:5, rep(1,5), xlim=c(-0.5,7), ylim=c(-0.5,2.8))
 par(tmp) # reset plotting region for future plots

The function definition is: 

squishplot - function(xlim,ylim,asp=1){
  if(length(xlim)  2) stop('xlim must be a vector of length 2')
  if(length(ylim)  2) stop('ylim must be a vector of length 2')

  tmp - par(c('plt','pin','xaxs','yaxs'))

  if( tmp$xaxs == 'i' ){ # not extended axis range
xlim - range(xlim)
  } else { # extended range
tmp.r - diff(range(xlim))
xlim - range(xlim) + c(-1,1)*0.04*tmp.r
  }

  if( tmp$yaxs == 'i' ){ # not extended axis range
ylim - range(ylim)
  } else { # extended range
tmp.r - diff(range(ylim))
ylim - range(ylim) + c(-1,1)*0.04*tmp.r
  }


  tmp2 - (ylim[2]-ylim[1])/(xlim[2]-xlim[1])

  tmp.y - tmp$pin[1] * tmp2 * asp

  if(tmp.y  tmp$pin[2]){ # squish vertically
par(pin=c(tmp$pin[1], tmp.y))
par(plt=c(tmp$plt[1:2], par('plt')[3:4]))
  } else { # squish horizontally
tmp.x - tmp$pin[2]/tmp2/asp
par(pin=c(tmp.x, tmp$pin[2]))
par(plt=c(par('plt')[1:2], tmp$plt[3:4]))
  }

  return(invisible(tmp['plt']))
}

Let me know of any improvments you can think of (including a better
name) and how it works.





-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111
 

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Antonio, Fabio Di
Narzo
Sent: Friday, May 26, 2006 1:29 AM
To: Peter Ehlers
Cc: R-help@stat.math.ethz.ch
Subject: Re: [R] missed ylim from plot.default

2006/5/26, Peter Ehlers [EMAIL PROTECTED]:


 Antonio, Fabio Di Narzo wrote:

  Hi all.
  On my R-2.3.0,  calling:
 
 
 plot(1:5, rep(1,5), xlim=c(-0.5,7), ylim=c(-0.5,2.8), asp=1)
 
 
  gives to me a plot (tried pdf and X11) whose y axis goes from about 
  -2
 to
  about 4.5. What I've missed? How can I show some pre-specified x and

  y ranges from a plot (keeping fixed the aspect ratio)?
 
  Tnx all,
  Antonio, Fabio Di Narzo.
 
[[alternative HTML version deleted]]

 Do you want a wide, but not very tall plot? If so, you could specify 
 the width/height in your X11() call.

 Peter Ehlers


Tnx for your indication, now I see: x and y plotting ranges are
influenced by the actual window size (at least when fixing aspect
ratio). I think this is obvious... What I missed/forgot above, was that
default plot width and height are fixed indipendently from the 'x-ylim'
and 'asp' arguments.

However now I have a problem: I can write a wrapper to set device width
and height depending on the (xlim,ylim,asp) triple. But there's often
extra-space in the device due to plot title, axis labels, etc. So, I
fear I can't simply set width=height*asp. Some suggestion?
(Anyway, by now I only have to produce few plots, so today I can go try
by try :-) ).

Bests,
Antonio.

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] R.oo question

2006-05-26 Thread Omar Lakkis
This is a simple R.oo question but I, thankfully, hope that someone
would explain it to me so I would better understand this work frame.
I create this class:

setConstructorS3(MyExample, function(param=0) {
  print(paste(called with param=, param))
  extend(Object(), MyExample,
.param = param
  );
})

From what is printed out, who made the second call to the class with
the default param?

 MyExample(1)
[1] called with param= 1
[1] called with param= 0# - this is line in question
[1] MyExample: 0x02831708

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] R newbie attempting to plot data

2006-05-26 Thread Rex Eastbourne
Hi,

I just started using R and am having trouble with the below error:

I type:

 df - read.csv(/home/rex/Desktop/mytable.csv)

which gives me what I want:

...
639 2006-05-26 16:46:54 4   16
640 2006-05-26 17:05:36 5   17
641 2006-05-26 17:30:48 6   17

But now I try:

 plot(df[4],df[3])
Error in pmatch(x, table, duplicates.ok) :
argument is not of mode character

Any ideas? I tried read.table and importing the data from MySQL too; same
error.

Thanks,

Rex

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] R newbie attempting to plot data

2006-05-26 Thread Chuck Cleland
Rex Eastbourne wrote:
 Hi,
 
 I just started using R and am having trouble with the below error:
 
 I type:
 
 df - read.csv(/home/rex/Desktop/mytable.csv)
 
 which gives me what I want:
 
 ...
 639 2006-05-26 16:46:54 4   16
 640 2006-05-26 17:05:36 5   17
 641 2006-05-26 17:30:48 6   17
 
 But now I try:
 
 plot(df[4],df[3])
 Error in pmatch(x, table, duplicates.ok) :
 argument is not of mode character

You probably intended:

plot(df[,4],df[,3])

 Any ideas? I tried read.table and importing the data from MySQL too; same
 error.
 
 Thanks,
 
 Rex
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
 

-- 
Chuck Cleland, Ph.D.
NDRI, Inc.
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] R.oo question

2006-05-26 Thread Henrik Bengtsson
On 5/26/06, Omar Lakkis [EMAIL PROTECTED] wrote:
 This is a simple R.oo question but I, thankfully, hope that someone
 would explain it to me so I would better understand this work frame.
 I create this class:

 setConstructorS3(MyExample, function(param=0) {
   print(paste(called with param=, param))
   extend(Object(), MyExample,
 .param = param
   );
 })

 From what is printed out, who made the second call to the class with
 the default param?

  MyExample(1)
 [1] called with param= 1
 [1] called with param= 0# - this is line in question
 [1] MyExample: 0x02831708

That is because of a new US government rule requiring that one
instance of every R.oo class is forwarded to the NSA.

Seriously, you will see that this happens once and only once for each
class defined this way and it normally happens when you create your
first instance (otherwise before that).  The first call to the
constructor creates a so called static instance of the class.  The
static instance is for instance used when you type 'MyExample' to get
the API coupled to class MyExample.  Another example:

 Object
Object
  public $(name)
  public $-(name, value)
  public [[(name)
  ...
  public objectSize(...)
  public print(...)
  public save(file=NULL, path=NULL, compress=TRUE, ...)
}

Technical details: The static instance is not always needed, but quite
often.  If needed, it has to be create at some stage and I found that
having extend() to do this is quite convenient.  The alternative would
be to let you do it explicitly, e.g. getStaticInstance(Object).  Note
that there is no central registry/database keeping track of the
classes defined by R.oo/Object, but all is just plain S3 classes
making use of the S3/UseMethod dispatching mechanisms - that's all.

There is a low-level way to figure out if the call to the constructor
is for the static instance or not, but normally it is easier not to
output anything in the constructor.  If you want to know the low-level
way, I'll have have to get back to, because I don't know it off the
top of my head.  I might provide a method for this if there is a need
for it, e.g. hasStaticInstance(MyExample).

Cheers

Henrik
(the author)

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html



__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Too many open files

2006-05-26 Thread Omar Lakkis
This may be more of an OS question ...
I have this call

r = get.hist.quote(symbol, start= format(start, %Y-%m-%d), end=
format(end, %Y-%m-%d))
which does a url request

in a loop and my program runs out of file handlers after few hundred
rotations. The error message is:  'Too many open files'. Other than
increasing the file handlers assigned to my process, is there a way to
cleanly release and reuse these connections?

I run on Mac OS X

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Too many open files

2006-05-26 Thread Seth Falcon
Omar Lakkis [EMAIL PROTECTED] writes:

 This may be more of an OS question ...
 I have this call

 r = get.hist.quote(symbol, start= format(start, %Y-%m-%d), end=
 format(end, %Y-%m-%d))
 which does a url request

 in a loop and my program runs out of file handlers after few hundred
 rotations. The error message is:  'Too many open files'. Other than
 increasing the file handlers assigned to my process, is there a way to
 cleanly release and reuse these connections?

Inside your loop you need to close the connection object created by
url().

for (i in 1:500) {
con - url(urls[i])
## ... stuff here ...
close(con)
}

R only allows you to have a fixed number of open connections at one
time, and they do not get closed automatically when they go out of
scope.  These commands may help make clear how things work...

 showConnections()
 description class mode text isopen can read can write
 f = url(http://www.r-project.org;, open=r)
 showConnections()
  descriptionclass mode text   isopen   can read can write
3 http://www.r-project.org; url r  text opened yesno 
 rm(f)
 showConnections()
  descriptionclass mode text   isopen   can read can write
3 http://www.r-project.org; url r  text opened yesno 
 f - getConnection(3)
 close(f)
 f
Error in summary.connection(x) : invalid connection
 showConnections()
 description class mode text isopen can read can write
 

+ seth

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Too many open files

2006-05-26 Thread Gabor Grothendieck
Try closeAllConnections()

On 5/26/06, Omar Lakkis [EMAIL PROTECTED] wrote:
 This may be more of an OS question ...
 I have this call

 r = get.hist.quote(symbol, start= format(start, %Y-%m-%d), end=
 format(end, %Y-%m-%d))
 which does a url request

 in a loop and my program runs out of file handlers after few hundred
 rotations. The error message is:  'Too many open files'. Other than
 increasing the file handlers assigned to my process, is there a way to
 cleanly release and reuse these connections?

 I run on Mac OS X

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Building V2.3.0 on Tru64 V5.1B

2006-05-26 Thread Ravi, Narendra, INFOT
Hi,
I am going to be attempting a build of R for the alphas (our
main need is to collate and present graphically, performance metrics we
gather on our applications).

We have:

CC=cc, CXX=cxx, Fortran V5.51 (f77, f95), GNU make 3.80, Perl 5.8.4

I am hesitant to rebuild GCC on this platform, and would like to know if
others have success in building R with native C, C++ and Fortran. I plan
on using the following configure string:

CPPFLAGS=-I/opt/gnu/include -I/opt/libjpeg/include
-I/opt/libpng/include -I/opt/ncurses/include CFLAGS=-std1
LDFLAGS=-L/opt/gnu/lib -L/opt/libjpeg/lib -L/opt/libpng/lib
-L/opt/ncurses/lib -lncurses -ljpeg -lpng -rpath
/opt/libpng/lib:/opt/libjpeg/lib:/opt/ncurses/lib:/opt/gnu/lib:/usr/shl
ib R_PAPERSIZE=letter CC=cc CXX=cxx ./configure prefix=/opt/R
--without-readline --without-blas

I expect to update the config.site file with the above values. If
someone has built R-2.3.0 on Tru64 V5.1B, I would appreciate hearing
how.

Thank you,

Narendra Ravi

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Trouble passing list or non-list to function using ...

2006-05-26 Thread Jason Barnhart
Hello,

Simply put, I'm trying to call a function testme with value age=NA.   I 
wish to use dotlist-list(...)  inside the function and have dotlist become:
$age
[1] NA

I'm modifying existing code and need to minimize changing that code so it's 
easiest to conform 
how I call the existing function.

My sample code fragment, results, and R.version information are listed below.   
I've been searching 
for prior questions regarding this phenomena, but haven't quite figured out how 
to do this.

Thank you in advance,
-jason



##BEGIN TEST CODE ###
testme-function(...){

dotlist-list(...)

dotlist

}

nm-age

testme(age=NA)

testme(nm=NA)

testme(age=NA,age2=NA)

testlist-list(age,age2)

for (i in 1:length(testlist)){

print(testme(testlist[i]))

}

##END TEST CODE #

##BEGIN RESULTS #
 testme-function(...){
+ 
+ dotlist-list(...)
+ 
+ dotlist
+ 
+ }
 
 nm-age
 
 testme(age=NA)   *This is what I'm really after
$age
[1] NA

 
 testme(nm=NA)
$nm
[1] NA

 #trying w/ 2 vars
 testme(age=NA,age2=NA)
$age
[1] NA

$age2
[1] NA

 
 testlist-list(age,age2)
 
 for (i in 1:length(testlist)){
+ 
+ print(testme(testlist[i]))
+ 
+ }
[[1]]
[[1]][[1]]
[1] age


[[1]]
[[1]][[1]]
[1] age2

##END RESULTS ###

platform   i386-pc-mingw32   
arch   i386  
os mingw32   
system i386, mingw32 
status   
major  2 
minor  3.0   
year   2006  
month  04
day24
svn rev37909 
language   R 
version.string Version 2.3.0 (2006-04-24)

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Trouble passing list or non-list to function using ...

2006-05-26 Thread Gabor Grothendieck
Its not clear from your post which of the tests does not give you
the desired output and what the desired output is in that case.

Just guessing but maybe you meant to use testlist[[i]] instead of testlist[i] in
the loop?


On 5/26/06, Jason Barnhart [EMAIL PROTECTED] wrote:
 Hello,

 Simply put, I'm trying to call a function testme with value age=NA.   I 
 wish to use dotlist-list(...)  inside the function and have dotlist become:
 $age
 [1] NA

 I'm modifying existing code and need to minimize changing that code so it's 
 easiest to conform
 how I call the existing function.

 My sample code fragment, results, and R.version information are listed below. 
   I've been searching
 for prior questions regarding this phenomena, but haven't quite figured out 
 how to do this.

 Thank you in advance,
 -jason



 ##BEGIN TEST CODE ###
 testme-function(...){

 dotlist-list(...)

dotlist

 }

 nm-age

 testme(age=NA)

 testme(nm=NA)

 testme(age=NA,age2=NA)

 testlist-list(age,age2)

 for (i in 1:length(testlist)){

print(testme(testlist[i]))

 }

 ##END TEST CODE #

 ##BEGIN RESULTS #
  testme-function(...){
 +
 + dotlist-list(...)
 +
 + dotlist
 +
 + }
 
  nm-age
 
  testme(age=NA)   *This is what I'm really after
 $age
 [1] NA

 
  testme(nm=NA)
 $nm
 [1] NA

  #trying w/ 2 vars
  testme(age=NA,age2=NA)
 $age
 [1] NA

 $age2
 [1] NA

 
  testlist-list(age,age2)
 
  for (i in 1:length(testlist)){
 +
 + print(testme(testlist[i]))
 +
 + }
 [[1]]
 [[1]][[1]]
 [1] age


 [[1]]
 [[1]][[1]]
 [1] age2

 ##END RESULTS ###

 platform   i386-pc-mingw32
 arch   i386
 os mingw32
 system i386, mingw32
 status
 major  2
 minor  3.0
 year   2006
 month  04
 day24
 svn rev37909
 language   R
 version.string Version 2.3.0 (2006-04-24)

[[alternative HTML version deleted]]

 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html