Résolution d'une équation intégrale simple par échantillonnage aléatoire

8

Soit une fonction non négative. Je suis intéressé à trouver tel que La mise en garde : tout ce que je peux faire, c'est échantillonner aux points de . Je peux cependant choisir les emplacements où j'échantillonne au hasard, si je le souhaite. fz[0,1]

0zf(x)dx=1201f(x)dx
f[0,1]f

Des questions:

  1. Est-il possible d'obtenir une estimation non biaisée de après un nombre fini d'échantillons? Si oui, quelle est la plus petite variance possible d'une telle estimation après échantillons?zk
  2. Sinon, quelles procédures sont disponibles pour estimer et quels sont les temps de convergence associés.z

Comme l'a souligné Douglas Zare dans les commentaires, cela peut être très difficile à faire si la fonction est proche de zéro ou très grande. Heureusement, la fonction pour laquelle j'utilise ceci est délimitée par le haut et le bas, donc supposons que . De plus, nous pouvons également supposer que est Lipschitz ou même différenciable si cela aide.1f(x)2f

robinson
la source
1
Si vous n'avez pas plus d'informations, vous pouvez avoir un très mauvais comportement. Imaginez que est entre et , etDe légères modifications de feront passer la médiane de moins de à plus de . f01/32/301/3f(x) dx1/2.f1/32/3
Douglas Zare
@robinson Pourriez-vous fournir plus d'informations sur ? Ou êtes-vous intéressé à résoudre le problème pour n'importe quelle densité ? ff
@DouglasZare - Merci pour le commentaire; voir mon montage.
robinson
@Procrastinator - J'ai édité la question avec un peu plus d'informations.
robinson
3
(+1) Pour la mise à jour. En divisant le côté gauche par le droit, on peut voir que cela revient à trouver la médiane d'une distribution de probabilité inconnue supportée sur . [0,1]
Cardinal

Réponses:

6

Comme l'a souligné le cardinal dans son commentaire, votre question peut être reformulée comme suit.

Par algèbre simple, l'équation intégrale peut être réécrite sous la forme dans laquelle est la fonction de densité de probabilité définie comme

0zg(x)dx=12,
g
g(x)=f(x)01f(t)dt.

Soit une variable aléatoire de densité . Par définition, , donc votre équation intégrale est équivalente à ce qui signifie que votre problème peut être indiqué comme suit:XgP{Xz}=0zg(x)dx

P{Xz}=12,

"Soit une variable aléatoire de densité . Trouvez la médiane de "XgX

Pour estimer la médiane de , utilisez n'importe quelle méthode de simulation pour dessiner un échantillon de valeurs de et prendre comme estimation la médiane de l'échantillon.XX

Une possibilité consiste à utiliser l'algorithme de Metropolis-Hastings pour obtenir un échantillon de points avec la distribution souhaitée. En raison de l'expression de la probabilité d'acceptation dans l'algorithme de Metropolis-Hastings, nous n'avons pas besoin de connaître la valeur de la constante de normalisation de la densité . Donc, nous n'avons pas à faire cette intégration.01f(t)dtg

Le code ci-dessous utilise une forme particulièrement simple de l'algorithme Metropolis-Hastings connu sous le nom d'Independent Sampler, qui utilise une proposition dont la distribution ne dépend pas de la valeur actuelle de la chaîne. J'ai utilisé des propositions uniformes indépendantes. A titre de comparaison, le script génère le minimum de Monte Carlo et le résultat trouvé avec l'optimisation standard. Les points d'échantillonnage sont stockés dans le vecteur chain, mais nous éliminons les premiers points qui forment la période dite de «rodage» de la simulation.10000

BURN_IN = 10000
DRAWS   = 100000

f = function(x) exp(sin(x))

chain = numeric(BURN_IN + DRAWS)

x = 1/2

for (i in 1:(BURN_IN + DRAWS)) {
    y = runif(1) # proposal
    if (runif(1) < min(1, f(y)/f(x))) x = y
    chain[i] = x
}

x_min = median(chain[BURN_IN : (BURN_IN + DRAWS)])

cat("Metropolis minimum found at", x_min, "\n\n")

# MONTE CARLO ENDS HERE. The integrations bellow are just to check the results.

A = integrate(f, 0, 1)$value

F = function(x) (abs(integrate(f, 0, x)$value - A/2))

cat("Optimize minimum found at", optimize(F, c(0, 1))$minimum, "\n")

Voici quelques résultats:

Metropolis minimum found at 0.6005409 
Optimize minimum found at 0.601365

Ce code est conçu comme un point de départ pour ce dont vous avez vraiment besoin. Par conséquent, utilisez-le avec précaution.

Zen
la source
Merci pour votre réponse. Je ne connais pas R, donc je ne sais pas comment analyser ce que vous faites. Pourriez-vous énoncer en mots / formules votre procédure? Je vous remercie. En particulier, je me demande si vous respectez la contrainte selon laquelle la seule chose que vous pouvez faire est d'évaluer f - vous n'êtes pas autorisé, par exemple, à intégrer , (bien que vous puissiez certainement former des approximations de Monte-Carlo aux intégrales basées sur évaluations aléatoires). f
robinson
Oui, j'évalue juste pour obtenir l'estimation de Monte Carlo. f
Zen
Le code n'est qu'un exemple. La syntaxe R est similaire à d'autres langages. Une déclaration particulière que vous ne comprenez pas? Consultez la page Wikipedia sur l'algorithme de Metropolis-Hastings. Bien sûr, l'idée générale est plus importante. Vous pouvez échantillonner à partir du utilisant n'importe quelle méthode dont vous disposez. f/f
le
Avez-vous suivi un cours d'introduction aux processus stochastiques, couvrant les chaînes de Markov à temps discret?
Zen
1
BTW: Procrastinateurs du monde, unissez-vous! Mais pas aujourd'hui ...
Zen
3

La qualité de l'approximation intégrale, au moins dans le cas aussi simple que 1D, est donnée par (Théorème 2.10 dans Niederreiter (1992) ): où est le module de continuité de la fonction (lié à la variation totale, et facilement exprimable pour les fonctions de Lipshitz), et est l'écart (extrême), ou la différence maximale entre la fraction de hits par la séquence

|1Nn=1Nf(xn)01f(u)du|ω(f;DN(x1,,xN))
ω(f;t)=sup{|f(u)f(v)|:u,v[0,1],|uv|t,t>0}
DN(x1,,xN)=supu|1Nn1{xn[0,u)}u|=12N+maxn|xn2n12N|
x1,,xNd'un intervalle semi-ouvert et sa mesure de Lebesgue . La première expression est la définition et la deuxième expression est la propriété des séquences 1D dans (Théorème 2.6 dans le même livre).[0,u)u[0,1]

Donc, évidemment, pour minimiser l'erreur dans l'approximation intégrale, au moins dans la RHS de votre équation, vous devez prendre . Vissez les évaluations aléatoires, elles courent le risque d'avoir un écart aléatoire sur une caractéristique importante de la fonction.xn=(2n1)/2N

Un gros inconvénient de cette approche est que vous devez vous engager sur une valeur de pour produire cette séquence uniformément distribuée. Si vous n'êtes pas satisfait de la qualité de l'approximation qu'il fournit, tout ce que vous pouvez faire est de doubler la valeur de et d'atteindre tous les points médians des intervalles créés précédemment.NN

Si vous voulez avoir une solution où vous pouvez augmenter le nombre de points plus progressivement, vous pouvez continuer à lire ce livre et en apprendre davantage sur les séquences de van der Corput et les inverses radicales. Voir Séquences à faible écart sur Wikipédia, il fournit tous les détails.

Mise à jour: pour résoudre , définissez la somme partielle Trouver tel que et interpoler pour trouver Cette interpolation suppose que est continue. Si en plus est deux fois différentiable, alors cette approximation en intégrant l'expansion du second ordre pour incorporer et , et en résolvant une équation cubique pour .z

Sk=1Nn=1kf(2n12N).
k
Sk12SN<Sk+1,
zN=2k12N+SN/2SkN(Sk+1Sk).
f()f()Sk1Sk+2z
StasK
la source
1
J'aime l'essentiel de cela. Je pense qu'il serait utile de rendre plus explicite la stratégie que vous proposez pour résoudre la question du PO. À l'heure actuelle, la réponse se lit (pour moi) principalement comme si elle traitait de la façon de calculer l'ERS de l'équation dans la question.
Cardinal
1
(+1) Belle mise à jour. peut simplement être considéré comme une approximation de somme de Riemann à l'intégrale où nous utilisons la valeur de au milieu de chaque intervalle défini par la partition, plutôt que le point final gauche ou droit. :-)SNf
Cardinal
Oui; il est intéressant cependant que cette somme de Riemann ait cette justification d'optimalité.
StasK