# R/update_amplitude.r In multiband: Period Estimation for Multiple Bands

```#' Update Amplitude parameter
#'
#' \code{update_amplitude} solves a simple linear system of equations (rank 1 perturbation on a diagonal matrix) to update the amplitude estimates.
#'
#' @param tms list of matrices whose rows are the triple (t,mu,sigma) for each band
#' @param beta vector of the current intercept estimates
#' @param rho vector of the current phase estimates
#' @param omega frequency
#' @param at prior vector
#' @param gamma nonnegative regularization parameter
#' @export
#' @examples
#' test_data <- synthetic_multiband()
#' B <- test_data\$B
#' tms <- test_data\$tms
#' beta <- test_data\$beta
#' rho <- test_data\$rho
#' omega <- test_data\$omega
#' gamma <- 1
#' at <- rnorm(B)
#' at <- as.matrix(at/sqrt(sum(at**2)),ncol=1)
#'
#' a_next <- update_amplitude(tms,beta,rho,omega,at,gamma)
#'
#' e <- double(B)
#' xi <- double(B)
#' for (b in 1:B) {
#'   nb <- length(tms[[b]][,1])
#'   s <- sin(omega*tms[[b]][,1] + rho[b])
#'   w <- 1/(tms[[b]][,3]**2)
#'   e[b] <- t(s)%*%(w*s) + gamma
#'   xi[b] <- t(s)%*%(w*(tms[[b]][,2]-beta[b]))
#' }
#' a_direct <- solve(diag(e)-gamma*at%*%t(at),xi)
#' norm(as.matrix(a_direct-a_next),'f')
update_amplitude <- function(tms,beta,rho,omega,at,gamma) {
B <- length(tms)
e <- double(B)
xi <- matrix(NA,B,1)
for (b in 1:B) {
band <- tms[[b]]
nb <- length(band[,1])
s <- sin(omega*band[,1] + rho[b])
w <- 1/(band[,3]**2)
iw <- band[,3]**2
xi[b] <- t(s)%*%((band[,2] - beta[b])/iw)
e[b] <- t(s)%*%(s/iw) + gamma
}
eixi <- xi/e
eiat <- at/e
a <- eixi + gamma*(sum(eixi*at)/(1 - gamma*sum(at*eiat)))*eiat
return(a)
}
```

## Try the multiband package in your browser

Any scripts or data that you put into this service are public.

multiband documentation built on May 2, 2019, 3:30 a.m.