Re: [R] Can DEoptim trace output be customized?
Dear David, The package doesn't have an option to customize the output of the trace. However, you can create a custom version of the package that doesn't print the parameters. Get the package source code, uncompress it, and find the file de4_0.c in the src/ directory. Then comment out the calls to Rprintf that print the parameter values, to read: Rprintf(Iteration: %d bestvalit: %f bestmemit:, i_iter, t_bestC); // for (j = 0; j i_D; j++) // Rprintf(%12.6f, gt_bestP[j]); Rprintf(\n); Then re-build and re-install the package (instructions are at http://cran.r-project.org/doc/manuals/r-release/R-exts.html#Checking-and-building-packages ). On Thu, May 9, 2013 at 11:27 AM, David Reiner david.rei...@xrtrading.comwrote: I'm running DEoptim - it works great! (A HUGE THANK YOU to David Ardia, Katharine Mullen, Brian Peterson, and Joshua Ulrich, and Kris Boudt!!!). Sometimes I set trace to a number so I can see a few intermediate points in the optimization. However, I have a large number of variables and would like to omit the bestmemit, which hides what I'm more interested in: bestvalit. Is there a way to customize the output of trace so it omits the best member? Thanks, David L. Reiner This e-mail and any materials attached hereto, including, without limitation, all content hereof and thereof (collectively, XR Content) are confidential and proprietary to XR Trading, LLC (XR) and/or its affiliates, and are protected by intellectual property laws. Without the prior written consent of XR, the XR Content may not (i) be disclosed to any third party or (ii) be reproduced or otherwise used by anyone other than current employees of XR or its affiliates, on behalf of XR or its affiliates. THE XR CONTENT IS PROVIDED AS IS, WITHOUT REPRESENTATIONS OR WARRANTIES OF ANY KIND. TO THE MAXIMUM EXTENT PERMISSIBLE UNDER APPLICABLE LAW, XR HEREBY DISCLAIMS ANY AND ALL WARRANTIES, EXPRESS AND IMPLIED, RELATING TO THE XR CONTENT, AND NEITHER XR NOR ANY OF ITS AFFILIATES SHALL IN ANY EVENT BE LIABLE FOR ANY DAMAGES OF ANY NATURE WHATSOEVER, INCLUDING, BUT NOT LIMITED TO, DIRECT, INDIRECT, CONSEQUENTIAL, SPECIAL AND PUNITIVE DAMAGES, LOSS OF PROFITS AND TRADING LOSSES, RESULTING FROM ANY PERSON'S USE OR RELIANCE UPON, OR INABILITY TO USE, ANY XR CONTENT, EVEN IF XR IS ADVISED OF THE POSSIBILITY OF SUCH DAMAGES OR IF SUCH DAMAGES WERE FORESEEABLE. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R crashes when loading rgl package before minqa package
I also cannot reproduce the crash.
sessionInfo()
R version 2.11.1 (2010-05-31)
x86_64-unknown-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] minqa_1.1.9 Rcpp_0.8.6 rgl_0.91
On Thu, 30 Sep 2010, Ravi Varadhan wrote:
No. It still does not crash in Windows.
library(rgl)
library(minqa)
Loading required package: Rcpp
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Thursday, September 30, 2010 11:43 am
Subject: Re: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
On Thu, 30 Sep 2010, Ravi Varadhan wrote:
I get this on Windows (it does not crash):
library(minqa)
library(rgl)
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
Does it crash when you load first rgl and then only minqa? Like this:
library(rgl)
library(minqa)
newuoa(initpar, optimft)
/Gaspard
This tells me that you should be constraining your parameter x[4]
(may be even x[5]) to be non-negative:
Here is what I get with `bobyqa':
bobyqa(initpar, optimft, lower=c(-Inf, -Inf, -Inf, 0, 0))
parameter estimates: -5.311767080681, -3861.89005072333,
979.239647766226, 0.268156271922112, 27.6418856936228
objective: 1457.20987728737
number of function evaluations: 78
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Wednesday, September 29, 2010 11:40 am
Subject: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
Calling newuoa (from the minqa package) makes R crash when the
package rgl is loaded first. This however only on certain selected
data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when
loading qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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-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,
[R] Registration deadline, useR! 2010
The final registration deadline for the R User Conference is June 20, 2010, one week away. Later registration will not be possible on site! Conference webpage: http://www.R-project.org/useR-2010 Conference program: http://www.R-project.org/useR-2010/program.html Registration: http://www.R-project.org/useR-2010/registration/registration.html The conference is scheduled for July 21-23, 2010, and will take place at the campus of the National Institute of Standards and Technology (NIST) in Gaithersburg, Maryland, USA. Following the successful useR! 2004, useR! 2006, useR! 2007, useR! 2008, and useR! 2009, conferences, the conference is focused on: 1. R as the `lingua franca' of data analysis and statistical computing, 2. providing a platform for R users to discuss and exchange ideas on how R can be used to do statistical computations, data analysis, visualization and exciting applications in various fields, 3. giving an overview of the new features of the rapidly evolving R project. As for the predecessor conferences, the program will consist of two parts: invited lectures and user-contributed sessions. Prior to the conference, there will be tutorials on R, descriptions of which are available at http://www.R-project.org/useR-2010/tutorials INVITED LECTURES Invited speakers will include Mark Handcock, Frank Harrell Jr, Friedrich Leisch, Michael Meyer, Richard Stallman, Luke Tierney, Diethelm Wuertz. USER-CONTRIBUTED SESSIONS The sessions will be a platform to bring together R users, contributors, package maintainers and developers in the S spirit that `users are developers'. People from different fields will show us how they solve problems with R in fascinating applications. The sessions are organized by members of the program committee, including Dirk Eddelbuettel, John Fox, Virgilio Gomez-Rubio, Richard Heiberger, Torsten Hothorn, Aaron King, Jan de Leeuw, Nicholas Lewin-Koh, Andy Liaw, Uwe Ligges, Martin Maechler, Katharine Mullen, Heather Turner, Ravi Varadhan, H. D. Vinod, John Verzani, Alan Zaslavsky, Achim Zeileis. The program will cover topics such as * Applied Statistics Biostatistics * Bayesian Statistics * Bioinformatics * Chemometrics and Computational Physics * Data Mining * Econometrics Finance * Environmetrics Ecological Modeling * High Performance Computing * Machine Learning * Marketing Business Analytics * Psychometrics * Robust Statistics * Social network analysis * Spatial Statistics * Statistics in the Social and Political Sciences * Teaching * Visualization Graphics * and many more. IMPORTANT DATES * ** 2010-06-20 registration deadline **(later registration NOT possible on site) * 2010-07-20 tutorials 2010-07-21 conference start 2010-07-23 conference end We hope to meet you in Gaithersburg! __ 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] Calling R-tists
Participants in the R User Conference, useR! 2010, July 21-23, (http://R-project.org/useR-2010) will each receive a t-shirt, thanks to the sponsorship of Mango Solutions (http://www.mango-solutions.com/). This email is a call for designs for the front of the t-shirt. The design should be made using a single color of your choice. The design should be in the form of a high-resolution (at least 300 dpi) jpeg, png or gif file, and will be printed on a 13'' x 13'' area on the shirt. There are no rules for the content. The shirts will be white. The back of the t-shirt will include the useR! logo (in blue and black), and the logo of Mango Solutions (in orange and black). The staff of Mango Solutions will decide which design to use. The creator of the chosen design will receive 5 free t-shirts plus a book token (i.e., a gift certificate for a book). You don't have to be registered for the conference to enter a design. Please email your design proposal by Sunday, June 6, to usert-sh...@mango-solutions.com. The designer that submits the chosen design will be notified by June 15th. Thanks in advance for enriching the conference and the wardrobes of useR! participants! Kate Mullen, for the useR! 2010 Organizing Committee __ 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] useR!, regular registration deadline today
This is a reminder that the regular registration deadline for the useR! 2010 conference, July 21-23 (http://www.R-project.org/useR-2010) is today. After today, late registration will be possible until June 20th. The link to submit registration information is: http://www.R-project.org/useR-2010/registration/registration.html __ 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] useR! 2010, Submission/early registration deadline: March 1
Dear R-help members, This message is to remind you that the submission deadline for abstracts for the R User Conference, useR! 2010, is one week away: Submission deadline: Monday, March 1, 2010 The deadline for early registration is also Monday, March 1, 2010. I encourage you all to submit abstracts (to give a presentation, or to present a poster). More information regarding the conference is pasted below. best regards, Kate Mullen, for the useR! 2010 organizing committee -- We are happy to announce that the R User Conference useR! 2010 http://www.R-project.org/useR-2010 is scheduled for July 21-23, 2010, and will take place at the headquarters of the National Institute of Standards and Technology (NIST) in Gaithersburg, Maryland, USA. Following the successful useR! 2004, useR! 2006, useR! 2007, useR! 2008, and useR! 2009, conferences, the conference is focused on: 1. R as the `lingua franca' of data analysis and statistical computing, 2. providing a platform for R users to discuss and exchange ideas on how R can be used to do statistical computations, data analysis, visualization and exciting applications in various fields, 3. giving an overview of the new features of the rapidly evolving R project. As for the predecessor conferences, the program will consist of two parts: invited lectures and user-contributed sessions. Prior to the conference, there will be tutorials on R, descriptions of which are available at http://www.R-project.org/useR-2010/tutorials All R users are invited to submit abstracts on exciting applications of R as specified in the call at http://www.R-project.org/useR-2010/#Call The deadline for abstract submission (and early registration) is March 1, 2010. INVITED LECTURES Invited speakers will include Mark Handcock, Frank Harrell Jr, Friedrich Leisch, Michael Meyer, Richard Stallman, Luke Tierney, Diethelm Wuertz. USER-CONTRIBUTED SESSIONS The sessions will be a platform to bring together R users, contributors, package maintainers and developers in the S spirit that `users are developers'. People from different fields will show us how they solve problems with R in fascinating applications. The sessions are organized by members of the program committee, including Dirk Eddelbuettel, John Fox, Virgilio Gomez-Rubio, Richard Heiberger, Torsten Hothorn, Aaron King, Jan de Leeuw, Nicholas Lewin-Koh, Andy Liaw, Uwe Ligges, Martin Maechler, Katharine Mullen, Heather Turner, Ravi Varadhan, H. D. Vinod, John Verzani, Alan Zaslavsky, Achim Zeileis. The program will cover topics such as * Applied Statistics Biostatistics * Bayesian Statistics * Bioinformatics * Chemometrics and Computational Physics * Data Mining * Econometrics Finance * Environmetrics Ecological Modeling * High Performance Computing * Machine Learning * Marketing Business Analytics * Psychometrics * Robust Statistics * Social network analysis * Spatial Statistics * Statistics in the Social and Political Sciences * Teaching * Visualization Graphics * and many more. IMPORTANT DATES 2009-10-01 open submission of abstracts 2009-10-01 open registration 2009-11-01 tutorial submission deadline ** 2010-03-01 early registration deadline ** 2010-03-01 submission deadline for abstracts Before 2010-03-15notification of acceptance 2010-06-20 registration deadline (later registration NOT possible on site) 2010-07-20 tutorials 2010-07-21 conference start 2010-07-23 conference end We hope to meet you in Gaithersburg! __ 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] nls.lm AIC
I will consider putting methods for AIC and logLik into the next version
of minpack.lm (contributions welcome).
For now, the following should work for logLik, where 'object' is the
return value of nls.lm.
logLik.nls.lm - function(object, REML = FALSE, ...)
{
res - object$fvec
N - length(res)
val - -N * (log(2 * pi) + 1 - log(N) + log(sum(res^2)))/2
## the formula here corresponds to estimating sigma^2.
attr(val, df) - 1L + length(coef(object))
attr(val, nobs) - attr(val, nall) - N
class(val) - logLik
val
}
On Tue, 16 Feb 2010, Baudron, Alan Ronan wrote:
Hi there,
I'm a PhD student investigating growth patterns in fish. I've been using
the minpack.lm package to fit extended von Bertalanffy growth models
that include explanatory covariates (temperature and density). I found
the nls.lm comand a powerful tool to fit models with a lot of
parameters. However, in order to select the best model over the possible
candidates (without covariates, with both covariates, or with only one
of them) I'd like to compare them based on their AIC criterion. However,
it seems that the nls.lm comand doesn't return an AIC, or a Log
Likelihood. Does someone have any idea of how I could proceed to get
such informations about my models?
Thanks for your help. Best regards,
Alan Baudron
The University of Aberdeen is a charity registered in Scotland, No SC013683.
[[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] Abstract submission deadline, useR! 2010
The abstract submission deadline for the R Users Conference, useR! 2010, is just one month away. Submission deadline: 2010-03-01 Early registration will also close on 2010-03-01. Now is the perfect time to prepare an abstract presenting innovations or exciting applications of R. The call for abstracts along with the link for submission is available at: http://www.R-project.org/useR-2010/#Call This meeting of the R user community will take place at the Gaithersburg, Maryland, USA campus of the National Institute of Standards and Technology (NIST), July 21-23, 2010. The conference schedule is comprised of invited lectures and user-contributed sessions as well as half-day tutorials presented by R experts on July 20, 2010, prior to the conference. Invited speakers will include Mark Handcock, Frank Harrell Jr, Friedrich Leisch, Michael Meyer, Richard Stallman, Luke Tierney, and Diethelm Wuertz. The full spectrum of conference information is available from the webpage: http://www.R-project.org/useR-2010/ We hope to meet you in Gaithersburg! Kate Mullen, for the organizing committee ___ r-annou...@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-announce __ 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] non-linear regression
The problem of estimation of parameters in R is that you have to know the
value of the initial estimates very accurately, otherwise it does not
converge.
This was discussed on R-help in the last 2 weeks; see the thread on
'Starting estimates for nls Exponential Fit'.
The example below could be resolved in Excel, however in does not converge.
How to solve the problem?
Can you consider a model with a different functional form or do you need
to fit this exact function?
If you want to minimize log(data) log(model) RSS, then:
tx.br - read.table('tx.br.H.txt',header=F,dec=',')
tx.br - tx.br[,1]
id - 1:100
qx.suav - function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
newD - log(tx.br)
HP - nls(newD~log(qx.suav(id,A,B,C,D,E,F,G,H,K)),
data=data.frame(id=id,newD=newD),
trace=TRUE,
nls.control(maxiter=5,warnOnly=TRUE),
algorithm='port',
start=list(A=0.000644,B=0.016761290,C=0.10927095582,D=0.00094877,
E=5.949082737,F=24.526811,G=0.46733960,H=1.0970550987,
K=0.771722501657),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
HP
par(mfrow=c(2,1))
matplot(cbind(fitted(HP), newD),type=l, main=model and fit)
matplot(cbind(exp(fitted(HP)), exp(newD)), type=l,
main=transformed back to original space)
I made the chart on a logarithmic scale to better visualize the differences.
Send the data file attached.
The commands are below:
tx.br - read.table('c:/tx.br.H.txt',header=F,dec=',')
tx.br -tx.br[,1]
id-1:100
qx.suav - function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
HP - nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5,warnOnly=TRUE,minFactor =
0.1),
algorithm='port',
start=list(A=0.000644,B=0.016761290,C=0.10927095582,D=0.00094877,
E=5.949082737,F=24.526811,G=0.46733960,H=1.0970550987,K=0.771722501657),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
HP
matplot(cbind(log(fitted(HP)), log(tx.br)),type=l)
- Original Message -
From: Katharine Mullen k...@few.vu.nl
To: AneSR citb...@terra.com.br
Cc: r-help@r-project.org
Sent: Thursday, December 10, 2009 9:55 PM
Subject: Re: [R] non-linear regression
You did not provide the data or a way of generating it.
I would guess that Excel finds the same solution (the same residual sum-of
squares) as nls but that it uses a weak convergence criterion and/or does
not give you information regarding why it terminates.
Regarding the step size: you can set the minimum step size factor via the
minFactor argument of control.
qx.suav - function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
## make noisy data from model
id - 1:1000
tx.br - qx.suav(id,A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638)
set.seed(1)
tx.br - tx.br + rnorm(length(tx.br),0,.0001)
HP - nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5000,warnOnly=TRUE),
algorithm='port',
start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
matplot(cbind(fitted(HP), tx.br),type=l)
On Thu, 10 Dec 2009, AneSR wrote:
I have a non-linear regression with 8 parameters to solve however it
does not converge ... easily solves the excel ... including the initial
estimates used in the R were found in the excel ... Another question is
how
to establish the increments of R by the parameters in the search ..
qx.suav-function(id,A,B,C,D,E,F,G,H,K){(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))}
HP-nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5000),algorithm='port',start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,K=0.382070638),lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
summary(HP)
How to solve this problem in R?
Ane
--
View this message in context:
http://n4.nabble.com/non-linear-regression-tp959977p959977.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
Re: [R] non-linear regression
You did not provide the data or a way of generating it.
I would guess that Excel finds the same solution (the same residual sum-of
squares) as nls but that it uses a weak convergence criterion and/or does
not give you information regarding why it terminates.
Regarding the step size: you can set the minimum step size factor via the
minFactor argument of control.
qx.suav - function(id,A,B,C,D,E,F,G,H,K)
(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
## make noisy data from model
id - 1:1000
tx.br - qx.suav(id,A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638)
set.seed(1)
tx.br - tx.br + rnorm(length(tx.br),0,.0001)
HP - nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5000,warnOnly=TRUE),
algorithm='port',
start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,
K=0.382070638),
lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
matplot(cbind(fitted(HP), tx.br),type=l)
On Thu, 10 Dec 2009, AneSR wrote:
I have a non-linear regression with 8 parameters to solve however it
does not converge ... easily solves the excel ... including the initial
estimates used in the R were found in the excel ... Another question is how
to establish the increments of R by the parameters in the search ..
qx.suav-function(id,A,B,C,D,E,F,G,H,K){(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))}
HP-nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),data=data.frame(id=id,tx.br=tx.br),
trace=TRUE,nls.control(maxiter=5000),algorithm='port',start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,E=0.0002127,F=38.5448420,G=0.115,H=1.1109286,K=0.382070638),lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
summary(HP)
How to solve this problem in R?
Ane
--
View this message in context:
http://n4.nabble.com/non-linear-regression-tp959977p959977.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-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] Starting estimates for nls Exponential Fit
You used starting values: pa - c(1,2,3) but both algorithms (port and Gauss-Newton) fail if you use the slightly different values: pa - c(1,2,3.5) Scaling does not fix the underlying sensitivity to starting values. pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both alg. also fail if there is insufficient spread (e.g., both fail for pa - c(1,1.25,1.5) ). For the record, DE uses (by default at least) a random start and for bad starts will sometimes fail to give good results; decrease the probability this happens by increasing itermax from the default itermax=200, as in: ss - DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10), control=list(trace=FALSE, itermax=1000)) On Wed, 2 Dec 2009, Prof. John C Nash wrote: Kate Mullen showed one approach to this problem by using DEOptim to get some good starting values. However, I believe that the real issue is scaling (Warning: well-ridden hobby-horse!). With appropriate scaling, as below, nls does fine. This isn't saying nls is perfect -- I've been trying to figure out how to do a nice job of helping folk to scale their problems. Ultimately, it would be nice to has an nls version that will do the scaling and also watch for some other situations that give trouble. Cheers, JN ## JN test rm(list=ls()) ExponValues - c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, 2468.32,2778.47) ExponCycles - c(17,18,19,20,21,22,23,24,25) mod - function(x) x[1] + x[2]*x[3]^ExponCycles modj - function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) fun - function(x) sum((ExponValues-mod(x))^2) pa-c(1,2,3) lo-c(0,0,0) up-c(20,20,20) names(pa) - c(Y0, a, E) ## fit w/port and GN Emodjn- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, algorithm='port', lower=lo, upper=up, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) Emodjn1 - nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) ## fit matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type=l) __ 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] Starting estimates for nls Exponential Fit
If you could reformulate your model as alpha * (y0 + E^t) then you could
use nls with alg=plinear (alpha then would be eliminated from the
nonlinear param and treated as conditionally linear) and this would help
with convergence.
Else you can try package DEoptim to get the starting values; the advantage
over optim is that you need then lower and upper bounds on the parameters,
not more starting values; the bounds however should be appropriate and as
tight as possible.
Also: by default nls uses max. 50 iter. Depending on where you start off
you may need more than this.
##
ExponValues - c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
2468.32,2778.47)
ExponCycles - c(17,18,19,20,21,22,23,24,25)
library(DEoptim)
mod - function(x) x[1] + x[2]*x[3]^ExponCycles
fun - function(x) sum((ExponValues-mod(x))^2)
ss - DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),
control=list(trace=FALSE))
pa - ss$optim$bestmem
## now have par est that give OK fit
matplot(cbind(ExponValues, mod(pa)),type=l)
names(pa) - c(Y0, a, E)
## fit w/port and GN
Emodel- nls(ExponValues ~ (Y0 + a*(E^ExponCycles)),
start=pa, algorithm=port,
control=list(maxiter=1000, warnOnly=TRUE))
Emodel1 - nls(ExponValues ~ (Y0 + a*(E^ExponCycles)), start=pa,
control=list(maxiter=1000, warnOnly=TRUE))
## fit
matplot(cbind(ExponValues, fitted(Emodel), fitted(Emodel1)),type=l)
On Tue, 1 Dec 2009, Anto wrote:
Hello everyone,
I have come across a bit of an odd problem: I am currently analysing data
PCR reaction data of which part is behaving exponential. I would like to fit
the exponential part to the following:
y0 + alpha * E^t
In which Y0 is the groundphase, alpha a fluorescence factor, E the
efficiency of the reaction t is time (in cycles)
I can get this to work for most of my reactions, but part fails the nls fit
(Convergence failure: iteration limit reached without convergence). This
mainly has to do with the starting estimates for the nls function. I have
used various 'smart' ways of calculating starting estimates (grid searches,
optim(), etc.) but however close the starting estimates are to the actual
values, nls still fails for many datasets.
The weirdest thing is, I have one set of starting estimates (alpha= 0.5 and
E=1.85) which are totaly arbitray and these DO work for, say 99% of the
datasets. I don't know why this is and I don't why my 'good estimates' do
not work. I am desperatly seeking a way of calculating working starting
estimates for my nls function. Can anybody give me a hand?
thanks,
Anto
R version 2.9.2
example dataset:
ExponValues
[1] 2018.34 2012.54 2018.85 2023.52 2054.58 2132.61 2247.17 2468.32 2778.47
ExponCycles
[1] 17 18 19 20 21 22 23 24 25
Example starting estimate calculation
Y0 - ExponValues[1]
E -
weighted.mean((ExponValues-eY0)[-1][-1]/(ExponValues-eY0)[-1][-(length(ExponValues)-1)],(1-abs(seq(-1,1,length=(length(ExponValues)-2)))^3)^3)
alpha -
weighted.mean((ExponValues[-1]-ExponValues[-length(ExponValues)])/((E^ExponCycles[-1])-(E^ExponCycles[-length(ExponCycles)])),(1-abs(seq(-1,1,length=(length(ExponCycles)-1)))^3)^3)
Optional parameter optimisation:
Esp -
optim(c(Y0=eY0,a=alpha,E=E),function(x){return(sum(abs(ExponValues-(x[1]+x[2]*x[3]^ExponCycles})$par
nls function:
Emodel-try(nls(ExponValues ~ (Y0 + a*(E^ExponCycles)) , start= Esp,
algorithm=port),silent=TRUE)
--
View this message in context:
http://n4.nabble.com/Starting-estimates-for-nls-Exponential-Fit-tp932230p932230.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] useR! 2010: submission registration started!
We are happy to inform you that abstract submission and registration for `useR! 2010' is now available online from http://www.R-project.org/useR-2010/ This meeting of the R user community will take place at the Gaithersburg, Maryland, USA campus of the National Institute of Standards and Technology (NIST), July 21-23, 2010. The conference schedule is comprised of invited lectures and user-contributed sessions as well as half-day tutorials presented by R experts on July 20, 2010, prior to the conference. Invited speakers will include Mark Handcock, Frank Harrell Jr, Friedrich Leisch, Michael Meyer, Richard Stallman, Luke Tierney, Diethelm Wuertz. We invite you to submit abstracts on topics presenting innovations or exciting applications of R. The call for papers along with the link for abstract submission is available at http://www.R-project.org/useR-2010/#Call ** Reminder: tutorial proposals are due by November 1, 2009, as specified in the call for tutorials available at http://www.R-project.org/useR-2010/tutorials/ ** We hope to meet you in Gaithersburg! The organizing committee, Kevin Coakley, Nathan Dodder, David Gil, William Guthrie, Walter Liggett, John Lu, Katharine Mullen, Jonathon Phillips, Antonio Possolo, Ravi Varadhan. ___ r-annou...@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-announce __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [R-pkgs] DEoptim 2.0-0
Dear All, We are happy to announce the release of the new version of DEoptim (version 2.0-0) which is now available from CRAN. The DEoptim package [3] performs Differential Evolution (DE) minimization, a genetic algorithm-based optimization technique [2,3]. This allows robust minimization over a continuous (bounded or not) domain. The new DEoptim function calls a C implementation of the DE algorithm similar to the MS Visual C++ v5.0 implementation distributed with [2]. More details on DE optimization can be found on the DE homepage [1]. We believe that the DE approach may be applicable in many fields of research and hope that the package DEoptim will be fruitful for many researchers. Kate Mullen and David Ardia MODIFICATIONS o The R-based implementation of Differential Evolution has been replaced with a C-based implementation similar to the MS Visual C++ v5.0 implementation accompanying the book `Differential Evolution - A Practical Approach to Global Optimization',downloaded from http://www.icsi.berkeley.edu/~storn/DeWin.zip. The new C implementation is significantly faster. o The S3 method for plotting has been enhanced. It allows now to plot the intermediate populations if provided. o The package maintainer has been changed to Katharine Mullen, katharine.mul...@nist.gov. o A NAMESPACE has been added. o Argument FUN for DEoptim is now called fn for compatibility with optim. o demo file has been removed o CITATION file modified REFERENCES [1] Differential Evolution homepage: http://www.icsi.berkeley.edu/~storn/code.html [2] Price, K.V., Storn, R.M., Lampinen J.A. (2005). Differential Evolution - A Practical Approach to Global Optimization. Springer-Verlag. ISBN 3540209506. [3] Ardia, D., Mullen, K. (2009). DEoptim: Differential Evolution Optimization in R. R package version 2.0-0. URL http://CRAN.R-project.org/package=DEoptim ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ 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] curve fitting
or use nls.lm as in install.packages(minpack.lm) library(minpack.lm) x - c(2, 8, 14, 20, 26, 32, 38, 44, 50, 56, 62, 68, 74) y - c(100, 99, 99, 98, 97, 94, 82, 66, 48, 38, 22, 10, 1) res - function(p, x, y) y - ff(p,x) ff - function(p, x) 100*exp(p[1]*(1-exp(p[2]*x))/p[2]) aa - nls.lm(par=c(.0001,.0001), fn=res, x=x, y=y) plot(x,y) lines(x, ff(coef(aa), x)) On Tue, 12 May 2009, Jorge Ivan Velez wrote: Dear Dmitry, Take a look at ?nls and its examples. HTH, Jorge On Tue, May 12, 2009 at 5:44 PM, Dmitry Gospodaryov gospodar...@rambler.ruwrote: I have the data: for x: 2, 8, 14, 20, 26, 32, 38, 44, 50, 56, 62, 68, 74, for y: 100, 99, 99, 98, 97, 94, 82, 66, 48, 38, 22, 10, 1. y depends on x by equation: y = 100*exp(b*(1-exp(c*x))/c), where b and c are coefficients. I need to find coefficients in this equation for given data. How can I do that by means of R? Thank you for advance. With regard, Dmitry. __ 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-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] Biexponential Fit
Using algorithm=plinear as shown by example below makes sum-of-exponentials fitting problems better conditioned. A1 - 1 A2 - 2 k1 - -.5 k2 - -2 x - seq(1,10,length=200) y - A1*exp(k1*x) + A2*exp(k2*x) + .001*rnorm(200) aa - nls(y~cbind(exp(k1*x), exp(k2*x)), algorithm=plinear, start=list(k1=-.1, k2=-6), data=list(y=y, x=x), trace=TRUE) On Thu, 9 Apr 2009, Peter Dalgaard wrote: Jonas Weickert wrote: Hi, I want to do a biexponential Fit, i.e. y ~ A1*exp(k1*x) + A2*exp(k2*x) Is this possible? I tried nls() but it stopped with several (different) errors. I'm using y and x as simple vectors and the formula for nls() exactly as mentioned above. Yes, it is possible, with some reservations. There is even a self-starting model for it (SSbiexp), at least if k1,k2 are negative. The reservations are that you need good data and good starting values for the fit. The model is theoretically unidentifiable because you can switch 1 and 2 and get the same model, which becomes an issue if you start too fra from the optimum. Worse, if k1 and k2 are similar, the whole system becomes ill-conditioned. There's a rule of thumb requiring a factor of 5 or more between k1 and k2, for pharmacokinetic applications. Thanks a lot! Jonas __ 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 ~~ - (p.dalga...@biostat.ku.dk) 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. __ 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] Gradient function for optim.
see the fmingr function in src/main/optim.c (https://svn.r-project.org/R/trunk/src/main/optim.c) On Wed, 25 Feb 2009 rkevinbur...@charter.net wrote: I have read that when the gradient function is not supplied (is null) then first order differencing is used to find the differential. I was trying to track down this for my own information but I run into .Internal(optim.). I was not sure where to look next to see the function that is automatically supplied for the gradient. Thank you. Kevin __ 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] Chromatogram deconvolution and peak matching
Hoi Bart, I think you're right that ALS should be applicable to this problem. Unfortunately in writing I see that there is a bug when the spectra are NOT constrained to nonnegative values (the package has been used to my knowledge only in fitting multiway mass spectra thus far, where this constraint is always applied); I'll have a patched version soon. I'd be interested in hearing about where the manual could be improved, too. #2D chromatogram generation time - seq(0,20,by=0.05) f - function(x,rt) dnorm((x-rt),mean=0,sd=rt/35) C1 - matrix(c( f(time,6.1), f(time,5.6), f(time,15)), nrow=401,ncol=3) C2 - matrix(c( f(time,2.1), f(time,4), f(time,8)), nrow=401,ncol=3) #spectrum generation spectra - function(x,a,b,c,d,e) a + b*(x-e) + c*((x-e)^2) + d*((x-e)^3) x - 220:300 s1 - spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),0) s2 - spectra(x,(-1.054e02),1.3,(-5.239e-03),(6.927e-06),-20) s3 - spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),20) spec - matrix(c(s1,s2,s3), nrow=81,ncol=3) ## data matrix 1 chrom1 - C1 %*% t(spec) ## data matrix 2 chrom2 - C2 %*% t(spec) require(ALS) ## you want to recover 2 (401 x 3) matrices C1 and C2 and a ## (82 x 3) matrix E such that ## chrom1 = C1 E^T, chrom2 = C2 E^T ## some starting guess for elution profiles ## these need to be pretty good Cstart1 - matrix(c( f(time,4), f(time,6), f(time,12)), nrow=401,ncol=3) Cstart2 - matrix(c( f(time,3), f(time,5), f(time,10)), nrow=401,ncol=3) xx - als( CList = list(Cstart1, Cstart2), PsiList = list(chrom1, chrom2), S=matrix(1, nrow=81, ncol=3), x=time, x2=x, uniC=TRUE, normS=TRUE, nonnegS=TRUE) par(mfrow=c(3,1)) matplot(time, xx$CList[[1]], col=2, type=l, main=elution profiles chrom1, lty=1, ylab=amplitude) matplot(time, C1, col=1, add=TRUE, type=l, lty=1) legend(10,2,legend=c(real, estimated), col=1:2, lty=1) matplot(time, xx$CList[[2]], col=2, type=l, main=elution profiles chrom2, lty=1, ylab=amplitude) matplot(time, C2, col=1, add=TRUE, type=l, lty=1) legend(10,6,legend=c(real, estimated), col=1:2, lty=1) matplot(x, xx$S, col=2, type=l, main=spectra, lty=1, ylab=amplitude) matplot(x, spec, col=1, add=TRUE, type=l, lty=1) legend(270,.95,legend=c(real, estimated), col=1:2, lty=1) On Tue, 17 Feb 2009, bartjoosen wrote: Hi, I'm trying to match peaks between chromatographic runs. I'm able to match peaks when they are chromatographed with the same method, but not when there are different methods are used and spectra comes in to play. While searching I found the ALS package which should be usefull for my application, but I couldn't figure it out. I made some dummy chroms with R, which mimic my actual datasets, to play with, but after looking at the manuals of ALS, I'm affraid I can't get the job done. Can someone put me on the right way? Here is my code to generate the dummy chroms, which also plots the 2 chroms and the spectra of the 3 peaks: #2D chromatogram generation par(mfrow=c(3,1)) time - seq(0,20,by=0.05) f - function(x,rt) dnorm((x-rt),mean=0,sd=rt/35) c1 - f(time,6.1) c2 - f(time,5.6) c3 - f(time,15) plot(c1+c2+c3~time,type=l,main=chrom1) #spectrum generation spectra - function(x,a,b,c,d,e) a + b*(x-e) + c*((x-e)^2) + d*((x-e)^3) x - 220:300 s1 - spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),0) s2 - spectra(x,(-1.054e02),1.3,(-5.239e-03),(6.927e-06),-20) s3 - spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),20) chrom1.tot - data.frame(time,outer(c1,s1,*)+outer(c2,s2,*)+outer(c2,s2,*)) names(chrom.tot)[-1] - x #generation of chromatogram 2 c1 - f(time,2.1) c2 - f(time,4) c3 - f(time,8) plot(c1+c2+c3~time,type=l,main=chrom2) chrom2.tot - data.frame(time,outer(c1,s1,*)+outer(c2,s2,*)+outer(c2,s2,*)) names(chrom.tot)[-1] - x plot(s1~x,type=l,main=spectra) lines(s2~x,col=2) lines(s3~x,col=3) Thanks for your time Kind Regards Bart -- View this message in context: http://www.nabble.com/Chromatogram-deconvolution-and-peak-matching-tp22057592p22057592.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] Task View for Chemometrics and Computational Physics
Dear All, A new task view ChemPhys on chemometrics and computational physics is available on CRAN (http://cran.r-project.org/web/views/ChemPhys.html). It describes packages and functions that are of use in modeling chemical/physical systems. Suggestions and comments regarding this task view are welcome. If you think a new category, package or function should be added, please mail. best regards, Kate Mullen Katharine Mullen mail: Department of Physics and Astronomy, Faculty of Sciences Vrije Universiteit Amsterdam, de Boelelaan 1081 1081 HV Amsterdam, The Netherlands room: T.1.06 tel: +31 205987870 fax: +31 205987992 e-mail: [EMAIL PROTECTED] homepage: http://www.nat.vu.nl/~kate/ __ 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 text and output to file with flexibility
you can combine write and write.table, using append=TRUE.
e.g.,
write1 - function(txt, mat, filename) {
for(i in 1:length(txt)) {
write(txt[i], file=filename, append=i!=1)
write.table(mat[[i]], file=filename, append=TRUE,sep =\t)
}
}
t - c(text1, text2, text3)
m - list(matrix(0,2,2), matrix(1,2,1), matrix(2,3,2))
write1(t, m, x.txt)
On Sat, 30 Aug 2008, [ISO-8859-1] Gon?alo Ferraz wrote:
Hi,
I want a function to write some of its output into a text file with
the following format:
'some text'
output matrix A
'some more text'
output matrix B
'some other text still'
output matrix C
...
The dimensions of matrices A, B, ... are different and the total
number of matrices
that I want to place in the text file depends on what I pass to the
function. (I would like to have values
from adjacent columns separated by tabs.)
'write' seems to place only one matrix into one file. 'sink' seems to
write only what I write, word by word.
Is there a standard solution for this type of problem?
Thank you,
Gon?alo
__
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] Optim stripping attributes from relistable objects
The following code is inspired by the help file for the relist()
function (see?relist), which explicitly details how you can use a
The example is in the 'Details' section and, indeed, it looks like it no
longer works.
relistable object in conjunction with optim to pass and reconstruct
complex parameter structures/groupings. The idea is that the optim()
function can only work with vectors, but in many cases you would like
to use a complex structure inside the objective function- relist is
one way to do that. The problem is that optim appears to be stripping
the attributes and therefore the example doesn't seem to run, giving
the error at the bottom.
You can get around this by specifying skeleton for relist:
rb.banana - function(params) {
params - relist(params, skeleton=list(x=NA,y=NA))
return( (1-params$x)^2 + 100*(params$y - params$x^2)^2)
}
ipar - as.relistable(list(x=5,y=0))
initial.params - unlist(ipar)
xx - optim(unlist(initial.params), rb.banana)
rb.banana - function(params) {
+ #Params is initially a vector
+ cat(Params initially has the attributes:\n)
+ print(names(attributes(params)))
+ #Relisting it turns it into a list...
+ params - relist(params)
+ cat(-\n)
+ #..which can then be called in the standard list manner
+ return( (1-params$x)^2 + 100*(params$y - params$x^2)^2)
+ }
ipar- as.relistable(list(x=5,y=0))
initial.params - unlist(ipar)
#Test to see if rb.banana works properly in the normal case
rb.banana(initial.params)
Params initially has the attributes:
[1] namesskeleton
-
[1] 62516
#OK, that's good. How about with optim though?
optim(initial.params,rb.banana)
Params initially has the attributes:
[1] names
Error in relist(params) :
The flesh argument does not contain a skeleton attribute.
Either ensure you unlist a relistable object, or specify the skeleton
separately.
What's going on here? Has this functionality been removed and the
documentation in relist() not updated? Or has the feature been broken?
Or have I misinterpreted something here (it wouldn't be the first
time!!)
Am running R version 2.7.1 (2008-06-23) under windows.
Cheers,
Mark
__
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] Port package
It is not an R package, but rather a collection of Fortran functions that R uses from netlib: http://www.netlib.org/port/ On Wed, 9 Jul 2008, Jos Kaefer wrote: Hi When I type: ?nls I come across this section: algorithm: character string specifying the algorithm to use. The default algorithm is a Gauss-Newton algorithm. Other possible values are 'plinear' for the Golub-Pereyra algorithm for partially linear least-squares models and 'port' for the 'nl2sol' algorithm from the Port package. The simple question is: where's the Port package? I can't find it on cran. Thanks, Jos __ 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] bug in nls?
Dear Petr, I think it's a feature. the formula interface also won't let you specify the slots of S4 objects in the model spec. How about coef(nls(y~a*x^b, data=list(x=DF[,1], y=DF[,2]), start=list(a=3, b=.7))) ? On Thu, 26 Jun 2008, Petr PIKAL wrote: Dear all Nobody responded to my previous post so far so I try with more offending subject. I just encountered a strange problem with nls formula. I tried to use nls in cycle but I was not successful. I traced the problem to some parse command. Here is an example DF-data.frame(x=1:10, y=3*(1:10)^.5+rnorm(10)) coef(lm(log(DF[,2])~log(DF[,1]))) (Intercept) log(DF[, 1]) 0.74373200.6831726 # this works coef(nls(y~a*x^b, data=DF, start=list(a=3, b=.7))) a b 2.6412881 0.5545907 # OK, this works too coef(nls(DF[,2]~a*DF[,1]^b, data=DF, start=list(a=3, b=.7))) Error in parse(text = x) : unexpected end of input in ~ coef(nls(DF[,2]~a*DF[,1]^b, start=list(a=3, b=.7))) Error in parse(text = x) : unexpected end of input in ~ # but this does not Browse[1] debug: mf$formula - as.formula(paste(~, paste(varNames[varIndex], collapse = +)), env = environment(formula)) Browse[1] Error in parse(text = x) : unexpected end of input in ~ Actually the problem is that with calling nls with DF[,n]~... varNames and varIndex is not correctly specified. I am not sure if this behaviour is a bug or feature. If it is a feature, please help me how to call variables from data frame when using nls inside cycle for (i in ) result[i,] - coef(nls( , )) Thank you Petr sessionInfo() R version 2.8.0 Under development (unstable) (2008-05-18 r45723) i386-pc-mingw32 locale: LC_COLLATE=Czech_Czech Republic.1250;LC_CTYPE=Czech_Czech Republic.1250;LC_MONETARY=Czech_Czech Republic.1250;LC_NUMERIC=C;LC_TIME=Czech_Czech Republic.1250 attached base packages: [1] stats grDevices datasets utils graphics methods base other attached packages: [1] nlme_3.1-88lattice_0.17-7 fun_0.1 loaded via a namespace (and not attached): [1] grid_2.8.0 tools_2.8.0 Petr Pikal [EMAIL PROTECTED] 724008364, 581252140, 581252257 __ 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] YAPMQ - Yet Another PlotMath Question
by example: slope - 45 plot(1:10) text(2, 7, labels = expression(45~degree)) ; text(2, 5, labels = substitute(paste(slope * degree), list(slope=45))) text(2, 3, labels = substitute(paste(lambda * = * mt * (nm)), list(mt = 33))) On Thu, 19 Jun 2008, Thompson, David (MNR) wrote: Hello, I'm having trouble finding (remembering) how to pass values into text functions in plots, as demonstrated by: slope - 45 ; plot(1:10) ; text(2, 7, labels = expression(45~degree)) ; text(2, 5, labels = paste(bquote(.(slope)), expression(degree))) Thanx, DaveT. * Silviculture Data Analyst Ontario Forest Research Institute Ontario Ministry of Natural Resources [EMAIL PROTECTED] http://ontario.ca/ofri __ 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] adding custom axis to image.plot() and strange clipping behavior
I also noticed that adding a custom axis with image.plot was a problem;
you can also do:
library(fields)
m - matrix(1:15,ncol=3)
par(mar=c(5,5,5,7))
image(m, axes=FALSE)
# add axis
axis(1,axTicks(1),lab=letters[1:length(axTicks(1))])
box()
## add legend
image.plot(m, legend.only=TRUE)
On Fri, 13 Jun 2008, Stephen Tucker wrote:
Hi list,
I wanted to plot an image with a colorbar to the right of the plot, but
set my own axis labels (text rather than numbers) to the image. I have
previously accomplished this with two calls to image(), but the package
'fields' has a wrapper function, image.plot(), which does this task
conveniently.
However, I could not add axes to the original image after a call to
image.plot(); I have found that I needed to set par(xpd=TRUE) within the
function to allow this to happen:
###=== begin code
library(fields)
## make data matrix
m - matrix(1:15,ncol=3)
## plot
image.plot(m,axes=FALSE)
axis(1) # doesn't work
par(xpd=TRUE)
axis(1) # still doesn't work
## replace the 28th element of the body of image.plot()
## and assign to new function called 'imp'
## here I just use the second condition of 'if' statement
## and set 'xpd = TRUE'
imp - `body-`(image.plot,value=`[[-`(body(image.plot),28,
quote({par(big.par)
par(plt = big.par$plt, xpd = TRUE)
par(mfg = mfg.save, new = FALSE)
invisible()})))
imp(m,axes=FALSE)
box()
axis(1,axTicks(1),lab=letters[1:length(axTicks(1))])
## clip to plotting region for additional
## graphical elements to be added:
par(xpd=FALSE)
abline(v=0.5)
###=== end code
I wonder if anyone has any insights into this behavior? Since in the axis()
documentation, it says:
Note that xpd is not accepted as clipping is always to the device region
I am surprised to find (1) that the par(xpd=TRUE) works in the case above,
and (2) that it must be called before the function call is terminated.
I wonder if anyone has any insights into this behavior. I have reproduced
this on both my Linux box (Ubuntu Gutsy Gibbon 64-bit, R 2.7.0, fields
package version 4.1) and Windows machine (32-bit XP Pro, R 2.7.0, fields
package version 4.1).
Thanks very much,
Stephen
__
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] optim, constrOptim: setting some parameters equal to each other
Here is an example w/optim where you have an objective function
F(x) where you have (possibly a few) constraints of form x_i=x_j, and you
can specify the constraints flexibly, which is what I _think_ you want.
An example from you would have been nice.
## objective function that depends on all parameters, here a, b, c
obfun - function(a,b,c,dd)
sum((dd - (exp(a * 1:100) + exp(b * 1:50) + exp(c * 1:25) ))^2)
## fn for optim
fr - function(x, eqspec, dd, obfun) {
## assign variables for parameter values given in x
for(i in 1:length(x))
assign(names(x)[i], x[[i]])
## use eqspec to assign parameter values determined w/equ. constr.
for(i in 1:length(eqspec))
assign( names(eqspec)[i], x[[ eqspec[[i]] ]])
## now have all parameter values, call objective function
obfun(a,b,c,d, dd)
}
dd - exp(-.05 * 1:100) + exp(-.05 * 1:50) + exp(.005 * 1:25)
## let par contain all parameters that are not eq. constrained and
## all parameters on right side of constraints of form a=b
## let eqspec give all dependent parameters on left side of
## the constraints of form a=b
## e.g.
## for constraint a=b
xx - optim(par=list(b=-4, c=-.1), fn=fr, eqspec=list(a=b),
dd=dd, obfun=obfun)
## for constraint b=a
x1 - optim(par=list(a=-4, c=-.1), fn=fr, eqspec=list(b=a),
dd=dd, obfun=obfun)
On Sun, 8 Jun 2008, Alex F. Bokov wrote:
Hello, and apologies for the upcoming naive questions. I am a biologist
who is trying to teach himself the appropriate areas of math and stats.
I welcome pointers to suggested background reading just as much as I do
direct answers to my question.
Let's say I have a function F() that takes variables (a,b,c,a1,b1,c1)
and returns x, and I want to find the values of these variables that
result in a minimum value of x. That's a straightforward application of
optim(). However, for the same function I also need to obtain values
that return the minimum value of x subject to the following constraints:
a=a1, b=b1, c=c1, a=a1 b=b1, a=a1 c=c1, ... and so on, for any
combination of these constraints including a=a1, b=b1, c=c1. The brute
force way to do this with optim() would be to write into F() an immense
switch statement anticipating every possible combination of constrained
variables. Obviously this is inefficient and unmaintainable. Therefore,
my question is:
Is the purpose of constrOptim() this exact type of problem? If so, how
does one express the constraint I described above in terms of the ui,
ci, and theta parameters? Are there any introductory texts I should have
read for this to be obvious to me from constrOptim's documentation?
If constrOptim() is not the answer to this problem, can anybody suggest
any more elegant alterntives to the switch-statement-from-hell approach?
Thank you very, very much in advance. I thought I understood R
reasonably well until I started banging my head against this problem!
__
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] optim, constrOptim: setting some parameters equal to each other
the example I just mailed had an error; it should have been:
## objective function that depends on all parameters, here a, b, c
obfun - function(a,b,c,dd)
sum((dd - (exp(a * 1:100) + exp(b * 1:50) + exp(c * 1:25) ))^2)
fr - function(x, eqspec, dd, obfun) {
## assign variables for parameter values given in x
for(i in 1:length(x))
assign(names(x)[i], x[[i]])
## use eqspec to assign parameter values determined w/equality constr.
for(i in 1:length(eqspec))
assign( names(eqspec)[i], x[[ eqspec[[i]] ]])
## now have all parameter values, call objective function
obfun(a,b,c,dd)
}
dd - exp(-.05 * 1:100) + exp(-.05 * 1:50) + exp(.005 * 1:25)
## let par contain all parameters that are not constrained and
## all parameters on right side of constraints of form a=b
## let eqspec give all dependent parameters on left side of
## the constraints of form a=b
## e.g.
## for constraint a=b
xx - optim(par=list(b=-4, c=-.1), fn=fr, eqspec=list(a=b),
dd=dd, obfun=obfun)
## for constraint b=a
x1 - optim(par=list(a=-3, c=-.1), fn=fr, eqspec=list(b=a),
dd=dd, obfun=obfun)
On Sun, 8 Jun 2008, Katharine Mullen wrote:
Here is an example w/optim where you have an objective function
F(x) where you have (possibly a few) constraints of form x_i=x_j, and you
can specify the constraints flexibly, which is what I _think_ you want.
An example from you would have been nice.
## objective function that depends on all parameters, here a, b, c
obfun - function(a,b,c,dd)
sum((dd - (exp(a * 1:100) + exp(b * 1:50) + exp(c * 1:25) ))^2)
## fn for optim
fr - function(x, eqspec, dd, obfun) {
## assign variables for parameter values given in x
for(i in 1:length(x))
assign(names(x)[i], x[[i]])
## use eqspec to assign parameter values determined w/equ. constr.
for(i in 1:length(eqspec))
assign( names(eqspec)[i], x[[ eqspec[[i]] ]])
## now have all parameter values, call objective function
obfun(a,b,c,d, dd)
}
dd - exp(-.05 * 1:100) + exp(-.05 * 1:50) + exp(.005 * 1:25)
## let par contain all parameters that are not eq. constrained and
## all parameters on right side of constraints of form a=b
## let eqspec give all dependent parameters on left side of
## the constraints of form a=b
## e.g.
## for constraint a=b
xx - optim(par=list(b=-4, c=-.1), fn=fr, eqspec=list(a=b),
dd=dd, obfun=obfun)
## for constraint b=a
x1 - optim(par=list(a=-4, c=-.1), fn=fr, eqspec=list(b=a),
dd=dd, obfun=obfun)
On Sun, 8 Jun 2008, Alex F. Bokov wrote:
Hello, and apologies for the upcoming naive questions. I am a biologist
who is trying to teach himself the appropriate areas of math and stats.
I welcome pointers to suggested background reading just as much as I do
direct answers to my question.
Let's say I have a function F() that takes variables (a,b,c,a1,b1,c1)
and returns x, and I want to find the values of these variables that
result in a minimum value of x. That's a straightforward application of
optim(). However, for the same function I also need to obtain values
that return the minimum value of x subject to the following constraints:
a=a1, b=b1, c=c1, a=a1 b=b1, a=a1 c=c1, ... and so on, for any
combination of these constraints including a=a1, b=b1, c=c1. The brute
force way to do this with optim() would be to write into F() an immense
switch statement anticipating every possible combination of constrained
variables. Obviously this is inefficient and unmaintainable. Therefore,
my question is:
Is the purpose of constrOptim() this exact type of problem? If so, how
does one express the constraint I described above in terms of the ui,
ci, and theta parameters? Are there any introductory texts I should have
read for this to be obvious to me from constrOptim's documentation?
If constrOptim() is not the answer to this problem, can anybody suggest
any more elegant alterntives to the switch-statement-from-hell approach?
Thank you very, very much in advance. I thought I understood R
reasonably well until I started banging my head against this problem!
__
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-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 can I display a characters table ?
Dear Maura, try the function textplot from the package gplots. you can say textplot(yourmatrix) and get a plot of a character matrix. On Fri, 6 Jun 2008, Maura E Monville wrote: I would like to generate a graphics text. I have a 67x2 table with 5-character string in col 1 and 2-character string in col 2. Is it possible to make such a table appear on a graphics or a message-box pop-up window ? Thank you so much. -- Maura E.M __ 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] nls() newbie convergence problem
dydata
x1 x2 y
1 9 27 248
2 9 22 213
3 4 23 190
4 11 16 183
5 -6 25 144
6 11 14 169
7 -4 13 72
8 2 8 73
9 10 13 156
10 8 30 263
11 12 10 147
12 -7 5 -2
13 0 10 75
14 12 0 77
15 9 8 115
16 12 24 245
17 34 23 370
18 12 1 84
19 10 37 324
20 26 30 371
weight - function(y,x1,x2,b0,b1,b2)
{
pred - b0+b1*x1 + b2*x2
parms - abs(b1*b2)^(1/3)
(y-pred)/parms
}
gmfit - nls(~weight(y,x1,x2,b0,b1,b2),data=dydata,trace=TRUE,
start=list(b0=1,b1=1,b2=1), control=list(warnOnly=TRUE))
library(minpack.lm)
weight2 - function(p, y,x1,x2)
{
pred - p[1]+p[2]*x1 + p[3]*x2
parms - abs(p[2]*p[3])^(1/3)
(y-pred)/parms
}
gmfit2 - nls.lm(par = c(1,1,1),
fn = weight2, y=dydata$y,x=dydata$x1,x2=dydata$x2,
control=list(nprint=1))
dydata is from Norman R. Draper, Yonghong (Fred) Yang, Generalization
of the geometric mean functional relationship, Computational
Statistics Data AnalysisVolume 23, Issue 3, 9 January 1997, Pages
355-372; I don't think the attachment got through.
nls uses an orthogonality convergence criterium. Bates talks about this
in the post at
https://stat.ethz.ch/pipermail/r-devel/2000-August/021059.html and it's
also described in the book by Bates and Watts (Nonlinear regression
analysis its applications). Apparently here the residuals are so small
they are effectively 0, and this criteria does not work; there is a
warning about using zero-residual data in the help page for nls.
nls.lm from the package minpack.lm stops if any of a few different
criteria are met; in this case it stops at the time you think is right.
On Wed, 4 Jun 2008, Bernard Leemon wrote:
I'm sure this must be a nls() newbie question, but I'm stumped.
I'm trying to do the example from Draper
and Yang (1997). They give this snippet of S-Plus code:
Specify the weight function:
weight - function(y,x1,x2,b0,b1,b2)
{
pred - b0+b1*x1 + b2*x2
parms - abs(b1*b2)^(1/3)
(y-pred)/parms
}
Fit the model
gmfit -nls(~weight(y,x1,x2,b0,b1,b2), observe,list(starting value))
in converting this to R, I left the weight function alone and replaced the
nls() with
gmfit -
nls(~weight(y,x1,x2,b0,b1,b2),data=dydata,trace=TRUE,start=list(b0=1,b1=1,b2=1))
where dydata is the appropriate data.frame.
nls() quickly (6 iterations) finds the exact values from Draper Yang for
b0, b1, and b2 but
despite reporting a discrepancy of only 3.760596e-29 by the 7th iteration,
it merrily goes on
to 50 iterations and thinks it never converged. how do I tell nls() that
I'm actually quite
happy with 3.760596e-29 and it need not work further?
thanks for any suggestions.
gary mcclelland (aka bernie)
univ of colorado
[[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.
Re: [R] optim error - repost
try:
vol - rep(c(0.03, 0.5, 2, 4, 8, 16, 32), 3)
time - rep(c(2,4,8),each=7)
p.mated - c(0.47, 0.48, 0.43, 0.43, 0.26, 0.23, NA, 0.68, 0.62, 0.64,
0.58, 0.53, 0.47, 0.24, 0.8, 0.79, 0.71, 0.56, 0.74, 0.8, 0.47)
eury - data.frame(vol=vol, time=time, p.mated=p.mated)
eury - na.omit(eury)
p0 - c(f=0.87, b=0.1, d=10)
eury.fit - function(p, time, vol, p.mated)
{
f - p[1]
b - p[2]
d - p[3]
xx - f*(1-exp(-b*time)-(b/(b-d/vol))*(exp(-d/vol*time)-exp(-b*time)))
## minimize sum of square differences
sum((p.mated - xx)^2)
}
eury.opt - optim(par=p0, fn=eury.fit, method = BFGS, hessian = TRUE,
time=eury$time, vol=eury$vol,
p.mated=eury$p.mated, control=list(trace=1))
## done with nls
eury.newfit1 - nls(p.mated ~ f * ( 1 -
exp(-b*time)-(b/(b-d/vol))*(exp(-d/vol*time)-exp(-b*time))),
data=eury, start=list(f=.87, b=0.1, d=10),
trace=TRUE)
coef(eury.newfit1)
eury.opt$par
On Sun, 1 Jun 2008 [EMAIL PROTECTED] wrote:
Here is a clean version. I did this with nls and it works (see below), but
I need to do it with optim. Keun-Hyung
# optim
vol-rep(c(0.03, 0.5, 2, 4, 8, 16, 32), 3)
time-rep(c(2,4,8),each=7)
p.mated-c(0.47, 0.48, 0.43, 0.43, 0.26, 0.23, NA, 0.68, 0.62, 0.64, 0.58,
0.53, 0.47,
0.24, 0.8, 0.79, 0.71, 0.56, 0.74, 0.8, 0.47)
eury-data.frame(vol=vol, time=time, p.mated=p.mated)
eury-na.omit(eury); eury
p0- c(f=0.87, b=0.1, d=10)
eury.fit - function (f, time, vol)
{
f-p[1]; b-p[2]; d-p[3]
p.mated = p[1] * ( (1 -
exp(-p[2]*time))-(p[2]/(p[2]-(p[3]/vol)))
* (exp(-p[3]/vol*time)-exp(-p[2]*time)))
}
eury.opt- optim(p0, fn=eury.fit, NULL, method = BFGS, hessian = TRUE)
# I received the following error message.
Error in fn(par, ...) : argument time is missing, with no default
## done with nls - this works
eury.newfit1 - nls(p.mated ~ f * ( 1 -
exp(-b*time)-(b/(b-d/vol))*(exp(-d/vol*time)-exp(-b*time))),
data=eury, start=list(f=.87,
b=0.1, d=10))
v - log10(range(eury$vol))
y - expand.grid(vol=10^seq(min(v), max(v), length=100), time=c(2,4,8))
y$pred.mate.new - predict(eury.newfit1,y)
plot (eury$vol, eury$p.mated, type=n, log=x, xlab=Volume L, ylab =
Fraction Mating)
for (i in c(2, 4, 8))
{
points(eury$vol[eury$time==i], eury$p.mated[eury$time==i], pch=16,
col=(i/2)+1)
lines(y$vol[y$time==i], y$pred.mate.new[y$time==i], lwd=3, col=i)
}
legend(.005,.2, c( 2h,4h,8h), col=c(2,4,8), lwd=3)
[[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.
Re: [R] nls diagnostics?
Dear Spencer,
I just saw your post.
If the singular gradient happens during or after iteration one (that is,
not at the initial estimates), then calling summary on the nls output
would give standard error estimates on the parameters useful for
diagnostics. You could also call chol2inv(xx$m$Rmat()) where xx is the
object returned by nls to get an estimate of the inverse of the hessian;
you could use this estimate to proceed with the diagnostics you were
discussing.
It would be possible (and in my opinion, desirable) to modify nls to
return an object that contains the information above even if the singular
gradient is at the intial estimates, but I don't think this can be
accomplished without quite a few changes.
You could also use nls.lm from the package minpack.lm to get a hessian
estimate.
By the way, I liked your idea from a while back
(https://stat.ethz.ch/pipermail/r-devel/2008-April/048984.html) to do some
diagnostics to identify problems with overparameterization automatically.
best,
Kate
On Fri, 23 May 2008, Spencer Graves wrote:
Hi, All:
What tools exist for diagnosing singular gradient problems with
'nls'? Consider the following toy example:
DF1 - data.frame(y=1:9, one=rep(1,9))
nlsToyProblem - nls(y~(a+2*b)*one, DF1, start=list(a=1, b=1),
control=nls.control(warnOnly=TRUE))
Error in nlsModel(formula, mf, start, wts) :
singular gradient matrix at initial parameter estimates
This example is obviously stupid, but other singular gradient
problems are not so obvious.
If we transfer this problem to 'optim', we can get diagnostics
from an eigen analysis of the hessian:
dumfun - function(x, y, one){
d - y-(x[1]+2*x[2])*one
sum(d^2)
}
optimToyProblem - optim(c(a=1, b=1), dumfun, hessian=TRUE,
y=DF1$y, one=DF1$one)
eigen(optimToyProblem$hessian, symmetric=TRUE)
$values
[1] 9.00e+01 -7.105427e-10
$vectors
[,1] [,2]
[1,] 0.4472136 -0.8944272
[2,] 0.8944272 0.4472136
The smallest eigenvalue is essentially numerically zero relative
to the largest,confirming the 'singular gradient' message. The
corresponding eigenvector helps diagnose the problem: Adding (-0.9,
0.45)*z to any solution gives another equally good solution, for any z.
I've used this technique to diagnose many subtle convergence
problems with 'optim'. Are tools of this nature available for 'nls'?
Thanks,
Spencer Graves
__
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] function in nls argument -- robust estimation
Dear Martin, Thanks for the ideas regarding the relation of what Fernando is doing with robust regression. Indeed, it's an important point that he can't consider the standard error estimates on his parameters correct. I know from discussion off-list that he's happy with the results he has now; nevertheless the robust regression route may be an interesting alternative. I'm posting a scipt to R-SIG-robust now that compares the 3 ways (nls, nlrob and nls.lm w/residuals above a certain quantile set to zero). best, Kate On Sat, 10 May 2008, Martin Maechler wrote: Hi Kate and Fernando, I'm late into this thread, but from reading it I get the impression that Fernando really wants to do *robust* (as opposed to least-squares) non-linear model fitting. His proposal to set residuals to zero when they are outside a given bound is a very special case of an M-estimator, namely (if I'm not mistaken) the so-called Huber skipped-mean, an M-estimator with psi-function psi - function(x, k) ifelse(abs(x) = k, x, 0) It is known that this can be far from optimal, and either using Huber-psi or a redescender such as Tukey's biweight can be considerably better. Also note that the standard inference (std.errors, P-values, ...) that you'd get from summary(nlsfit) or anova(nls1, nl2) is *invalid* here, since you are effectively using *random* weighting. The nlrob() function in package 'robustbase' implements M-estimation of nonlinear models directly. Unfortunately, how to do correct inference in this situation is a hard problem, probably even an open research question in parts. I would expect that the bootstrap should work if you only have a few outliers. I don't have time at the moment to look at the example data and the model, and show you how to use it for nlrob(); if you find a way to you it for nls() , then the same should work for nlrob(). I'm CCing this to the specialists for Robust Stats with R mailing list, R-SIG-robust. Best regards, Martin Maechler ETH Zurich KateM == Katharine Mullen [EMAIL PROTECTED] on Fri, 9 May 2008 15:50:08 +0200 (CEST) writes: KateM You can take minpack.lm_1.1-0 (source code and MS Windows build, KateM respectively) from here: KateM http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.tar.gz KateM http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.zip KateM The bug that occurs when nprint = 0 is fixed. Also fixed is another KateM problem suggested your example: when the argument par is a list, calling KateM summary on the output of nls.lm was not working. KateM I'll submit the new version to CRAN soon. KateM This disscusion has been fruitful - thanks for it. KateM On Fri, 9 May 2008, Katharine Mullen wrote: You indeed found a bug. I can reproduce it (which I should have tried to do on other examples in the first place!). Thanks for finding it. It will be fixed in version 1.1-0 which I will submit to CRAN soon. On Fri, 9 May 2008, elnano wrote: Find the data (data_nls.lm_moyano.txt) here: ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano Katharine Mullen wrote: Thanks for the details - it sounds like a bug. You can either send me the data in an email off-list or make it available on-line somewhere, so that I and other people can download it. __ 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://www.nabble.com/function-in-nls-argument-tp17108100p17146812.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. KateM __ KateM R-help@r-project.org mailing list KateM https://stat.ethz.ch/mailman/listinfo/r-help KateM PLEASE do read the posting guide http://www.R-project.org/posting-guide.html KateM and provide commented, minimal, self-contained, reproducible code. __ R-help@r
Re: [R] function in nls argument
Thanks for the details - it sounds like a bug. You can either send me the
data in an email off-list or make it available on-line somewhere, so that
I and other people can download it.
On Fri, 9 May 2008, elnano wrote:
Thank you Katharine. I am certain nprint is affecting my solution. Let me
know how I can send the data (~300Kb). The script I used it:
ST1 - ST04
SM1 - SM08
SR1 - SRch2
ST - ST1 [!is.na(SR1)]
SM - SM1 [!is.na(SR1)]
SR - SR1 [!is.na(SR1)]
q - 0.90
p - c(a=-0.003, b=0.13, c=0.50, E=400)
SRf - function(ST, SM, a, b, c, E)
{
expr -
expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)
eval(expr)
}
optim.f - function(p, ST, SM, SR, SRfcall)
{
res - SR - do.call(SRfcall, c(list(ST = ST, SM = SM), as.list(p)))
abserr - abs(res)
qnum - quantile(abserr, probs=q, na.rm=TRUE)
res[abserr qnum]=0
res
}
SRmodel- nls.lm(par = p, fn = optim.f, SRfcall = SRf, ST = ST, SM = SM,
SR = SR, control = nls.lm.control(nprint=1))
summary(SRmodel)
SRpred - do.call(SRf, c(list(ST = ST1, SM = SM1), SRmodel$par))
plot(SR1~SRpred)
a=c(0,2,4,6)
b=c(0,2,4,6)
lines(a,b,col=4)
Changing the nprint argument to 0 (or removing nprint) results in different
and completely wrong solutions. The same is true for this equivalent
simplified script from your example.
ST1 - ST04
SM1 - SM08
SR1 - SRch2
ST - ST1 [!is.na(SR1)]
SM - SM1 [!is.na(SR1)]
SR - SR1 [!is.na(SR1)]
q - 0.9
p - list(a=-0.003, b=0.13, c=0.50, E=400)
optim.f - function(p, ST, SM, SR, q) {
res - SR -
(p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13
abserr - abs(res)
qnum - quantile(abserr, probs=q, na.rm=TRUE)
res[abserr qnum] - 0
res
}
SRmodel - nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, SR = SR, q = q,
control = nls.lm.control(nprint=1))
summary(SRmodel)
SRfun - function(ST, SM, a, b, c, E)
{(a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13}
SRpred - do.call(SRfun, c(list(ST = ST1, SM = SM1), SRmodel$par))
plot(SR1~SRpred)
a=seq(0,6,1)
b=seq(0,6,1)
lines(a,b,col=4)
Here, however, I get the following additional error after summary(SRmodel):
Error in param/se : non-numeric argument to binary operator
although the solution is same as for the first script.
Note that a q of 95 is OK, a q of 90 is not, but a q of 50 is OK again...
To make your example reproducible you have to provide the data somehow;
I am pretty sure nprint doesn't effect the solution, but if it does this
would be a bug and I would appreciate a reproducible report.
The example in nls.lm is a little complicated in order to show how to use
an analytical expression for the gradient (possible to provide in the
argument jac); since you don't need this, note that your residual function
can be simplified into something along the lines of
p - list(a=-0.003, b=0.13, c=0.50, E=400)
optim.f - function(p, ST, SM, R, q) {
res - R
-(p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))
abserr - abs(res)
qnum - quantile(abserr, probs=q, na.rm=T)
res[abserr qnum] - 0
res
}
Rmodel - nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, R = R, q = q)
On Thu, 8 May 2008, Fernando Moyano wrote:
I solved the problem arising from using certain quantile values simply
by printing the iterations with the argument nprint. This gave me
correct estimates. I have no idea why.
--
View this message in context:
http://www.nabble.com/function-in-nls-argument-tp17108100p17145801.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.
Re: [R] function in nls argument
You indeed found a bug. I can reproduce it (which I should have tried to do on other examples in the first place!). Thanks for finding it. It will be fixed in version 1.1-0 which I will submit to CRAN soon. On Fri, 9 May 2008, elnano wrote: Find the data (data_nls.lm_moyano.txt) here: ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano Katharine Mullen wrote: Thanks for the details - it sounds like a bug. You can either send me the data in an email off-list or make it available on-line somewhere, so that I and other people can download it. __ 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://www.nabble.com/function-in-nls-argument-tp17108100p17146812.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.
Re: [R] function in nls argument
You can take minpack.lm_1.1-0 (source code and MS Windows build, respectively) from here: http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.tar.gz http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.zip The bug that occurs when nprint = 0 is fixed. Also fixed is another problem suggested your example: when the argument par is a list, calling summary on the output of nls.lm was not working. I'll submit the new version to CRAN soon. This disscusion has been fruitful - thanks for it. On Fri, 9 May 2008, Katharine Mullen wrote: You indeed found a bug. I can reproduce it (which I should have tried to do on other examples in the first place!). Thanks for finding it. It will be fixed in version 1.1-0 which I will submit to CRAN soon. On Fri, 9 May 2008, elnano wrote: Find the data (data_nls.lm_moyano.txt) here: ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano Katharine Mullen wrote: Thanks for the details - it sounds like a bug. You can either send me the data in an email off-list or make it available on-line somewhere, so that I and other people can download it. __ 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://www.nabble.com/function-in-nls-argument-tp17108100p17146812.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-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] function in nls argument
To make your example reproducible you have to provide the data somehow;
I am pretty sure nprint doesn't effect the solution, but if it does this
would be a bug and I would appreciate a reproducible report.
The example in nls.lm is a little complicated in order to show how to use
an analytical expression for the gradient (possible to provide in the
argument jac); since you don't need this, note that your residual function
can be simplified into something along the lines of
p - list(a=-0.003, b=0.13, c=0.50, E=400)
optim.f - function(p, ST, SM, R, q) {
res - R
-(p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))
abserr - abs(res)
qnum - quantile(abserr, probs=q, na.rm=T)
res[abserr qnum] - 0
res
}
Rmodel - nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, R = R, q = q)
On Thu, 8 May 2008, Fernando Moyano wrote:
I solved the problem arising from using certain quantile values simply
by printing the iterations with the argument nprint. This gave me
correct estimates. I have no idea why.
- Original Message
From: elnano [EMAIL PROTECTED]
To: r-help@r-project.org
Sent: Thursday, May 8, 2008 5:43:31 PM
Subject: Re: [R] function in nls argument
I've basically solved the problem using the nls.lm function from the
minpack.lm (thanks Katharine) with some modifications for ignoring residuals
above a given percentile. This is to avoid the strong influence of points
which push my modeled vs. measured values away from the 1:1 line.
I based it on the example given for nls.lm. Here it is:
R # soil respiration data
ST - ST [!is.na(R)] # soil temeprature data. Had to remove na to make
nls.lm work
SM - SM [!is.na(R)]# soil moisture data
R - R [!is.na(R)]
q - 0.95# quantile
p - c(a=-0.003, b=0.13, c=0.50, E=400)# model parameters with
estimated values
Rf - function(ST, SM, a, b, c, E)
{
expr -
expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)
eval(expr)
}
optim.f - function(p, ST, SM, R, Rfcall)
{
res - R - do.call(Rfcall, c(list(ST = ST, SM = SM), as.list(p))) #
the nls.lm example divides this by sqrt(R), I don't know why. I removed
that.
abserr - abs(res)
qnum - quantile(abserr, probs=q, na.rm=T)# calculate the q
quantile of the absolute errors
res[abserr qnum]=0 # convert
residuals above qnum to 0
res
}
Rmodel- nls.lm(par = p, fn = optim.f, Rfcall = Rf, ST = ST, SM = SM, R
= R)
summary(Rmodel)
The only error I still get is when using lower q values. A q of around 0.95
or less (depending on the dataset) gives me completely wrong parameter
estimates resulting in negative predicted values. Maybe someone has a
suggestion here.
Fernando
Katharine Mullen wrote:
The error message means that the gradient (first derivative of residual
vector with respect to the parameter vector) is not possible to work with;
calling the function qr on the gradient multiplied by the square root of
the weight vector .swts (in your case all 1's) fails.
If you want concrete advice it would be helpful to provide the commented,
minimal, self-contained, reproducible code that the posting guide
requests. what are the values of ST04, SM08b, ch2no, and tower?
General comments: If your goal is to minimize sum( (observed -
predicted)^2), the function you give nls to minimize (optim.fun in your
case) should return the vector (observed - predicted), not the scalar sum(
(observed - predicted)^2). You may want to see the nls.lm function in
package minpack.lm for another way to deal with the problem.
On Wed, 7 May 2008, Fernando Moyano wrote:
Greetings R users, maybe there is someone who can help
me with this problem:
I define a function optim.fun and want as output the
sum of squared errors between predicted and measured
values, as follows:
optim.fun - function (ST04, SM08b, ch2no, a, b, d, E)
{
predR -
(a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13
abserr - abs(ch2no-predR)
qnum - quantile(abserr, probs=0.95, na.rm=T)
is.na(abserr) - (abserr qnum)
errsq - sum(abserr^2, na.rm=T)
errsq
}
Then I want to optimize parameters a,b,d and E as to
minimize the function output with the following:
optim.model-nls(nulldat ~ optim.fun(ST04, SM08b,
ch2no, a, b, d, E), data=tower,
start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action =
na.exclude )
were nulldat=0
At this point I get the following error message:
Error in qr.default(.swts * attr(rhs, gradient)) :
NA/NaN/Inf in foreign function call (arg 1)
Warning messages:
1: In if (na.rm) x - x[!is.na(x)] else if
(any(is.na(x))) stop(missing values and NaN's not
allowed if 'na.rm' is FALSE) ... :
the condition has length 1
Re: [R] NLS plinear question
recall that 0 ^{-.2} = 1/0^{.2}, and that dividing by 0 gives Inf.
so when 0 is in trl, part of your model for RT is Inf:
trl - 0:14
p - -.2
cbind(1,trl, trl^p)
trl
[1,] 1 0 Inf
[2,] 1 1 1.000
[3,] 1 2 0.8705506
[4,] 1 3 0.8027416
[5,] 1 4 0.7578583
[6,] 1 5 0.7247797
[7,] 1 6 0.6988271
[8,] 1 7 0.6776109
[9,] 1 8 0.6597540
[10,] 1 9 0.6443940
[11,] 1 10 0.6309573
[12,] 1 11 0.6190439
[13,] 1 12 0.6083643
[14,] 1 13 0.5987029
[15,] 1 14 0.5898946
On Tue, 6 May 2008, Rick DeShon wrote:
Hi All.
I've run into a problem with the plinear algorithm in nls that is confusing
me.
Assume the following reaction time data over 15 trials for a single unit.
Trials are coded from 0-14 so that the intercept represents reaction time in
the first trial.
trl RT
01132.0
1 630.5
21371.5
3 704.0
4 488.5
5 575.5
6 613.0
7 824.5
8 509.0
9 791.0
10 492.5
11 515.5
12 467.0
13 556.5
14 456.0
Now fit a power function to this data using nls with the plinear algorithm
fit.pw -nls(RT ~ cbind(1,trl, trl^p), start = c(p = -.2), algorithm =
plinear, data=df.one)
Yields the following error message
Error in numericDeriv(form[[3]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
Now, recode trial from 1-15 and run the same model.
fit.pw -nls(RT ~ cbind(1,trl, trl^p), start = c(p = -.2), algorithm =
plinear, data=df.one)
Seems to work fine now...
Nonlinear regression model
model: RT ~ cbind(1, trl, trl^p)
data: df.one
p .lin1.lin.trl .lin3
-0.2845 200.3230-8.9467 904.7582
residual sum-of-squares: 555915
Number of iterations to convergence: 11
Any idea why having a zero for the first value of X causes this problem?
Thanks in advance,
Rick DeShon
[[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.
Re: [R] rzinb (VGAM) and dnbinom in optim - Solved.
Did you notice that the residual function I gave flipped what your
original example problem used as starting values for munb and size? so you
can sometimes afford to give a bad starting value for these param.
in trying to calculate the gradient L-BFGS-B may try evaluate fn with
parameter values outside of the box constraints; to categorically avoid
problems you'd have to know the max. step size for each parameter. it
seems better to use a try statement to catch problems and impose a big
penalty if the try fails, as in
library(VGAM)
phi - 0
munb - 72
size - 0.7
x - rzinb(200, phi=phi, size=size, munb=munb)
## x should have some non-zero elements
fit - vglm(x~1, zinegbinomial,trace=TRUE)
fncZiNB - function(xx, x){
phi - xx[1]
munb - xx[2]
size - xx[3]
yy - try(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))
if(class(yy) == try-error || !all(is.finite(yy))) 1e10 else -sum(yy)
}
# size, the 3rd param must be strictly 0
solut - optim(fn=fncZiNB,par=c(0, 72, .7),x=x,method =
L-BFGS-B, lower=c(0,0,1e-10),upper=c(1,Inf,Inf),
control=list(trace=TRUE),hessian=TRUE)
##optim
solut$par
##vglm
Coef(fit)
On Sat, 19 Apr 2008, Tim Howard wrote:
Thank you very much for the help and the great suggestion on cleaning up
the residual function.
Setting the bounds and being very careful with the starting position (on
my real data) do seem to do the trick, which I wouldn't have guessed
given the error message that was getting thrown.
Best,
Tim
Thomas Yee [EMAIL PROTECTED] 4/18/2008 6:02 PM
Hi,
using something like
lower=c(1e-5,1e-5,1e-5),upper=c(1-1e-5,Inf,Inf)
is even more safer since it sometime evaluates the function slightly
outside the lower and upper limits.
cheers
Thomas
Katharine Mullen wrote:
I'm not going to look into what's happening in-depth but it looks like the
bounds for your parameters need to be set with care; the below (with
slight re-def. of your residual function) gives results that seem to match
vglm approximately, with 1e-10 standing in for a bound of 0 and 1-e-10
standing in for a bound of 1.
library(VGAM)
phi - 0.09
size - 0.7
munb - 72
x - rzinb(200, phi=phi, size=size, munb=munb)
fit - vglm(x~1, zinegbinomial,trace=TRUE)
fncZiNB - function(xx, x){
phi - xx[1]
size - xx[2]
munb - xx[3]
-sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))
}
solut - optim(fn=fncZiNB,par=c(0.09, 72, 0.7),x=x,method =
L-BFGS-B, lower=c(1e-10,1e-10,1e-10),upper=c(1-1e-10,Inf,Inf),
control=list(trace=TRUE),hessian=TRUE)
##optim
solut$par
##vglm
Coef(fit)
On Fri, 18 Apr 2008, Tim Howard wrote:
Dear R-help gurus (and T.Yee, the VGAM maintainer) -
I've been banging my head against the keyboard for too long now, hopefully
someone can pick up on the errors of my ways...
I am trying to use optim to fit a zero-inflated negative binomial
distribution. No matter what I try I can't get optim to recognize my
initial parameters. I think the problem is that dnbinom allows either
'prob' or 'mu' and I am trying to pass it 'mu'. I've tried a wrapper
function but it still hangs. An example:
#set up dummy data,load VGAM, which is where I am getting the zero-inflated
neg binom
#I get the same problem without VGAM, using dnbinom in the wrapper instead.
library(VGAM)
phi - 0.09
size - 0.7
munb - 72
x - rzinb(200, phi=phi, size=size, munb=munb)
#VGAM can fit using vglm... but I'd like to try optim
fit - vglm(x~1, zinegbinomial,trace=TRUE)
VGLMlinear loop 1 : loglikelihood = -964.1915
VGLMlinear loop 2 : loglikelihood = -964.1392
VGLMlinear loop 3 : loglikelihood = -964.1392
Coef(fit)
phi munb k
0.1200317 62.8882874 0.8179183
#build my wrapper function for dzinb
fncZiNB - function(phi, size, munb){
+ -sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))}
#test the wrapper, seems to work.
fncZiNB(phi=0.09, size=0.7, munb=72)
[1] 969.1435
#try it in optim
solut - optim(fn=fncZiNB,par=c(phi=0.09, munb=72, size=0.7),method =
L-BFGS-B, hessian=TRUE)
Error in dnbinom(x = x, size = size, prob = prob, log = log.arg) :
argument size is missing, with no default
I have the same problem with dnbinom.
I'm sure I'm doing something obvious any suggestions?
Session info below. Thanks in advance.
Tim Howard
New York Natural Heritage Program
sessionInfo()
R version 2.6.2 (2008-02-08)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats4splines stats graphics grDevices utils datasets
methods
[9] base
other attached packages:
[1] VGAM_0.7-6 randomForest_4.5-23
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman
Re: [R] rzinb (VGAM) and dnbinom in optim
I'm not going to look into what's happening in-depth but it looks like the
bounds for your parameters need to be set with care; the below (with
slight re-def. of your residual function) gives results that seem to match
vglm approximately, with 1e-10 standing in for a bound of 0 and 1-e-10
standing in for a bound of 1.
library(VGAM)
phi - 0.09
size - 0.7
munb - 72
x - rzinb(200, phi=phi, size=size, munb=munb)
fit - vglm(x~1, zinegbinomial,trace=TRUE)
fncZiNB - function(xx, x){
phi - xx[1]
size - xx[2]
munb - xx[3]
-sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))
}
solut - optim(fn=fncZiNB,par=c(0.09, 72, 0.7),x=x,method =
L-BFGS-B, lower=c(1e-10,1e-10,1e-10),upper=c(1-1e-10,Inf,Inf),
control=list(trace=TRUE),hessian=TRUE)
##optim
solut$par
##vglm
Coef(fit)
On Fri, 18 Apr 2008, Tim Howard wrote:
Dear R-help gurus (and T.Yee, the VGAM maintainer) -
I've been banging my head against the keyboard for too long now, hopefully
someone can pick up on the errors of my ways...
I am trying to use optim to fit a zero-inflated negative binomial
distribution. No matter what I try I can't get optim to recognize my initial
parameters. I think the problem is that dnbinom allows either 'prob' or 'mu'
and I am trying to pass it 'mu'. I've tried a wrapper function but it still
hangs. An example:
#set up dummy data,load VGAM, which is where I am getting the zero-inflated
neg binom
#I get the same problem without VGAM, using dnbinom in the wrapper instead.
library(VGAM)
phi - 0.09
size - 0.7
munb - 72
x - rzinb(200, phi=phi, size=size, munb=munb)
#VGAM can fit using vglm... but I'd like to try optim
fit - vglm(x~1, zinegbinomial,trace=TRUE)
VGLMlinear loop 1 : loglikelihood = -964.1915
VGLMlinear loop 2 : loglikelihood = -964.1392
VGLMlinear loop 3 : loglikelihood = -964.1392
Coef(fit)
phi munb k
0.1200317 62.8882874 0.8179183
#build my wrapper function for dzinb
fncZiNB - function(phi, size, munb){
+ -sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))}
#test the wrapper, seems to work.
fncZiNB(phi=0.09, size=0.7, munb=72)
[1] 969.1435
#try it in optim
solut - optim(fn=fncZiNB,par=c(phi=0.09, munb=72, size=0.7),method =
L-BFGS-B, hessian=TRUE)
Error in dnbinom(x = x, size = size, prob = prob, log = log.arg) :
argument size is missing, with no default
I have the same problem with dnbinom.
I'm sure I'm doing something obvious any suggestions?
Session info below. Thanks in advance.
Tim Howard
New York Natural Heritage Program
sessionInfo()
R version 2.6.2 (2008-02-08)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats4splines stats graphics grDevices utils datasets
methods
[9] base
other attached packages:
[1] VGAM_0.7-6 randomForest_4.5-23
__
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] using rbind() on multiple objects at once
do.call(rbind, list.foo) On Fri, 18 Apr 2008, Andrew Yee wrote: Is there an efficient way to use rbind() with the five dataframes described in the following example: a - c(1:5) list.foo - lapply(a, function(x) data.frame(beta=a*rnorm(10), deta=a*rnorm(10))) big.data.frame - rbind(list.foo[[1]], list.foo[[2]], list.foo[[3]], list.foo[[4]], list.foo[[5]]) #is there an easier method? For example, I naively thought you could do something like rbind(list.foo[[1:5]]) #gives an error message but that results in an error message. Thanks, Andrew [[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.
Re: [R] Suggestions: Terminology Pkgs for following spectra over time
Dear Bryan, For work in that general direction, look at the R News that was on 'R in chemistry' (Volume 6/3, August 2006, http://www.r-project.org/doc/Rnews/Rnews_2006-3.pdf) and the special vol. of the Journal of Statistical Software on 'Spectroscopy and Chemometrics in R' (vol 18, http://www.jstatsoft.org/v18) For ideas about modeling spectra resolved in time and frequency, the papers of Sabine van Huffel at K U Leuven may be of interest (http://www.esat.kuleuven.be/~sistawww/cgi-bin/newsearch.pl?Name=Van+Huffel+S) If you decide that you want to describe the time domain with a parametric model, then look at the package TIMP (and perhaps contact me off list). On Wed, 16 Apr 2008, Bryan Hanson wrote: Hi Folks... No code to troubleshoot here. I need some suggestions about the right terminology to use in further searching, and any suggestions about R pkgs that might be appropriate. I am in the planning stages of a project in which IR, NMR and other spectra (I'm a chemist) would be collected on various samples, and individual samples would be followed over time. The spectra will be feature rich/complex, so one can't see the changes by visual inspection. The spectra are basically 2D matrices: peaks as a function of frequencies. So the data set is in the form of spectra of a single sample over time, for multiple samples. I am wondering about methods R pkgs that can be used to analyze changes in the spectra over time. For instance, I would like to find specific peaks that are changing over time, sets of peaks that are changing in a correlated way over time etc. I'd like to do this in an efficient and statistically valid way. What I am thinking of is somewhat like a time series, somewhat like image analysis (but only 2D), but it's not quite either of those and I need to know what it's really called to investigate further. Any suggestions as to R pkgs and key words/phrases will be appreciated. TIA, Bryan * Bryan Hanson Professor of Chemistry Biochemistry DePauw University, Greencastle Indiana USA __ 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] replace NULL with NA in matrix
if your matrix is stored in a text file then xx - read.table(x, na.strings=NULL) where x is the name of the text file. On Thu, 17 Apr 2008, Tim Smith wrote: Hi, I had a matrix with NULL values, which I wanted to replace with NA. Is there an efficient way to do this? Small sample input matrix: A B C D E 1 5222.18 6355.10 4392.68 2607.41 4524.09 2NULL 257.33NULL 161.51 119.44 3NULL 274.80 305.28 443.27NULL 4 1759.76 1556.45 1224.06 1558.73 1837.09 Tim Be a better friend, newshound, and [[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.
Re: [R] Singular Gradient in nls
sorry, the below should read: A QR decomposition is done on the weighted gradient matrix; if the estimate of the rank that results is less than min( the number of columns in the gradient (the number of nonlinear parameters), the number of rows (the number of observations) ) , nls stops. so you look at the number of observations, and the number of parameters, and take the smaller of the two. if the estimated rank of the weighted gradient is small than this number, you stop. On Fri, 28 Mar 2008, Katharine Mullen wrote: A QR decomposition is done on the weighted gradient matrix; if the estimate of the rank that results is less than the number of columns in the gradient (the number of nonlinear parameters), or less than the number of rows (the number of observations), nls stops. You can see the calls in the source code of nlsModel (https://svn.r-project.org/R/trunk/src/library/stats/R/nls.R). On Fri, 28 Mar 2008, glenn andrews wrote: //Referring to the response posted many years ago, copied below, what is the specific criterium used for singularity of the gradient matrix? Is a Singular Value Decomposition used to determine the singular values? Is it the gradient matrix condition number or some other criterion for determining singularity? // //Glenn // / / / What does the error 'singular gradient' mean during a nonlinear regression? / The gradient matrix to which the message refers is the derivative of the vector of predicted values with respect to the vector of parameters at the current parameter estimates. If you have 20 observations and three parameters, this will be a matrix with 20 rows and three columns. For the model to be estimable in a region of the current estimates, this matrix must have full column rank. When it fails to have full column rank the singular gradient message is given and the iterations stop. Generally this indicates that the model is overparameterized or that the starting estimates were poorly chosen. Try using trace = TRUE in the call to nls and watching the progress of the iterations. This will often show that the estimates are wandering into unreasonable regions of the parameter space. -- Douglas Bates[EMAIL PROTECTED] Statistics Department608/262-2598 University of Wisconsin - Madisonhttp://www.stat.wisc.edu/~bates/ http://www.stat.wisc.edu/%7Ebates/ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html http://www.ci.tuwien.ac.at/%7Ehornik/R/R-FAQ.html Send info, help, or [un]subscribe (in the body, not the subject !) To: [EMAIL PROTECTED] _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ [[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-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] Singular Gradient in nls
A QR decomposition is done on the weighted gradient matrix; if the estimate of the rank that results is less than the number of columns in the gradient (the number of nonlinear parameters), or less than the number of rows (the number of observations), nls stops. You can see the calls in the source code of nlsModel (https://svn.r-project.org/R/trunk/src/library/stats/R/nls.R). On Fri, 28 Mar 2008, glenn andrews wrote: //Referring to the response posted many years ago, copied below, what is the specific criterium used for singularity of the gradient matrix? Is a Singular Value Decomposition used to determine the singular values? Is it the gradient matrix condition number or some other criterion for determining singularity? // //Glenn // / / / What does the error 'singular gradient' mean during a nonlinear regression? / The gradient matrix to which the message refers is the derivative of the vector of predicted values with respect to the vector of parameters at the current parameter estimates. If you have 20 observations and three parameters, this will be a matrix with 20 rows and three columns. For the model to be estimable in a region of the current estimates, this matrix must have full column rank. When it fails to have full column rank the singular gradient message is given and the iterations stop. Generally this indicates that the model is overparameterized or that the starting estimates were poorly chosen. Try using trace = TRUE in the call to nls and watching the progress of the iterations. This will often show that the estimates are wandering into unreasonable regions of the parameter space. -- Douglas Bates[EMAIL PROTECTED] Statistics Department608/262-2598 University of Wisconsin - Madisonhttp://www.stat.wisc.edu/~bates/ http://www.stat.wisc.edu/%7Ebates/ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html http://www.ci.tuwien.ac.at/%7Ehornik/R/R-FAQ.html Send info, help, or [un]subscribe (in the body, not the subject !) To: [EMAIL PROTECTED] _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ [[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] [R-pkgs] new version of minpack.lm
The package minpack.lm allows nonlinear regression problems to be addressed with a modification of the Levenberg-Marquardt algorithm based on the implementation of 'lmder' and 'lmdif' in MINPACK. Version 1.0-8 of the package is now available on CRAN. Changes in version 1.0-8 include: o possibility to obtain standard error estimates on the parameters via new methods for the generic functions 'summary' and 'vcov' o possibility to extract other information via new methods for the generic functions 'coef', 'deviance', 'df.residual', 'print', and 'residuals' o the argument 'control' of 'nls.lm' now defaults to 'nls.lm.control()'; 'nls.control.lm' allows a maximum number of iterations to be specified; when the element 'nprint' of the 'control' argument of a call to 'nls.lm' is an integer greater than 0, the residual sum of squares is now included in the information printed every 'nprint' iterations ` o the list returned by 'nls.lm' includes elements 'niter' and 'deviance' that represent the number of iterations performed and the residual sum of squares, respectively side-note on Levenberg-Marquardt (LM) versus Gauss-Newton (GN): There was some discussion (http://finzi.psych.upenn.edu/R/Rhelp02a/archive/108758.html) on Rhelp regarding whether one comes across real-world problems in which LM performs better than GN. I have been seeing such problems recently in some applications where GN as implemented in 'nls' reduces the step to a very small value, resulting in little change in the residual sum of squares from the starting values, whereas both NL2SOL applied via 'nls' called with 'algorithm=port' or LM as implemented in 'minpack.lm::nls.lm' significantly reduce the RSS. The implementation of NL2SOL is slower by a significant factor on these problems as compared to either the GN or LM implementations, making use of 'minpack.lm::nls.lm' attractive. Note that these problems may be considered pathological; there are issues with near collinearity of columns of the Jacobian and with the assumption that the residuals are Gaussian. Kate Mullen Timur Elzhov ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ 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
It is in the VR bundle (collection of packages), which is 'recommended', so it is included in all binary distributions of R; if you have the source code you find VR in R-2.6.2/src/library/Recommended On Sat, 8 Mar 2008, Svyatkovskiy Alexey wrote: Hello, Tell me please, can I get a sourcse code for 'lda' from MASS? Is MASS package still present in the latest R-2.6.2 version (I couldn't find it)? Thankful in advance, Alexey. __ 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] [OT] normal (as in Guassian)
There is some information and references regarding the name 'normal' in the internet article 'Earliest Known Uses of Some of the Words of Mathematics (N)', http://members.aol.com/jeff570/n.html, by John Aldrich. It contains the comment, Galton does not explain why he uses the term normal but the sense of conforming to a norm ( = 'A standard, model, pattern, type.' (OED)) seems implied. On Sun, 2 Mar 2008 [EMAIL PROTECTED] wrote: Hi Folks, Apologies to anyone who'd prefer not to see this query on this list; but I'm asking because it is probably the forum where I'm most likely to get a good answer! I'm interested in the provenance of the name normal distribution (for what I'd really prefer to call the Gaussian distribution). According to Wikipedia, The name normal distribution was coined independently by Charles S. Peirce, Francis Galton and Wilhelm Lexis around 1875. So be it, if that was the case -- but I would like to know why they chose the name normal: what did they intend to convey? As background: I'm reflecting a bit on the usage in statistics of everyday language as techincal terms, as in significantly different. This, for instance, is likely to be misunderstood by the general publidc when they encounter statements in the media. Likewise, normally distributed would probably be interpreted as distributed in the way one would normally expect or, perhaps, there was nothing unusual about the distribution. Comments welcome! With thanks, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 02-Mar-08 Time: 13:04:17 -- XFMail -- __ 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] install ncdf package on a 64-bit machine
I am running 64-bit Ubuntu 7.10 and unfortunately remember seeing that error message but not how I got it to go away. I would first try compiling netcdf-3.6.2 from source, without changing the default directories for installation. Looking at the error message it seems like the 3 shared libraries get made but then linking them gives the error - that step seems to be where it says to compile with -fPIC. My output of install.packages(ncdf) is below with that of sessionInfo() - I still see a lot of warnings. * Installing *source* package 'ncdf' ... checking for gcc... gcc -std=gnu99 checking for C compiler default output file name... a.out checking whether the C compiler works... yes checking whether we are cross compiling... no checking for suffix of executables... checking for suffix of object files... o checking whether we are using the GNU C compiler... yes checking whether gcc -std=gnu99 accepts -g... yes checking for gcc -std=gnu99 option to accept ANSI C... none needed checking how to run the C preprocessor... gcc -std=gnu99 -E checking for egrep... grep -E checking for ANSI C header files... yes checking for sys/types.h... yes checking for sys/stat.h... yes checking for stdlib.h... yes checking for string.h... yes checking for memory.h... yes checking for strings.h... yes checking for inttypes.h... yes checking for stdint.h... yes checking for unistd.h... yes checking netcdf.h usability... yes checking netcdf.h presence... yes checking for netcdf.h... yes Found netcdf.h in: . checking for nc_open in -lnetcdf... yes Found netcdf library file libnetcdf.a in directory . configure: creating ./config.status config.status: creating R/load.R config.status: creating src/Makevars ** libs gcc -std=gnu99 -I/usr/local/lib64/R/include -I/usr/local/lib64/R/include -I. -I/usr/local/include-fpic -g -O2 -c ncdf2.c -o ncdf2.o gcc -std=gnu99 -I/usr/local/lib64/R/include -I/usr/local/lib64/R/include -I. -I/usr/local/include-fpic -g -O2 -c ncdf3.c -o ncdf3.o ncdf3.c: In function R_nc_get_vara_charvarid: ncdf3.c:221: warning: assignment discards qualifiers from pointer target type ncdf3.c: In function R_nc_get_vara_numvarid: ncdf3.c:267: warning: assignment discards qualifiers from pointer target type gcc -std=gnu99 -I/usr/local/lib64/R/include -I/usr/local/lib64/R/include -I. -I/usr/local/include-fpic -g -O2 -c ncdf.c -o ncdf.o ncdf.c: In function R_nc_ttc_to_nctype: ncdf.c:424: warning: implicit declaration of function exit ncdf.c:424: warning: incompatible implicit declaration of built-in function exit gcc -std=gnu99 -shared -L/usr/local/lib64 -o ncdf.so ncdf2.o ncdf3.o ncdf.o -L. -lnetcdf ** R ** preparing package for lazy loading ** help Building/Updating help pages for package 'ncdf' Formats: text html latex example ancdf texthtmllatex aput.var.ncdf texthtmllatex example Note: unmatched right brace in 'att.get.ncdf' on or after line 6 att.get.ncdf texthtmllatex example att.put.ncdf texthtmllatex example close.ncdftexthtmllatex example create.ncdf texthtmllatex example dim.def.ncdf texthtmllatex example enddef.ncdf texthtmllatex example get.var.ncdf texthtmllatex example ncdf-internal texthtmllatex open.ncdf texthtmllatex example print.ncdftexthtmllatex example redef.ncdftexthtmllatex example set.missval.ncdf texthtmllatex example sync.ncdf texthtmllatex example var.add.ncdf texthtmllatex example var.def.ncdf texthtmllatex example Note: unmatched right brace in 'version.ncdf' on or after line 6 version.ncdf texthtmllatex ** building package indices ... * DONE (ncdf) The downloaded packages are in /tmp/RtmpcJ0Af9/downloaded_packages Updating HTML index of packages in '.Library' sessionInfo() R version 2.6.1 (2007-11-26) x86_64-unknown-linux-gnu locale: LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C attached base packages: [1] splines grid tcltk stats graphics grDevices utils [8] datasets methods base other attached packages: [1] ncdf_1.6 xtable_1.5-2 TIMP_1.5 minpack.lm_1.0-5 [5] quadprog_1.4-11 nnls_1.2 gclus_1.2gplots_2.3.2 [9] gdata_2.3.1 gtools_2.4.0 fields_3.5 vcd_1.0-6 [13]
Re: [R] Non linear regression with 2 explanatory variables
Dear Rolf,
One thing that sometimes makes nls easier to apply is using the 'formula'
argument like you would use the 'fn' argument of optim. That is, if you
have a residual function that has arguments x, y, a, b and you need to
optimize a and b, you would make a call like
nls(~resid(x,y,a=astart, b=bstart), control = nls.control(warnOnly =
TRUE, printEval = TRUE), start = list(a=astart, b=bstart))
This did not work easily before R-2.6.0, but does now. The Puromycin
analysis from the help files is an example of this useage and below is
another.
Or do you already use nls this way and still have problems?
# get data as a sum of exponentials
dataSumOfExp - function(rates = seq(.05, .005, length=3),
times = 1:100,
amps = rep(1, length(rates))) {
tfun - function(t,r) exp(-r*t)
## get C with tfun
C - mapply(tfun, r=rates, MoreArgs=list(t=times))
## add the columns of C with relative amplitudes 1, and add noise
C %*% amps + rnorm( nrow(C) ) * max(C) * .1
}
# residual function
resFun - function(rates, amps, measured, times = 1:100) {
tfun - function(t,r) exp(-r*t)
CEst - mapply(tfun, r=rates, MoreArgs=list(t=times))
measured - CEst %*% amps
}
# get data
measured - dataSumOfExp()
# optimize rates of exponentials and their relative amplitudes
res - nls(~resFun(rates = rates, measured = measured, amps = amps),
control = nls.control(warnOnly = TRUE, printEval = TRUE),
start = list(rates = c(.04, .1, .001),
amps = rep(1,3)), trace = TRUE)
summary(res)
On Thu, 17 Jan 2008, Rolf Turner wrote:
I have never had much success in using nls(). If you scan the archives
you will find one or two postings from me on this topic. I have
received
no useful responses to these postings.
I have found that anything that I tried (and failed) to do using nls()
could be done quite easily using optim().
cheers,
Rolf Turner
On 17/01/2008, at 3:56 AM, Gavin Simpson wrote:
hits=-2.6 tests=BAYES_00
X-USF-Spam-Flag: NO
On Wed, 2008-01-16 at 11:02 +0100, Janice Kielbassa wrote:
Hello!
I want to do a non-linear regression with 2 explanatory variables
(something like : length ~ a * time * exp( b* temperature)), having a
data set (length, time, temperature). Which function could I use (I
tried nls but I think it doesn't work)
Janice, I'll start by saying I can't help you as I have never used
nls()
myself and I am not familiar with this type of analysis.
Why do you think that nls() doesn't work? It is a widely used
part of
R and thus probably very well tested.
My understanding of these things is that nls is a sophisticated tool
that requires some effort on the part of the user, such as selecting
appropriate starting values.
##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
__
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] [R-pkgs] new package 'bvls', update of 'nnls'
A new package 'bvls' is available on CRAN. The package provides an R interface to the Stark-Parker algorithm for bounded-variable least squares (BVLS) that solves A x = b with the constraint l = x = u under least squares criteria, where l,x,u \in R^n, b \in R^m and A is an m \times n matrix. The Stark-Parker BVLS algorithm was published in Stark PB, Parker RL (1995). Bounded-variable least-squares: an algorithm and applications, Computational Statistics, 10, 129-141. The packages interfaces the Fortran77 code distributed via the statlib on-line software repository at Carnegie Mellon University (http://lib.stat.cmu.edu/general/bvls), modified very slightly for compatibility with the gfortran compiler. Stark and Parker have agreed to distribution under GPL version 2 or newer. The function 'bvls::bvls' returns an object of (S3) class 'bvls', which has methods for 'coefficients', 'fitted.values', 'deviance' and 'residuals'. Version 1.1 of the package 'nnls' is available on CRAN. Changes between Version 1.0 and 1.1: o The function 'nnls::nnls' returns an object of (S3) class 'nnls', which has methods for 'coefficients', 'fitted.values', 'deviance' and 'residuals' o The function 'nnnpls::nnnpls' allows each element of x to be constrained to either a non-positive or a non-negative value Katharine Mullen mail: Department of Physics and Astronomy, Faculty of Sciences Vrije Universiteit Amsterdam, de Boelelaan 1081 1081 HV Amsterdam, The Netherlands room: T.1.06 tel: +31 205987870 fax: +31 205987992 e-mail: [EMAIL PROTECTED] homepage: http://www.nat.vu.nl/~kate/ ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ 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] Organising tick-marks in plot()
something like this sounds like what you want:
labs-NA
x-200
for(i in 1:5){ labs-append(labs, c(rep(NA,3),x)); x-x+200 }
plot(1:900, xaxt=n, xlim=c(0, 1000))
axis(side = 1, at = seq(0,1000,by=50), labels=labs)
On Fri, 30 Nov 2007 [EMAIL PROTECTED] wrote:
Hi Folks,
I'm advising someone who's a beginner with R,
and therefore wants the simplest answer possible.
The issue is to produce a plot using
plot(x,y,...)
where, on the X-axis, the tick-marks should be
on the lines of:
-- Range of X-axis: 0:1000
-- tick-marks labelled 0,200,...,800,1000
-- unlabelled tick-marks every 50 from 0 to 1000.
regardless of the actual range of x-values (which
however would normally range over most of 0:1000).
I've looked at the Q-Q Plot example in the MASS
book (Section 3.4), which does achieve this kind
of effect, but the code is too complicated for
the present context. (Whatever code is used needs
to be readily changeable for different plots of
the same general kind).
Is there a way of supplying straightforward
arguments to plot() which would achieve this?
I've been reading and site-searching for a while,
and the best I can find is on the lines of
plot((x,y,pch=+,col=blue,xlim=c(0,1000),xaxt=n)
par(xaxp=c(0, 1000, 50))
axis(1)
which is already complicated enough in the present
context; but it doesn't do exactly what is required
since (for example) if the x-values range from
0 to 900 the labelled tick-marks are at
0 60 140 220 300 380 460 540 620 700 780 860 940
and have therefore been computed from the range of
the data, and not from the specified range of the x-axis.
Help please!
Best wishes to all,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 30-Nov-07 Time: 11:55:16
-- XFMail --
__
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] How is the Gauss-Newton method compared to Levenberg-Marquardt for curve-fitting?
Almost every problem I have seen in the context of fitting mixtures of exponentials/Gaussians has near-equal performance under LM or GN in terms of iterations/evaluations to convergence and the SSE of the solution found. (However, very recently I have encountered a problem with a large number of nonlinear parameters (~100) for which LM finds a better optimum than GN; I am in the process of looking into why this is the case.) Regarding profiling: I have never looked for such literature; if you find it, realize that the results may not generalize to the SSE surface of the residual functions you will consider. Therefore you might be better off doing the profiling yourself for the residual functions that are typical in your problem domain (you can use the R package minpack.lm for LM). For general discussion of LM and GN, standard referencs are Bates and Watts, Nonlinear regression analysis and its applications, and Seber and Wild, Nonlinear regression. I believe that one of the main reasons that GN gets more use than LM by R users is that the former is in base R, and for many simple problems, writing the function to minimize using a formula is handy, which only nls, but not minpack.lm, allows. I don't think it's the case that LM is hardly competitive in general. On Tue, 20 Nov 2007, Ali - wrote: Hi, It seems to me that the most suitable method in R for curve-fitting is the use of nls, which uses a Gauss-Newton (GN) algorithm, while the use of the Levenberg-Marquardt (LM) algorithm does not seem to be very stressed in R. According to this [1] by Ripley, 'Levenberg-Marquardt is hardly competitive these days' which could imply the low emphasize on LM in R. The position of LM is, to some extend, confusing. Bonnans et al [2] introduce the trust-region-based method of LM like this: 'This chapter is mostly devoted to methods which, although less universal than the preceding, are useful in a good number of cases. The frst one (trust-region) is actually extremely important, and might supersede line-searches, sooner or later.' The above should demonstrate the contradiction. Since some R developers are indeed the pioneers in the optimisation theory, I would like to ask for references involving profiling of various methods, including more modern techniques, with an application in general model-fitting. [1] http://tolstoy.newcastle.edu.au/R/help/00b/2492.html [2] Numerical Optimization, 2nd ed _ 100?s of Music vouchers to be won with MSN Music __ 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] specifying decimal places
Your calculations are done with a precision that depends on your machine. I guess that you want to know how to _print_ more decimal digits. See how many digits you print now with getOption(digits) Change this with, e.g., options(digits=12) On Mon, 12 Nov 2007, raymond chiruka wrote: hie how do I specify the number of decuimal places for my calulations thanks __ [[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.
Re: [R] strange `nls' behaviour
I can confirm this behavior on R-2.6.0 but don't have time to look into it
further at the moment.
On Mon, 12 Nov 2007, Joerg van den Hoff wrote:
I initially thought, this should better be posted to r-devel
but alas! no response. so I try it here. sory for the
lengthy explanation but it seems unavoidable. to quickly see
the problem simply copy the litte example below and execute
f(n=5)
which crashes. called with n != 5 (and of course n3 since
there are 3 parameters in the model...) everything is as it
should be.
in detail:
I stumbled over the follwing _very_ strange behaviour/error
when using `nls' which I'm tempted (despite the implied
dangers) to call a bug:
I've written a driver for `nls' which allows specifying the
model and the data vectors using arbitrary symbols. these
are internally mapped to consistent names, which poses a
slight complication when using `deriv' to provide analytic
derivatives. the following fragment gives the idea:
#-
f - function(n = 4) {
x - seq(0, 5, length = n)
y - 2 * exp(-1*x) + 2;
y - rnorm(y,y, 0.01*y)
model - y ~ a * exp (-b*x) + c
fitfunc - deriv(model[[3]], c(a, b, c), c(a, b, c, x))
#standard call of nls:
res1 - nls(y ~ fitfunc(a, b, c, x), start = c(a=1, b=1, c=1))
call.fitfunc -
c(list(fitfunc), as.name(a), as.name(b), as.name(c), as.name(x))
call.fitfunc - as.call(call.fitfunc)
frml - as.formula(y ~ eval(call.fitfunc))
#computed call of nls:
res2 - nls(frml, start = c(a=1, b=1, c=1))
list(res1 = res1, res2 = res2)
}
#-
the argument `n' defines the number of (simulated) data
points x/y which are going to be fitted by some model ( here
y ~ a*exp(-b*x)+c )
the first call to `nls' is the standard way of calling `nls'
when knowing all the variable and parameter names.
the second call (yielding `res2') uses a constructed formula
in `frml' (which in this example is of course not necessary,
but in the general case 'a,b,c,x,y' are not a priori known
names).
now, here is the problem: the call
f(4)
runs fine/consistently, as does every call with n 5.
BUT: for n = 5 (i.e. issuing f(5))
the second fit leads to the error message:
Error in model.frame(formula, rownames, variables, varnames, extras,
extranames, :
invalid type (language) for variable 'call.fitfunc'
I cornered this to a spot in `nls' where a model frame is
constructed in variable `mf'. the parsing/constructing here
seems simply to be messed up for n = 5: `call.fitfunc' is
interpreted as variable.
I, moreover, empirically noted that the problem occurs when
the total number of parameters plus dependent/independent
variables equals the number of data points (in the present
example a,b,c,x,y).
so it is not the 'magic' number of 5 but rather the identity
of data vector length and number of parameters+variables in
the model which leads to the problem.
this is with 2.5.0 (which hopefully is not considered
ancient) and MacOSX 10.4.10.
any ideas?
thanks
joerg
__
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] strange `nls' behaviour
On about line 547 in nls.R there is
mf$formula - # replace by one-sided linear model formula
as.formula(paste(~, paste(varNames[varIndex], collapse = +)),
env = environment(formula))
If this is replaced with
mf$formula - # replace by one-sided linear model formula
as.formula(paste(~, paste(varNames[1], collapse = +)),
env = environment(formula))
then the behviour seems to be OK at least for Joerg's examples.
It remains to be determined whether this is a _good_ solution (that is, if
it's a general solution to get the desired behavior).
On Mon, 12 Nov 2007, Martin Maechler wrote:
DM == Duncan Murdoch [EMAIL PROTECTED]
on Mon, 12 Nov 2007 07:36:34 -0500 writes:
DM On 11/12/2007 6:51 AM, Joerg van den Hoff wrote:
I initially thought, this should better be posted to r-devel
but alas! no response.
DM I think the reason there was no response is that your example is too
DM complicated. You're doing a lot of strange things (fitfunc as a
result
DM of deriv, using as.name, as.call, as.formula, etc.) You should
simplify
DM it down to isolate the bug. Thats a lot of work, but you're the one
in
DM the best position to do it. I'd say there's at least an even chance
DM that the bug is in your code rather than in nls().
yes.. and.. no :
- His code is quite peculiar, but I think only slightly too complicated
- one could argue that the bug is in Joerg's thinking that
something like
nls(y ~ eval(fitfunc), )
should be working at all.
But then he had found by experiment that it (accidentally I d'say)
does work in many cases.
DM And 2.5.0 *is* ancient; please confirm the bug exists in R-patched if
it
DM turns out to be an R bug.
You are right, but indeed (as has Kate just said) it *does*
exist in current R versions.
I agree that the behavior of nls() here is sub-optimal.
It *should* be consistent, i.e. work the same for n=4,5,6,..
I had spent about an hour after Joerg's R-devel posting,
and found to be too busy with more urgent matters --
unfortunately forgetting to give *some* feedback about my findings.
It may well be that we find that nls() should give an
(intelligible) error message in such eval() cases - rather than
only in one case...
Martin Maechler
DM Duncan Murdoch
DM so I try it here. sory for the
lengthy explanation but it seems unavoidable. to quickly see
the problem simply copy the litte example below and execute
f(n=5)
which crashes. called with n != 5 (and of course n3 since
there are 3 parameters in the model...) everything is as it
should be.
in detail:
I stumbled over the follwing _very_ strange behaviour/error
when using `nls' which I'm tempted (despite the implied
dangers) to call a bug:
I've written a driver for `nls' which allows specifying the
model and the data vectors using arbitrary symbols. these
are internally mapped to consistent names, which poses a
slight complication when using `deriv' to provide analytic
derivatives. the following fragment gives the idea:
#-
f - function(n = 4) {
x - seq(0, 5, length = n)
y - 2 * exp(-1*x) + 2;
y - rnorm(y,y, 0.01*y)
model - y ~ a * exp (-b*x) + c
fitfunc - deriv(model[[3]], c(a, b, c), c(a, b, c, x))
#standard call of nls:
res1 - nls(y ~ fitfunc(a, b, c, x), start = c(a=1, b=1, c=1))
call.fitfunc -
c(list(fitfunc), as.name(a), as.name(b), as.name(c),
as.name(x))
call.fitfunc - as.call(call.fitfunc)
frml - as.formula(y ~ eval(call.fitfunc))
#computed call of nls:
res2 - nls(frml, start = c(a=1, b=1, c=1))
list(res1 = res1, res2 = res2)
}
#-
the argument `n' defines the number of (simulated) data
points x/y which are going to be fitted by some model ( here
y ~ a*exp(-b*x)+c )
the first call to `nls' is the standard way of calling `nls'
when knowing all the variable and parameter names.
the second call (yielding `res2') uses a constructed formula
in `frml' (which in this example is of course not necessary,
but in the general case 'a,b,c,x,y' are not a priori known
names).
now, here is the problem: the call
f(4)
runs fine/consistently, as does every call with n 5.
BUT: for n = 5 (i.e. issuing f(5))
the second fit leads to the error message:
Error in model.frame(formula, rownames, variables, varnames, extras,
extranames, :
invalid type (language) for variable 'call.fitfunc'
I cornered this to a spot in
Re: [R] help needed: taking a function out of a package and it can not find some funtions
bw.SJ calls C code from the file /usr/local/bin/R-2.6.0/src/library/stats/src/bandwidths.c You have to make this code available. You can do the following: 1. copy bandwidths.c and the definition of bw.SJ (the latter renamed) to a directory 2. compile bandwidths.c into a shared library with the command R CMD SHLIB bandwidths.c (see the manual writing R extensions for details) 3. in your definition of bw.SJ, in the 3 places .C is called, remove R_ from the first argument and quote it; e.g., C(R_band_phi4_bin, ... becomes C(band_phi4_bin, 4. in R, load the shared library you built with the dyn.load function; now your definition of the (renamed) function bw.SJ can be modified as you like. On Tue, 6 Nov 2007, Jason Liao wrote: I tried to modify the R function bw.SJ (from the stats package) for my own use. But if I just get the source code and run directly (without any modification), it can no longer find some key functions called within it. I understand this has something to do with searching path which I do not understand well. Can anyone tell me how to modify the source code of an R function and still make it part of an existing package? Thanks. Jason Liao, http://www.geocities.com/jg_liao Associate Professor of Biostatistics Drexel University School of Public Health 245 N. 15th Street, Mail Stop 660 Philadelphia, PA 19102-1192 phone 215-762-3934 __ 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] help needed: taking a function out of a package and it can not find some funtions
OK, well, the way I suggested will allow you to modify the C code too -- but reading your question I realize you probably don't want to do this. If you want to use the definitions of the C functions from stats, you can skip 2. below and add another modification to the calls to .C in your def. of bw.SJ; add the argument PACKAGE=stats. On Tue, 6 Nov 2007, Katharine Mullen wrote: bw.SJ calls C code from the file /usr/local/bin/R-2.6.0/src/library/stats/src/bandwidths.c You have to make this code available. You can do the following: 1. copy bandwidths.c and the definition of bw.SJ (the latter renamed) to a directory 2. compile bandwidths.c into a shared library with the command R CMD SHLIB bandwidths.c (see the manual writing R extensions for details) 3. in your definition of bw.SJ, in the 3 places .C is called, remove R_ from the first argument and quote it; e.g., C(R_band_phi4_bin, ... becomes C(band_phi4_bin, 4. in R, load the shared library you built with the dyn.load function; now your definition of the (renamed) function bw.SJ can be modified as you like. On Tue, 6 Nov 2007, Jason Liao wrote: I tried to modify the R function bw.SJ (from the stats package) for my own use. But if I just get the source code and run directly (without any modification), it can no longer find some key functions called within it. I understand this has something to do with searching path which I do not understand well. Can anyone tell me how to modify the source code of an R function and still make it part of an existing package? Thanks. Jason Liao, http://www.geocities.com/jg_liao Associate Professor of Biostatistics Drexel University School of Public Health 245 N. 15th Street, Mail Stop 660 Philadelphia, PA 19102-1192 phone 215-762-3934 __ 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-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] correct wording and notation for R stuff in LaTex
If you are free to format your references to packages and functions as you
please, you might follow the convention used by the Journal of Statistical
Software.
you can do this by adding the following to your latex file:
\newcommand{\pkg}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\let\proglang=\textsf
\let\code=\texttt
then you can refer to R as \proglang{R}, to packages with e.g.,
\pkg{packagename} and to functions/other code as \code{functionname}.
On Thu, 1 Nov 2007, Edna Bell wrote:
Hi R Gurus:
I'm putting together an article about some R stuff in Latex.
I refer to packages and functions.
I think that I use {\em} for packages and {\tt} for functions.
Is that correct, please?
Thank you in advance!
Sincerely,
Edna Bell
__
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] Find A, given B where B=A'A
B is symmetric by definition; if it's also real positive-definite then A is the upper triangular factor of the Choleski decomposition, and you can use chol(B) to get A. On Wed, 31 Oct 2007, Michael Gormley wrote: Given a matrix B, where B=A'A, how can I find A? In other words, if I have a matrix B which I know is another matrix A times its transpose, can I find matrix A? Thanks, Mike __ 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] How to make own function load automatically on startup
Regarding your first issue: you could make a package to contain your function (see ?package.skeleton and the manual Writing R extensions) - then you could use library(yourpackagename). Or you could save the definition of your function(s) in a text file, and use source(textfilename) to load the definition(s) after you start R. Regarding your third issue: see help(::) On Sat, 27 Oct 2007, Jonas Malmros wrote: Dear list members, I have written a function, called say Analysis. I supply an Excel file name as an argument, it does analysis on this file and returns a pdf file with specific plots, and a text file with several statistical models' output (I extract certain values from the output and create my own custom dataframe with output). As of now I have to open the script file and load the function by hand every time I start R. I want to make it load automatically into the environment when I start R. How can I make this happen? Also, is there a simple way to combine plots and my custom dataframe with model output in one pdf file? I do use LaTeX, but my collegues who want to run this function as well do not know what LaTeX is. So I wonder if I can avoid LaTeX when creating this combined PDF file? (Avoid in the sense that one does not have to fiddle with LaTeX itself). Finally, another small question. My code contains following snippet: xlsReadWrite::read.xls, this is because another package I load masks read.xls function and :: help me to select correct read.xls. I guessed this use of :: but I would really like to know what :: is and where can I read more about this function? Thank you very much in advance for your time and help! -- Jonas Malmros Stockholm University Stockholm, Sweden __ 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] flipping vector and matrix
sorry, just do rev on the columns (no t()): x [,1] [,2] [,3] [1,]123 [2,]456 [3,]789 apply(x,2,rev) [,1] [,2] [,3] [1,]789 [2,]456 [3,]123 On Tue, 23 Oct 2007, Katharine Mullen wrote: One way: x-1:9 rev(x) [1] 9 8 7 6 5 4 3 2 1 for x as the matrix you gave: x [,1] [,2] [,3] [1,]123 [2,]456 [3,]789 apply(t(x),1,rev) [,1] [,2] [,3] [1,]789 [2,]456 [3,]123 On Tue, 23 Oct 2007, Rainer M Krug wrote: Hi I have a vector x - c(1, 2, 3, 4, 5) and I want to flip it around, i.e. I need 5, 4, 3, 2, 1 Is there a ssolution apart from y - x[length(x):1] I am also looking for the same for a matrix M, i.e. 1 2 3 4 5 6 7 8 9 should become 7 8 9 4 5 6 1 2 3 again, I am using M[1:dim(M)[1], dim(M)[2]:1] Thanks Rainer __ 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-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] export R-data to VisIt
Below is an example R - netCDF - R for rows of a dataframe that are numeric vectors --note however that your dataframe includes character vectors. I can't look into that case at the moment - maybe it's easy to solve, or maybe you have to do some hashing. ## begin ex. library(ncdf) dat - matrix(rnorm(20),10,2) c1 - dim.def.ncdf( c1, c1units, 1:nrow(dat) ) c2 - dim.def.ncdf( c2, c2units, 1:nrow(dat) ) x1 - var.def.ncdf(name=v1, units=c1units,dim = c1, missval=0) x2 - var.def.ncdf(name=v2, units=c2units,dim = c2, missval=0) ## for some reason, when vars is the vector c(x1,x2) this does't work ## it may be a bug; get around my adding other vars later ncnew - create.ncdf(filename=file.cdf, vars = x1) ncnew - var.add.ncdf(ncnew, x2) put.var.ncdf(ncnew, v1, dat[,1]) put.var.ncdf(ncnew, v2, dat[,2]) close.ncdf(ncnew) ofile - open.ncdf(file.cdf) c_x1 - get.var.ncdf(ofile, v1) c_x2 - get.var.ncdf(ofile, v2) Peter, what a quick response! But unfortunately, yes I tried the ncdf package, I looked at the examples, but after 2 hours trying and many, many errors, I gave up. Bart - Original Message - From: Peter Dalgaard [EMAIL PROTECTED] To: Bart Joosen [EMAIL PROTECTED] Cc: [EMAIL PROTECTED] Sent: Friday, October 19, 2007 4:56 PM Subject: Re: [R] export R-data to VisIt Bart Joosen wrote: Hello, Is there anyone porting R data to VisIt (http://www.llnl.gov/visit/)? Altough VisIt accepts 5 dozen of data formats, I can't get my data into VisIt. I currently ran a simulation which gave me a data frame, which I wanted to import into VisIt to further explore the dataframe. Let's say I have a data frame as follows: dat - data.frame(cbind( 1, 1:10),X3= sample(LETTERS[1:3], 10, repl=TRUE)) My currently data export is write.table(dat,C:/filename.csv) But I can't import this kind of data in Visit. An option is to export my dataframe as a .CDF file, but I couldn't get the right output of my dataframe with netcdf. So how do I put my dataframe in a netCDF format, or is there anyone who knows the easiest way to transport my data to VisIt? Have you tried the ncdf package? -- 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. __ 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] Problems with paste and blank
use sep=, e.g., paste(Laufwerk,Stadtwerksname,.csv,sep=) On Wed, 17 Oct 2007 [EMAIL PROTECTED] wrote: Hi there, I've got the following problem under Windows XP, R 2.5.1: When I'm pasting some objects: Stadtwerksname-Mannheim Laufwerk-C:\\ paste(Laufwerk,Stadtwerksname,.csv) I get the result: [1] C:\\ Mannheim .csv The problem's are the superfluous gaps/blanks between the three parts. Is there a way to get rid off this problem? Thx, Thomas __ Thomas Schwander MVV Energie Konzern-Risikocontrolling Telefon 0621 - 290-3115 Telefax 0621 - 290-3664 E-Mail: [EMAIL PROTECTED] . Internet: www.mvv.de MVV Energie AG . Luisenring 49 . 68159 Mannheim Handelsregister-Nr. HRB 1780, Amtsgericht Mannheim Vorsitzender des Aufsichtsrates: Herr Dr. Peter Kurz Vorstand: Dr. Rudolf Schulten (Vorsitzender) . Matthias Br?ckmann . Dr. Werner Dub . Hans-J?rgen Farrenkopf [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [R-pkgs] new package 'nnls'
A new package 'nnls' is now available on CRAN. The package provides an R interface to the Lawson-Hanson NNLS algorithm for non-negative least squares that solves the least squares problem A x = b with the constraint x = 0. The Lawson-Hanson NNLS algorithm was published in Lawson CL, Hanson RJ (1974). Solving Least Squares Problems. Prentice Hall, Englewood Cliffs, NJ. Lawson CL, Hanson RJ (1995). Solving Least Squares Problems. Classics in Applied Mathematics. SIAM, Philadelphia. and is available as Fortran77 code on Netlib (file lawson-hanson/all). The 'nnls' package interfaces to this code. Included in the examples section of the function 'nnls' is a test problem comparing NNLS to the L-BFGS-B method of 'optim' and to the 'solve.QP' function of the package 'quadprog'. NNLS is shown to be faster and slightly more accurate than these more general purpose algorithms for the test problem examined. I do not have access to S-PLUS or the S-PLUS source, but the help page for the S-SPLUS function 'nnls.fit' references Lawson and Hanson (1974), and is probably close to the port here. I have not written any methods for printing or summaries, and only return the solution x. On request, I can modify this behavior. I would be interested in suggestions, bug reports, and other comments. Kate Mullen Katharine Mullen mail: Department of Physics and Astronomy, Faculty of Sciences Vrije Universiteit Amsterdam, de Boelelaan 1081 1081 HV Amsterdam, The Netherlands room: T.1.06 tel: +31 205987870 fax: +31 205987992 e-mail: [EMAIL PROTECTED] homepage: http://www.nat.vu.nl/~kate/ ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] the use of the .C function
Let me also comment that you are trying to interface to a function ported
to C++ from the Cephes library, which is in C. You have not explained why
you are trying to interface to the function (is this an exercise, or do
you suspect a problem with digamma()?), but if you just want access to it
as implemented in Cephes, you may as well use the original C:
http://www.netlib.org/cephes/
On Sun, 14 Oct 2007, Duncan Temple Lang wrote:
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1
The code is C++ and so is compiled
using the C++ compiler, not the C compiler.
This is because the name of the file is .cpp
(and you include iostream.h, but don't seem to make any use of it.)
As a result, the names of the routines are mangled
and so psi is not the actual name of the compiled routine,
but something more like _Z3psid which encodes the
types of the parameters in the name (to allow overloaded routines).
To avoid the mangling, use
#ifdef __cplusplus
extern C
#endif
double psi(double)
Alternatively, you can use the registration mechanism to
associate a name with a routine directly for R. See
Writing R Extensions which would be a good thing to read
as trying to get this sort of thing working by trial and
error can be both very time consuming, and also very confusing
to the point that you unlearn things.
After that you will still need to make some changes
to use the routine with the .C() interface.
That expects inputs as pointers, i.e. the double
should be a double * and you fetch the value from
that. And a routine called by .C() should return nothing
directly, but provide the results in space provided
by the inputs
So
void R_psi(double *x)
{
*x = psi(*x);
}
and make certain R_psi is not mangled.
Also, you might want to upgrade to the latest version of R.
Bernardo Lagos Alvarez wrote:
Dear All,
could someone please shed some light on the use of the .C or .Fortran
function:
I am trying load and running on R the following function
// psi.cpp -- psi function for real arguments.
// Algorithms and coefficient values from Computation of Special
// Functions, Zhang and Jin, John Wiley and Sons, 1996.
//
// (C) 2003, C. Bond. All rights reserved.
//
// Returns psi function for real argument 'x'.
// NOTE: Returns 1e308 for argument 0 or a negative integer.
//
#include iostream.h // or stdio.h?
#include math.h
#define el 0.5772156649015329
double psi(double x)
{
double s,ps,xa,x2;
int n,k;
static double a[] = {
-0.8e-01,
0.8e-02,
-0.39682539682539683e-02,
0.41667e-02,
-0.75757575757575758e-02,
0.21092796092796093e-01,
-0.8e-01,
0.4432598039215686};
xa = fabs(x);
s = 0.0;
if ((x == (int)x) (x = 0.0)) {
ps = 1e308;
return ps;
}
if (xa == (int)xa) {
n = xa;
for (k=1;kn;k++) {
s += 1.0/k;
}
ps = s-el;
}
else if ((xa+0.5) == ((int)(xa+0.5))) {
n = xa-0.5;
for (k=1;k=n;k++) {
s += 1.0/(2.0*k-1.0);
}
ps = 2.0*s-el-1.386294361119891;
}
else {
if (xa 10.0) {
n = 10-(int)xa;
for (k=0;kn;k++) {
s += 1.0/(xa+k);
}
xa += n;
}
x2 = 1.0/(xa*xa);
ps = log(xa)-0.5/xa+x2*(((a[7]*x2+a[6])*x2+a[5])*x2+
a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]);
ps -= s;
}
if (x 0.0)
ps = ps - M_PI*cos(M_PI*x)/sin(M_PI*x)-1.0/x;
return ps;
}
However, when applicated the codes
digamma(-0.9)
[1] -9.312644 OK!
But when
dyn.load(psi.so)
out-.C(psi, as.double(-0.9))
Erro en .C(psi, as.double(-0.9)) : C symbol name psi not in load table
out
Error: objeto out no encontrado
More information on OS:
version
_
platform i486-pc-linux-gnu
arch i486
os linux-gnu
system i486, linux-gnu
status
major 2
minor 4.1
year 2006
month 12
day18
svn rev40228
language R
version.string R version 2.4.1 (2006-12-18)
Many thanks for any insight or comments.
Bernardo Lagos A.
__
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.
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.7 (Darwin)
Comment: Using GnuPG with Mozilla -
Re: [R] Anybody has ever met the problem to add a legend to a figure generated by image()
have you tried the function image.plot in the package fields? On Fri, 12 Oct 2007, zhijie zhang wrote: Dear friends, Anybody has ever met the problem to add a legend to a figure generated by image()? I have three variables,x,y and z. x and y are the coordinates, and z is the third values. we can use image(x, y, z,...) to generate a figure according to the z-values, but the problem is the figure legend. How can the legend be added to a figure generated by image()? Note that filled.contour() can add the figure legend automatically, but there are some problems sometime. I want to know the specific method for adding the legend to the figure generated by image(). Thanks very much. -- With Kind Regards, oooO: (..): :\.(:::Oooo:: ::\_)::(..):: :::)./::: ::(_/ : [***] Zhi Jie,Zhang ,PHD Tel:86-21-54237149 Dept. of Epidemiology,School of Public Health,Fudan University Address:No. 138 Yi Xue Yuan Road,Shanghai,China Postcode:200032 Email:[EMAIL PROTECTED] Website: www.statABC.com [***] oooO: (..): :\.(:::Oooo:: ::\_)::(..):: :::)./::: ::(_/ : [[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.
Re: [R] psi using .C with R
in case you are not aware: the function psi is available without a port; see help(digamma). On Wed, 10 Oct 2007, Bernardo Lagos Alvarez wrote: Hi All, Anybody know as run the function psi on http://www.alglib.net/translator/dl/specialfunctions.psi.csharp.zip or http://www.alglib.net/translator/dl/specialfunctions.psi.cpp.zip using .C or .Fortran ? Thank for your attention. Bernardo. __ 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] convert a tar.gz to a Windows zip
Uwe Ligges seems to be providing a simple way to make the conversion (though I have not used it): http://win-builder.r-project.org/ On Mon, 1 Oct 2007, Edna Bell wrote: Dear R Gurus; Is there a simple way to convert a Linux produced tar.gz file (a package) to a Windows binary zip package, please? Thanks in advance __ 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] R routines vs. MATLAB/SPSS Routines
I would stress the advantages of the free and open source nature of R over the proprietary programs you mention. Because R is free (as in beer), your student will have access to it even when they are free of the university that I presume buys a MATLAB/SPSS license for them. And because R is open source, when your student has an idea regarding how to make R work faster/more efficiently/more accurately, they can modify the source code at whatever level they like. They can also examine at whatever level they like what R is doing, so there are no black boxes involved. On Tue, 2 Oct 2007, Matthew Dubins wrote: Hi all, I've become quite enamored of R lately, and have decided to try to teach some of its basics (reading in data, manipulation and classical stats analyses) to my fellow grad students at the University of Toronto. I sent out a mass email and have already received some positive responses. One student, however, wanted to know what differentiates the routines that R uses, from those that MATLAB and SPSS use. In other words, in what respects do R routines work faster/more efficiently/more accurately than those of MATLAB/SPSS. I thank you in advance for any answer you can give me (or rather, the inquiring student). Cheers, Matthew Dubins __ 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] Creating nice looking lists: how?
The source code for print.summary.lm and summary.lm is in the file /src/library/stats/R/lm.R of the R source code, which you can download from CRAN if you don't have it already. The file lm.R is also at https://svn.r-project.org/R/trunk/src/library/stats/R/lm.R On Fri, 28 Sep 2007, Sergey Goriatchev wrote: Hello, For my functions I want to create output similar in appearance to that of what you get when you print a summary of lm model: Residuals: Min1Q Median3Q Max -0.209209 -0.043133 0.001793 0.044105 0.234750 Coefficients: Estimate Std. Error t valuePr(|t|) (Intercept) 0.981762 0.004089 240.103 2e-16 *** Factor 1-0.009581 0.006381 -1.501 0.134296 Factor 2-0.008993 0.009182 -0.979 0.328163 Factor 3 0.029960 0.009547 3.138 0.001866 ** Factor 4-0.026575 0.007370 -3.606 0.000363 *** Factor 5-0.004847 0.006382 -0.760 0.448138 Factor 6 0.005099 0.006483 0.786 0.432202 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 I want: 1) no $ before the list component names 2) component names that take values from outside variables (ex.: number - 10 There are 'number' factors in the model:) There is not much information on how to create nice output in terms of lists, so I was looking for core to write the summary(lm) output, but could not find much. Obviously, I can type summary.lm, but it does not show how to create the name Coefficients: Could someone give me pointers on how to create nice lists? Thanks Sergey __ 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] extracting data using strings as delimiters
have you seen help(strsplit)? On Tue, 25 Sep 2007, lucy b wrote: Dear List, I have an ascii text file with data I'd like to extract. Example: Year Built: 1873 Gross Building Area: 578 sq ft Total Rooms: 6 Living Area: 578 sq ft There is a lot of data I'd like to ignore in each record, so I'm hoping there is a way to use strings as delimiters to get the data I want (e.g. tell R to take data between Built: and Gross - incidentally, not always numeric). I think an ugly way would be to start at the end of each record and use a substitution expression to chip away at it, but I'm afraid it will take forever to run. Is there a way to use strings as delimiters in an expression? Thanks in advance for ideas. LB __ 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] nls fits by groups
It is not clear from your post what changes per-group. If only the
starting values change (but the data and the model structure are the
same), then you can just store the starting values you want to use for
each group in a list, and then index into this list in your call to nls.
e.g., modifying an example in the help page for nls:
x - 1:10
y - 2*x + 3# perfect fit
yeps - y + rnorm(length(y), sd = 0.01) # added noise
startlist - list(
list(a = 0.12345, b = 0.54321), ##group 1 start val
list(a = 0.12, b = 0.54) ## group 2 start val.
)
reslist - list() ## filling this with results from different start val
for(i in 1:length(startlist)) {
reslist[[i]] - nls(yeps ~ a + b*x, start = startlist[[i]],
trace = TRUE)
}
On Sun, 23 Sep 2007, Aleksi Lehtonen wrote:
Dear Colleagues,
I am trying to estimate several non-linear models simultaneously. I don't
want to use non-linear mixed model, but non-linear model with same form, but
it should be estimated separately according to variable group (I have lots
of groups that have lots of observations). I would like to have unique
parameters for each group.
e.g. something like this
mod - nls(y ~ a*x^b, start=c(a=1, b=1), group=group)
but knowing that group option does not work. If someone has an idea (or has
done it already) how to implement this either using just nls statement or by
building a simple function in R, I would be very grateful for hints
regards, Aleksi Lehtonen
[[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.
Re: [R] Smooth line in graph
how about: require(splines) x-1:5 y - c(0.31, 0.45, 0.84, 0.43, 0.25) yy -predict(interpSpline(x, y)) plot(x, y) lines(yy) Katharine Mullen mail: Department of Physics and Astronomy, Faculty of Sciences Vrije Universiteit Amsterdam, de Boelelaan 1081 1081 HV Amsterdam, The Netherlands room: T.1.06 tel: +31 205987870 fax: +31 205987992 e-mail: [EMAIL PROTECTED] homepage: http://www.nat.vu.nl/~kate/ On Wed, 19 Sep 2007, Nestor Fernandez wrote: Hi, I’m trying to get smooth curves connecting points in a plot using spline but I don’t get what I whant. Eg.: x-1:5 y - c(0.31, 0.45, 0.84, 0.43, 0.25) plot(x,y) lines(spline(x,y)) Creates a valley between the first and second points, then peaks at 3rd, and another valley between 4th and 5th. I’m trying to get a consistently growing curve up to the 3rth point and then a decrease like with SigmaPlot spline curves or with Excel. I tried with different spline arguments and also lowess and loess, with no success. Any ideas? Thanks. __ 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] help with unit indication when plotting
see ?plotmath you want to use mu for micro, as in plot(1:10, ylab = expression(paste(mu, M))) On Thu, 13 Sep 2007, Zheng Lu wrote: Dear all: Is there any one could tell me how I can represent Micro-molar as an unit of concentration when I plot with R(S-plus), I don't want write 'uM' from keyboard, I am thinking to write it like in word, in word, people insert symbol for 'u' for uM. Am I clear? Thank you very much for your consideration and help. Lu __ 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.
