[R] LogLik of nls
Hello there, Can anyone point me to the code for logLik of an nls object? I found the code for logLik of an lm but could not find exactly what function is used for calculating the logLik of nls function? I am using the nls to fit the following model to data - Model 1: y ~ Ae^(-mx) + Be^(-nx) +c and want to understand what is the likelihood function used by nls. Presumably it is using - N(Ae^(-mx) + Be^(-nx) +c, var = residuals from nls) Thank you in advance, Diviya [[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] Legend symbols
Hi there, I was wondering if there is any R package that one can use for plotting that has more legend symbols - the standard pch has 18 symbols but I need ~30 for my application- and just using different colors is not an option. Thank you in advance, Diviya [[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] nls for sum of exponentials
Hi there, I am trying to fit the following model with a sum of exponentials - y ~ Ae^(-md) + B e^(-nd) + c the model has 5 parameters A, b, m, n, c I am using nls to fit the data and I am using DEoptim package to pick the most optimal start values - fm4 - function(x) x[1] + x[2]*exp(x[3] * -dist) + x[4]*exp(x[5] * -dist) fm5 - function(x) sum((wcorr-fm4(x))^2) fm6 - DEoptim(fm5, lower=c(0,0.1,1,0.1,1), upper=c(10e7, 100, 300,100, 300 ), control=list(trace=FALSE)) par2 - fm6$optim$bestmem names(par2) - c(c, A, n1, B, n2) fit2 - nls(wcorr ~ (c + A*exp(n1 * -dist) + B*exp(n2 * -dist)), start=par2,control =list(maxiter=1000, warnOnly=TRUE)) However, I keep getting the following error - Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model OR In nls(wcorr ~ (c + A * exp(n1 * -dist) + B * exp(n2 * -dist)), : singular gradient Can anyone suggest how to resolve the problem? Thanks, Priya [[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] Rscript installing packages
Hi there
I have an Rscript and I am looking for a way to install a package
non-interactively. In Rscript {Utils}, I saw an example which does
something like this, however this does not seem to work for my particular
example. I am trying to install the following package in an Rscript
(without switching to interactive mode).
res - try(install.packages(DEoptim))
if(inherits(res, try-error)) q(status=1) else q()
Error:
Installing package(s) into /home/xx/R/x86_64-pc-linux-gnu-library/2.13
(as lib is unspecified)
Error in contrib.url(repos, type) :
trying to use CRAN without setting a mirror
Any idea how I can fix this or specify lib? Any help would be most
appreciated.
Thank you in advance,
Diviya
[[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] Solving quadratic equation in R
Hi there, I would like to solve a simple equation in R a^2 - a = 8.313 There is no real solution to this problem but I would like to get an approximate numerical solution. Can someone suggest how I can set this up? Thanks in advance, Diviya [[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] Solving quadratic equation in R
Sorry it is important for me to constrain the value of 'a' between c(0,1) On Thu, Jul 26, 2012 at 4:48 PM, Diviya Smith diviya.sm...@gmail.comwrote: Hi there, I would like to solve a simple equation in R a^2 - a = 8.313 There is no real solution to this problem but I would like to get an approximate numerical solution. Can someone suggest how I can set this up? Thanks in advance, Diviya [[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] Solving quadratic equation in R
On Thu, Jul 26, 2012 at 5:16 PM, Diviya Smith diviya.sm...@gmail.comwrote:
Thank you for pointing me to the uniroot function?
Is there a way to constrain this solution so that it only gives me values
of 'a' between c(0,1)?
I tried using nlminb and for some reason it always estimates a = 0, even
when I change the objective function.
f - function(vals) {
x - vals[1]
sum(c(x^2-x-8.013)^2)
}
nlminb(c(0.05),f, lower=0.00, upper=1.00)
On Thu, Jul 26, 2012 at 5:00 PM, è¾è´¨æ ggg...@gmail.com wrote:
I think it's impossible for 'a' between c(0,1) to satisfy that
a^2 - a = 8.313
Without the restraint, you can try following lines to get roots:
fun=function(a) a^2-a-8.313
root1=uniroot(fun,c(3,4))
root2=uniroot(fun,c(-3,-2))
Eric
On Thu, Jul 26, 2012 at 4:49 PM, Diviya Smith diviya.sm...@gmail.comwrote:
Sorry it is important for me to constrain the value of 'a' between c(0,1)
On Thu, Jul 26, 2012 at 4:48 PM, Diviya Smith diviya.sm...@gmail.com
wrote:
Hi there,
I would like to solve a simple equation in R
a^2 - a = 8.313
There is no real solution to this problem but I would like to get an
approximate numerical solution. Can someone suggest how I can set this
up?
Thanks in advance,
Diviya
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Solving equations in R
Hi there, I would like to solve the following equation in R to estimate 'a'. I have the amp, d, x and y. amp*y^2 = 2*a*(1-a)*(-a*d+(1-a)*x)^2 test data: amp = 0.2370 y= 0.0233 d= 0.002 x= 0.091 Can anyone suggest how I can set this up? Thanks, Diviya [[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] Linear regression
Hello there, I am new to using regression in R. I wanted to solve a simple regression problem where I have 2 equations and 2 unknowns. So lets say - y1 = alpha1*A + beta1*B y2 = alpha2*A + beta2*B y1 - runif(10, 0,1) y2 - runif(10,0,1) alpha1 - 0.6 alpha2 - 0.75 beta1 - 1-alpha1 beta2 - 1-apha2 I now want this equation to estimate the values of A and B. Both A and B are constrained to be between (0,1). I would like to use lm with these constraints and I am having a little trouble in defining the equations correctly. Any help would be most appreciated. Thank you, Diviya [[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] Linear regression
Hello Bert, This is definitely not for a homework problem. I am trying to estimate frequencies of mutations in different groups. The mutation frequencies can be modeled as a linear relation in cases of mixtures. So I have a lot of populations that follow the relationship - y = alpha*A + beta*B and I want to estimate A and B; given y, alpha and beta. A and B are both vectors of the same size as y. Can you suggest where I can find some information about your suggestion #4...that is exactly what I was hoping to do. Thanks, Diviya On Mon, Mar 19, 2012 at 3:02 PM, Bert Gunter gunter.ber...@gene.com wrote: 1. Homework assignment? We don't do homework here. 2. If not, a mixture model of some sort? I suggest you state the context of the problem more fully. R has several packages to do mixture modeling, if that's what you're trying to do. 3. In any case, this cannot be done with lm() (at least without tricks). 4. In your notation below, the separate regressions can be stacked into a single constrained regression model. 5. You might do better to find local statistical help, as you may have bitten off more than you can chew. -- Bert On Mon, Mar 19, 2012 at 11:29 AM, Diviya Smith diviya.sm...@gmail.com wrote: Hello there, I am new to using regression in R. I wanted to solve a simple regression problem where I have 2 equations and 2 unknowns. So lets say - y1 = alpha1*A + beta1*B y2 = alpha2*A + beta2*B y1 - runif(10, 0,1) y2 - runif(10,0,1) alpha1 - 0.6 alpha2 - 0.75 beta1 - 1-alpha1 beta2 - 1-apha2 I now want this equation to estimate the values of A and B. Both A and B are constrained to be between (0,1). I would like to use lm with these constraints and I am having a little trouble in defining the equations correctly. Any help would be most appreciated. Thank you, Diviya [[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. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm [[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] Strange plotting error
Hi there, I am trying to compute the autocorrelation in a dataset using R's acf function. ACF automatically plots the results. This works well except in some cases xlim doesnt work data - rnorm(2000,0,1) acf(data,xlim=c(1,10)) # works - the plot starts at 1 acf(data,lag=100,xlim=c(1,100)) # this does not work and the plot still starts at 0 Is there another way to specify the xlim or the starting value for x? Any help would be most appreciated. Thanks, Diviya [[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 function implementation
Hello there, I recently wrote some code to perform pairwise correlations between all samples in a large dataset. So we are talking about performing pairwise correlations between 400K vectors. Since R has a very rich library of functions, it was very easy to code this in R. However, R was probably not the best choice as it is super slow for this large job. So my plan is to recode it in C. I was wondering if it is possible to see the code between some of the functions. I am looking for the raw code for simple data handling functions such as split, cut, etc so that I can directly use those in C. I am fairly new to programming in C and so this will be a big help. Thanks for you help in advance. Diviya [[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] Binning the data based on a value
Hello there, I have a matrix with some data and I want to split this matrix based on the values in one column. Is there a quick way of doing this? I have looked at cut but I am not sure how to exactly use it? for example: I would like to split the matrix a based on the spending such that the data is binned groups [0..99],[100..199]...and so on. a - data.frame(patient=1:7, charges=c(100,500,200,90,400,500,600), age=c(0,3,5,7,10,16,19), spending=c(10, 60, 110, 200, 250, 400, 450)) Expected output - bin[1] - c(10, 60) bin[2] - c(110, 200, 250) bin[3] - c(400, 450) NOTE that the number of data points in each bin is not the same and the empty bins are removed (since there are no points between [199..299], bin[3] starts at 400. Any help would be most appreciated. Thank you in advance. Diviya [[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] Binning the data based on a value
Thank you very much Michael. This is very helpful. However, if there is any way to exclude zero length bins. Lets imagine that the matrix was as follows - a - data.frame(patient=1:7, charges=c(100,500,200,90,400,500,600), age=c(0,3,5,7,10,16,19), spending=c(10, 60, 110, 200, 400, 450, 500)) bins - split(a, cut(a$spending, breaks = (0:5)*100) then bins[3] = $`(200,300]` [1] patient charges age spending 0 rows (or 0-length row.names) Is there a way to exclude this? Priya On Mon, Dec 5, 2011 at 3:43 PM, R. Michael Weylandt michael.weyla...@gmail.com wrote: Just a clarification: I can't get round to work as I first expected so if you want to do bins by 100's you'd probably want: split(a, cut(a$spending, breaks = (0:5)*100)) Michael On Mon, Dec 5, 2011 at 3:41 PM, R. Michael Weylandt michael.weyla...@gmail.com wrote: I'd so something like split(a, a$spending) and you can include a round(a$spending, -2) or something similar if you want to group by the 100's. Michael On Mon, Dec 5, 2011 at 3:37 PM, Diviya Smith diviya.sm...@gmail.com wrote: Hello there, I have a matrix with some data and I want to split this matrix based on the values in one column. Is there a quick way of doing this? I have looked at cut but I am not sure how to exactly use it? for example: I would like to split the matrix a based on the spending such that the data is binned groups [0..99],[100..199]...and so on. a - data.frame(patient=1:7, charges=c(100,500,200,90,400,500,600), age=c(0,3,5,7,10,16,19), spending=c(10, 60, 110, 200, 250, 400, 450)) Expected output - bin[1] - c(10, 60) bin[2] - c(110, 200, 250) bin[3] - c(400, 450) NOTE that the number of data points in each bin is not the same and the empty bins are removed (since there are no points between [199..299], bin[3] starts at 400. Any help would be most appreciated. Thank you in advance. Diviya [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] barplot names.arg
Hello there, I have a question regarding bar plots. I am trying to plot the data from the following matrix as a barplot - # input data mdat - matrix(c(0.1,0.9,0.9,0.1,0.5,0.5,0.45,1-0.45,0.6,0.4,0.8,0.2), nrow = 6, ncol=2, byrow=TRUE, +dimnames = list(c(Mon, Mon, Tues, Tues, Thurs, Friday), +c(C.1, C.2))) # plot barplot(t(as.matrix(mdat)), col=rainbow(2), ylab=Sales, names.arg = rownames(mdat), border=NA, cex.names=0.5) I am using names.arg to print the label for the bars. However, I would like to group all the entries for the Mon and print the label only once. Is there a way to do this? The documentation for barplots suggests that this can be done but I was not able to figure it out. Please help. Thanks in advance, Diviya [[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] barplot names.arg
The documentation for bar plots includes - names.arga vector of names to be plotted below each bar or *group of bars*. If this argument is omitted, then the names are taken from the names attribute of height if this is a vector, or the column names if it is a matrix Thanks for your suggestion. However, I am still not sure how to group all the entries for Mon and include one label. Your solution sort of does that but it actually does not print the individual values . I would like to show that there is a lot of variation in C.1 for Mon. On Thu, Nov 10, 2011 at 11:32 PM, R. Michael Weylandt michael.weyla...@gmail.com wrote: I'm not sure you can do that data aggregation in barplot directly (a quick skim doesn't reveal anything in the documentation that suggests it to me, thought I might have missed it) though I think this does what you are talking about: barplot(sapply(unique(rownames(mdat)), function(n) colSums(mdat[n,,drop=F])), col = rainbow(2)) Michael On Thu, Nov 10, 2011 at 9:21 PM, Diviya Smith diviya.sm...@gmail.com wrote: Hello there, I have a question regarding bar plots. I am trying to plot the data from the following matrix as a barplot - # input data mdat - matrix(c(0.1,0.9,0.9,0.1,0.5,0.5,0.45,1-0.45,0.6,0.4,0.8,0.2), nrow = 6, ncol=2, byrow=TRUE, +dimnames = list(c(Mon, Mon, Tues, Tues, Thurs, Friday), +c(C.1, C.2))) # plot barplot(t(as.matrix(mdat)), col=rainbow(2), ylab=Sales, names.arg = rownames(mdat), border=NA, cex.names=0.5) I am using names.arg to print the label for the bars. However, I would like to group all the entries for the Mon and print the label only once. Is there a way to do this? The documentation for barplots suggests that this can be done but I was not able to figure it out. Please help. Thanks in advance, Diviya [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] NLS- check convergence
Hi there, I am having some trouble with NLS convergence for my function. I was wondering if there is a way to check the value of the indicator function for all the iterations to look at the surface of the likelihood function. Any help would be most appreciated. Thanks, Diviya [[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] NLS error
Hello there, I am using NLS for fitting a complex model to some data to estimate a couple of the missing parameters. The model is - y ~ (C+((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*d+d*B+log(1-r) where A, B and C are unknown. In order to test the model, I generate data by setting values for all parameters and add some noise (C). A - 20 B - 500 r - 0.6881 d - (1:1000)/1000 y ~ (rnorm(d)*0.005)+(((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*d+d*B+log(1-r) I use Deoptim package to pick the optimum starting values. The model works fine in most cases, but every now and again, I get the following error - Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model Any suggestions on how I can resolve this? Can you suggest a better way for picking the starting parameters? Thanks, Diviya [[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] NLS error
I dont think *r* is related to the problem. I am not trying to estimate *r* and so basically I am giving the model the correct value of *r* and so log(1-r) should not go to infinity. For test data, I generate data from the same model and add noise (using *r norm*), with the following parameters - A - 20 B - 500 r - 0.6881 d - (1:1000)/1000 y - (rnorm(d)*0.005)+(((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*+d*B+log(1-r) And then I try estimating A, B and C (Note: C is capturing the noise term, and I am specifying r and d) but I get the following error - Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model Any suggestions how to fix this? Diviya On Tue, Sep 20, 2011 at 2:37 PM, Jean V Adams jvad...@usgs.gov wrote: Diviya Smith wrote on 09/20/2011 01:03:22 PM: Hello there, I am using NLS for fitting a complex model to some data to estimate a couple of the missing parameters. The model is - y ~ (C+((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*d+d*B+log(1-r) where A, B and C are unknown. In order to test the model, I generate data by setting values for all parameters and add some noise (C). A - 20 B - 500 r - 0.6881 d - (1:1000)/1000 y ~ (rnorm(d)*0.005)+(((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*d +d*B+log(1-r) I use Deoptim package to pick the optimum starting values. The model works fine in most cases, but every now and again, I get the following error - Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model Any suggestions on how I can resolve this? Can you suggest a better way for picking the starting parameters? Thanks, Diviya I'm not sure if this is the problem, but if r grows greater than 1, log(1-r) will be undefined, and you'll get an error. You can impose a constraint on r by rewriting your formula in terms of a variable that can take on any real value: R = log(1-r) So, replace log(1-r) in your formula with R, replace (-1 + r) with -exp(R), and replace r with 1 - exp(R): y ~ (C+(R*exp(-A*d)*(-exp(R)+exp(d*(A-B)))/((1 - exp(R))*(-A*d+d*B+R If that doesn't fix the problem, then you are likely getting infinite values as result of large numbers in your exponents. Without example data to work through, I can only speculate. Jean [[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] NLS error
Also I have this problem that many times the method converges to a local
maximum...if I run it a few times, I often get the expected answer.
Is there any way to fix this problem by changing the convergence conditions?
I have tried changing the nls.control {stats} but hasn't helped much. Any
suggestions?
On Tue, Sep 20, 2011 at 5:59 PM, Diviya Smith diviya.sm...@gmail.comwrote:
I dont think *r* is related to the problem. I am not trying to estimate *r
* and so basically I am giving the model the correct value of *r* and so
log(1-r) should not go to infinity.
For test data, I generate data from the same model and add noise (using *r
norm*), with the following parameters -
A - 20
B - 500
r - 0.6881
d - (1:1000)/1000
y
-
(rnorm(d)*0.005)+(((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*+d*B+log(1-r)
And then I try estimating A, B and C (Note: C is capturing the noise term,
and I am specifying r and d) but I get the following error -
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
Any suggestions how to fix this?
Diviya
On Tue, Sep 20, 2011 at 2:37 PM, Jean V Adams jvad...@usgs.gov wrote:
Diviya Smith wrote on 09/20/2011 01:03:22 PM:
Hello there,
I am using NLS for fitting a complex model to some data to estimate a
couple
of the missing parameters. The model is -
y ~
(C+((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*d+d*B+log(1-r)
where A, B and C are unknown.
In order to test the model, I generate data by setting values for all
parameters and add some noise (C).
A - 20
B - 500
r - 0.6881
d - (1:1000)/1000
y ~
(rnorm(d)*0.005)+(((log(1-r))*exp(-A*d)*(-1+r+exp(d*(A-B)))/(r*(-A*d
+d*B+log(1-r)
I use Deoptim package to pick the optimum starting values. The model
works
fine in most cases, but every now and again, I get the following error -
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
Any suggestions on how I can resolve this? Can you suggest a better way
for
picking the starting parameters?
Thanks,
Diviya
I'm not sure if this is the problem, but if r grows greater than 1,
log(1-r) will be undefined, and you'll get an error. You can impose a
constraint on r by rewriting your formula in terms of a variable that can
take on any real value:
R = log(1-r)
So,
replace log(1-r) in your formula with R,
replace (-1 + r) with -exp(R), and
replace r with 1 - exp(R):
y ~ (C+(R*exp(-A*d)*(-exp(R)+exp(d*(A-B)))/((1 - exp(R))*(-A*d+d*B+R
If that doesn't fix the problem, then you are likely getting infinite
values as result of large numbers in your exponents. Without example data
to work through, I can only speculate.
Jean
[[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] Optimization package
Hi there, I have a complex math equation which does not have a closed form solution. It is - y - (p*exp(-a*d)*(1-exp((d-p)*(a-x[1]/((p-d)*(1-exp(-p*(a-x[1] For this equation, I have all the values except for x[1]. So I need to solve this problem numerically. Can anyone suggest an optimization package that I can use to estimate the value for x[1]? Thanks in advance, Diviya [[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] Print the summary of a model to file
Hi there, I am having a strange problem. I am running nls on some data. #data x - -(1:100)/10 y - 100 + 10 * (exp(-x / 2) Using nls I fit an exponential model to this data and get a great fit summary(fit) Formula: wcorr ~ (Y0 + a * exp(m1 * -dist/100)) Parameters: Estimate Std. Error t value Pr(|t|) Y0 -0.0001821 0.0002886 -0.6310.528 a 0.1669675 0.0015223 109.678 2e-16 *** m1 10.8045840 0.1575957 68.559 2e-16 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 0.007244 on 997 degrees of freedom Number of iterations to convergence: 7 Achieved convergence tolerance: 4.214e-06 I would like to print this summary to a file but when I use cat it adds a lot of the intermediate values and does not give me the Pvalues of the t-test. Can you suggest a good way to print this to a file? I have used- cat(paste(a), file=testing.txt, append = FALSE) output in file: wcorr ~ (Y0 + a * exp(m1 * -dist/100)) c(0.168587967690081, 0.0410369885074373, 0.015604934198253, 0.00651400975192817, -0.00198358204504023, -0.00297164058538482, 0.00200403258235102, -0.00585536622301963, -0.00677164279240891, -0.00663060500378226, -0.00429206279972985, -0.00752682816527953, -0.00772771510594769, -0.0140235396260263, -0.00905111970710343, -0.00527727528681396, -0.00637982823781946, -0.0101696023470134, -0.0101074232949501, -0.00715411863549439, -0.00662351777569056, -0.00284945195584713, -0.00401775422983583, -0.00833925944560227, 0.00724377249911614 c(3, 997) c(0.001587401493019, 3.32626693693668e-05, 0.412938194338572, 3.32626693693668e-05, 0.0441665508678859, 2.86654588531843, 0.412938194338572, 2.86654588531843, 473.324480704723) nls(formula = wcorr ~ (Y0 + a * exp(m1 * -dist/100)), start = pa1, control = list(maxiter = 1000, tol = 1e-05, minFactor = 0.0009765625, printEval = FALSE, warnOnly = TRUE), algorithm = default, trace = FALSE) list(isConv = TRUE, finIter = 7, finTol = 4.21396379219722e-06, stopCode = 0, stopMessage = converged) list(maxiter = 1000, warnOnly = TRUE) NULL c(-0.00018212632256334, 0.166967461792613, 10.8045839935834, 0.000288607886496774, 0.00152233960007251, 0.157595671762849, -0.631051094181987, 109.678196497457, 68.5588878978997, 0.528151752824769, 0, 0) c(-0.00018212632256334, 0.166967461792613, 10.8045839935834, 0.000288607886496774, 0.00152233960007251, 0.157595671762849, -0.631051094181987, 109.678196497457, 68.5588878978997, 0.528151752824769, 0, 0) Thanks, Diviya [[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] NLS fit for exponential distribution
Hello there, I am trying to fit an exponential fit using Least squares to some data. #data x - c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) y - c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, -0.06, -0.004626, -0.004626, -0.004626, -0.004626) sub - data.frame(x,y) #If model is y = a*exp(-x) + b then fit - nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0), trace = TRUE) This works well. However, if I want to fit the model : y = a*exp(-mx)+c then I try - fit - nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m= 0), trace = TRUE) It fails and I get the following error - Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates Any suggestions how I can fix this? Also next I want to try to fit a sum of 2 exponentials to this data. So the new model would be y = a*exp[(-m1+ m2)*x]+c . Any suggestion how I can do this... Any help would be most appreciated. Thanks in advance. Diviya [[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] Likelihood ratio test
Hello there, I want to perform a likelihood ratio test to check if a single exponential or a sum of 2 exponentials provides the best fit to my data. I am new to R programming and I am not sure if there is a direct function for doing this and whats the best way to go about it? #data x - c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100) y - c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436, -0.06, -0.004626, -0.004626, -0.004626, -0.004626) data - data.frame(x,y) Specifically, I would like to test if the model1 or model2 provides the best fit to the data- model 1: y = a*exp(-m*x) + c model 2: y = a*exp(-(m1+m2)*x) + c Likelihood ratio test = L(data| model1)/ L(data | model2) Any help would be most appreciated. Thanks in advance. Diviya [[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] Error in NLS example in the documentation
Hello there, I am trying to use R function NLS to analyze my data and one of the examples in the documentation is - ## the nls() internal cheap guess for starting values can be sufficient: x - -(1:100)/10 y - 100 + 10 * exp(x / 2) + rnorm(x)/10 nlmod - nls(y ~ Const + A * exp(B * x), trace=TRUE) plot(x,y, main = nls(*), data, true function and fit, n=100) curve(100 + 10 * exp(x / 2), col=4, add = TRUE) lines(x, predict(nlmod), col=2) However, when I try to execute this - I get the following error- Error in nls(y ~ Const + A * exp(B * x), trace = TRUE) : step factor 0.000488281 reduced below 'minFactor' of 0.000976562 In addition: Warning message: In nls(y ~ Const + A * exp(B * x), trace = TRUE) : No starting values specified for some parameters Can someone propose how I can fix this error? Thanks. Diviya [[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] NLS fit for exponential distribution
Thank you for very informative answers. This is great help. Peter, thanks
for correcting the model. That was exactly what I meant - apologies for the
typo.
I was able to run the analysis on some simulated data and the subset of data
that I had posted earlier. However, when I apply the analysis to full
dataset, I usually get the following errors-
fm2 - nls(y ~ b + a * (exp(m1 * -x) + exp(m2 * -x)), start = list(a = 1, b
= 0, m1 = 0.1, m2 = 0.5), trace = TRUE)
84.42045 : 1.0 0.0 0.1 0.5
0.1593364 : 0.072399953 -0.000226868 0.101135423 0.571436093
0.09032837 : 7.897930e-02 7.905359e-06 1.271435e-01 1.872206e+00
0.06834495 : 0.0808842841 0.0005815865 0.1725302473 4.9048478439
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
OR
Error in nls(y ~ b + a * (exp(m1 * -x) + exp(m2 * -x)), start = list(a = 1,
:
number of iterations exceeded maximum of 50
Is there a way to optimize the start values? I was unable to try the RStudio
solution as RStudio keeps crashing on my computer.
Also is there any way to set the start value in NLS - I mean if I want to
fit the model to only a subset of data starting at x = 10, how can I specify
that?
Thanks,
Diviya
On Sun, Jun 12, 2011 at 7:19 PM, Dennis Murphy djmu...@gmail.com wrote:
Hi:
If you use RStudio, then you can use its manipulate package to figure
out starting values for the model visually through sliders. Walmes
Zeviani posted the template of how to do this last week; I've just
adapted it for the models under consideration. It's interesting to
play with this because it shows how strongly some of the parameters
are correlated with one another...and it's fun :) Thanks to Walmes for
the excellent tutorial example.
Firstly, x and y (or the inputs in general) need to be in the global
environment as vectors of the same length.
[RStudio is a GUI, so I'm assuming that you run things from its
command line.] Load the manipulate package (which comes with RStudio)
and use the manipulate() function to create a plot of the data and fit
a curve to it. The sliders adjust the parameter values and you play
with them until the curve fits closely to the data. The range of the
sliders should match the range of plausible values you think they
might take. The output of the manipulate() function is a vector of
starting values that you can pass into nls(), which is done by closing
the window containing the sliders.
# Model 1:
x - c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
y - c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436,
-0.06, -0.004626, -0.004626, -0.004626, -0.004626)
plot(x, y) # Initial plot of the data
start - list() # Initialize an empty list for the starting values
library(manipulate)
manipulate(
{
plot(x, y)
k - kk; b0 - b00; b1 - b10
curve(k*exp(-b1*x) + b0, add=TRUE)
start - list(k=k, b0=b0, b1=b1)
},
kk=slider(0.01, 0.2, step = 0.01, initial = 0.1),
b10=slider(0.01, 0.4, step = 0.01, initial = 0.3),
b00=slider(-0.01, -0.001, step=0.001,initial= -0.005))
# When done, start() is a list of named parameters in
# the global environment
# Model fit:
fit1 - nls(y ~ k*exp(-b1*x) + b0, start = start)
summary(fit1)
### Model 2: [following Peter Dalgaard's suggested model]
### Use the estimates from fit1 and shrink b1 to anticipate
### potential effect of b2; the initial estimate in the slider for
### b1 will be too big, so it needs to be shrunk (a lot :)
start - list()
manipulate(
{
plot(x, y)
k - kk; b0 - b00; b1 - b10; b2 - b20
curve(k*(exp(-b1*x) + exp(-b2*x)) + b0, add=TRUE)
start - list(k=k, b0=b0, b1=b1, b2 = b2)
},
kk=slider(0.01, 0.1, step = 0.01, initial = 0.04),
b10=slider(0.01, 0.4, step = 0.01, initial = 0.3),
b20 = slider(0.01, 0.1, step = 0.005, initial = 0.05),
b00=slider(-0.01, -0.001, step=0.001,initial= -0.004))
fit2 - nls(y ~ k*(exp(-b1*x) + exp(-b2*x)) + b0, start = start)
summary(fit2)
anova(fit1, fit2)
HTH,
Dennis
On Sun, Jun 12, 2011 at 9:57 AM, Diviya Smith diviya.sm...@gmail.com
wrote:
Hello there,
I am trying to fit an exponential fit using Least squares to some data.
#data
x - c(1 ,10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
y - c(0.033823, 0.014779, 0.004698, 0.001584, -0.002017, -0.003436,
-0.06, -0.004626, -0.004626, -0.004626, -0.004626)
sub - data.frame(x,y)
#If model is y = a*exp(-x) + b then
fit - nls(y ~ a*exp(-x) + b, data = sub, start = list(a = 0, b = 0),
trace
= TRUE)
This works well. However, if I want to fit the model : y = a*exp(-mx)+c
then
I try -
fit - nls(y ~ a*exp(-m*x) + b, data = sub, start = list(a = 0, b = 0, m=
0), trace = TRUE)
It fails and I get the following error -
Error
