salut

il me semble par integration

En fait je dois faire une simulation de valeurs suivant la loi .

Et pour cela, il me faut une formule fixe de la fonction de répartition de .

Où est la difficulté ? Le calcul est un peu casse-bonbons mais il s'agit juste de couper une intégrale en deux et de faire des changements de variable

Si , alors .

Pour calculer la première intégrale, pour tout , .
Tu fais le changement de variable , tu trouves .


Pour l'autre c'est le même principe mais tu remplaces par , ça reste un calcul de primitive de fonction exponentielle

Merci Ulmiere

Finalement je trouve .

Et comme je dois simuler la loi de pour , j'ai alors calculé la bijection réciproque de :

Alors pour tout

Tu vas un peu trop vite là !
Le calcul me semble correct dans le cas où x>x0, mais quand x<=x0 le second terme n'existe pas. Ton expression ne tient pas compte des deux cas possibles.
Il y aura forcément une indicative en fonction de x0 dans l'un des deux termes.

L'autre problème c'est que f* n'est pas forcément une densité. Il faut calculer la limite de F en l'infini pour trouver le facteur de normalisation et alors tu pourras donner l'expression finale de F_r

Ensuite il faudra trouver l'inverse... S'il existe !
Ce que tu cherches est l'inverse généralisée à droite de la fonction F_r

Salut Ulmiere

Si ,



si

Un peu pessimiste à l'idée de poursuivre l'étude des différents cas de choix de dans  ..

L'inversion de la fonction de répartition devient alors assez délicate il me semble

Pour le cas , je t'ai donné le résultat
Pour le cas , j'ai cassé l'intégrale en deux. Le premier terme se calcule grâce au premier cas, puisque , ce qui donne .
Le terme qui reste est un calcul du même genre que le premier


Donc si on résume la situation jusqu'ici, ton calcul était correct et


Autrement dit,

Maintenant, regarde la définition de l'inverse généralisé à droite et applique

Je t'aide pour le début.
Calculons , où et est un réel que je te laisse calculer.

On constate déjà que F est continue à droite, ce qui est une bonne nouvelle pour une fonction de répartition. Ca veut dire que la loi est diffuse. Est-ce que c'est vrai aussi à gauche ? F étant croissante, est-elle inversible ?

Comme définition de l'inverse, on prend .
En toutes lettres : on trace une droite horizontale par dessus la courbe représentative de , on récupère la plus petite des abscisses des points d'intersection des deux courbes.

La chance qu'on a, c'est que F est croissante. Par conséquent,
* si , il s'agit de résoudre une certaine équation.
* si , même chose avec une autre équation.

Tu réunis ensuite ces deux cas sous une seule expression avec indicatrice comme nous l'avons fait pour et tu regardes si la fonction obtenue est continue.

Bonsoir Ulmiere

Si ,



Donc,

Si ,



Donc,

Finalement :

Oui, mais

il faut faire attention à bien expliquer pourquoi tu as le droit de passer tel réel à la puissance 1/(eta - 1)
il manque une indicatrice
et le facteur K a été oublié dans tes calculs

Il manque aussi une petite étude de la continuité de F^{-1}, bien que ce ne soit pas nécessaire pour appliquer une méthode de Monte Carlo

Bonsoir Ulmiere, merci beaucoup.

Concernant la simulation de la loi de , je code en langage R.

Voici un début de code :

# Fonction pour simuler n valeurs suivant la loi P_r
rremb <- function(n, params) {
  
  # Extraire les paramètres depuis la liste params
  alpha <- params$alpha
  x0 <- params$x0
  eta <- params$eta
  
  # Vérifier si les paramètres sont corrects
  if (alpha <= 0 | x0 <= 0 | eta <= 3) {
    stop("Mauvais paramètres : alpha doit être > 0, x0 > 0 et eta > 3")
  }
  
  # Fonction inverse de la fonction de répartition F^(-1)
  F_inv <- function(y) {
    if (y <= 0.5) {
      # Cas où x <= x0
      return(alpha + x0 - (y * (eta - 1))^(1 / (1 - eta)))
    } else {
      # Cas où x > x0
      return(x0 - alpha + (2 * alpha^(1 - eta) - y * (eta - 1))^(1 / (eta - 1)))
    }
  }
  
  # Générer n valeurs uniformes sur (0, 1)
  u <- runif(n)
  
  # Utiliser F_inv pour générer des valeurs selon P_r
  x <- sapply(u, F_inv)
  
  return(x)
}

# Paramètres
params <- list(alpha = 2.0, x0 = 50.0, eta = 4.0)

# Ajustement de n pour un temps d'exécution proche de 10 secondes
n <- 10000  # Exemple de valeur, à ajuster selon les performances

# Simulation des valeurs
valeurs_simulees <- rremb(n, params)

# Afficher les résultats
print(valeurs_simulees)

# Optionnel : Tracer l'histogramme des valeurs simulées
hist(valeurs_simulees, breaks = 50, main = "Histogramme des valeurs simulées", xlab = "Valeurs", ylab = "Fréquence")

Il manque toujours le facteur K de normalisation dans tous tes calculs
Quand tu l'auras pris en compte, tu pourras superposer à l'histogramme le graphe de et regarder si tout coincide

Salut Ulmiere, peux tu me donner des indications sur comment calculer ce facteur de normalisation ?

J'ai ajouté, cette ligne dans mon algorithme pour la normalisation de ,

1/integrate(f_star, 0, Inf, alpha = alpha, x0 = x0, eta = eta)$value

Il ne faut pas faire le calcul approché, sinon la complexité explose !
Le facteur de normalisation est simplement la limite en l'infini de , dont tu connais l'expression explicitement. Il suffit de calculer cette limite à la main et de modifier le calcul en prenant en compte cette constante !

Je voulais dire la limite de F* bien sûr, pas celle de F, qui vaut 1 par construction

On a

Donc

En multipliant cette valeur par dans mon code, j'ai que des 0 et des NaN..

F = (F*) / K donc

F(x) = y
ssi (F*)(x) = Ky
ssi x = (F*)^(-1)(Ky)

Ca veut donc dire que F^(-1)(y) = (F*)^(-1)(Ky)
Tu peux réutiliser les expressions que tu as déjà trouvé pour l'inverse de F* pour calculer F_inv dans ton code

L'équivalent en python (sans introduire numpy, donc très lent) serait en gros

K = 2 * alpha ** (1-eta) / float(eta - 1)

def F_star_inv(y):
  ..
  return ...


def F_inv(y):
  return F_start_inv(K * y)


def monte_carlo(n,  s):
  # B est l'évènement R>S, de probabilité P(B) = 1 - F(s)
  # alors E(R-s | B) = E((R-s)1_B) / P(B) = E((R-s)_+) / P(B)  

  U = generer_unif(n)
  V = filter(lambda r: r>s, map(F_inv, U))
  W = map(lambda r: r-s, V) # v.a iid de même loi que (R-s)_+

  m = sum(V) / float(n) # monte carlo / estimateur de la moyenne de W
  p = 1 - F(s) # P(B)

  return m / float(p)

Changeons de notation alors. Soient K est une constante et f et g deux fonctions inversibles sur R+ telles que f = g / K.
On suppose g^(-1), la réciproque de g, connue et on cherche f^(-1), la réciproque de f.

Soient x et y deux réels. Pour trouver f^(-1)(y), on utilise comme d'habitude une série d'équivalences
f(x) = y
ssi machin
ssi bidule
ssi x = truc
et on en déduit alors que f^(-1)(y) = truc, parce que f(x) = y équivaut aussi à f^(-1)(y) = x

Dans notre cas, f(x) = y
ssi g(x) / K = y
ssi g(x) = Ky
ssi x = g^(-1)(Ky)

Ce qui veut dire que f^(-1)(y) = g^(-1)(Ky), pour tout y où cette égalité a un sens.

Reste plus qu'à appliquer ceci avec f = F et g = F*

rremb <- function(n, params) {
  # Extraire les paramètres depuis la liste params
  alpha <- params$alpha
  x0 <- params$x0
  eta <- params$eta
  
  # Vérifier si les paramètres sont corrects
  if (alpha <= 0 | x0 <= 0 | eta <= 3) {
    stop("Mauvais paramètres : alpha doit être > 0, x0 > 0 et eta > 3")
  }
  
  # Facteur de normalisation K
  K <- 2 * alpha^(1 - eta) / (eta - 1)
  
  # Fonction inverse de F* (sans facteur de normalisation)
  F_star_inv <- function(z) {
    if (z <= 0.5) {
      # Cas où x <= x0
      return((alpha + x0 - (z * (eta - 1))^(1 / (1 - eta))))
    } else {
      # Cas où x > x0
      return((x0 - alpha + ((1 - z) * (eta - 1))^(1 / (eta - 1))))
    }
  }
  
  # Fonction inverse de F en tenant compte de K
  F_inv <- function(y) {
    return(F_star_inv(K * y))
  }
  
  # Générer n valeurs uniformes sur (0, 1)
  u <- runif(n)
  
  # Utiliser F_inv pour générer des valeurs selon P_r
  x <- sapply(u, F_inv)
  
  return(x)
}

# Paramètres
params <- list(alpha = 2.0, x0 = 50.0, eta = 4.0)

# Simulation des valeurs
n <- 100000  # Ajuster selon vos besoins
valeurs_simulees <- rremb(n, params)

# Afficher les résultats
print(head(valeurs_simulees, 10))

# Optionnel : Tracer l'histogramme des valeurs simulées
hist(valeurs_simulees, breaks = 50, main = "Histogramme des valeurs simulées", 
     xlab = "Valeurs", ylab = "Fréquence", col = "skyblue")

Alors j'ai utilisé des packages R pour essayer d'optimiser le code et j'arrive à 6 secondes pour 1 million de simulations

À vue de nez ton code passe l'essentiel de son temps à vérifier des conditions à cause des if et des else.

Tu peux probablement déjà gagner quelques millions en écrivant directement F sous forme explicite sans passer par F*. Ça évitera déjà de construire une call stack (pour te la faire simple une fonction qui en appelle une autre doit se souvenir d'où elle est pour pouvoir y revenir une fois que la fonction appelée est terminée).

Ensuite il faut te débarrasser des if/else en les remplaçant par exemple par des multiplications par des fonctions qui approchent tes indicatrices. Des trucs du genre exp(-1/x^2)  avec une pente très abrupte.

R n'est pas un langage très rapide, de base. C'est pratique et vectorialisé mais ça reste un langage de script qui fait appel à des librairies écrites en C++, comme Python.

Après en dernier recours il reste des techniques qui dépassent d'assez loin ton niveau actuel. Utiliser des suites pseudo aléatoires à faible discrépance, suites de Halton, quasi Monte Carlo, etc.

rremb_optimized <- function(n, params) {
  alpha <- params$alpha
  x0 <- params$x0
  eta <- params$eta
  
  # Vérification des paramètres une seule fois
  if (alpha <= 0 || x0 <= 0 || eta <= 3) {
    stop("Mauvais paramètres : alpha doit être > 0, x0 > 0 et eta > 3")
  }
  
  # Calcul du facteur de normalisation
  K <- 2 * alpha^(1 - eta) / (eta - 1)
  
  # Générer n valeurs uniformes
  u <- runif(n)

  # Approche vectorisée sans `if/else`
  z <- (eta - 1) * u*K
  x <- x0 + alpha * (1 - exp(-1 / z^2)) - 
       alpha * (1 - exp(-1 / (1 - z)^2)) 
  
  return(x)
}

# Paramètres
params <- list(alpha = 2.0, x0 = 50.0, eta = 4.0)

# Simulation optimisée
n <- 63000000
valeurs_simulees_optim <- rremb_optimized(n, params)

# Afficher les résultats
print(valeurs_simulees_optim)

# Résultat
hist(valeurs_simulees_optim, breaks = 50, main = "Histogramme des valeurs simulées (Optimisé)",
     xlab = "Valeurs", ylab = "Fréquence", col = "skyblue")

Ben non c'est bigrement incorrect. Ma suggestion d'exponentielle portait sur la forme de la fonction de régularisation de , ce n'était pas à prendre au mot près!

Je voulais dire qu'au lieu d'avoir un saut, tu aurais une pente raide mais continue du genre arctan ou expression de la tension dans un dipôle RC, ...

rremb <- function(n, params) {
  # Extraire les paramètres depuis la liste params
  alpha <- params$alpha
  x0 <- params$x0
  eta <- params$eta
  
  # Vérifier si les paramètres sont corrects
  if (alpha <= 0 | x0 <= 0 | eta <= 3) {
    stop("Mauvais paramètres : alpha doit être > 0, x0 > 0 et eta > 3")
  }
  
  # Fonction sigmoïde logistique pour la régularisation
  logistic <- function(x, y0, k = 100) {
    1 / (1 + exp(-k * (x - y0)))
  }
  
  # Fonction inverse modifiée avec transitions douces
  F_inv <- function(x) {
    # Terme pour x <= 0.5
    term1 <- alpha + x0 - (2 * x^(1 - eta))^(1 / (1 - eta))
    # Terme pour x > 0.5
    term2 <- x0 - alpha + (2 * alpha^(1 - eta) - 2 * x * alpha^(1 - eta))^(1 / (eta - 1))
    
    # Transition continue autour de 0.5
    transition <- logistic(x, 0.5)

    # Combinaison pondérée des deux termes
    return(transition * term2 + (1 - transition) * term1)
  }
  
  # Générer n valeurs uniformes sur (0, 1)
  u <- runif(n)
  
  # Utiliser F_inv pour générer des valeurs selon P_r
  x <- sapply(u, F_inv)
  
  return(x)
}

# Paramètres
params <- list(alpha = 2.0, x0 = 50.0, eta = 4.0)

# Simulation des valeurs
n <- 10000  # Ajuster selon vos besoins
valeurs_simulees <- rremb(n, params)

# Afficher les résultats
print(valeurs_simulees)

# Optionnel : Tracer l'histogramme des valeurs simulées
hist(valeurs_simulees, breaks = 50, main = "Histogramme des valeurs simulées", 
     xlab = "Valeurs", ylab = "Fréquence", col = "skyblue")

halton <- function(n, base) {
  sequence <- numeric(n)
  for (i in 1:n) {
    f <- 1
    value <- 0
    index <- i
    while (index > 0) {
      f <- f / base
      value <- value + f * (index %% base)
      index <- floor(index / base)
    }
    sequence[i] <- value
  }
  return(sequence)
}