Re: [R] bptest

2010-09-25 Thread Andrew Miles

First load the package lmtest.  Then run the bptest.

library(lmtest)
bptest(modelCH)

You don't need to tell the function which variables are explanatory or  
dependent.  Just give it the fitted model object, and it will sort all  
of that out and return the statistic.


Andrew Miles
Department of Sociology
Duke University

On Sep 24, 2010, at 9:53 AM, robm wrote:



Hi

I'm very new to R but have plenty of experience with statistics and  
other

packages like SPSS, SAS etc.

I have a dataset of around 20 columns and 200 rows.   I'm trying to  
fit a
very simple linear model between two variables.   Having done so, I  
want to

test the model for heteroscedasticity using the Breusch-Pagan test.
Apparently this is easy in R by simply doing

bptest(modelCH, data=KP)

I've tried this but I'm told it cannot find function bptest.   It's  
here
where I'm struggling.   I'm probably wrong but as far as I can see,  
bptest

is part of the lm package which, as far as I know, I have installed.

Irrespective of the fact I'm not sure how to tell bptest which is the
dependent and explanatory variables - there's a more fundamental  
problem if

it can't find the bptest function.

I have searched the documentation - albeit briefly so if anyone  
could help

I'd be very grateful

Rob

QBE Management
--
View this message in context: 
http://r.789695.n4.nabble.com/bptest-tp2553506p2553506.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] How to store regex expression in a variable

2010-09-25 Thread Yong Wang
dear list

I know how to store a regex expression in perl and ruby, no clue on R.
I do read R regex manual , archives, and searched on line,
still I need somebody help me out on how to store a regular expression
in a variable.

Thank you very much

yong

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Multiple graph in one graph window

2010-09-25 Thread Andrew Miles
Another option is to use the par() command before executing a plot  
command, like so:

par(mfrow=c(2,2))
plot(...)

This will put the next four plots all in the same window.

Andrew Miles
Department of Sociology
Duke University

On Sep 24, 2010, at 10:28 PM, Dennis Murphy wrote:

 Which one do you want?

 library(lattice)
 d - data.frame(x = rep(1:3, 4), g = factor(rep(1:4, each = 3)),
 y = c(1, 2, 3, 3, 1, 7, 4, 2, 11, 5, 5, 9))
 xyplot(y ~ x | g, data = d, pch = 16)
 xyplot(y ~ x, data = d, groups = g, type = c('p', 'l'), pch = 16)

 Both provide 'multiple plots' and both are on the 'same page'.

 HTH,
 Dennis


 On Fri, Sep 24, 2010 at 2:37 PM, mrdaw...@abo.fi wrote:

 Dear R users;
 Could you tell me how to let R create multiple graphs in a graph  
 window.
 Please note that I am talking about creating multiple plots in the  
 same page
 of the graph window. For example:
 g1=c(1,2,3)
 g2=c(3,1,7)
 g3=c(4,2,11)
 g4=c(5,5,9)
 ...
 all in one graph window.
 Best Regards
 Reza

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


   [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] why I could not reproduce the Mandelbrot plot demonstrated on R wiki

2010-09-25 Thread Peter Dalgaard
On 09/24/2010 10:41 PM, xin wei wrote:
 
 I am trying to reproduce the nice looking of Mandelbrot demonstrated by R
 wiki page by the following code:
 
 library(caTools)# external package providing write.gif function
 jet.colors = colorRampPalette(c(#7F, blue, #007FFF, cyan,
 #7FFF7F, 
 yellow, #FF7F00, red, #7F)) 
 m = 600 # define size
 C = complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ), 
  imag=rep(seq(-1.2,1.2, length.out=m), m ) ) 
 C = matrix(C,m,m)   # reshape as square matrix of complex numbers
 Z = 0   # initialize Z to zero
 X = array(0, c(m,m,20)) # initialize output 3D array
 for (k in 1:20) {   # loop with 20 iterations
   Z = Z^2+C # the central difference equation  
   X[,,k] = exp(-abs(Z)) # capture results
 } 
 write.gif(X, Mandelbrot.gif, col=jet.colors, delay=100)

Hmm, I couldn't be bothered with the caTools, but it looks fine for me with

image(X[,,20],col=jet.colors(100))

Perhaps you need jet.colors(n) as well?



 however, the gif file created by this looks much worse than what is shown on
 R wiki page, see the comparison as follows (left one is what i created)
 
 http://r.789695.n4.nabble.com/file/n2591429/Picture1.png 

Save for the odd color scheme, the one on the left looks like a
Mandelbrot set, the one on the right appears to be iteration 4. I
couldn't find your original source for this on wiki.r-project.org?


-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to store regex expression in a variable

2010-09-25 Thread Peter Dalgaard
On 09/25/2010 03:45 AM, Yong Wang wrote:
 dear list
 
 I know how to store a regex expression in perl and ruby, no clue on R.
 I do read R regex manual , archives, and searched on line,
 still I need somebody help me out on how to store a regular expression
 in a variable.

A regex is just a character string in R. So, e.g.

  str - 'Now is the time  '
  sub(' +$', '', str)  ## spaces only
[1] Now is the time
 re - ' +$'
  sub(re, '', str)  ## spaces only
[1] Now is the time

-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Reading in .aux (ESRI raster files) into R

2010-09-25 Thread Barry Rowlingson
On Fri, Sep 24, 2010 at 3:21 PM, JoH jh...@york.ac.uk wrote:


 Error in
 readAsciiGrid(F:/GIS.LandcoverEuropeForRisk/Sept10kmmaps/Sp10KPointID.aux)
 :
  object 'cellsize' not found

 My original data in Arc GIS is have a cell size an i'mm curious as to how to
 make sure all the details are included.

 Well at this point absolutely none of the file has been included.

 I also tried to use

 Spain10km-getRasterData(F://RMap//sp10kpointid1, ,band=NULL,

 This is a different file to the one you tried above.

  Spain10km - system.file(F://RMap//sp10kpointid1, package=rgdal)

 That won't do anything useful.

 #These lines were used to then try to get it to recognise sp10kpointid1 as a
 member og the GDALReadOnly Dataset.
 x - new(GDALReadOnlyDataset, Spain10km)
 Error in .local(.Object, ...) : empty file name
 x - new(sp10kpointid1, Spain10km)


 Which is the best way to read this file into R? Why didn't it include my
 classes?

 Assuming it *is* a valid raster file, then readGDAL from the rgdal
package should handle it:

 require(rgdal)
 data = readGDAL(file.choose())

Using file.choose() here makes sure you dont muck up the path because
it will prompt you to browse for the file. readGDAL is an interface to
all the raster data types that rgdal knows about.

 If this fails, then I am wondering if the file is a valid raster
file. .aux is an odd file extension - how did you create it exactly?

Barry

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Power Function

2010-09-25 Thread jethi


hi, the part with the power function does not work. this is the   gute  
function in my programm. because i´m not shure on what dependet my power 
function gute. and ofcourse the graphik does not work. thnx a lot:)

regards
Kaja

Date: Fri, 24 Sep 2010 19:21:17 -0700
From: ml-node+2653667-1721823656-167...@n4.nabble.com
To: kart...@hotmail.com
Subject: Re: Power Function



Hi,


What part of your code does not work or does not do what you want?  I

am guessing that you want gute to return more than a single value

(which is what it does in your code).  Can you show us an example of

the results you would like?


Best regards,


Josh




I reworked some of your code to simplify, clean, and wrap it all in a function:


myfun - function(n, m, alpha = .05, seeder = 1000) {

  set.seed(seeder)

  x - matrix(rnorm(n, 0, 0.5), ncol = m)

  y - matrix(rnorm(n, 0, 0.8), ncol = m)

  l - diag(cor(x, y))

  cat(Correlations between two random variables \n, l, fill = TRUE)

  gute - function(x, m, alpha) {

q_1 - qnorm(alpha, 0, 0.05)

q_2 - qnorm(1 - alpha, 0, 0.05)

p - (x^2)/sum(x^2)

H - log(m) - sum(p * log(p), na.rm = TRUE)

1 - mean(q_1 = H  H = q_2)

  }

  dat - seq(0, 1, length.out = 10)

  output - gute(x = dat, m = m, alpha = alpha)

  return(output)

}


This results in:


 myfun(100, 5)

Correlations between two random variables

 -0.5887829 0.07042411 -0.06082638

0.2395193 -0.1038213

[1] 1



On Fri, Sep 24, 2010 at 5:26 PM, jethi [hidden email] wrote:



 Hi, at first, i´m from germany, so sorry for my bad english. but i need ur

 help in R to programm a power function  and to make at last a graphik of it.



 i have already tried my best. but it doesn´t work.the topic is: the

 homogeneity test of correlation based entropy.



 so it means, that i have to check if all correlations of a bivariate random

 vectors are same or not. for that i saperate the  n bivariate random vectors

 (x1,y1),...,(xn,yn)  in blocks  m so, so that i at first calculate  the

 correlation in a  block. n=m*k. the numbers of the blocks m are

 user-defined. the test value is the entropy. pls help me!!!





 set.seed(1000)

 n=100

 m=5

 k=n/m



 x=rnorm(n,0,0.5)

 y=rnorm(n,0,0.8)





  #alpha/2 Quantil

  q_1=qnorm(0.05,0,0.05)



  #1-alpha/2 Quantile

  q_2=qnorm(0.95,0,0.05)





  l=matrix(0,nrow=m,ncol=1)

 for(i in 1:m){





 l[i]=print(cor((x[(((i-1)*k)+1):(((i-1)*k)+k)]),

 (y[(((i-1)*k)+1):(((i-1)*k)+k)])))





 }

 güte=function(l){

 p=matrix(0,nrow=m,ncol=1)

 for(i in 1:m){

 p[i]=l[i]^2/sum(l^2)



 }

 H=log(m)-sum(p*log(p))

 1-mean(q_1=H  H =q_2)

 }





 l=seq(0,1,len=10)











 plot(l,güte, type=o,pch=20,ylim=c(0,1),col=red)

 --

 View this message in context: 
 http://r.789695.n4.nabble.com/Power-Function-tp2631929p2631929.html
 Sent from the R help mailing list archive at Nabble.com.



 __

 [hidden email] mailing list

 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



--

Joshua Wiley

Ph.D. Student, Health Psychology

University of California, Los Angeles

http://www.joshuawiley.com/

__

[hidden email] mailing list


PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.







View message @ 
http://r.789695.n4.nabble.com/Power-Function-tp2631929p2653667.html


To unsubscribe from Power Function, click here.


  
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Power-Function-tp2631929p2694493.html
Sent from the R help mailing list archive at Nabble.com.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Measure Difference Between Two Distributions

2010-09-25 Thread Lorenzo Isella

Dear All,
Suppose you are given two distributions (or better: two equally-sized 
lists of data); how can you evaluate the difference between them?
I need something like an overlap measure of the two (let us say 0 == no 
overlap and 1== complete overlap). I should add that there is a 1-1 
correspondence of the data in the two distributions (they are ordered 
lists and e.g. the third element in the first distribution matches the 
third element in the second distribution).
The two distributions are not analytical (I mean they do not belong to 
any well known family).

I wonder if mutual information could be what I am looking for.
Cheers

Lorenzo

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] labels in (box)plot

2010-09-25 Thread Peter Ehlers

On 2010-09-21 8:23, Ivan Calandra wrote:

   Dear users,

I would like all the ticks on a boxplot (x and y) to be labeled
I have checked all the par() arguments but couldn't find what I'm
looking for

Here is an example to show it:
df- structure(list(SPECSHOR = structure(c(1L, 1L, 1L, 3L, 3L, 3L, 3L,
3L, 4L, 4L), .Label = c(cotau, dibic, eqgre, gicam), class =
factor), Sq122.median = c(2.335835, 1.76091, 1.64717, 1.285505,
1.572405, 1.86761, 1.82541, 1.62458, 0.157813, 0.864523)), .Names =
c(SPECSHOR, Sq122.median), class = data.frame, row.names = c(9L,
16L, 23L, 74L, 83L, 90L, 98L, 109L, 121L, 139L))

par(mfrow=c(2,2))  ## not necessary in that example, but I do need it
with my real data
boxplot(Sq122.median~SPECSHOR, data=df, horizontal=TRUE)  ## the upper
label on the y-axis (gicam) is not shown

I know I can decrease the size of the labels by setting cex.axis=0.8,
but I would prefer that all labels are always there, independent of
their size, without me setting their size explicitly, and even if they
then overlap.

Am I clear? Is it possible, and how?

Thanks in advance
Ivan



Two suggestions:

1) use staxlab() in the plotrix package (but the staggered
labels might look a bit strange if you have short labels).

 library(plotrix)
 boxplot(Sq122.median~SPECSHOR, data=df, horizontal=TRUE, axes=FALSE)
 box(); axis(1)
 staxlab(side=2, at=1:4, lab=levels(df$SPECSHOR), cex=.8)


2) use mtext() to control your labelling.

 boxplot(Sq122.median~SPECSHOR, data=df, horizontal=TRUE, axes=FALSE)
 box(); axis(1)
 axis(2, at=1:4, lab=NA)
 mtext(levels(df$SPECSHOR), side=2, at=1:4, line=1, cex=0.8)


  -Peter Ehlers

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] (no subject)

2010-09-25 Thread Kartija Oompragash

hi how can i plot now this function??? have to be m= 2???  because of the 
dimensions?thanks for ur help


myfun - function(n, m, alpha = .05, seeder = 1000) {
 
  set.seed(seeder)
 
  x - matrix(rnorm(n, 0, 0.5), ncol = m)
 
  y - matrix(rnorm(n, 0, 0.8), ncol = m)
 
  l - diag(cor(x, y))
 
  cat(Correlations between two random variables \n, l, fill = TRUE)
 
  gute - function(x, m, alpha) {
 
q_1 - qnorm(alpha, 0, 0.05)
 
q_2 - qnorm(1 - alpha, 0, 0.05)
 
p - (x^2)/sum(x^2)
 
H - log(m) - sum(p * log(p), na.rm = TRUE)
 
1 - mean(q_1 = H  H = q_2)
 
  }
 
  dat - seq(0, 1, length.out = 10)
 
  output - gute(x = dat, m = m, alpha = alpha)
 
  return(output)
 
}


  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] plot.ts versus plot.zoo

2010-09-25 Thread Achim Zeileis

On Wed, 22 Sep 2010, Alex van der Spek wrote:


plot.ts has an argument yax.flip, plot.zoo does not.


Thanks for the info. I've just added an yax.flip argument to plot.zoo in 
the devel version of zoo on R-Forge.


hth,
Z


Is there a way to make the yaxis flip in plot.zoo?

I tried using a custom panel function:

panel.yaxis-function(...) {
   npnl-parent.frame$panel.number
   if (npnl %% 2 == 0) {
   axis(side=3)
   } else {
   axis(side=2)
   }
}

This leads to a blank window. I am stuck with the intricacies of the plotting 
and axis generation here.


Thank you in advance,
Alex van der Spek

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Measure Difference Between Two Distributions

2010-09-25 Thread Dennis Murphy
Hi:

On Sat, Sep 25, 2010 at 3:47 AM, Lorenzo Isella lorenzo.ise...@gmail.comwrote:

 Dear All,
 Suppose you are given two distributions (or better: two equally-sized lists
 of data); how can you evaluate the difference between them?
 I need something like an overlap measure of the two (let us say 0 == no
 overlap and 1== complete overlap). I should add that there is a 1-1
 correspondence of the data in the two distributions (they are ordered lists
 and e.g. the third element in the first distribution matches the third
 element in the second distribution).


To visualize the two distributions, you could do an empirical Q-Q plot
(qqplot(x, y)); if the distributions are identical, they should lie on a 45
degree line - location shifts are indicated by level shifts (parallel lines)
to the 45 degree line, differences in scale by slope differences away from
1.

There are many types of tests to test equality of distribution, the simplest
one being the two-sample Kolmogorov-Smirnov test. It is sensitive to changes
in both location and scale of the two edf's (empirical distribution
functions). An alternative
is the two-sample Cramer-von Mises test. I have no doubt there are others...


HTH,
Dennis

The two distributions are not analytical (I mean they do not belong to any
 well known family).
 I wonder if mutual information could be what I am looking for.
 Cheers

 Lorenzo

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Measure Difference Between Two Distributions

2010-09-25 Thread Lorenzo Isella

On 09/25/2010 03:23 PM, Rainer M Krug wrote:



Evaluate, for me, does not necessary mean test if they are
significantly different, but rather to quantify the difference. If that
is what you are looking for, you could look at the Earth Movers
Distance, where a package is available at R-forge
(https://r-forge.r-project.org/projects/earthmovdist/) which I co-wrote
and used before.

Cheers,

Rainer



Thanks Rainer. I had a quick look at wikipedia and the package you 
mention, and it seems what I am looking for.
Just a question about normalization of the distance calculated by the 
algorithm.
Let us say that I have 4 distributions A,B,C,D coupled this way (A,B) 
and (C,D).
The length of data in A is equal to the length of data in B, same 
applies to C and D but length(A)!=length(C).
Now, the argument I would like to make is that A and B are more similar 
than C and D and show a couple of numbers to prove this.
Bottom line: provided my data lists are long enough, does this distance 
scale with the number of data? and if they do, how should I normalize 
this distance to compare the results?

Cheers

Lorenzo

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Uncertainty propagation

2010-09-25 Thread Maayt

I have a small model running under R. This is basically running various
power-law relations on a variable (in this case water level in a river)
changing spatially and through time. I'd like to include some kind of error
propagation to this.
My first intention was to use a kind of monte carlo routine and run the
model many times by changing the power law parameters. These power laws were
obtained by fitting data points under R. I thus have std error associated to
them: alpha (±da) * WaterHight ^ beta (±db). Is it statistically correct to
sample alpha and beta for each run by picking them from a normal
distribution centered on alpha (resp. beta) with a standard deviation of da
(resp. db) and to perform my statistics (mean and standrad edviation of the
model result) on the model output?
It seems to me that da and db are correlated in some way and by doing what I
entended to, I would overestimate the final error of my model...
My statistical skills are rather weak, is there a way people usually deal
with this kind of problems?

Thanks

-- 
View this message in context: 
http://r.789695.n4.nabble.com/Uncertainty-propagation-tp2713499p2713499.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Measure Difference Between Two Distributions

2010-09-25 Thread Rainer M Krug
On Sat, Sep 25, 2010 at 3:53 PM, Lorenzo Isella lorenzo.ise...@gmail.comwrote:

 On 09/25/2010 03:23 PM, Rainer M Krug wrote:


 Evaluate, for me, does not necessary mean test if they are
 significantly different, but rather to quantify the difference. If that
 is what you are looking for, you could look at the Earth Movers
 Distance, where a package is available at R-forge
 (https://r-forge.r-project.org/projects/earthmovdist/) which I co-wrote
 and used before.

 Cheers,

 Rainer


 Thanks Rainer. I had a quick look at wikipedia and the package you mention,
 and it seems what I am looking for.


Great - could you please give me some feedback after using the package, if
something could be improved? Thanks.


 Just a question about normalization of the distance calculated by the
 algorithm.
 Let us say that I have 4 distributions A,B,C,D coupled this way (A,B) and
 (C,D).
 The length of data in A is equal to the length of data in B, same applies
 to C and D but length(A)!=length(C).
 Now, the argument I would like to make is that A and B are more similar
 than C and D and show a couple of numbers to prove this.
 Bottom line: provided my data lists are long enough, does this distance
 scale with the number of data? and if they do, how should I normalize this
 distance to compare the results?


You could represent the distance as the proportion of maximum possible
distance, i.e. scaling it to be between 0 and 1.

An example:
A and B have the same length (x), and you calculate the emd(A, B), which is
d.
Now you have to determine the maximum distance between these two:
remembering the analogy of moving earth, the biggest distance between the
two distributions would be if in A, all elements would be in A(1) and all
other would be zero, and in B all elements would be zero, except of B(x).
Now you can calculate the difference between these two, and you get dmax
The last step is to divide d/dmax, i.e. scaling to a value between 0 and 1.

this value then can be compared with the same ratio obtained from C and D
with length y.

One important point to keep in mind when using the emd: if the sum(A) is not
the same as sum(B), emd(A,B) is NOT EQUAL to emd(B,A). If this applies to
your case, you have to decide what to do, but one option is to standardise A
and B so that their sum is the same (effectively comparing the SHAPES and
not the actual values.


If you need a reference where we used this approach (for comparison of
different maps from different areas), see :

@ARTICLE{Roura-Pascual2009_rmkc,
  author = {Roura-Pascual, N\'{u}ria and Krug, Rainer M. and Richardson,
David
M. and Hui, Cang},
  title = {Spatially-explicit sensitivity analysis for conservation
management:
exploring the influence of decisions in invasive alien plant management},
  journal = {Diversity and Distributions},
  year = {2010},
  volume = {16},
  pages = {426--438},
  doi = {10./j.1472-4642.2010.00659.x},
  file = {Article:Roura-Pascual2009_rmkc.pdf:PDF},
  owner = {rkrug},
  timestamp = {2009.03.11}
}

Please feel free to contact me if you have further questions,

Cheers,

Rainer



 Cheers

 Lorenzo




-- 
NEW GERMAN FAX NUMBER!!!

Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology,
UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Cell:   +27 - (0)83 9479 042
Fax:+27 - (0)86 516 2782
Fax:+49 - (0)321 2125 2244
email:  rai...@krugs.de

Skype:  RMkrug
Google: r.m.k...@gmail.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] OT: What distribution is this?

2010-09-25 Thread Rainer M Krug
Hi

This is OT, but I need it for my simulation in R.

I have a special case for sampling with replacement: instead of sampling
once and replacing it immediately, I sample n times, and then replace all n
items.


So:

N entities
x samples with replacement
each sample consists of n sub-samples WITHOUT replacement, which are all
replaced before the next sample is drawn

My question is: which distribution can I use to describe how often each
entity of the N has been sampled?

Thanks for your help,

Rainer

-- 
NEW GERMAN FAX NUMBER!!!

Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology,
UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Cell:   +27 - (0)83 9479 042
Fax:+27 - (0)86 516 2782
Fax:+49 - (0)321 2125 2244
email:  rai...@krugs.de

Skype:  RMkrug
Google: r.m.k...@gmail.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] margin control in lattice package

2010-09-25 Thread Jonathan Flowers
Hi all,

I am difficulty with simple layout of plots in the lattice package

I have created a series of levelplots and would like to plot them to a
single device, but need to reduce the margin areas.  This is easily
accomplished with par(oma) and par(mar) in the base graphics package but I
am having problems finding the equivalent features in the lattice package.
Ideally, I would like to reduce the amount of white space among plots in the
following example. Thanks in advance.

library(lattice)
p1 - levelplot( matrix(c(1:25),nr=5,nc=5),row.
values=1:5,column.values=1:5)
p2 - levelplot(matrix
(rnorm(25),nr=5,nc=5),row.values=1:5,column.values=1:5)
p3 - levelplot( matrix(c(1:25),nr=5,nc=5),row.values=1:5,column.values=1:5)
p4 - levelplot(matrix
(rnorm(25),nr=5,nc=5),row.values=1:5,column.values=1:5)

print(p1,split=c(1,1,2,2),more=T)
print(p2,split=c(2,1,2,2),more=T)
print(p3,split=c(1,2,2,2),more=T)
print(p4,split=c(2,2,2,2))

Thanks in advance,

Jonathan

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] OT: What distribution is this?

2010-09-25 Thread Berwin A Turlach
G'day Rainer,

On Sat, 25 Sep 2010 16:24:17 +0200
Rainer M Krug r.m.k...@gmail.com wrote:

 This is OT, but I need it for my simulation in R.
 
 I have a special case for sampling with replacement: instead of
 sampling once and replacing it immediately, I sample n times, and
 then replace all n items.
 
 
 So:
 
 N entities
 x samples with replacement
 each sample consists of n sub-samples WITHOUT replacement, which are
 all replaced before the next sample is drawn
 
 My question is: which distribution can I use to describe how often
 each entity of the N has been sampled?

Surely, unless I am missing something, any given entity would have
(marginally) a binomial distribution:

A sub-sample of size n either contains the entity or it does not.  The
probability that a sub-sample contains the entity is a function of N
and n alone.

x sub-samples are drawn (with replacement), so the number of times that
an entity has been sampled is the number of sub-samples in which it
appears.  This is given by the binomial distribution with parameters x
and p, where p is the probability determined in the previous paragraph.

I guess the fun starts if you try to determine the joint distribution
of two (or more) entities.

HTH.

Cheers,

Berwin 

== Full address 
Berwin A Turlach  Tel.: +61 (8) 6488 3338 (secr)
School of Maths and Stats (M019)+61 (8) 6488 3383 (self)
The University of Western Australia   FAX : +61 (8) 6488 1028
35 Stirling Highway   
Crawley WA 6009e-mail: ber...@maths.uwa.edu.au
Australiahttp://www.maths.uwa.edu.au/~berwin

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] OT: What distribution is this?

2010-09-25 Thread Peter Dalgaard
On 09/25/2010 04:24 PM, Rainer M Krug wrote:
 Hi
 
 This is OT, but I need it for my simulation in R.
 
 I have a special case for sampling with replacement: instead of sampling
 once and replacing it immediately, I sample n times, and then replace all n
 items.
 
 
 So:
 
 N entities
 x samples with replacement
 each sample consists of n sub-samples WITHOUT replacement, which are all
 replaced before the next sample is drawn
 
 My question is: which distribution can I use to describe how often each
 entity of the N has been sampled?
 
 Thanks for your help,
 
 Rainer
 

How did you know I was in the middle of preparing lectures on the
variance of the hypergeometric distribution and such? ;-)

If you look at a single item, the answer is of course that you have a
binomial with size=x and prob=n/N. The problem is that these binomials
are correlated between items.

If you can make do with a 2nd order approximation, then the covariances
between the indicators for two items being selected is easily found from
the symmetry and the fact that if you sum all N indicators you get the
constant n. I.e. the variance is p(1-p) and the covariance is
-p(1-p)/(N-1). For sums over repeated samples, just multiply everything
by the number x of samples.

If you intend to just count the frequency of a particular feature in
each of your n-samples, i.e., you have x replications of a
hypergeometric experiment, then you can do somewhat better by computing
the explicit convolution of x hypergeometrics (convolve(x, rev(y),
type=o) and Reduce() are your friends). I'm not sure this is actually
worth the trouble, but it should be doable for decent-sized N and x.



-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Measure Difference Between Two Distributions

2010-09-25 Thread Lorenzo Isella

ld represent the distance as the proportion of maximum possible

distance, i.e. scaling it to be between 0 and 1.

An example:
A and B have the same length (x), and you calculate the emd(A, B), which
is d.
Now you have to determine the maximum distance between these two:
remembering the analogy of moving earth, the biggest distance between
the two distributions would be if in A, all elements would be in A(1)
and all other would be zero, and in B all elements would be zero, except
of B(x). Now you can calculate the difference between these two, and you
get dmax
The last step is to divide d/dmax, i.e. scaling to a value between 0 and 1.

this value then can be compared with the same ratio obtained from C and
D with length y.

One important point to keep in mind when using the emd: if the sum(A) is
not the same as sum(B), emd(A,B) is NOT EQUAL to emd(B,A). If this
applies to your case, you have to decide what to do, but one option is
to standardise A and B so that their sum is the same (effectively
comparing the SHAPES and not the actual values.


OK, I see. The standardization part is not a terrible problem, I guess.
The other bit is less clear (to me). What are A(1) and B(x)? Am I piling 
up all the elements in A and B in a single bin?

Cheers

Lorenzo

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Uncertainty propagation

2010-09-25 Thread Ben Bolker
Maayt m.lupker at hotmail.com writes:

[snip]
 My first intention was to use a kind of monte carlo routine and run the
 model many times by changing the power law parameters. These power laws were
 obtained by fitting data points under R. I thus have std error associated to
 them: alpha (±da) * WaterHight ^ beta (±db). Is it statistically correct to
 sample alpha and beta for each run by picking them from a normal
 distribution centered on alpha (resp. beta) with a standard deviation of da
 (resp. db) and to perform my statistics (mean and standrad edviation of the
 model result) on the model output?
 It seems to me that da and db are correlated in some way and by doing what I
 entended to, I would overestimate the final error of my model...

  How have you fitted the models?
  Many of the fitting procedures in R give you access not just to
the standard errors of the parameters, but also to their correlations/
covariances.  If you have this information, you can sample the pairs
of parameters from an appropriate multivariate normal distribution.
Typically you could do something like ...

params - MASS::mvrnorm(1000,mu=coef(modelfit),Sigma=vcov(modelfit))
predictions - apply(params,1,predictfun)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Measure Difference Between Two Distributions

2010-09-25 Thread Rainer M Krug
On Sat, Sep 25, 2010 at 6:24 PM, Lorenzo Isella lorenzo.ise...@gmail.comwrote:

 ld represent the distance as the proportion of maximum possible

 distance, i.e. scaling it to be between 0 and 1.

 An example:
 A and B have the same length (x), and you calculate the emd(A, B), which
 is d.
 Now you have to determine the maximum distance between these two:
 remembering the analogy of moving earth, the biggest distance between
 the two distributions would be if in A, all elements would be in A(1)
 and all other would be zero, and in B all elements would be zero, except
 of B(x). Now you can calculate the difference between these two, and you
 get dmax
 The last step is to divide d/dmax, i.e. scaling to a value between 0 and
 1.

 this value then can be compared with the same ratio obtained from C and
 D with length y.

 One important point to keep in mind when using the emd: if the sum(A) is
 not the same as sum(B), emd(A,B) is NOT EQUAL to emd(B,A). If this
 applies to your case, you have to decide what to do, but one option is
 to standardise A and B so that their sum is the same (effectively
 comparing the SHAPES and not the actual values.


 OK, I see. The standardization part is not a terrible problem, I guess.
 The other bit is less clear (to me). What are A(1) and B(x)? Am I piling up
 all the elements in A and B in a single bin?
 Cheers


OK. Some code:

 set.seed(13)
 B - sample(1:10, 10)
 B
 [1]  8  3  4  1  6  7  9 10  2  5
 set.seed(13)
 A - sample(1:10, 10)
 B - sample(1:10, 10)
 A
 [1]  8  3  4  1  6  7  9 10  2  5
 B
 [1]  7  8  9  4 10  2  5  6  3  1
 A[1] - sum(A)
 A[-1] - 0
 B[length(B)] - sum(B)
 B[-length(B)] - 0
 A
 [1] 55  0  0  0  0  0  0  0  0  0
 B
 [1]  0  0  0  0  0  0  0  0  0 55

And now you can calculate the emd(A, B), which then is the maximum distance
between A and B. Imagine: the distance is the work you have to do to convert
A into B. Work equals distance times mass you have to move. Therefore you
have to maximise the distance you have to carry the earth and the amount you
have to carry. Therefore, in A, piling everything up in the first element,
and in B, piling everything up in the last element, gives you the most work
you have to du, which equals the largest distance.

Even though it is rather straight forward, I should probably integrate a
function in the package which gives you the largest distance between two
distributions - I'll think about it.

Hope this helps,

Cheers,

Rainer



 Lorenzo




-- 
NEW GERMAN FAX NUMBER!!!

Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology,
UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Cell:   +27 - (0)83 9479 042
Fax:+27 - (0)86 516 2782
Fax:+49 - (0)321 2125 2244
email:  rai...@krugs.de

Skype:  RMkrug
Google: r.m.k...@gmail.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Help required

2010-09-25 Thread Malik Shahzad

Is it possible to read jpeg files into R?

If yes please guide, Thanks.. I tried to search many time but failed to do.

Thankis in advance..

with Best Regards,
 
Malik Shahzad
Visiting Researcher
National Institute of Informatics (NII)
Tokyo, Japan
 
Doctoral Student
Asian Institute of Technology (AIT)
Bangkok, Thailand
+66-8-7676-5616



  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] apply union function vectorially

2010-09-25 Thread statquant2

Yes
Thank you very much, NICE ONE.
Do you have some interesting references about R (I do know some common books
but may be you have some nice hidden books I have never heard about...
)
-- 
View this message in context: 
http://r.789695.n4.nabble.com/apply-union-function-vectorially-tp2550162p2704852.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Measure Difference Between Two Distributions

2010-09-25 Thread Rainer M Krug
On Sat, Sep 25, 2010 at 3:16 PM, Dennis Murphy djmu...@gmail.com wrote:

 Hi:

 On Sat, Sep 25, 2010 at 3:47 AM, Lorenzo Isella lorenzo.ise...@gmail.com
 wrote:

  Dear All,
  Suppose you are given two distributions (or better: two equally-sized
 lists
  of data); how can you evaluate the difference between them?


Evaluate, for me, does not necessary mean test if they are significantly
different, but rather to quantify the difference. If that is what you are
looking for, you could look at the Earth Movers Distance, where a package
is available at R-forge (https://r-forge.r-project.org/projects/earthmovdist/)
which I co-wrote and used before.

Cheers,

Rainer


 I need something like an overlap measure of the two (let us say 0 == no
  overlap and 1== complete overlap). I should add that there is a 1-1
  correspondence of the data in the two distributions (they are ordered
 lists
  and e.g. the third element in the first distribution matches the third
  element in the second distribution).
 

 To visualize the two distributions, you could do an empirical Q-Q plot
 (qqplot(x, y)); if the distributions are identical, they should lie on a 45
 degree line - location shifts are indicated by level shifts (parallel
 lines)
 to the 45 degree line, differences in scale by slope differences away from
 1.

 There are many types of tests to test equality of distribution, the
 simplest
 one being the two-sample Kolmogorov-Smirnov test. It is sensitive to
 changes
 in both location and scale of the two edf's (empirical distribution
 functions). An alternative
 is the two-sample Cramer-von Mises test. I have no doubt there are
 others...


 HTH,
 Dennis

 The two distributions are not analytical (I mean they do not belong to any
  well known family).
  I wonder if mutual information could be what I am looking for.
  Cheers
 
  Lorenzo
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
NEW GERMAN FAX NUMBER!!!

Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology,
UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Cell:   +27 - (0)83 9479 042
Fax:+27 - (0)86 516 2782
Fax:+49 - (0)321 2125 2244
email:  rai...@krugs.de

Skype:  RMkrug
Google: r.m.k...@gmail.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Help required

2010-09-25 Thread Ben Bolker
Malik Shahzad geomalik at live.com writes:

 
 
 Is it possible to read jpeg files into R?
 
 If yes please guide, Thanks.. I tried to search many time but failed to do.
 

install.packages(sos)
library(sos)
findFn(read jpeg)

(I initially tried findFn(import jpeg) and didn't
get any hits, then tried findFn(jpeg) and scrolled
down to find these results)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Fitting GLMM models with glmer

2010-09-25 Thread Ben Bolker
Clodio Almeida clodioalm at gmail.com writes:

 
 Hi everybody:
 
[snip]

 
 I obtained the same results for parameters estimates and
 standarderrors, by using:
 
 glmer(cens ~ treat + heart + offset(log(SURVT) + (1 | INST),
 data=liver, family=poisson)

  It's nice that this check works.  I would like to see a lot
more NLMIXED/glmer comparison ...
 
 After that, the authors present the same model, but now bi = log(gi)
 where gi has the following gamma distribution:
 
 f(gi | theta1) = gi ^ (1/theta1 - 1) exp(-gi/theta1)/
 GAMMA(1/theta1)(theta1^(1/theta1)
 
 They used a set of transformations proposed in the paper, and with the
 following rotine fitted the model achieving
 estimates for the same fixed parameters and for theta1:
 
 proc nlmixed data=liver;
 parms theta1 = 1 b0 = -2:8 btrt = -0.54 bhrt =0:18;
 pi = CDF('NORMAL',ai);
 if(pi  0:99) then pi =0:99;
 gi2 = quantile('GAMMA',pi, 1/theta1);
 gi = theta1 * gi2;
 xb= log_s + b0 + btrt * treat + bhrt * heart+ log(gi);
 lambda=exp(xb);
 model cens ~ poisson(lambda);
 random ai ~ normal(0,1) subject=INST;
 run;
 
 This time I' m lost. Could anyone please give me a hint?
 
 The data set is:
 liver - as.data.frame(matrix(c(1,1,1,1,2,2,2,3,3,3,3,4,4,4,4,
 + 5,5,6,6,7,7,23.286,6.429,26.857,11.143,3.857,9.000,8.714,
 + 1.143,23.143,2.571,2.857,76.429,35.857,25.857,52.286,25.143,
 + 25.143,29.286,28.857,1.857,14.286,1,0,1,0,0,1,0,0,1,1,0,1,
 + 1,0,0,1,0,1,0,0,1,1,0,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0),
 + ncol=4,dimnames=list(NULL,c(INST,SURVT,treat,heart

  It's not immediately obvious that you can do this in glmer;
in fact, I suspect you can't. I don't know of an R package that
can fit non-normal random effects 'out of the box'; the only solutions
I know of other than SAS PROC NLMIXED are alternative programs
like WinBUGS and AD Model Builder (which do have R interfaces,
but have their own fairly significant learning curves).

  You might try asking this question on r-sig-mixed-models
instead.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Power Function

2010-09-25 Thread jethi

hey, how can i plot this function?

n=1000
m=2
k=n/m
N=100
myfun - function(n, m, alpha = .05, seeder = 1000) {
l=matrix(0,nrow=m,ncol=N)
for(i in 1:N){
set.seed(i)
for(j in 1:m){
x=rnorm(n,0,0.5)
y=rnorm(n,0,0.8)
l[j,i]=cor((x[(((j-1)*k)+1):(((j-1)*k)+k)]),
(y[(((j-1)*k)+1):(((j-1)*k)+k)]))
}
}

gute - function(x) {
q_1 - qnorm(alpha, 0, 0.05)
 
q_2 - qnorm(1 - alpha, 0, 0.05)
 
p=matrix(0,nrow=m,ncol=N)
H=matrix(0,nrow=N,ncol=1)
for(i in 1:N){
for (j in 1:m){
p[j,i]=x[j]^2/sum(x[,i]^2)
}
H[i]=log(m)-sum(p[,i]*log(p[,i]))
}
   1 - mean(q_1 = H  H = q_2)
}
output - gute(x = l)
 return(output)
}

regards
jethi
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Power-Function-tp2631929p2713799.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] 3D plot

2010-09-25 Thread jethi

hey, how can i plot this function??? thanks for ur help

n=1000
m=2
k=n/m
N=100
myfun - function(n, m, alpha = .05, seeder = 1000) {
l=matrix(0,nrow=m,ncol=N)
for(i in 1:N){
set.seed(i)
for(j in 1:m){
x=rnorm(n,0,0.5)
y=rnorm(n,0,0.8)
l[j,i]=cor((x[(((j-1)*k)+1):(((j-1)*k)+k)]),
(y[(((j-1)*k)+1):(((j-1)*k)+k)]))
}
}

for(i in 1:N){
for (j in 1:m){
gute - function() {
q_1 - qnorm(alpha, 0, 0.05)
 
q_2 - qnorm(1 - alpha, 0, 0.05)
 
p=matrix(0,nrow=m,ncol=N)
H=matrix(0,nrow=N,ncol=1)

p[j,i]=x[j]^2/sum(x[,i]^2)
}
H[i]=log(m)-sum(p[,i]*log(p[,i]))
}
   1 - mean(q_1 = H  H = q_2)
}
output - gute(a = l[,i])
 return(output)
}

regards 
jethi
-- 
View this message in context: 
http://r.789695.n4.nabble.com/3D-plot-tp2713818p2713818.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] 3D plot

2010-09-25 Thread Steve Lianoglou
Hi Jethi,

On Sat, Sep 25, 2010 at 3:42 PM, jethi kart...@hotmail.com wrote:

 hey, how can i plot this function??? thanks for ur help

I don't really have the time (or inclination, sorry) to look through
your code, but if it relates to your subject line, then there are
several ways to plot 3d objects in R.

You can look at the scatterplot3d package, or more interestingly (I
think) the rgl package, and its plot3d, lines3d, etc. functions.

-steve


 n=1000
 m=2
 k=n/m
 N=100
 myfun - function(n, m, alpha = .05, seeder = 1000) {
 l=matrix(0,nrow=m,ncol=N)
 for(i in 1:N){
 set.seed(i)
 for(j in 1:m){
 x=rnorm(n,0,0.5)
 y=rnorm(n,0,0.8)
 l[j,i]=cor((x[(((j-1)*k)+1):(((j-1)*k)+k)]),
 (y[(((j-1)*k)+1):(((j-1)*k)+k)]))
 }
 }

 for(i in 1:N){
 for (j in 1:m){
 gute - function() {
    q_1 - qnorm(alpha, 0, 0.05)

    q_2 - qnorm(1 - alpha, 0, 0.05)

 p=matrix(0,nrow=m,ncol=N)
 H=matrix(0,nrow=N,ncol=1)

 p[j,i]=x[j]^2/sum(x[,i]^2)
 }
 H[i]=log(m)-sum(p[,i]*log(p[,i]))
 }
   1 - mean(q_1 = H  H = q_2)
 }
 output - gute(a = l[,i])
  return(output)
 }

 regards
 jethi
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/3D-plot-tp2713818p2713818.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Steve Lianoglou
Graduate Student: Computational Systems Biology
 | Memorial Sloan-Kettering Cancer Center
 | Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] 3D plot

2010-09-25 Thread jethi

hi steve, thnx for ur reply. but my problem is, that i have a matrix as a
functions value and as a result i get a single number. so  how can i plot
this? and it would be nice if u read over my function, because i´m not sure
if it correct. anyway thnx
regards 
jethi
-- 
View this message in context: 
http://r.789695.n4.nabble.com/3D-plot-tp2713818p2713856.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Help required

2010-09-25 Thread Prof Brian Ripley

On Sat, 25 Sep 2010, Malik Shahzad wrote:



Is it possible to read jpeg files into R?

If yes please guide, Thanks.. I tried to search many time but failed to do.


On my system ??jpeg gave

ReadImages::read.jpeg   Read JPEG file
biOps::readJpeg Read jpeg file
rimage::read.jpeg   Read JPEG file

There are a number of other possibilities, including
EBImage::readImage from Bioconductor which uses ImageMagick to convert 
the file format (and that or many other image convertors could be used 
in a simple R wrapper).




Thankis in advance..

with Best Regards,

Malik Shahzad
Visiting Researcher
National Institute of Informatics (NII)
Tokyo, Japan

Doctoral Student
Asian Institute of Technology (AIT)
Bangkok, Thailand
+66-8-7676-5616




[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] (no subject)

2010-09-25 Thread Tal Galili
Did you look at:

?curve




Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--




On Sat, Sep 25, 2010 at 2:08 PM, Kartija Oompragash kart...@hotmail.comwrote:


 hi how can i plot now this function??? have to be m= 2???  because of the
 dimensions?thanks for ur help


 myfun - function(n, m, alpha = .05, seeder = 1000) {

  set.seed(seeder)

  x - matrix(rnorm(n, 0, 0.5), ncol = m)

  y - matrix(rnorm(n, 0, 0.8), ncol = m)

  l - diag(cor(x, y))

  cat(Correlations between two random variables \n, l, fill = TRUE)

  gute - function(x, m, alpha) {

q_1 - qnorm(alpha, 0, 0.05)

q_2 - qnorm(1 - alpha, 0, 0.05)

p - (x^2)/sum(x^2)

H - log(m) - sum(p * log(p), na.rm = TRUE)

1 - mean(q_1 = H  H = q_2)

  }

  dat - seq(0, 1, length.out = 10)

  output - gute(x = dat, m = m, alpha = alpha)

  return(output)

 }



[[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] margin control in lattice package

2010-09-25 Thread Peter Ehlers

On 2010-09-25 8:59, Jonathan Flowers wrote:

Hi all,

I am difficulty with simple layout of plots in the lattice package

I have created a series of levelplots and would like to plot them to a
single device, but need to reduce the margin areas.  This is easily
accomplished with par(oma) and par(mar) in the base graphics package but I
am having problems finding the equivalent features in the lattice package.
Ideally, I would like to reduce the amount of white space among plots in the
following example. Thanks in advance.

library(lattice)
p1- levelplot( matrix(c(1:25),nr=5,nc=5),row.
values=1:5,column.values=1:5)
p2- levelplot(matrix
(rnorm(25),nr=5,nc=5),row.values=1:5,column.values=1:5)
p3- levelplot( matrix(c(1:25),nr=5,nc=5),row.values=1:5,column.values=1:5)
p4- levelplot(matrix
(rnorm(25),nr=5,nc=5),row.values=1:5,column.values=1:5)

print(p1,split=c(1,1,2,2),more=T)
print(p2,split=c(2,1,2,2),more=T)
print(p3,split=c(1,2,2,2),more=T)
print(p4,split=c(2,2,2,2))

Thanks in advance,

Jonathan


Here are a couple of things you can play with.

First, the default for the matrix method of levelplot() is to
produce square plots using the argument aspect=iso. So
unless you set aspect=some other value or fill, you can
only reduce the outer white space at the expense of
increasing the inner white space. To reduce the outer white
space but keep the aspect=iso, you can set the size of
the graphics device and then fiddle with the top.padding
and/or bottom.padding components of the layout.heights
parameter. Something like this:

## some simple data
 m - matrix(1:25, nr=5)

## create 4 (identical for illustration only) plots
 p1 - p2 - p3 - p4 - levelplot(m, aspect=iso,
   par.settings=list(layout.heights=list(top.padding=-2)))

## open a trellis device (I'm on Windows)
 trellis.device(windows, height=6, width=7)

## print the plots
 print(p1, split=c(1,1,2,2), more=TRUE)
 print(p2, split=c(2,1,2,2), more=TRUE)
 print(p3, split=c(1,2,2,2), more=TRUE)
 print(p4, split=c(2,2,2,2))


## alternatively, use aspect=fill and adjust size in the
## print() calls
 p1 - p2 - p3 - p4 - levelplot(m, aspect=fill)

 trellis.device(windows, height=6, width=7)
 print(p1, split=c(1,1,2,2),
panel.height=list(x=2, units=in), more=TRUE)
 print(p2, split=c(2,1,2,2),
panel.height=list(x=2, units=in), more=TRUE)
 print(p3, split=c(1,2,2,2),
panel.height=list(x=2, units=in), more=TRUE)
 print(p4, split=c(2,2,2,2),
panel.height=list(x=2, units=in))


Probably the best way, if your levels are roughly the same
for all plots, is to convert your data to a data frame,
define a 4-level factor, and create a standard 4-panel plot
instead of using the 'split' argument.

  -Peter Ehlers

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Question on levels function and extracting the associated level number

2010-09-25 Thread Jeffrey Cexun Cai

 Dear Sir/Madam,

 I have a quick question, which I hope someone can help.

 I am using the levels function in R, which helps to summarize the number of
 factors that I have in a vector in an ordered manner.

 Using an already existing example:

  state - c(tas, sa,  qld, nsw, nsw, nt,  wa,  wa,
   qld, vic, nsw, vic, qld, qld, sa,  tas,
   sa,  nt,  wa,  vic, qld, nsw, nsw, wa,
   sa,  act, nsw, vic, vic, act)

  statef - factor(state)

  levels(statef)

  [1] act nsw nt  qld sa  tas vic wa


 I will now like to go through the entire vector state, and generate the
 corresponding level. Using the example, there are 8 levels, and so I will
 like to generate the corresponding vector (6, 5, 4, 2, 2, 3, ...1) which
 is the level number corresponding to the vector state.

 Any guidance/assistance will be much appreciated. Thank you so much!!

 sincerely,
 Jeff


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] 3D plot

2010-09-25 Thread jethi

hey,  is there anybody who can help me? its very urgent because i have to
send my bachelor thesis on monday. pls help me
-- 
View this message in context: 
http://r.789695.n4.nabble.com/3D-plot-tp2713818p2713909.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Question on levels function and extracting the associated level number

2010-09-25 Thread Jorge Ivan Velez
Jeffrey,

You can try

 as.numeric(state)
 [1] 6 5 4 2 2 3 8 8 4 7 2 7 4 4 5 6 5 3 8 7 4 2 2 8 5 1 2 7 7 1

HTH,
Jorge


On Sat, Sep 25, 2010 at 5:51 PM, Jeffrey Cexun Cai  wrote:

 
  Dear Sir/Madam,
 
  I have a quick question, which I hope someone can help.
 
  I am using the levels function in R, which helps to summarize the number
 of
  factors that I have in a vector in an ordered manner.
 
  Using an already existing example:
 
   state - c(tas, sa,  qld, nsw, nsw, nt,  wa,  wa,
qld, vic, nsw, vic, qld, qld, sa,  tas,
sa,  nt,  wa,  vic, qld, nsw, nsw, wa,
sa,  act, nsw, vic, vic, act)
 
   statef - factor(state)
 
   levels(statef)
 
   [1] act nsw nt  qld sa  tas vic wa
 
 
  I will now like to go through the entire vector state, and generate the
  corresponding level. Using the example, there are 8 levels, and so I will
  like to generate the corresponding vector (6, 5, 4, 2, 2, 3, ...1)
 which
  is the level number corresponding to the vector state.
 
  Any guidance/assistance will be much appreciated. Thank you so much!!
 
  sincerely,
  Jeff
 

[[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Finding Zeros of a Function

2010-09-25 Thread Lorenzo Isella

Dear All,
I need to find the (possible multiple) zeros of a function f within an 
interval. I gave uniroot a try, but it just returns one zero and I need 
to provide it with an interval [a,b] such that f(a)f(b)0.
Is there any function to find the multiple zeros of f in (a,b) without 
constraints on the sign of f(a) and f(b)?

Many thanks

Lorenzo

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] A little complex puzzle

2010-09-25 Thread Bill.Venables
I was looking for an example of complex variables in R.  This one is trivial, 
but rather cute (though World War II aficionados may 'come over all funny').

See if you can guess the image before you try the function.  It's not difficult.

jif - function(res = 100) {
  z - sample(do.call(complex, subset(expand.grid(real =
  seq(-3, 4, len = 7*res + 1), 
  imaginary =
  seq(-2, 2, len = 4*res + 1)),
  real  -2.439  real  3.717)))
  del - 2*base::pi/32
  plot(z, type = n, asp = 1, ann = FALSE, axes = FALSE)
  points(z[((Arg(z) + del/2) %/% del) %% 2 == 0
   | Mod(z)  1.15], col = red, pch = .)
}

Njoi.

Bill Venables.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Finding Zeros of a Function

2010-09-25 Thread Jorge Ivan Velez
Dear Lorenzo,

You could try the BB package:

# function we need the roots for

fn1 - function(x, a) - x^2 + x + a


 # plot

curve(fn1(x, a = 5), -4, 4)

abline(h=0, col = 2)


 # searching the roots in (-4, 4)

# install.packages('BB')

require(BB)

BBsolve(par = c(-1, 1), fn = fn1, a = 5)

values - BBsolve(par = c(-1, 1), fn = fn1, a = 5)$par

values


 # adding the roots to the plot

points(values, c(0, 0), pch = 8, cex = 1.1, col = 4)

HTH,
Jorge


On Sat, Sep 25, 2010 at 8:24 PM, Lorenzo Isella  wrote:

 Dear All,
 I need to find the (possible multiple) zeros of a function f within an
 interval. I gave uniroot a try, but it just returns one zero and I need to
 provide it with an interval [a,b] such that f(a)f(b)0.
 Is there any function to find the multiple zeros of f in (a,b) without
 constraints on the sign of f(a) and f(b)?
 Many thanks

 Lorenzo

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] why I could not reproduce the Mandelbrot plot demonstrated on R wiki

2010-09-25 Thread xin wei

hi, peter:
thank you for your attention. adding the line you suggested did display the
static Mandelbrot plot with good resolution on R graphics device. However,
the resulting gif file still come out ugly. the R wiki page I was referring
to is the following:
http://en.wikipedia.org/wiki/R_(programming_language)
where the nice Mandelbrot plot and sample codes are provided.

i would appreciate your help if you can provide further hint. 

thanks





Peter Dalgaard-2 wrote:
 
 On 09/24/2010 10:41 PM, xin wei wrote:
 
 I am trying to reproduce the nice looking of Mandelbrot demonstrated by R
 wiki page by the following code:
 
 library(caTools)# external package providing write.gif function
 jet.colors = colorRampPalette(c(#7F, blue, #007FFF, cyan,
 #7FFF7F, 
 yellow, #FF7F00, red, #7F)) 
 m = 600 # define size
 C = complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ), 
  imag=rep(seq(-1.2,1.2, length.out=m), m ) ) 
 C = matrix(C,m,m)   # reshape as square matrix of complex numbers
 Z = 0   # initialize Z to zero
 X = array(0, c(m,m,20)) # initialize output 3D array
 for (k in 1:20) {   # loop with 20 iterations
   Z = Z^2+C # the central difference equation  
   X[,,k] = exp(-abs(Z)) # capture results
 } 
 write.gif(X, Mandelbrot.gif, col=jet.colors, delay=100)
 
 Hmm, I couldn't be bothered with the caTools, but it looks fine for me
 with
 
 image(X[,,20],col=jet.colors(100))
 
 Perhaps you need jet.colors(n) as well?
 
 
 
 however, the gif file created by this looks much worse than what is shown
 on
 R wiki page, see the comparison as follows (left one is what i created)
 
 http://r.789695.n4.nabble.com/file/n2591429/Picture1.png 
 
 Save for the odd color scheme, the one on the left looks like a
 Mandelbrot set, the one on the right appears to be iteration 4. I
 couldn't find your original source for this on wiki.r-project.org?
 
 
 -- 
 Peter Dalgaard
 Center for Statistics, Copenhagen Business School
 Phone: (+45)38153501
 Email: pd@cbs.dk  Priv: pda...@gmail.com
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 

-- 
View this message in context: 
http://r.789695.n4.nabble.com/why-I-could-not-reproduce-the-Mandelbrot-plot-demonstrated-on-R-wiki-tp2591429p2714024.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] formatting data for predict()

2010-09-25 Thread Andrew Miles
I'm trying to get predicted probabilities out of a regression model,  
but am having trouble with the newdata option in the predict()  
function.  Suppose I have a model with two independent variables, like  
this:


y=rbinom(100, 1, .3)
x1=rbinom(100, 1, .5)
x2=rnorm(100, 3, 2)
mod=glm(y ~ x1 + x2, family=binomial)

I can then get the predicted probabilities for the two values of x1,  
holding x2 constant at 0 like this:


p2=predict(mod, type=response, newdata=as.data.frame(cbind(x1, x2=0)))
unique(p2)

However, I am running regressions as part of a function I wrote, which  
feeds in the independent variables to the regression in matrix form,  
like this:


dat=cbind(x1, x2)
mod2=glm(y ~ dat, family=binomial)

The results are the same as in mod.  Yet I cannot figure out how to  
input information into the newdata option of predict() in order to  
generate the same predicted probabilities as above.  The same code as  
above does not work:


p2a=predict(mod2, type=response, newdata=as.data.frame(cbind(x1,  
x2=0)))

unique(p2a)

Nor does creating a data frame that has the names datx1 and datx2,  
which is how the variables appear if you run a summary() on mod2.   
Looking at the model matrix of mod2 shows that the fitted model only  
shows two variables, the dependent variable y and one independent  
variable called dat.  It is as if my two variables x1 and x2 have  
become two levels in a factor variable called dat.


names(mod2$model)

My question is this:  if I have a fitted model like mod2, how do I use  
the newdata option in the predict function so that I can get the  
predicted values I am after?  I.E. how do I recreate a data frame with  
one variable called dat that contains two levels which represent my  
(modified) variables x1 and x2?


Thanks in advance!

Andrew Miles
Department of Sociology
Duke University

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Finding Zeros of a Function

2010-09-25 Thread Berend Hasselman

You may also try package nleqslv, which uses  Newton or Broyden to solve  a
system of nonlinear equations. It is an alternative to BB.
Without further information from your side, no additional information can be
provided.

/Berend
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Finding-Zeros-of-a-Function-tp2713980p2714073.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.