Re: [R] Writing list object to a file
Arun Kumar Saha [EMAIL PROTECTED] wrote:
Hi all,
I am wondering how to write a 'list' object to a file. I already gone
through some threads like
http://mail.python.org/pipermail/python-list/2001-April/080639.html, however
could not trace out any reliable solution. I tried following :
write.table(calc, file=c:/data.csv)
Error in data.frame(200501 = c(-0.000387071806652095,
-0.000387221689252648, :
arguments imply differing number of rows: 25, 24, 16, 17, 18, 26, 27, 19,
23, 20, 11
Anybody can help?
Remember that a list is not a rectangular table. It may have
components of various types and shapes. Therefore, it's unlikely
you'll be able to write a file that can be read easily by
another program.
That said, we do this occasionally for users who want to use
spreadsheets (yuck!) to analyze part of a list. We use the
following function, which simply prints the list to an ASCII
file with a long line length:
##
# File: Robj2txt.r
# Language: R
# Programmer:Michael H. Prager
# Date: July 7, 2004
# Synopsis:
# Function to write a complex R object to an ASCII file.
# Main use is to write objects created from .rdat files;
# however, could be used to save other objects as well.
# This can be used to create files for those who want to use
# a spreadsheet or other program on the data.
###
Robj2txt - function(x, file = paste(deparse(substitute(x)),
.txt, sep = )) {
#
# ARGUMENTS:
# x: R data object to save as ascii
# filename: Name of file to save. Default is name of x
with .txt extension
#
tmp.wid = getOption(width) # save current width
options(width = 1)# increase output width
sink(file)# redirect output to file
print(x) # print the object
sink()# cancel redirection
options(width = tmp.wid) # restore linewidth
return(invisible(NULL)) # return (nothing) from
function
}
##
--
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Can I get rid of this for loop using apply?
Hey all,
The code below creates a partial dependence plot for the variable x1 in the
linear model y ~ x1 + x1^2 + x2.
I have noticed that the for loop in the code takes a long time to run if the
size of the data is increased. Is there a way to change the for loop into
an apply statement? The tricky part is that I need to change the values of
x1 in each step of the loop to give me the appropriate dataset to make
predictions on cbind(m[,-match(x1,names(m))],x1=a[1,i+1]) . If I try
and add the 1112 columns to the dataset a priori and use apply, the code
won't work because the predict function needs the column labeled x1. I
realize I could just grab the form of the linear function and use that
instead of predict(), but I don't want to do that because I want to make
this code applicable to generic model fits.
#create fake data and fit a simple linear regression model
x1 - rep(c(1,3,4),100)
x2 - rep(c(1:6),50)
y - 40+2*x1^.5 - 6*x2 + rnorm(100,0,2)
m - as.data.frame(cbind(y,x1,x2))
lm1 - lm(y~x1+I(x1^2)+x2,data=m)
#super small version of R code for partial dependence plot
a - rbind(c(0:)*(max(m$x1)-min(m$x1))/ +
min(m$x1),c(0:)*0-9)
for(i in c(0:))
{
a[2,i+1] -
mean(predict(lm1,cbind(m[,-match(x1,names(m))],x1=a[1,i+1])))
}
plot(a[1,],a[2,],xlab=x1,ylab=Response,type=l,main=Partial Dependence
Plot)
Many thanks,
Mike
[[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] Can I get rid of this for loop using apply?
The answer to my post is yes (which I just figured out).
Solution:
#super small version of R code for pd plot using apply
a - rbind(c(0:)*(max(m$x1)-min(m$x1))/ +
min(m$x1),c(0:)*0-9)
b - matrix(rep(c(0:)*(max(m$x1)-min(m$x1))/ + min(m$x1), nrow(m)),
nrow(m), 1112, byrow=T)
a[2,] - apply(b,2,FUN=function(x)
{mean(predict(lm1,cbind(m[,-match(x1,names(m))],x1=x))) })
plot(a[1,],a[2,],xlab=x1,ylab=Response,type=l,main=Partial Dependence
Plot)
Mike Dugas
[[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] Can I get rid of this for loop using apply?
Thanks for the help. That explains why my time testing showed no
difference. Is there any way to speed up the program? It is unbearably
slow if I increase the number of loops.
Mike
On Wed, Apr 23, 2008 at 6:23 PM, hadley wickham [EMAIL PROTECTED] wrote:
On Wed, Apr 23, 2008 at 4:23 PM, Mike Dugas [EMAIL PROTECTED] wrote:
The answer to my post is yes (which I just figured out).
Switching from for to apply isn't going to speed up your code. If you
carefully read the source code of apply, you'll see the guts of the
work is done by:
for (i in 1:d2) {
tmp - FUN(array(newX[, i], d.call, dn.call), ...)
if (!is.null(tmp))
ans[[i]] - tmp
}
i.e apply uses for internally. The reason to use apply instead of a
for loop is so that you can better express the intent of your
algorithm.
Hadley
--
http://had.co.nz/
[[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] Can I get rid of this for loop using apply?
Sure, I am creating a partial dependence plot (reference Friedman's stochastic gradient paper from, I want to say, 2001). The idea is to find the relationship between one of the predictors, say x1, and y by creating the following plot: take a random sample of actual data points, hold other predictors fixed (x2-xp), vary x1 across its range, create a string of predictions for each value of x1, repeat for all observations in sample, and finally average all the predictions for each value of x1. If you think about it, this plot solves Simpson's paradox under fairly mild conditions. The code I wrote does this using predict() which is useful for modeling approaches like GAMs. Mike On Wed, Apr 23, 2008 at 8:47 PM, hadley wickham [EMAIL PROTECTED] wrote: On Wed, Apr 23, 2008 at 7:31 PM, Mike Dugas [EMAIL PROTECTED] wrote: Thanks for the help. That explains why my time testing showed no difference. Is there any way to speed up the program? It is unbearably slow if I increase the number of loops. Could you explain exactly what you're trying to do with your code? It's a little hard to understand. Hadley -- http://had.co.nz/ [[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] graphics history
Duncan Murdoch [EMAIL PROTECTED] wrote: One thing I'd like to do, but didn't have time to implement before 2.7.0, is to have history set to some finite size, e.g. a default might be the last 3 or 10 plots. The problem with record=TRUE is that it keeps a record of all the plots, so memory use just increases and increases. Why not just startup another device with record=FALSE? I'd like to have recording always on, but I don't need an infinite history. But this isn't urgent enough to have prodded me into writing it before now. A finite size would be nice. I've been using this code in scripts: graphics.off() windows(record = TRUE) .SavedPlots - NULL Not exactly the same thing, but it limits memory use. Are there side effects that could bite me? -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 data with a colour scale
merca duria [EMAIL PROTECTED] wrote: I have a two dimensional data (X,Y) that I want to represent with different colours, I want to make a plot with a graduate color scale at right, or below Take a look at the levelplot function in the lattice package. require(lattice) ?levelplot -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Adding average into a matplot?
-Halcyon- [EMAIL PROTECTED] wrote: I have a matrix which is filled with simulation results for several years. Example of an output (7 years, 4 simulations): [...] My matplot gives me 4 lines, but I would like to add a line with the averages of each year for all simulations. Can anyone help me with this? Take a look at help pages of the following functions: mean, apply, abline. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 add background color of a 2D chart by quadrant
tom soyer [EMAIL PROTECTED] wrote: I have a 2D chart that is divided into four quadrants, I, II, III, IV: plot(1:10,ylim=c(0,10),xlim=c(0,10),type=n) abline(v=5,h=5) text(x=c(7.5,7.5,2.5,2.5),y=c(2.5,7.5,7.5,2.5),labels=c(I,II,III,IV)) I would like to fill each quadrant with a background color unique to the quadrant. Does anyone know how to do this in R? In response to a similar question no more than two weeks ago, I posted detailed code as an example. I expect it would be easy to modify my example to fit your question. If you search the group archives, you should be able to find it. The thread title was background color in scatterplots. MHP -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] background color in scatterplots
Georg Ehret [EMAIL PROTECTED] wrote:
Dear R community,I have a scatterplot with multiple vertical ablines. I
wish to color each interval between two ablines in a different color...
Could you please indicate me how to do this efficiently?
Thank you!
Georg.
Dear Georg,
Try this example:
# Demonstrate scatterplot w/ different color bands
# M.H.P. - March 2008
graphics.off()
# Generate dummy data make a blank plot:
a - runif(100)
b - 2 * runif(100)
plot(a, b, type = n, main = Demonstration)
# Set the values of the x-axis references:
vrefs - c(0.2, 0.4, 0.6)
# Get coordinates of plot
pdims - par()$usr
# Concatenate them to the vector of x_axis references:
vrefs2 - c(pdims[1], vrefs, pdims[2])
# Generate some weak colors:
cols - rainbow(n = length(vrefs2), s = 0.15)
# Add the colors to the plot:
for (i in 2:length(vrefs2)){
polygon(c(vrefs2[(i-1):i], vrefs2[i:(i-1)]),
c(rep(pdims[3],2), rep(pdims[4],2)), col = cols[i],
border = NA)
}
# Draw the reference lines, points, and box:
abline(v = vrefs)
points(a, b, lwd = 1.5)
box()
--
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Extracting single output data into a vector or matrix
Ayman, It is difficult to say without seeing some code, but your output seems to be not a list in the R sense but a collection of vectors, each of length 1. The best way to put the values into a vector probably is to assign them to the elements of the vector during your computations. Mike Prager Ayman Oweida [EMAIL PROTECTED] wrote: I've been struggling to do the following: After a lengthy computation, I receive an output along the lines of the list below. This list has 41 values and is not the end of my computations. I have another computation to do on the list below, but in this final computation the list is supposed to be a vector. I've tried to assign the list below to a data frame and then extract it, but not luck! Cleary, this is because each of the outputs in the list represents an individual data point that is not regarded as part of a matrix. Any help? I desperately need to be able to extract all the output data into a Vector so I can perform the final step of my computation. Thanks in advance. [1] 1.573233e-10 [1] 2.939187e-10 [1] 5.491124e-10 [1] 1.025877e-09 [1] 1.916591e-09 [1] 3.580663e-09 [1] 6.689559e-09 [1] 1.249774e-08 [1] 2.334885e-08 [1] 4.36214e-08 [1] 8.149551e-08 [1] 1.522537e-07 [1] 2.844473e-07 [1] 5.314175e-07 [1] 9.928186e-07 [1] 1.854829e-06 [1] 3.465277e-06 [1] 6.47399e-06 [1] 1.209501e-05 [1] 2.259645e-05 [1] 4.221572e-05 [1] 7.886935e-05 [1] 0.0001473473 [1] 0.0002752811 [1] 0.0005142927 [1] 0.0009608253 [1] 0.001795058 [1] 0.00335361 [1] 0.006265368 [1] 0.01171493 [1] 0.02188637 [1] 0.04088913 [1] 0.07639071 [1] 0.1427168 [1] 0.2666308 [1] 0.4981322 [1] 0.9306848 [1] 1.738748 [1] 3.248409 [1] 6.068827 [1] 11.33806 - [[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. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Constrained regression
Dear Carlos, One approach is to use structural equation modeling (SEM). Some SEM packages, such as LISREL, Mplus and Mx, allow inequality and nonlinear constraints. Phantom variables (Rindskopf, 1984) may be used to impose inequality constraints. Your model is basically: y = b0 + b1*b1*x1 + b2*b2*x2 +...+ bp*bp*xp + e 1 = b1*b1 + b2*b2 +...+ bp*bp Alternatively, you can set some condition bounds on the parameter estimates. Then you only have to impose the second constraint. Rindskopf, D. (1984). Using phantom and imaginary latent variables to parameterize constraints in linear structural models. Psychometrika, 49, 37-47. Regards, Mike -- - Mike W.L. Cheung Phone: (65) 6516-3702 Department of Psychology Fax: (65) 6773-1843 National University of Singapore http://courses.nus.edu.sg/course/psycwlm/internet/ - On Mon, Mar 3, 2008 at 11:52 AM, Carlos Alzola [EMAIL PROTECTED] wrote: Dear list members, I am trying to get information on how to fit a linear regression with constrained parameters. Specifically, I have 8 predictors , their coeffiecients should all be non-negative and add up to 1. I understand it is a quadratic programming problem but I have no experience in the subject. I searched the archives but the results were inconclusive. Could someone provide suggestions and references to the literature, please? Thank you very much. Carlos Carlos Alzola [EMAIL PROTECTED] (703) 242-6747 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] new to latex to pdf
Dear All,
I'm trying to teach myself latex along with the latex function in Hmisc
and have hit a roadblock that I can't seem to get around. I'd greatly
appreciate any pointers.
I'm running R 2.6.0 on Windows XP and have Miktex 2.7 installed.
I've reproduced the code below, taken from Frank Harrell's latexsummary
introduction. My question relates to getting a pdf version of the table
from the following code. The pdfs of the graphics (f1a and f1b)
generated by setpdf are fine. However, after a number of attempts using
different methods, I don't seem to be able to get a pdf of the table
from the s1 object (I see the right table in my previewer just fine).
I've tried texi2dvi (texfilename, pdf=T) but get only a series of errors
and an unformatted table in the pdf.
I'm sure I'm missing some fundamental concept here, but I'm afraid I'm
not seeing it. I'd appreciate if anyone could point me in the right
direction. ( I have no trouble writing my own simple latex code and
converting it to pdf using pdftex in miktex).
Thanks,
Mike Babyak
Department of Psychiatry and Behavioral Science
Duke University Medical Center
**
Here's the code:
library(Hmisc)
library(survival)
getHdata(pbc)
pbc-upData(pbc, moveUnits=TRUE,
labels=c(stage='Histologic Stage\nLudwig Criteria'))
kmsurv - function(S, times) {
f - survfit.km(factor(rep(1,nrow(S))), S)
tt - c(0, f$time)
ss - c(1, f$surv) # add first point to survival curve
approx(tt, ss, xout=times, method='constant', f=0)$y
}
describe.survival - function(y) {
km - kmsurv(y, c(2,5))
c('2 Year'=km[1], '5 Year'=km[2], 'Mean, y'=sum(y[,1])/sum(y[,2]))
}
S - with(pbc, Surv(fu.days/365.25, status))
s1 - summary(S ~ age + albumin + ascites + bili + drug + edema + chol,
fun=describe.survival, data=pbc)
for(w in 1:2) {
if(w==1) setpdf(f1a,sublines=1,h=5.25) else
setpdf(f1b,sublines=1,h=5)
plot(s1, which=if(w==1)1:2 else 3,
cex.labels=.7, cex.group.labels=.7*1.15, subtitles=T, main='',
pch=if(w==2) 16 else c('2','5'), # 16=solid circle
xlab=if(w==2)'Survival Time' else 'Survival Probability')
dev.off()
}
w - latex(s1, cdec=c(2,2,1), ctable=TRUE, caption='Survival')
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 a file created with Fortran
Ben Bolker [EMAIL PROTECTED] wrote: Dennis Fisher fisher at plessthan.com writes: I am trying to read a file written by Fortran. Several lines of the file are pasted below: Perhaps read.fwf is what you want? (fwf stands for fixed width format). You would have to work out the field widths, but it would seem to be pretty straightforward). A couple of points. First, since you know the format statement, perhaps you control the Fortran program. Then, it might be nicer to introduce whitespace between the data items, which would serve two purposes: making read.table() work on the data set and making it easier for humans to check the data file more easily. Second, you could look at read.fortran() -- a function that takes a lightly modified Fortran format specification as an argument. That seems even better for your purposes than read.fwf. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] numeric format
Rolf Turner [EMAIL PROTECTED] wrote: I have often wanted to suppress these row numbers and for that purpose wrote the following version of print.data.frame() [...] The ``srn'' argument means ``suppress row numbers''; [...] I once suggested to an R Core person that my version of print.data.frame() be adopted as the system version, but was politely declined. Rolf-- Clearly, and appropriately, R development is not a democratic process. Still, if a vote were held, I would support your version. I have also needed to suppress row names from time to time. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Font analysis?
Hi all, My initial searches on this yielded bupkis; I'm looking to access the pixel-by-pixel data (preferably as a matrix, 0 for black, 1 for white) of rendered letters(given a specific font, size, and dpi. Any suggestions as to how to achieve this? Mike -- Mike Lawrence Graduate Student, Department of Psychology, Dalhousie University Website: http://memetic.ca Public calendar subscribe link for iCal users: webcal://icalx.com/public/informavore/Public.ics The road to wisdom? Well, it's plain and simple to express: Err and err and err again, but less and less and less. - Piet Hein __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Conditionally incrementing a loop counter
Hi, I am trying a for loop from 1 to 100 by 1. However, if a condition does not get met, I want to throw away that iteration. So if my loop is for (i in 1:100) and i is say, 25 and the condition is not met then I don't want i to go up to 26. Is there a way to do that? I can't seem to manually adjust i because from what I understand, R creates 100 long vector and uses that to loops thru and I'm not sure how to get at the index of that vector. Any suggestions? Thanks in advance. Mike Jones Westat 1650 Research Blvd. RE401 Rockville, MD 20850 Ph: 240.314.2312 [[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] Conditionally incrementing a loop counter: Take 2
My apologies for not including a working example.
Here it is:
for (i in 1:10){
cat(initial i = ,i,\n)
x - runif(1)
if (x 0.7){
i - i-1
}
cat(second i = ,i,\n)
}
When I ran this i got what follows, so there were four cases where I
wanted the i not to increment.
initial i = 1
second i = 1
initial i = 2
second i = 1
initial i = 3
second i = 3
initial i = 4
second i = 3
initial i = 5
second i = 4
initial i = 6
second i = 6
initial i = 7
second i = 7
initial i = 8
second i = 7
initial i = 9
second i = 9
initial i = 10
second i = 10
-Original Message-
From: Mike Jones
Sent: Thursday, December 27, 2007 4:35 PM
To: '[EMAIL PROTECTED]'
Subject: Conditionally incrementing a loop counter
Hi,
I am trying a for loop from 1 to 10 by 1. However, if a condition does
not get met, I want to throw away that iteration. So if my loop is
for (i in 1:10)
and i is say, 4 and the condition is not met then I don't want i to go
up to 5. Is there a way to do that? I can't seem to manually adjust i
because from what I understand, R creates 10 long vector and uses that
to loops thru and I'm not sure how to get at the index of that
vector. Any suggestions? Thanks in advance.
Mike Jones
Westat
1650 Research Blvd. RE401
Rockville, MD 20850
Ph: 240.314.2312
[[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] Conditionally incrementing a loop counter: Take 2
Since I didn't want the i to increment in the loop when the condition is not
met, then in my example I wanted the loop to actually run 14 times instead of
the 10 since I wanted 4 of the iterations to be thrown away, or ignored. I
still haven't been able to figure this out. Going the while route doesn't
seem to work for me either.
nums - numeric(10)
i - 1
garbage - 0
while (i = 10){
x - runif(1)
cat(x = ,x,\n)
if (x 0.1){
nums[i] - x
i - i + 1
}
else{
garbage - garbage+1
}
cat(i = ,i,garbage = ,garbage,\n)
}
-Original Message-
From: Peter Dalgaard [mailto:[EMAIL PROTECTED]
Sent: Thursday, December 27, 2007 5:36 PM
To: Mike Jones
Cc: [EMAIL PROTECTED]
Subject: Re: [R] Conditionally incrementing a loop counter: Take 2
Mike Jones wrote:
My apologies for not including a working example.
Here it is:
for (i in 1:10){
cat(initial i = ,i,\n)
x - runif(1)
if (x 0.7){
i - i-1
}
cat(second i = ,i,\n)
}
When I ran this i got what follows, so there were four cases where I
wanted the i not to increment.
initial i = 1
second i = 1
initial i = 2
second i = 1
initial i = 3
second i = 3
initial i = 4
second i = 3
initial i = 5
second i = 4
initial i = 6
second i = 6
initial i = 7
second i = 7
initial i = 8
second i = 7
initial i = 9
second i = 9
initial i = 10
second i = 10
Is this the kind of effect you want?
x - runif(10)
cbind(x, 1:10, cumsum(x .7))
x
[1,] 0.384165631 1 1
[2,] 0.392715845 2 2
[3,] 0.895936431 3 2
[4,] 0.910242185 4 2
[5,] 0.689987301 5 3
[6,] 0.237071326 6 4
[7,] 0.225032680 7 5
[8,] 0.001856286 8 6
[9,] 0.392034868 9 7
[10,] 0.655076045 10 8
If you insist on using a loop, you need to separate the loop control
from the manipulation of i, as in (e.g.)
i - 0
for (j in 1:10){
i - i + 1
cat(initial i = ,i,\n)
x - runif(1)
if (x 0.7){
i - i-1
}
cat(second i = ,i,\n)
}
-Original Message-
From:Mike Jones
Sent:Thursday, December 27, 2007 4:35 PM
To: '[EMAIL PROTECTED]'
Subject: Conditionally incrementing a loop counter
Hi,
I am trying a for loop from 1 to 10 by 1. However, if a condition
does not get met, I want to throw away that iteration. So if my
loop is for (i in 1:10) and i is say, 4 and the condition is not met
then I don't want i to go up to 5. Is there a way to do that? I
can't seem to manually adjust i because from what I understand, R
creates 10 long vector and uses that to loops thru and I'm not sure
how to get at the index of that vector. Any suggestions? Thanks in
advance.
Mike Jones
Westat
1650 Research Blvd. RE401
Rockville, MD 20850
Ph: 240.314.2312
[[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.
--
O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Rating R Helpers
John Sorkin [EMAIL PROTECTED] wrote: I believe we need to know the following about packages: (1) Does the package do what it purports to do, i.e. are the results valid? (2) Have the results generated by the package been validate against some other statistical package, or hand-worked example? (3) Are the methods used in the soundly based? (4) Does the package documentation refer to referred papers or textbooks? (5) In addition to the principle result, does the package return ancillary values that allow for proper interpretation of the main result, (e.g. lm gives estimates of the betas and their SEs, but also generates residuals)?. (6) Is the package easy to use, i.e. do the parameters used when invoking the package chosen so as to allow the package to be flexible? (7) Are the error messages produced by the package helpful? (8) Does the package conform to standards of R coding and good programming principles in general? (9) Does the package interact will with the larger R environment, e.g. does it have a plot method etc.? (10) Is the package well documented internally, i.e. is the code easy to follow, are the comments in the code adequate? (11) Is the package well documented externally, i.e. through man pages and perhaps other documentation (e.g. MASS and its associated textbook)? In addition to package evaluation and reviews, we also need some plan for the future of R. Who will maintain, modify, and extend packages after the principle author, or authors, retire? Software is never done. Errors need to be corrected, programs need to be modified to accommodate changes in software and hardware. I have reasonable certainty that commercial software (e.g. SAS) will be available in 10-years (and that PROC MIXED will still be a part of SAS). I am far less sanguine about any number of R packages. John Interesting questions. Re, the future : LaTeX provides an example. The more complex packages tend to stop developing when the original programmer loses interest. Sometimes another person picks one up, but not frequently. I think, for example, of the many slide-preparation packages, each more complex than the next, that have come and gone during my relatively short (15 yr) professional use of LaTeX. At its root, this is a rather deep question: how open-source, largely volunteer-developed software can survive over the long term, while continuing to improve and maintain high standards. We are rather early in the history of free software development to know the answer. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] putting URLs in Latex
Edna Bell [EMAIL PROTECTED] wrote: Hi R Gurus! This is definitely off topic, but I thought I'd try: what is the way to put in url's into a Latex file, please? In future, you might want to post such questions in group comp.text.tex Mike -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R graph window
Duncan Murdoch [EMAIL PROTECTED] wrote: I am using R 2.6.0 on Windows XP. I am wondering if R can show multiple graph at the same graph window with different tabs at the bottom of the window, like in S-plus. Does anybody have experience on this? No, it can't. If you turn history recording on (see the History menu entry when the graphics window has the focus), then R will save the plots and PgUp/PgDn will scroll through them. To enable that, I usually begin an R script: graphics.off() .SavedPlots - NULL # Deletes any existing plot history windows(record = TRUE, ...) You might also be interested in the savePlot function, well described in the R help files. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] writing summary() to a text file
Federico Calboli [EMAIL PROTECTED] wrote: I would like to output the results of a function into a text file, legible as a such. The function produces a summary quite like: Take a look at the sink() function. Does that do what you need? -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Packages - a great resource, but hard to find the right one
hadley wickham [EMAIL PROTECTED] wrote: Which moves somewhat back towards my original suggestion of review articles. To me, an article which compared and contrasted four or five packages on a given topic would be much more useful than an article which reviewed only a single package. I think basing reviews around a specific topic/methodology would be more useful than basing them around a single package. I agree: Such articles would be welcome resources if published either in JSS or in R-News. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Packages - a great resource, but hard to find the right one.
John Sorkin [EMAIL PROTECTED] wrote: The multitude of packages is one of the great strengths of R. Unfortunately there is no (or at least I am not aware of) any single source that lists all available packages and gives a synopsis of what each package does. One can install and load packages one-by-one and look at the help pages to see what each package does, but this is at best an inefficient and a worst a very frustrating task. Might there be a way to put together a searchable database that will allow a user to easily search for a given function or technique in all contributed packages? Besides the excellent answers already given, don't overlook Google. Searching on r statistics box-cox transform turns up a reference to MASS as the third entry. When programming in any language, I now find it quicker to search for syntax (and other) help by Googling than to pull the reference manual off the shelf or start up an online help file. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Aggregate with non-scalar function
R-Helpers,
I'm sorry to have to ask this -- I've not used R very much in the last
8 or 10 months, and I've gotten rusty.
I have the following (ff2 is a subset of a much, much larger dataset):
ff2
hostName user sys idle obsTime
10142 fred 0.4 0.5 98.0 2007-11-01 02:02:18
16886 barney 0.5 0.2 94.6 2007-10-25 19:12:12
8795 fred 0.0 0.1 99.8 2007-10-30 05:08:22
5261 fred 0.1 0.2 99.7 2007-10-25 07:20:32
12427 barney 0.1 0.2 93.2 2007-10-19 14:34:10
18067 barney 0.1 0.2 99.4 2007-10-27 10:34:08
973 fred 0.0 0.2 99.8 2007-10-19 08:24:22
5426 fred 0.2 0.3 99.5 2007-10-25 12:50:33
7067 fred 0.1 0.2 99.4 2007-10-27 19:32:27
13159 barney 0.1 0.4 84.3 2007-10-20 14:58:11
17481 barney 1.2 2.0 92.6 2007-10-26 15:02:11
21632 barney 0.1 0.1 99.6 2007-11-01 09:24:09
206 fred 19.4 4.8 53.7 2007-10-18 06:50:34
18151 barney 0.1 0.2 94.9 2007-10-27 13:22:09
10662 fred 0.1 0.2 99.6 2007-11-01 19:22:27
10376 fred 0.0 0.2 99.7 2007-11-01 09:50:24
3630 fred 43.7 7.0 33.0 2007-10-23 00:58:27
1118 fred 0.6 0.4 98.9 2007-10-19 13:14:23
5122 fred 0.1 0.2 99.6 2007-10-25 02:42:21
22117 barney 0.0 0.2 99.4 2007-11-02 01:34:04
doit(ff2)
hostName hour user.mean sys.mean idle.mean user.max sys.max idle.max
1barney 01 0.00 0.20 99.40 0.0 0.2 99.4
2barney 09 0.10 0.10 99.60 0.1 0.1 99.6
3barney 10 0.10 0.20 99.40 0.1 0.2 99.4
4barney 13 0.10 0.20 94.90 0.1 0.2 94.9
5barney 14 0.10 0.30 88.75 0.1 0.4 93.2
6barney 15 1.20 2.00 92.60 1.2 2.0 92.6
7barney 19 0.50 0.20 94.60 0.5 0.2 94.6
8 fred 00 43.70 7.00 33.00 43.7 7.0 33.0
9 fred 02 0.25 0.35 98.80 0.4 0.5 99.6
10 fred 05 0.00 0.10 99.80 0.0 0.1 99.8
11 fred 06 19.40 4.80 53.70 19.4 4.8 53.7
12 fred 07 0.10 0.20 99.70 0.1 0.2 99.7
13 fred 08 0.00 0.20 99.80 0.0 0.2 99.8
14 fred 09 0.00 0.20 99.70 0.0 0.2 99.7
15 fred 12 0.20 0.30 99.50 0.2 0.3 99.5
16 fred 13 0.60 0.40 98.90 0.6 0.4 98.9
17 fred 19 0.10 0.20 99.50 0.1 0.2 99.6
doit
function(x){
x.mean-aggregate(x[,c(user,sys,idle)],
by=list(hostName=x$hostName,
hour=strftime(as.POSIXlt(x$obsTime),%H)),
mean)
x.max-aggregate(x[,c(user,sys,idle)],
by=list(hostName=x$hostName,
hour=strftime(as.POSIXlt(x$obsTime),%H)),
max)
t1-merge(x.mean,x.max,by=c(hostName,hour),suffixes=c(.mean,.max))
return(t1)
}
The point of the doit function is to make a new dataframe in which
the columns are summary statistics of certain columns in the argument.
Is there a function similar to:
magic.function(ff2[,c(user,system,idle)],
by=list(hostName=ff2$hostName,hour=strftime(as.POSIXlt(ff2$obsTime),%H)),
function(x){c(mean.user=mean(x$user),
mean.system=mean(x$system),
mean.idle=mean(x$idle),
max.user=max(x$user),
max.system=max(x$system),
max.idle=max(x$idle))})
ie. an aggregate that can cope with a non-scalar function and do
what I mean? My doit function gets more and more ugly the more
summary statistics I add, and I worry about the merge with hundreds
of thousands of rows.
I'm almost sure I've seen a solution to what I know is a simple
problem, but I guess my search skills are as bad as my R: I've
rummaged around the r-help archives and came up with nothing to show
for it.
Pointers would be gratefully received.
Many thanks.
--
Regards,
Mike Nielsen
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 for Beginner!!
Leandre Bassole [EMAIL PROTECTED] wrote: I am a new user of R. I am very familar to Stata, but few days ago I have decided to switch to R. But R langage is very difficult.I really want to know the best way to learn this famous and interesting software. I agree with the suggestions made so far. Also, consider buying or borrowing a copy of Introductory Statistics with R by Peter Dalgaard. It is well written and brief. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Tricky vectorization problem
Had a realization this morning that I think achieves ceiling
performance and thought I'd share with the list. I was previously
generating sample data per participant and then calculating a mean,
but of course this could be simplified by simply getting a single
value from the sampling distribution of the mean for that
participant. This speeds things up immensely (from 10.5 seconds to .5
seconds) and should be sufficient for my needs, but I'd still welcome
further improvement suggestions.
#start a timer
start = proc.time()[1]
#set the number of monte carlo experiments
mce = 1e3
#set the true correllation
rho = .5
#set the number of Subjects
Ns = 100
#for each subject, set a number of observations in A
a.No = 1:100
#for each subject, set a number of observations in B
b.No = 1:100
#set the between Ss variance in condition A
v.a = 1
#set the between Ss variance in condition B
v.b = 2
#for each subject, set a standard deviation in A
s.a.w = 1:100
#for each subject, set a standard deviation in B
s.b.w = 1:100
#set up a collector for the simulated correlations
sim.r = vector(length=mce)
#define a covariance matrix for use in generating the correlated data
Sigma=matrix(c(v.a,sqrt(v.a*v.b)*rho,sqrt(v.a*v.b)*rho,v.b),2,2)
eS - eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
ev - eS$values
fact - eS$vectors %*% diag(sqrt(pmax(ev, 0)), 2)
#set up a collector for subject means in A
a = vector(length=Ns)
#set up a collector for subject means in B
b = vector(length=Ns)
#generate correlated ideal means for for each subject in each condition
X - rnorm(2 * Ns * mce)
dim(X) - c(2, Ns * mce)
sub.means - t(fact %*% X)
#generate observed means for each Subject in A
a=sub.means[,1]+rnorm(Ns*mce,0,s.a.w/sqrt(a.No))
#generate observed means for each Subject in B
b=sub.means[,2]+rnorm(Ns*mce,0,s.b.w/sqrt(b.No))
#Get the observed correlation for each monte carlo experiment
for(i in 1:mce){
sim.r[i] = cor(
a[(i*Ns-Ns+1):(i*Ns)]
,b[(i*Ns-Ns+1):(i*Ns)]
)
}
#show the total time this took
print(proc.time()[1]-start)
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Tricky vectorization problem
Posting this for posterity and to demonstrate that a speed-up of 2
orders of magnitude is indeed possible! I can now run 1e5 monte carlo
experiments in the time the old code could only run 1e3. The final
change was to replace my use of cor() with .Internal(cor()), which
avoids some checks that are unnecessary in this case:
#start a timer
start = proc.time()[1]
#set the number of monte carlo experiments
mce = 1e5
#set the true correllation
rho = 1
#set the number of Subjects
Ns = 100
#for each subject, set a number of observations in A
a.No = 1:100
#for each subject, set a number of observations in B
b.No = 1:100
#set the between Ss variance in condition A
v.a = 1
#set the between Ss variance in condition B
v.b = 2
#for each subject, set a standard deviation in A
s.a.w = 1:100
#for each subject, set a standard deviation in B
s.b.w = 1:100
#set up a collector for the simulated correlations
sim.r = vector(length=mce)
#define a covariance matrix for use in generating the correlated data
Sigma=matrix(c(v.a,sqrt(v.a*v.b)*rho,sqrt(v.a*v.b)*rho,v.b),2,2)
eS - eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
ev - eS$values
fact - eS$vectors %*% diag(sqrt(pmax(ev, 0)), 2)
#set up a collector for subject means in A
a = vector(length=Ns)
#set up a collector for subject means in B
b = vector(length=Ns)
#generate correlated ideal means for for each subject in each condition
X - rnorm(2 * Ns * mce)
dim(X) - c(2, Ns * mce)
sub.means - t(fact %*% X)
a=sub.means[,1]+rnorm(Ns*mce,0,s.a.w/sqrt(a.No))
b=sub.means[,2]+rnorm(Ns*mce,0,s.b.w/sqrt(b.No))
for(i in 1:mce){
end=i*Ns
sim.r[i] = .Internal(
cor(
a[(end-Ns+1):end]
,b[(end-Ns+1):end]
,TRUE
,FALSE
)
)
}
#show the total time this took
print(proc.time()[1]-start)
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Tricky vectorization problem
Hi all,
I'm using the code below within a loop that I run thousands of times
and even with the super-computing resources at my disposal this is
just too slow. The snippet below takes about 10s on my machines,
which is an order of magnitude or two slower than would be
preferable; in the end I'd like to set the number of monte carlo
experiments to 1e4 or even 1e5 to ensure stable results. I've tried
my hand at vectorizing it myself but I'm finding it tricky. Any help
would be greatly appreciated!
This code attempts to find the distribution of observed correlation
values between A B when each are measured imperfectly:
start = proc.time()[1] #start a timer
rho = .5 #set the true correllation
Ns = 100 #set the number of Subjects
a.No = 1:100 #for each subject, set a number of observations in A
(1:100 chosen arbitrarily for demonstration)
b.No = 1:100 #for each subject, set a number of observations in B
v.a = 1 #set the between Ss variance in condition A
v.b = 2 #set the between Ss variance in condition B
s.a.w = 1:100 #for each subject, set a standard deviation in A (1:100
chosen arbitrarily for demonstration)
s.b.w = 1:100 #for each subject, set a standard deviation in B
mce = 1e3 #set the number of monte carlo experiments
sim.r = vector(length=mce) #set up a collector for the simulated
correlations
Sigma=matrix(c(v.a,sqrt(v.a*v.b)*rho,sqrt(v.a*v.b)*rho,v.b),2,2)
#define a covariance matrix for use in generating the correlated data
a = vector(length=Ns) #set up a collector for subject means in A
b = vector(length=Ns) #set up a collector for subject means in B
for(i in 1:mce){ #start a monte carlo experiment
sub.means=mvrnorm(Ns,rep(0,2),Sigma) #generate correlated ideal
means for for each subject in each condition
for(s in 1:Ns){ #for each subject
a[s] = mean(rnorm(a.No[s],sub.means[s,1],s.a.w[s])) #generate
some
data for A and grab the mean
b[s] = mean(rnorm(b.No[s],sub.means[s,2],s.b.w[s])) #generate
some
data for B and grab the mean
}
sim.r[i] = cor(a,b) #store the observed correlation between A and B
}
print(proc.time()[1]-start) #show the total time this took
--
Mike Lawrence
Graduate Student, Department of Psychology, Dalhousie University
Website: http://memetic.ca
Public calendar: http://icalx.com/public/informavore/Public
The road to wisdom? Well, it's plain and simple to express:
Err and err and err again, but less and less and less.
- Piet Hein
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Tricky vectorization problem
Seems there were some line break issues when pasting the code, trying
again with a different commenting style:
#start a timer
start = proc.time()[1]
#set the true correllation
rho = .5
#set the number of Subjects
Ns = 100
#for each subject, set a number of observations in A
a.No = 1:100
#for each subject, set a number of observations in B
b.No = 1:100
#set the between Ss variance in condition A
v.a = 1
#set the between Ss variance in condition B
v.b = 2
#for each subject, set a standard deviation in A
s.a.w = 1:100
#for each subject, set a standard deviation in B
s.b.w = 1:100
#set the number of monte carlo experiments
mce = 1e3
#set up a collector for the simulated correlations
sim.r = vector(length=mce)
#define a covariance matrix for use in generating the correlated data
Sigma=matrix(c(v.a,sqrt(v.a*v.b)*rho,sqrt(v.a*v.b)*rho,v.b),2,2)
#set up a collector for subject means in A
a = vector(length=Ns)
#set up a collector for subject means in B
b = vector(length=Ns)
#start a monte carlo experiment
for(i in 1:mce){
#generate correlated ideal means for for each subject in each condition
sub.means=mvrnorm(Ns,rep(0,2),Sigma)
#for each subject
for(s in 1:Ns){
#generate some data for A and grab the mean
a[s] = mean(rnorm(a.No[s],sub.means[s,1],s.a.w[s]))
#generate some data for B and grab the mean
b[s] = mean(rnorm(b.No[s],sub.means[s,2],s.b.w[s]))
}
#store the observed correlation between A and B
sim.r[i] = cor(a,b)
}
#show the total time this took
print(proc.time()[1]-start)
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] optimize() stuck in local plateau ?
I may be misunderstanding, but my example seems to violate the local
minimum inside [x_1,x_2] will be found as solution, even when f is
constant in there rule since the search in my example continues on
towards +1.
On 1-Oct-07, at 10:40 AM, Ravi Varadhan wrote:
Please read the help for optimize() carefully. The following
excerpted from
there should help explain your problem:
The first evaluation of f is always at x_1 = a + (1-phi)(b-a)
where (a,b) =
(lower, upper) and phi = (sqrt 5 - 1)/2 = 0.61803.. is the golden
section
ratio. Almost always, the second evaluation is at x_2 = a + phi(b-
a). Note
that a local minimum inside [x_1,x_2] will be found as solution,
even when f
is constant in there, see the last example.
In your case,
x_1 = a + (1-phi)*(b-a)
x_1
[1] -0.236068
x_2 = a + phi*(b-a)
x_2
[1] 0.236068
Since your function is constant in (x_1, x_2), you get your
solution in that
interval.
Try a different interval, and you'll get your answer:
optimize(my.func,interval=c(-1,0),maximum=TRUE)
[1] -0.618034 0.618034
[1] -0.381966 0.00
[1] -0.763932 0.763932
[1] -0.8090170 0.4045085
[1] -0.7016261 0.7016261
[1] -0.7401333 0.7401333
[1] -0.781153 0.781153
[1] -0.791796 0.791796
[1] -0.7983739 0.7983739
[1] -0.8024392 0.4012196
[1] -0.7958614 0.7958614
[1] -0.7999267 0.7999267
[1] -0.8008864 0.4004432
[1] -0.7993336 0.7993336
[1] -0.8002933 0.4001466
[1] -0.7997001 0.7997001
[1] -0.8000667 0.4000334
[1] -0.7998402 0.7998402
[1] -0.7999802 0.7999802
[1] -0.8000209 0.4000104
[1] -0.7999802 0.7999802
$maximum
[1] -0.7999802
$objective
[1] 0.7999802
Best,
Ravi.
--
--
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/
Varadhan.html
--
--
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
project.org] On
Behalf Of Mike Lawrence
Sent: Monday, October 01, 2007 1:29 AM
To: Rhelp
Subject: [R] optimize() stuck in local plateau ?
Hi all,
Consider the following function:
my.func = function(x){
y=ifelse(x-.5,0,ifelse(x -.8,abs(x)/2,abs(x)))
print(c(x,y)) #print what was tested and what the result is
return(y)
}
curve(my.func,from=-1,1)
When I attempt to find the maximum of this function, which should be
-.8, I find that optimize gets stuck in the plateau area and doesn't
bother testing the more interesting bits of the function:
optimize(my.func,interval=c(-1,1),maximum=TRUE)
I really don't understand why the search moves to the positive/
constant area of the function and neglects the more negative area of
the function. On step #4, after finding that there is no difference
between tests at -.23, .23 .52, shouldn't the algorithm try -.52?
In fact, it seems to me that it would make sense to try -.52 on step
3, so that we've tested one negative, one positive (found no
difference), now one negative again. Thoughts?
Of course I could define my interval more reasonably for this
particular function, but this is in fact simply one of a class of
functions I'm exploring, none of which have known formal descriptions
as above (I'm exploring a large number of 'black boxes'). I do know
that the maximum must occur between -1 and 1 for all however. Please
advise on how I might use optimize more usefully.
Mike
--
Mike Lawrence
Graduate Student, Department of Psychology, Dalhousie University
Website: http://memetic.ca
Public calendar: http://icalx.com/public/informavore/Public
The road to wisdom? Well, it's plain and simple to express:
Err and err and err again, but less and less and less.
- Piet Hein
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-
guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Mike Lawrence
Graduate Student, Department of Psychology, Dalhousie University
Website: http://memetic.ca
Public calendar: http://icalx.com/public/informavore/Public
The road to wisdom? Well, it's plain and simple to express:
Err and err and err again, but less and less and less.
- Piet Hein
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] optimize() stuck in local plateau ?
Hi all,
Consider the following function:
my.func = function(x){
y=ifelse(x-.5,0,ifelse(x -.8,abs(x)/2,abs(x)))
print(c(x,y)) #print what was tested and what the result is
return(y)
}
curve(my.func,from=-1,1)
When I attempt to find the maximum of this function, which should be
-.8, I find that optimize gets stuck in the plateau area and doesn't
bother testing the more interesting bits of the function:
optimize(my.func,interval=c(-1,1),maximum=TRUE)
I really don't understand why the search moves to the positive/
constant area of the function and neglects the more negative area of
the function. On step #4, after finding that there is no difference
between tests at -.23, .23 .52, shouldn't the algorithm try -.52?
In fact, it seems to me that it would make sense to try -.52 on step
3, so that we've tested one negative, one positive (found no
difference), now one negative again. Thoughts?
Of course I could define my interval more reasonably for this
particular function, but this is in fact simply one of a class of
functions I'm exploring, none of which have known formal descriptions
as above (I'm exploring a large number of 'black boxes'). I do know
that the maximum must occur between -1 and 1 for all however. Please
advise on how I might use optimize more usefully.
Mike
--
Mike Lawrence
Graduate Student, Department of Psychology, Dalhousie University
Website: http://memetic.ca
Public calendar: http://icalx.com/public/informavore/Public
The road to wisdom? Well, it's plain and simple to express:
Err and err and err again, but less and less and less.
- Piet Hein
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 find p in binomial(n,p)
I think the function you need is 'help.search'; try: help.search(binomial) and look for something obvious in the 'stats' package. A good deal quicker and easier than posting to an internet forum! Cheers, Mike. cathelf wrote: Dear all, I am trying a find the value p in binomial. X ~ Bin(n,p) I want to find the value p, so that Pr(X = k) = alpha Here, n, k and alpha are known. n, k are integers. alpha is between (0,1). Thanks a lot! Catherine -- View this message in context: http://www.nabble.com/how-to-find-%22p%22-in-binomial%28n%2Cp%29-tf4484227.html#a12787900 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.
