Re: [R] bptest
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
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
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
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
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
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
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
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
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)
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
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
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
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
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
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?
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
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?
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?
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
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
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
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
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
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
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
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
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
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
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
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
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
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)
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
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
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
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
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
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
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
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
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()
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
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.