Fonction qui retourne "NULL" en guise de résultat

[Taleb Mandelbrot]

Bonjour,

En exécutant une fonction (whittST, elle me revoit un résultat "NULL" au lieu d'un nombre calculé préalablement. Ci dessous le code de la fonction, ainsi que ceux de deux autres méthodes appelées. Merci d'avance pour votre aide.

Code : Tout sélectionner

whittST <- function(fBm, Hprel)
{
         if (missing(fBm)) fBm <- circFBM(n=500, H=0.6)
         if (missing(Hprel))
         {
            Db4 <- c(0.4829629, -0.8365163, 0.22414386, 0.12940952)
            Hprel <- VaPkolST(fBm=fBm, k=2, a=Db4, M=5, 0)$Hols
         }
         fGn <- c(fBm[1], diff(fBm))
         n <- length(fGn)
         nstar <- trunc((n - 1)/2)
         perio <- (Mod(fft(fGn, inverse = F)))^2/(2 * pi * n)
         perio <- perio[1:nstar]
         global <- vector("list", 2)
         global[[1]] <- perio
         global[[2]] <- n
         assign ("global", global)  ##
         
## ---------------------
## Criterion to minimize
## ---------------------

         Q <- function(Htry)
         {
            perio <- global[[1]]
            n <- global[[2]]
            spd <- spdFGN(Htry, n)
            result <- (4 * pi)/n * sum(perio/spd)
            drop(result)
         }
         epsilon <- 1.0000000000000001e-05
         Hest <- nlminb(objective = Q, start = Hprel, lower = epsilon, upper = 1 - epsilon)$parameters   
         drop(Hest)
}
whittST(fBm=circFBM(n=1000, H=0.9), Hprel=0.9)

Code : Tout sélectionner

spdFGN <- function(Htry, n)
{
if (missing(Htry)) Htry <- 0.6
      if (missing(n)) n <- 500
      alpha <- 2 * Htry + 1
      nstar <- trunc((n - 1)/2)
      clambda <- (sin(pi * Htry) * gamma(alpha))/pi
 
  ## -----------------------------------------
  ## Computation of spd at Fourier frequencies
  ## -----------------------------------------
 
      j <- 2 * pi * ((-300):300)
      spd <- rep(0, nstar)
      for (k in (1:nstar))
      {
        lambda <- (2 * pi * k)/n
        stocksum <- sum(abs(j + lambda)^( - alpha))
        spd[k] <- clambda * (1 - cos(lambda)) * stocksum
      }
 
  ## ----------------------
  ## Renormalization of spd
  ## ----------------------
 
      theta <- exp(2/n * sum(log(spd)))
      spd <- spd/theta
      drop(spd)
}

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 <- Re(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 (Wood and Chan) ---- parameters", NHleg), collapse = " : ")
               plot(time, fBm, type = "l", main = leg)
            }
             drop(fBm)            
            }
            


}

[Facundo Muñoz]

Bonjour Taleb,

Quand quelque chose ne donne les résultats attendus, vous devez apprendre à isoler le problème pour déterminer ce qui ne va pas.

Dans votre situation, j'ai ajouté une ligne

Code : Tout sélectionner

browser()

juste à la fin de la fonction `whittST` et exécuté ensuite l'appel à la fonction. Cela arrête l'exécution juste avant la fin de la fonction, ce qui m'a permis d'explorer les valeurs de toutes les variables.
J'ai appris ainsi que l'objet `Hest` est effectivement NULL, c'est pourquoi la fonction retourne NULL.

Fait l'excercise avec `browser()` et regardez comment l'objet Hest est calculé. Fait d'abord appel à la fonction `nlminb()` (sans aller chercher les `$parameters`) et regardez qu'est-ce que ça donne. Je vous laisse trouver la solution. Regardez aussi la fiche d'aide de la fonction `nlminb()` pour voir la structure attendue de l'objet retourné.

Si vous lisez l'anglais, vous pouvez creuser sur le sujet de débogage ici [1].

[1] http://adv-r.had.co.nz/Exceptions-Debugging.html

Cordialement,