Robust CI Computation Issue

[ghost]

Stéphane pointed out one possible bug in function wvar(): Roberto, could you please take a look at the method that we used to calculate the variance? James thinks you might be familiar with the method.

[robertomolinari]

J-dog, I can’t access the cpp code to see how the WV CIs are computed. Could you send me the code or the formulas used? From what I recall, we used the classic bounds: multiplying the lower bound by eff and dividing the upper bound by eff. In this case we could use sqrt(eff) or sqrt(sqrt(eff)).

On 25 Jul 2015, at 01:28, wyang40 notifications@github.com<mailto:notifications@github.com> wrote:

Assigned #67https://github.com/TimeSeriesLab/GMWM/issues/67 to @robertomolinarihttps://github.com/robertomolinari.

— Reply to this email directly or view it on GitHubhttps://github.com/TimeSeriesLab/GMWM/issues/67#event-365029309.

[coatless]

See wave_variance.cpp

Specifics:

//' @title Generate eta3 confidence interval
//' @description Computes the eta3 CI
//' @param y A \code{vec} that computes the brickwalled modwt dot product of each wavelet coefficient divided by their length.
//' @param dims A \code{String} indicating the confidence interval being calculated.
//' @param alpha_ov_2 A \code{double} that indicates the \eqn{\left(1-p\right)*\alpha}{(1-p)*alpha} confidence level 
//' @return A \code{matrix} with the structure:
//' \itemize{
//'  \item{Column 1}{Wavelet Variance}
//'  \item{Column 2}{Chi-squared Lower Bounds}
//'  \item{Column 3}{Chi-squared Upper Bounds}
//' }
//' @examples
//' x = rnorm(100)
//' # Uses the internal MODWT function not associated with an S3 class.
//' decomp = modwt_cpp(x, "haar", 4, "periodic")
//' signal_modwt_bw = brick_wall(decomp, select_filter("haar"), "modwt")
//' y = wave_variance(signal_modwt_bw)
//' ci_wave_variance(signal_modwt_bw, y, type = "eta3", p = 0.025)
// [[Rcpp::export]]
arma::mat ci_eta3(arma::vec y,  arma::vec dims, double alpha_ov_2) {

    unsigned int num_elem = dims.n_elem;

    arma::mat out(num_elem, 3);

    for(unsigned int i = 0; i<num_elem;i++){
      double eta3 = std::max(dims(i)/pow(2,i+1),1.0);
      out(i,1) = eta3 * y(i)/R::qchisq(1-alpha_ov_2, eta3, 1, 0); // Lower CI
      out(i,2) = eta3 * y(i)/R::qchisq(alpha_ov_2, eta3, 1, 0); // Upper CI
    }

    out.col(0) = y;

    return out;
}

//' @title Generate eta3 robust confidence interval
//' @description Computes the eta3 robust CI
//' @param y A \code{vec} that computes the brickwalled modwt dot product of each wavelet coefficient divided by their length.
//' @param dims A \code{String} indicating the confidence interval being calculated.
//' @param alpha_ov_2 A \code{double} that indicates the \eqn{\left(1-p\right)*\alpha}{(1-p)*alpha} confidence level
//' @param eff A \code{double} that indicates the efficiency.
//' @return A \code{matrix} with the structure:
//' \itemize{
//'  \item{Column 1}{Robust Wavelet Variance}
//'  \item{Column 2}{Chi-squared Lower Bounds}
//'  \item{Column 3}{Chi-squared Upper Bounds}
//' }
//' @examples
//' x = rnorm(100)
//' # Uses the internal MODWT function not associated with an S3 class.
//' decomp = modwt_cpp(x, "haar", 4, boundary="periodic")
//' signal_modwt_bw = brick_wall(decomp, select_filter("haar"), "modwt")
//' y = wave_variance(signal_modwt_bw, robust = TRUE,  eff = 0.6)
//' ci_wave_variance(signal_modwt_bw, y, type = "eta3", alpha_ov_2 = 0.025, robust = TRUE, eff = 0.6)
// [[Rcpp::export]]
arma::mat ci_eta3_robust(arma::vec y, arma::vec dims, double alpha_ov_2, double eff) {

    unsigned int num_elem = dims.n_elem;

    arma::mat out(num_elem, 3);
    eff = sqrt(eff);
    for(unsigned int i = 0; i<num_elem;i++){
      double eta3 = std::max(dims(i)/pow(2,i+1),1.0);
      out(i,1) = eff * eta3 * y(i)/(R::qchisq(1-alpha_ov_2, eta3, 1, 0)); // Lower CI
      out(i,2) = eta3 * y(i)/(eff*R::qchisq(alpha_ov_2, eta3, 1, 0)); // Upper CI
    }

    out.col(0) = y;

    return out;
}

//' @title Generate a Confidence intervval for a Univariate Time Series
//' @description Computes an estimate of the multiscale variance and a chi-squared confidence interval
//' @param signal_modwt_bw A \code{field<vec>} that contains the brick walled modwt or dwt decomposition
//' @param y A \code{vec} that contains the wave variance.
//' @param type A \code{String} indicating the confidence interval being calculated.
//' @param alpha_ov_2 A \code{double} that indicates the \eqn{\left(1-p\right)*\alpha}{(1-p)*alpha} confidence level.
//' @param robust A \code{boolean} to determine the type of wave estimation.
//' @param eff A \code{double} that indicates the efficiency.
//' @return A \code{matrix} with the structure:
//' \itemize{
//'  \item{Column 1}{Wavelet Variance}
//'  \item{Column 2}{Chi-squared Lower Bounds}
//'  \item{Column 3}{Chi-squared Upper Bounds}
//' }
//' @details 
//' This function can be expanded to allow for other confidence interval calculations.
//' @examples
//' set.seed(1337)
//' x = rnorm(100)
//' # Uses the internal MODWT function not associated with an S3 class.
//' decomp = modwt_cpp(x, "haar", 4, boundary="periodic")
//' signal_modwt_bw = brick_wall(decomp, select_filter("haar"), "modwt")
//' y = wave_variance(signal_modwt_bw)
//' ci_wave_variance(signal_modwt_bw, y, type = "eta3", p = 0.025)
// [[Rcpp::export]]
arma::mat ci_wave_variance(const arma::field<arma::vec>& signal_modwt_bw, const arma::vec& y,
                            std::string type = "eta3", double alpha_ov_2 = 0.025, bool robust = false, double eff = 0.6){

  unsigned int nb_level = y.n_elem;
  arma::vec dims(nb_level);

  for(unsigned int i = 0; i < nb_level; i++){
    dims(i) = signal_modwt_bw(i).n_elem;
  }

  arma::mat out(y.n_elem , 3);

  if(type == "eta3"){
      if(robust){
        out = ci_eta3_robust(y, dims, alpha_ov_2, eff);
      }
      else{
        out = ci_eta3(y, dims, alpha_ov_2);  
      }
  }
  else{
      stop("The wave variance type supplied is not supported. Please use: eta3");
  }

  return out;
}

//' @title Generate a Wave Variance for a Univariate Time Series
//' @description Computes an estimate of the wave variance
//' @param signal_modwt_bw A \code{field<vec>} that contains the brick walled modwt or dwt decomposition
//' @param robust A \code{boolean} to determine the type of wave estimation.
//' @param eff A \code{double} that indicates the efficiency.
//' @return A \code{vec} that contains the wave variance.
//' @examples
//' set.seed(1337)
//' x = rnorm(100)
//' signal_modwt_bw = brick_wall(modwt_cpp(x), haar_filter())
//' wave_variance(signal_modwt_bw)
//' 
//' wave_variance(signal_modwt_bw, robust = TRUE, eff = 0.6)
// [[Rcpp::export]]
arma::vec wave_variance(const arma::field<arma::vec>& signal_modwt_bw, bool robust = false, double eff = 0.6){

  unsigned int nb_level = signal_modwt_bw.n_elem;
  arma::vec y(nb_level);

  if(robust){
    // Robust wavelet variance estimation
    for(unsigned int i=0; i < nb_level; i++){
      arma::vec wav_coef = sort(signal_modwt_bw(i));
      y(i) = sig_rob_bw(wav_coef, eff);
    }
  }else{
    // Classical wavelet variance estimation
    for(unsigned int i=0; i < nb_level;i++){
      arma::vec temp = signal_modwt_bw(i);
      y(i) = dot(temp,temp)/temp.n_elem;
    }
  }

  return y;
}

//' @title Computes the (MODWT) wavelet variance
//' @description Calculates the (MODWT) wavelet variance
//' @param signal_modwt A \code{field<vec>} that contains the modwt decomposition.
//' @param robust A \code{boolean} that triggers the use of the robust estimate.
//' @param eff A \code{double} that indicates the efficiency as it relates to an MLE.
//' @param alpha A \code{double} that indicates the \eqn{\left(1-p\right)*\alpha}{(1-p)*alpha} confidence level 
//' @param ci_type A \code{String} indicating the confidence interval being calculated. Valid value: "eta3"
//' @param strWavelet A \code{String} indicating the type of wave filter to be applied. Must be "haar"
//' @return A \code{mat} with the structure:
//' \itemize{
//'   \item{"variance"}{Wavelet Variance}
//'   \item{"low"}{Lower CI}
//'   \item{"high"}{Upper CI}
//' }
//' @details 
//' This function powers the wvar object. It is also extendable...
//' @examples
//' x=rnorm(100)
//' decomp = modwt(x)
//' wvar_cpp(decomp$data, robust = FALSE)
// [[Rcpp::export]]
arma::mat wvar_cpp(const arma::field<arma::vec>& signal_modwt,
                   bool robust=false, double eff=0.6, double alpha = 0.05, 
                   std::string ci_type="eta3", std::string strWavelet="haar") {
  double alpha_ov_2 = alpha/2.0;
  // MODWT transform
  arma::field<arma::vec> signal_modwt_bw = brick_wall(signal_modwt, select_filter(strWavelet), "modwt");

  // Wavelet Variance
  arma::vec y = wave_variance(signal_modwt_bw, robust, eff);

  // Confidence Interval
  return ci_wave_variance(signal_modwt_bw, y, ci_type, alpha_ov_2, robust, eff);
}

[robertomolinari]

Try implementing the change I put in red… (i.e. sqrt(eff))

On 25 Jul 2015, at 20:07, James J Balamuta notifications@github.com<mailto:notifications@github.com> wrote:

See wave_variance.cpphttps://github.com/TimeSeriesLab/GMWM/blob/master/src/wave_variance.cpp

Specifics:

//' @title Generate eta3 confidence interval //' @description Computes the eta3 CI //' @param y A \code{vec} that computes the brickwalled modwt dot product of each wavelet coefficient divided by their length. //' @param dims A \code{String} indicating the confidence interval being calculated. //' @param alpha_ov2 A \code{double} that indicates the \eqn{\left(1-p\right)\alpha}{(1-p)_alpha} confidence level //' @return A \code{matrix} with the structure: //' \itemize{ //' \item{Column 1}{Wavelet Variance} //' \item{Column 2}{Chi-squared Lower Bounds} //' \item{Column 3}{Chi-squared Upper Bounds} //' } //' @examples //' x = rnorm(100) //' # Uses the internal MODWT function not associated with an S3 class. //' decomp = modwt_cpp(x, "haar", 4, "periodic") //' signal_modwt_bw = brick_wall(decomp, select_filter("haar"), "modwt") //' y = wave_variance(signal_modwt_bw) //' ci_wave_variance(signal_modwt_bw, y, type = "eta3", p = 0.025) // [[Rcpp::export]] arma::mat ci_eta3(arma::vec y, arma::vec dims, double alpha_ov_2) {

unsigned int num_elem = dims.n_elem;

arma::mat out(num_elem, 3);

for(unsigned int i = 0; i<num_elem;i++){
  double eta3 = std::max(dims(i)/pow(2,i+1),1.0);
  out(i,1) = eta3 * y(i)/R::qchisq(1-alpha_ov_2, eta3, 1, 0); // Lower CI
  out(i,2) = eta3 * y(i)/R::qchisq(alpha_ov_2, eta3, 1, 0); // Upper CI
}

out.col(0) = y;

return out;

}

//' @title Generate eta3 robust confidence interval //' @description Computes the eta3 robust CI //' @param y A \code{vec} that computes the brickwalled modwt dot product of each wavelet coefficient divided by their length. //' @param dims A \code{String} indicating the confidence interval being calculated. //' @param alpha_ov2 A \code{double} that indicates the \eqn{\left(1-p\right)\alpha}{(1-p)_alpha} confidence level //' @param eff A \code{double} that indicates the efficiency. //' @return A \code{matrix} with the structure: //' \itemize{ //' \item{Column 1}{Robust Wavelet Variance} //' \item{Column 2}{Chi-squared Lower Bounds} //' \item{Column 3}{Chi-squared Upper Bounds} //' } //' @examples //' x = rnorm(100) //' # Uses the internal MODWT function not associated with an S3 class. //' decomp = modwt_cpp(x, "haar", 4, boundary="periodic") //' signal_modwt_bw = brick_wall(decomp, select_filter("haar"), "modwt") //' y = wave_variance(signal_modwt_bw, robust = TRUE, eff = 0.6) //' ci_wave_variance(signal_modwt_bw, y, type = "eta3", alpha_ov_2 = 0.025, robust = TRUE, eff = 0.6) // [[Rcpp::export]] arma::mat ci_eta3_robust(arma::vec y, arma::vec dims, double alpha_ov_2, double eff) {

unsigned int num_elem = dims.n_elem;

arma::mat out(num_elem, 3);
eff = sqrt(eff);
for(unsigned int i = 0; i<num_elem;i++){
  double eta3 = std::max(dims(i)/pow(2,i+1),1.0);
  out(i,1) = sqrt(eff) * eta3 * y(i)/(R::qchisq(1-alpha_ov_2, eta3, 1, 0)); // Lower CI
  out(i,2) = eta3 * y(i)/(sqrt(eff)*R::qchisq(alpha_ov_2, eta3, 1, 0)); // Upper CI
}

out.col(0) = y;

return out;

}

//' @title Generate a Confidence intervval for a Univariate Time Series //' @description Computes an estimate of the multiscale variance and a chi-squared confidence interval //' @param signal_modwt_bw A \code{field} that contains the brick walled modwt or dwt decomposition //' @param y A \code{vec} that contains the wave variance. //' @param type A \code{String} indicating the confidence interval being calculated. //' @param alpha_ov2 A \code{double} that indicates the \eqn{\left(1-p\right)\alpha}{(1-p)_alpha} confidence level. //' @param robust A \code{boolean} to determine the type of wave estimation. //' @param eff A \code{double} that indicates the efficiency. //' @return A \code{matrix} with the structure: //' \itemize{ //' \item{Column 1}{Wavelet Variance} //' \item{Column 2}{Chi-squared Lower Bounds} //' \item{Column 3}{Chi-squared Upper Bounds} //' } //' @details //' This function can be expanded to allow for other confidence interval calculations. //' @examples //' set.seed(1337) //' x = rnorm(100) //' # Uses the internal MODWT function not associated with an S3 class. //' decomp = modwt_cpp(x, "haar", 4, boundary="periodic") //' signal_modwt_bw = brick_wall(decomp, select_filter("haar"), "modwt") //' y = wave_variance(signal_modwt_bw) //' ci_wave_variance(signal_modwt_bw, y, type = "eta3", p = 0.025) // [[Rcpp::export]] arma::mat ci_wave_variance(const arma::fieldarma::vec& signal_modwt_bw, const arma::vec& y, std::string type = "eta3", double alpha_ov_2 = 0.025, bool robust = false, double eff = 0.6){

unsigned int nb_level = y.n_elem; arma::vec dims(nb_level);

for(unsigned int i = 0; i < nb_level; i++){ dims(i) = signal_modwt_bw(i).n_elem; }

arma::mat out(y.n_elem , 3);

if(type == "eta3"){ if(robust){ out = ci_eta3_robust(y, dims, alpha_ov_2, eff); } else{ out = ci_eta3(y, dims, alpha_ov_2); } } else{ stop("The wave variance type supplied is not supported. Please use: eta3"); }

return out; }

//' @title Generate a Wave Variance for a Univariate Time Series //' @description Computes an estimate of the wave variance //' @param signal_modwt_bw A \code{field} that contains the brick walled modwt or dwt decomposition //' @param robust A \code{boolean} to determine the type of wave estimation. //' @param eff A \code{double} that indicates the efficiency. //' @return A \code{vec} that contains the wave variance. //' @examples //' set.seed(1337) //' x = rnorm(100) //' signal_modwt_bw = brick_wall(modwt_cpp(x), haar_filter()) //' wave_variance(signal_modwt_bw) //' //' wave_variance(signal_modwt_bw, robust = TRUE, eff = 0.6) // [[Rcpp::export]] arma::vec wave_variance(const arma::fieldarma::vec& signal_modwt_bw, bool robust = false, double eff = 0.6){

unsigned int nb_level = signal_modwt_bw.n_elem; arma::vec y(nb_level);

if(robust){ // Robust wavelet variance estimation for(unsigned int i=0; i < nb_level; i++){ arma::vec wav_coef = sort(signal_modwt_bw(i)); y(i) = sig_rob_bw(wav_coef, eff); } }else{ // Classical wavelet variance estimation for(unsigned int i=0; i < nb_level;i++){ arma::vec temp = signal_modwt_bw(i); y(i) = dot(temp,temp)/temp.n_elem; } }

return y; }

//' @title Computes the (MODWT) wavelet variance //' @description Calculates the (MODWT) wavelet variance //' @param signalmodwt A \code{field} that contains the modwt decomposition. //' @param robust A \code{boolean} that triggers the use of the robust estimate. //' @param eff A \code{double} that indicates the efficiency as it relates to an MLE. //' @param alpha A \code{double} that indicates the \eqn{\left(1-p\right)\alpha}{(1-p)_alpha} confidence level //' @param ci_type A \code{String} indicating the confidence interval being calculated. Valid value: "eta3" //' @param strWavelet A \code{String} indicating the type of wave filter to be applied. Must be "haar" //' @return A \code{mat} with the structure: //' \itemize{ //' \item{"variance"}{Wavelet Variance} //' \item{"low"}{Lower CI} //' \item{"high"}{Upper CI} //' } //' @details //' This function powers the wvar object. It is also extendable... //' @examples //' x=rnorm(100) //' decomp = modwt(x) //' wvar_cpp(decomp$data, robust = FALSE) // [[Rcpp::export]] arma::mat wvar_cpp(const arma::fieldarma::vec& signal_modwt, bool robust=false, double eff=0.6, double alpha = 0.05, std::string ci_type="eta3", std::string strWavelet="haar") { double alpha_ov_2 = alpha/2.0; // MODWT transform arma::fieldarma::vec signal_modwt_bw = brick_wall(signal_modwt, select_filter(strWavelet), "modwt");

// Wavelet Variance arma::vec y = wave_variance(signal_modwt_bw, robust, eff);

// Confidence Interval return ci_wave_variance(signal_modwt_bw, y, ci_type, alpha_ov_2, robust, eff); }

— Reply to this email directly or view it on GitHubhttps://github.com/TimeSeriesLab/GMWM/issues/67#issuecomment-124868515.

[coatless]

Switched eff from sqrt(eff) to be sqrt(sqrt(eff)) in the robust WV.

Let me know if it's better.

[coatless]

@stephaneguerrier + @robertomolinari

Better or worse with new change?

[coatless]

See 519307b for fix.