[vsnet-chat 7803] Re: Period Analysis using the Lasso

[Taichi Kato]

   Lasso may be too new to astronomy and may not be readily understandable.

> For the simplest case p=0, s=1,
> x(t)=a_0+(C_1*cos(w*t)+S_1*sin(w*t))
> with 3(!) parameters a_0, C_1, S_1, instead of two C_1, S_1,
> x(t)-xm=(C_1*cos(w*t)+S_1*sin(w*t))
> as generally the sample mean xm does not coincide with the best fit parameter a_0 
> (when the observations are irregularly distributed over phases).
>
> Generally, the trend a_0+a_1*(t-tm)+a_2*(t-tm)^2+...+a_p*(t-tm)^p

   Although this is not directly related to the lasso analysis,
if the number of parameters is fixed, your problem can be very easily 
(and gracefully) solved by using Markov-Chain Monte Carlo (MCMC) method.
The error estimates of parameters can be easily obtained.
We may then evaluate the fit by other method (AIC might be another
choice).

See Appendix in:
http://adsabs.harvard.edu/abs/2010PASJ...62.1525K

   A sample R code.  This is a four-parameter problem applied to BW Scl
(0.0032 and other parameters are dependent on the data).
Let's see how simple R code is (x is observation data).

psin <- function(p,z) {
    sum(-(z$V2-p[1]*sin((2*pi/p[2])*z$V1+p[3])-p[4])^2/0.0032)
}

n <- 1000

zz <- matrix(0,n,4)
pp <- numeric(n)

zp <- c(0.05814, 0.05497, 1.4109, 0.0011)
parameter <- c(0.001, 0.000003, 0.018, 0.0004)

for (i in 1:n) {
    zpnew <- zp + rnorm(4)*parameter
    p1 <- psin(zp,x)
    p2 <- psin(zpnew,x)
    r <- p2-p1
    if (r > log(runif(1))) {
        zp <- zpnew
        pp[i] <- p2
    } else {
        pp[i] <- p1
    }
    zz[i,] <- zp
}