R/SynthesizeData.R

In boris109able/ChangePointCalc: Estimate change-points via Sparse Group Lasso and Dynamic Programming.

SynthesizeData <- function(ChangePoints, sigma1, sigma2, n=1000, p=20, s=4, coefficients=NULL)
{
  #This function generates synthesized data, where each x_i is generated by Gaussian distribution
  #ChangePoints is a list of change-points exclulding 1
  #sigma1 is the standard deviation of data x
  #sigma2 is the standard deviation of noise
  #n is the number of samples
  #p is the dimension of coefficients
  #coefficients is a matrix of size K*+1 by p of real coefficients provided by user
  #coefficients by default is not provided and it is all 2 or -2
  
  np <- n*p
  KReal <- length(ChangePoints) #number of change-points, excluding 1.
  
  #TmpPoints is a variant of ChangePoints and just append the 1 and n+1 from the beginning and end
  #This is just for convenience
  #The length of TmpPoints is KReal+2
  TmpPoints <- c(1, ChangePoints, n+1)
  
  x <- matrix(rnorm(n*p)*sigma1,ncol=p,nrow=n)
  #construct XTilde
  #if(type=="SGL") {
  #XTilde <- matrix(0, n, np)
  #for (i in 1:n)
  #{
  #  for (j in 1:i)
  #  {
  #    XTilde[i, ((j-1)*p+1):(j*p)] <- x[i,]
  #  }
  #}}
  
  #construct alpha and beta
  #coefficients is the real value of each interval i.e. alpha
  if(is.null(coefficients))
  {
    coefficients <- matrix(0, KReal+1, p)
    for (i in 1:(KReal+1))
    {
      if (i %% 2 == 1)
      {
        coefficients[i,1:s] <- 2 * rep.int(1, s)
      }
      else
      {
        coefficients[i,1:s] <- -2 * rep.int(1, s)
      }
    }
  }
  beta <- matrix(NA, n, p)
  for (i in 1:(KReal+1))
  {
    beta[TmpPoints[i]:(TmpPoints[i+1]-1), ] <- t(matrix(rep(coefficients[i,], TmpPoints[i+1]-TmpPoints[i]), p, TmpPoints[i+1]-TmpPoints[i]))  
  }

  #construct y
  y <- matrix(NA, n, 1)
  for (i in 1:n)
  {
    y[i] <- (x[i, ]%*%beta[i, ])[1]+matrix(rnorm(1)*sigma2,ncol=1,nrow=1)
  }
  #if(type=="SGL")
  #{data <- list(x=XTilde, y=y, coefficients=coefficients, beta=beta)}
  #if(type=="DP")
  #{data <- list(x=x, y=y, coefficients=coefficients, beta=beta)}
  data <- list(x=x, y=y, coefficients=coefficients, beta=beta)
  return(data)
}