Re: [R] fitting cosine curve

2017-06-21 Thread Charles C. Berry

On Wed, 21 Jun 2017, J C Nash wrote:


Using a more stable nonlinear modeling tool will also help, but key is to get
the periodicity right.



The model is linear up to omega after transformation as Don and I noted.

Taking a guess that 2*pi/240 = 0.0262 is about right for omega:


rsq <- function(x) {t2<-t*x;summary(lm(y~cos(t2)+sin(t2)))$r.squared}
vrsq <- Vectorize(rsq)
optimise(rsq, c(0.8,1.2)*2*pi/240,maximum=TRUE)

$maximum
[1] 0.02794878

$objective
[1] 0.8127072


curve(vrsq,0.025,0.03)



Isn't this stable enough?

And as you note

plot(lm(y~cos(t*0.0279)+sin(t*0.0279)))

reveals lack-of-fit.

Of course there are some other issues not addressed here such as 
possible autoregression.


HTH,

Chuck


y=c(16.82, 16.72, 16.63, 16.47, 16.84, 16.25, 16.15, 16.83, 17.41, 17.67,
17.62, 17.81, 17.91, 17.85, 17.70, 17.67, 17.45, 17.58, 16.99, 17.10)
t=c(7,  37,  58,  79,  96, 110, 114, 127, 146, 156, 161, 169, 176, 182,
190, 197, 209, 218, 232, 240)

lidata <- data.frame(y=y, t=t)

#I use the method to fit a curve, but it is different from the real curve,
#which can be seen in the figure.
linFit  <- lm(y ~ cos(t))
library(nlsr)
#fullFit <- nls(y ~ A*cos(omega*t+C) + B,
#start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.4))
#omega cannot be set to 1, don't know why.
fullFit <- nlxb(y ~ A*cos(omega*t+C) + B, data=lidata,
start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.04), trace=TRUE)
co <- coef(fullFit)
fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
plot(x=t, y=y)
curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
,lwd=2, col="steelblue")
jstart <- list(A=20, B=100, C=0, omega=0.01)
jfit <- nlxb(y ~ A*cos(omega*t+C) + B, data=lidata,
   start=jstart, trace=TRUE)
co <- coef(jfit)
fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
plot(x=t, y=y)
curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
 ,lwd=2, col="steelblue")

JN

On 2017-06-21 12:06 AM, lily li wrote:

I'm trying the different parameters, but don't know what the error is:
Error in nlsModel(formula, mf, start, wts) :
   singular gradient matrix at initial parameter estimates

Thanks for any suggestions.

On Tue, Jun 20, 2017 at 7:37 PM, Don Cohen 
wrote:



If you know the period and want to fit phase and amplitude, this is
equivalent to fitting a * sin + b * cos

  > >>> > I don't know how to set the approximate starting values.

I'm not sure what you meant by that, but I suspect it's related to
phase and amplitude.

  > >>> > Besides, does the method work for sine curve as well?

sin is the same as cos with a different phase
Any combination of a and b above = c * sin (theta + d) for
some value of c and d and = e * cos (theta + f) for some value
of e and f.
Also for any c,d and for any e,f there is an a,b.
the c and e are what I'm calling amplitude, the d and f are what
I'm calling phase.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

and provide commented, minimal, self-contained, reproducible code.





Charles C. Berry Dept of Family Medicine & Public Health
cberry at ucsd edu   UC San Diego / La Jolla, CA 92093-0901
http://biostat.ucsd.edu/ccberry.htm

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


Re: [R] fitting cosine curve

2017-06-21 Thread J C Nash

Using a more stable nonlinear modeling tool will also help, but key is to get
the periodicity right.

y=c(16.82, 16.72, 16.63, 16.47, 16.84, 16.25, 16.15, 16.83, 17.41, 17.67,
17.62, 17.81, 17.91, 17.85, 17.70, 17.67, 17.45, 17.58, 16.99, 17.10)
t=c(7,  37,  58,  79,  96, 110, 114, 127, 146, 156, 161, 169, 176, 182,
190, 197, 209, 218, 232, 240)

lidata <- data.frame(y=y, t=t)

#I use the method to fit a curve, but it is different from the real curve,
#which can be seen in the figure.
linFit  <- lm(y ~ cos(t))
library(nlsr)
#fullFit <- nls(y ~ A*cos(omega*t+C) + B,
#start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.4))
#omega cannot be set to 1, don't know why.
fullFit <- nlxb(y ~ A*cos(omega*t+C) + B, data=lidata,
 start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.04), trace=TRUE)
co <- coef(fullFit)
fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
plot(x=t, y=y)
curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
,lwd=2, col="steelblue")
jstart <- list(A=20, B=100, C=0, omega=0.01)
jfit <- nlxb(y ~ A*cos(omega*t+C) + B, data=lidata,
start=jstart, trace=TRUE)
co <- coef(jfit)
fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
plot(x=t, y=y)
curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
  ,lwd=2, col="steelblue")

JN

On 2017-06-21 12:06 AM, lily li wrote:

I'm trying the different parameters, but don't know what the error is:
Error in nlsModel(formula, mf, start, wts) :
   singular gradient matrix at initial parameter estimates

Thanks for any suggestions.

On Tue, Jun 20, 2017 at 7:37 PM, Don Cohen 
wrote:



If you know the period and want to fit phase and amplitude, this is
equivalent to fitting a * sin + b * cos

  > >>> > I don't know how to set the approximate starting values.

I'm not sure what you meant by that, but I suspect it's related to
phase and amplitude.

  > >>> > Besides, does the method work for sine curve as well?

sin is the same as cos with a different phase
Any combination of a and b above = c * sin (theta + d) for
some value of c and d and = e * cos (theta + f) for some value
of e and f.
Also for any c,d and for any e,f there is an a,b.
the c and e are what I'm calling amplitude, the d and f are what
I'm calling phase.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] fitting cosine curve

2017-06-20 Thread Don Cohen

If you know the period and want to fit phase and amplitude, this is
equivalent to fitting a * sin + b * cos

 > >>> > I don't know how to set the approximate starting values.

I'm not sure what you meant by that, but I suspect it's related to
phase and amplitude.

 > >>> > Besides, does the method work for sine curve as well? 

sin is the same as cos with a different phase
Any combination of a and b above = c * sin (theta + d) for
some value of c and d and = e * cos (theta + f) for some value
of e and f.
Also for any c,d and for any e,f there is an a,b.
the c and e are what I'm calling amplitude, the d and f are what
I'm calling phase.

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


Re: [R] fitting cosine curve

2017-06-20 Thread Charles C. Berry

On Tue, 20 Jun 2017, lily li wrote:


Hi R users,

I have a question about fitting a cosine curve. I don't know how to set the
approximate starting values.


See

Y.L. Tong (1976) Biometrics 32:85-94

The method is known as `cosinor' analysis.  It takes advantage of the 
*intrinsic* linearity of y = a + b * cos( omega*t - c ), when omega is 
given.


It you are scratching your head saying 'that thing is not linear', you 
need to go back to your linear models text and review what `linearity' 
actually refers to.  Also, reading the Tong paper is recomended as you 
will need the transformations given there in any case.


What you end up doing is fitting

fit <- lm(y~cos(t.times.omega)+sin(t.times.omega))

and then transforming coef(fit) to get back a, b, and c. So, you only need 
to have omega. If it is not obvious what value to use,  then that will be 
more of a challenge.


The paper gives asymptotics for the dispersion matrix of (a, b, c), I 
recall.



Besides, does the method work for sine curve as well?


Seriously? See

https://en.wikipedia.org/wiki/List_of_trigonometric_identities#Shifts_and_periodicity

HTH,

Chuck

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


Re: [R] fitting cosine curve

2017-06-20 Thread lily li
I'm trying the different parameters, but don't know what the error is:
Error in nlsModel(formula, mf, start, wts) :
  singular gradient matrix at initial parameter estimates

Thanks for any suggestions.

On Tue, Jun 20, 2017 at 7:37 PM, Don Cohen 
wrote:

>
> If you know the period and want to fit phase and amplitude, this is
> equivalent to fitting a * sin + b * cos
>
>  > >>> > I don't know how to set the approximate starting values.
>
> I'm not sure what you meant by that, but I suspect it's related to
> phase and amplitude.
>
>  > >>> > Besides, does the method work for sine curve as well?
>
> sin is the same as cos with a different phase
> Any combination of a and b above = c * sin (theta + d) for
> some value of c and d and = e * cos (theta + f) for some value
> of e and f.
> Also for any c,d and for any e,f there is an a,b.
> the c and e are what I'm calling amplitude, the d and f are what
> I'm calling phase.
>

[[alternative HTML version deleted]]

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


Re: [R] fitting cosine curve

2017-06-20 Thread lily li
Thanks. I will do a trial first. Also, is it okay to have the datasets that
have only part of the cycle, or better to have equal or more than one
cycle? That is to say, I cannot have the complete datasets sometimes.

On Tue, Jun 20, 2017 at 7:37 PM, Don Cohen 
wrote:

>
> If you know the period and want to fit phase and amplitude, this is
> equivalent to fitting a * sin + b * cos
>
>  > >>> > I don't know how to set the approximate starting values.
>
> I'm not sure what you meant by that, but I suspect it's related to
> phase and amplitude.
>
>  > >>> > Besides, does the method work for sine curve as well?
>
> sin is the same as cos with a different phase
> Any combination of a and b above = c * sin (theta + d) for
> some value of c and d and = e * cos (theta + f) for some value
> of e and f.
> Also for any c,d and for any e,f there is an a,b.
> the c and e are what I'm calling amplitude, the d and f are what
> I'm calling phase.
>

[[alternative HTML version deleted]]

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


Re: [R] fitting cosine curve

2017-06-20 Thread Jim Lemon
What I did was to plot your initial values, then plot the smoothed
values and guess the constants. That is, I got an "eyeball" fit to the
smoothed values. As I have described this as "gross cheating" in the
past, you should either split your data, estimate on one subset and
then test on another, or estimate on your data and test on a
replication. If you get pretty much the same values, you can be
reasonably confident that they are reliable.

Jim

On Wed, Jun 21, 2017 at 9:52 AM, lily li  wrote:
> For example, how do you know the value 0.6, and the frequency within cos?
>
> On Tue, Jun 20, 2017 at 5:49 PM, lily li  wrote:
>>
>> Thanks, that is cool. But would there be a way that can approximate the
>> curve by trying more starting values automatically?
>>
>> On Tue, Jun 20, 2017 at 5:45 PM, Jim Lemon  wrote:
>>>
>>> Hi lily,
>>> You can get fairly good starting values just by eyeballing the curves:
>>>
>>> plot(y)
>>> lines(supsmu(1:20,y))
>>> lines(0.6*cos((1:20)/3+0.6*pi)+17.2)
>>>
>>> Jim
>>>
>>>
>>> On Wed, Jun 21, 2017 at 9:17 AM, lily li  wrote:
>>> > Hi R users,
>>> >
>>> > I have a question about fitting a cosine curve. I don't know how to set
>>> > the
>>> > approximate starting values. Besides, does the method work for sine
>>> > curve
>>> > as well? Thanks.
>>> >
>>> > Part of the dataset is in the following:
>>> > y=c(16.82, 16.72, 16.63, 16.47, 16.84, 16.25, 16.15, 16.83, 17.41,
>>> > 17.67,
>>> > 17.62, 17.81, 17.91, 17.85, 17.70, 17.67, 17.45, 17.58, 16.99, 17.10)
>>> > t=c(7,  37,  58,  79,  96, 110, 114, 127, 146, 156, 161, 169, 176, 182,
>>> > 190, 197, 209, 218, 232, 240)
>>> >
>>> > I use the method to fit a curve, but it is different from the real
>>> > curve,
>>> > which can be seen in the figure.
>>> > linFit  <- lm(y ~ cos(t))
>>> > fullFit <- nls(y ~ A*cos(omega*t+C) + B,
>>> > start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.4)) #omega
>>> > cannot
>>> > be set to 1, don't know why.
>>> > co <- coef(fullFit)
>>> > fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
>>> > plot(x=t, y=y)
>>> > curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
>>> > ,lwd=2, col="steelblue")
>>> >
>>> > __
>>> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> > https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] fitting cosine curve

2017-06-20 Thread lily li
Thanks, that is cool. But would there be a way that can approximate the
curve by trying more starting values automatically?

On Tue, Jun 20, 2017 at 5:45 PM, Jim Lemon  wrote:

> Hi lily,
> You can get fairly good starting values just by eyeballing the curves:
>
> plot(y)
> lines(supsmu(1:20,y))
> lines(0.6*cos((1:20)/3+0.6*pi)+17.2)
>
> Jim
>
>
> On Wed, Jun 21, 2017 at 9:17 AM, lily li  wrote:
> > Hi R users,
> >
> > I have a question about fitting a cosine curve. I don't know how to set
> the
> > approximate starting values. Besides, does the method work for sine curve
> > as well? Thanks.
> >
> > Part of the dataset is in the following:
> > y=c(16.82, 16.72, 16.63, 16.47, 16.84, 16.25, 16.15, 16.83, 17.41, 17.67,
> > 17.62, 17.81, 17.91, 17.85, 17.70, 17.67, 17.45, 17.58, 16.99, 17.10)
> > t=c(7,  37,  58,  79,  96, 110, 114, 127, 146, 156, 161, 169, 176, 182,
> > 190, 197, 209, 218, 232, 240)
> >
> > I use the method to fit a curve, but it is different from the real curve,
> > which can be seen in the figure.
> > linFit  <- lm(y ~ cos(t))
> > fullFit <- nls(y ~ A*cos(omega*t+C) + B,
> > start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.4)) #omega
> cannot
> > be set to 1, don't know why.
> > co <- coef(fullFit)
> > fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
> > plot(x=t, y=y)
> > curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
> > ,lwd=2, col="steelblue")
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] fitting cosine curve

2017-06-20 Thread Jim Lemon
Hi lily,
You can get fairly good starting values just by eyeballing the curves:

plot(y)
lines(supsmu(1:20,y))
lines(0.6*cos((1:20)/3+0.6*pi)+17.2)

Jim


On Wed, Jun 21, 2017 at 9:17 AM, lily li  wrote:
> Hi R users,
>
> I have a question about fitting a cosine curve. I don't know how to set the
> approximate starting values. Besides, does the method work for sine curve
> as well? Thanks.
>
> Part of the dataset is in the following:
> y=c(16.82, 16.72, 16.63, 16.47, 16.84, 16.25, 16.15, 16.83, 17.41, 17.67,
> 17.62, 17.81, 17.91, 17.85, 17.70, 17.67, 17.45, 17.58, 16.99, 17.10)
> t=c(7,  37,  58,  79,  96, 110, 114, 127, 146, 156, 161, 169, 176, 182,
> 190, 197, 209, 218, 232, 240)
>
> I use the method to fit a curve, but it is different from the real curve,
> which can be seen in the figure.
> linFit  <- lm(y ~ cos(t))
> fullFit <- nls(y ~ A*cos(omega*t+C) + B,
> start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.4)) #omega cannot
> be set to 1, don't know why.
> co <- coef(fullFit)
> fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
> plot(x=t, y=y)
> curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
> ,lwd=2, col="steelblue")
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] fitting cosine curve

2017-06-20 Thread lily li
Hi R users,

I have a question about fitting a cosine curve. I don't know how to set the
approximate starting values. Besides, does the method work for sine curve
as well? Thanks.

Part of the dataset is in the following:
y=c(16.82, 16.72, 16.63, 16.47, 16.84, 16.25, 16.15, 16.83, 17.41, 17.67,
17.62, 17.81, 17.91, 17.85, 17.70, 17.67, 17.45, 17.58, 16.99, 17.10)
t=c(7,  37,  58,  79,  96, 110, 114, 127, 146, 156, 161, 169, 176, 182,
190, 197, 209, 218, 232, 240)

I use the method to fit a curve, but it is different from the real curve,
which can be seen in the figure.
linFit  <- lm(y ~ cos(t))
fullFit <- nls(y ~ A*cos(omega*t+C) + B,
start=list(A=coef(linFit)[1],B=coef(linFit)[2],C=0,omega=.4)) #omega cannot
be set to 1, don't know why.
co <- coef(fullFit)
fit <- function(x, a, b, c, d) {a*cos(b*x+c)+d}
plot(x=t, y=y)
curve(fit(x, a=co['A'], b=co['omega'], c=co['C'],d=co['B']), add=TRUE
,lwd=2, col="steelblue")


curve1.pdf
Description: Adobe PDF document
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.