Error in eigenvalC > 0 : invalid comparison with complex values

[Taleb Mandelbrot]

Bonjour,

En exécutant mon code R, j'ai le message d'erreur suivant : "invalid comparison with complex values". Pouvez-vous m'aider s'il vous plaît ? Merci d'avance.

[Facundo Muñoz]

Bonjour,

C'est impossible de vous aider sans voir votre code. Aidez nous à vous aider.

• Qu'est-ce qu'un code reproductible ?
• Comment insérer des données dans un message ?
• FAQ: What's a reproducible example (`reprex`) and how do I do one?

Cordialement,

[Taleb Mandelbrot]

Bonjour,
Ci-dessous le code concerné. Merci d'avance pour votre aide.

Code : Tout sélectionner

circFBM <- function(n, H, plotfBm)
 
{
 
  if(missing(n)) n <- 500
  if(missing(H)) H <- 0.6
  if(missing(plotfBm)) plotfBm <- 1
 
  ## -------------------------------------------------------------------
  ## First line of the circulant matrix, C, built via covariances of fGn
  ## -------------------------------------------------------------------
 
  lineC <- function(n, H, m)
  {
    k <- 0:(m - 1)
    H2 <- 2 * H
    v <- (abs((k - 1)/n)^H2 - 2 * (k/n)^H2 + ((k + 1)/n)^H2)/2
    ind <- c(0:(m/2 - 1), m/2, (m/2 - 1):1)
    v <- v[ind + 1]
    drop(v)
  }
 
 
  ## ---------------------
  ## Next power of two > n
  ## ---------------------
 
  m <- 2
  repeat
  {
    m <- 2 * m
    if (m >= (n - 1))
      break
  }
  stockm <- m ##
 
  ## ---------------------------------------------------------------------
  ## Research of the power of two (<2^18) such that C is definite positive
  ## ---------------------------------------------------------------------
 
  repeat
  {
    m <- 2 * m
    eigenvalC <- lineC(n, H, m)
    eigenvalC <- fft(c(eigenvalC), inverse = F)
    if ((all(eigenvalC > 0)) | (m > 2^17))
      break
  }
 
  if(m > 2^17)
  {
    cat("----> exact method, impossible!!", fill = T)
    cat("----> can't find m such that C is definite positive", fill = T)
    break
  }
  else
  {
   
    ## ----------------------------------------------------
    ## Simulation of W=(Q)^t Z, where Z leads N(0, I_m) and
    ## (Q)_{jk} = m^(-1/2) exp(-2i pi jk/m)
    ## ----------------------------------------------------
   
    ar <- rnorm(m/2 + 1)
    ai <- rnorm(m/2 + 1)
    ar[1] <- sqrt(2) * ar[1]
    ar[(m/2 + 1)] <- sqrt(2) * ar[(m/2 + 1)]
    ai[1] <- 0
    ai[(m/2 + 1)] <- 0
    ar <- c(ar[c(1:(m/2 + 1))], ar[c((m/2):2)])
    aic <- - ai
    ai <- c(ai[c(1:(m/2 + 1))], aic[c((m/2):2)])
    W <- complex(real = ar, imaginary = ai) ##
   
    ## -------------------------
    ## Reconstruction of the fGn
    ## -------------------------
   
    W <- (sqrt(eigenvalC)) * W
    fGn <- fft(W, inverse = F)
    fGn <- (1/(sqrt(2 * m))) * fGn
    fGn <- Re(fGn[c(1:n)])
    fBm <- cumsum(fGn)
    fBm[1] <- 0 ##
   
    ## -----------
    ## Plot of fBm
    ## -----------
   
    if(plotfBm == 1)
    {
      par(mfrow = c(1, 1))
      time <- (0:(n - 1)) / n
      Nchar <- as.character(n)
      Nleg <- paste(c("N=", Nchar), collapse = " ")
      Hchar <- as.character (round(H, 3))
      Hleg <- paste(c(", H=", Hchar), collapse = " ")
      NHleg <- paste(c(Nleg, Hleg), collapse = " ")
      leg <- paste(c("Path of a fractional Browian motion ---- parameters", NHleg), collapse = " : ")
      plot(time, fBm, type = "1", main = leg)
    }
    drop(fBm)            
  }
}
circFBM(n=500,H=0.7)

[Pierre-Yves Berrard]

Bonjour,

En utilisant les outils de débuggage comme browser(), j'ai repéré que l'erreur se situe ligne 46

Code : Tout sélectionner

eigenvalC > 0

Comme le message l'indique, eigenvalC est un nombre complexe et donc ça n'aurait pas de sens de le comparer à 0.

[Facundo Muñoz]

Bonjour,

J'ai réduit votre problème à une expression minimale.
Vous calculez une transformée de Fourier rapide (FFT) ce qui donne un vecteur à valeurs complexes. Ensuite vous comparez ces valeurs contre zéro. Ce qui cause l'erreur.

À l'occurrence, les valeurs complexes ont une composante imaginaire presque nulle (en dessous de la précision de l'ordinateur). Si vous êtes sûr que cette FFT doit donner des valeurs réels, alors il suffit de prendre la partie réelle systématiquement.

Code : Tout sélectionner

lineC <- function(n, H, m) {
  k <- 0:(m - 1)
  H2 <- 2 * H
  v <- (abs((k - 1)/n)^H2 - 2 * (k/n)^H2 + ((k + 1)/n)^H2)/2
  ind <- c(0:(m/2 - 1), m/2, (m/2 - 1):1)
  v <- v[ind + 1]
  drop(v)
}

eigenvalC <- lineC(500, .7, 1024)  # a numeric vector
eigenvalC <- fft(c(eigenvalC), inverse = F)  # a complex vector
all(eigenvalC > 0)  # Erreur: comparer des valeurs complexes contre 0
#> Error in eigenvalC > 0: invalid comparison with complex values

## Cependant, les composantes imaginaires sont presque nulles:
summary(Im(eigenvalC)) 
#>       Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
#> -8.047e-20 -1.444e-20 -2.303e-21  0.000e+00  1.002e-20  1.838e-19

## Si c'est le cas en général, on peut définir:
eigenvalC <- Re(fft(c(eigenvalC), inverse = F))  # a real vector
all(eigenvalC > 0)  # TRUE
#> [1] TRUE


Cordialement,

[Taleb Mandelbrot]

Ah oui vous avez raison, cela fonctionne mieux en prenant un vecteur de réels. Merci beaucoup pour votre aide.

Cordialement,