Resolver una ecuación integral simple mediante muestreo aleatorio

Dejar $f$ser una función no negativa Estoy interesado en encontrar tal que La advertencia : todo lo que puedo hacer es muestrear en puntos en . Sin embargo, puedo elegir los lugares donde muestreo al azar, si así lo deseo. $z \in [0,1]$

$$\int_0^{z} f(x)\,dx = \frac{1}{2}\int_0^1 f(x)\,dx$$ $f$$[0,1]$$f$

Preguntas:

1. ¿Es posible obtener una estimación imparcial de después de muchas muestras? Si es así, ¿cuál es la varianza más pequeña posible de tal estimación después de muestras?$z$$k$
2. Si no, qué procedimientos están disponibles para estimar y cuáles son sus tiempos de convergencia asociados.$z$

Como señaló Douglas Zare en los comentarios, esto puede ser muy difícil de hacer si la función es cercana a cero o muy grande. Afortunadamente, la función para la que necesito usar esto está limitada desde arriba y desde abajo, así que supongamos que . Además, también podemos suponer que es Lipschitz o incluso diferenciable si eso ayuda.$1 \leq f(x) \leq 2$$f$

Como lo señaló el cardenal en su comentario, su pregunta puede reformularse como sigue.

Por álgebra simple, la ecuación integral se puede reescribir como

$$\int_0^z g(x)\,dx = \frac{1}{2} \, ,$$ en el cual $g$ es la función de densidad de probabilidad definida como $$g(x)=\frac{f(x)}{\int_0^1 f(t)\,dt} \, .$$

Dejar $X$ ser una variable aleatoria con densidad $g$. Por definición,$P\{X\leq z\}=\int_0^z g(x)\,dx$, entonces tu ecuación integral es equivalente a

$$P\{X\leq z\}=\frac{1}{2} \, ,$$ lo que significa que su problema puede expresarse como:

"Dejar $X$ ser una variable aleatoria con densidad $g$. Encuentra la mediana de$X$".

Para estimar la mediana de $X$, use cualquier método de simulación para dibujar una muestra de valores de $X$ y tome como estimación la mediana de la muestra.

Una posibilidad es utilizar el algoritmo Metropolis-Hastings para obtener una muestra de puntos con la distribución deseada. Debido a la expresión de la probabilidad de aceptación en el algoritmo de Metropolis-Hastings, no necesitamos saber el valor de la constante de normalización$\int_0^1 f(t)\,dt$ de la densidad $g$. Entonces, no tenemos que hacer esta integración.

El siguiente código usa una forma particularmente simple del algoritmo Metropolis-Hastings conocido como Indepence Sampler, que usa una propuesta cuya distribución no depende del valor actual de la cadena. He usado propuestas uniformes independientes. A modo de comparación, el script genera el mínimo de Monte Carlo y el resultado encontrado con la optimización estándar. Los puntos de muestra se almacenan en el vector chain, pero descartamos el primer$10000$ puntos que forman el llamado período de "quemado" de la simulación.

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")

Aquí hay algunos resultados:

Metropolis minimum found at 0.6005409 
Optimize minimum found at 0.601365

Este código es solo un punto de partida para lo que realmente necesita. Por lo tanto, use con cuidado.

La calidad de la aproximación integral, al menos en el caso tan simple como 1D, viene dada por (Teorema 2.10 en Niederreiter (1992) ):

$$\Bigl|\frac 1N \sum_{n=1}^N f(x_n) - \int_0^1 f(u) \, {\rm d}u \Bigr| \le \omega (f; D_N^*(x_1, \ldots, x_N) )$$ dónde $$\omega(f;t) = \sup \{ |f(u)-f(v)| : u, v \in [0,1], |u-v|\le t , t>0\}$$ es el módulo de continuidad de la función (relacionado con la variación total, y fácilmente expresable para las funciones de Lipshitz), y $$D_N^*(x_1,\ldots,x_N) = \sup_u \Bigl| \frac1N \sum_n 1\bigl\{ x_n \in [0,u) \bigr\} - u \Bigr| = \frac1{2N} + \max_n \Bigl|x_n - \frac{2n-1}{2N}\Bigr|$$ es la discrepancia (extrema), o la diferencia máxima entre la fracción de golpes por la secuencia $x_1, \ldots, x_N$ de un intervalo semiabierto $[0,u)$ y su medida de Lebesgue $u$. La primera expresión es la definición, y la segunda expresión es propiedad de las secuencias 1D en$[0,1]$ (Teorema 2.6 en el mismo libro).

Obviamente, para minimizar el error en la aproximación integral, al menos en el RHS de su ecuación, debe tomar $x_n = (2n-1)/2N$. Atornille las evaluaciones aleatorias, corren el riesgo de tener una brecha aleatoria en una característica importante de la función.

Una gran desventaja de este enfoque es que debe comprometerse con un valor de $N$para producir esta secuencia uniformemente distribuida. Si no está satisfecho con la calidad de aproximación que proporciona, todo lo que puede hacer es duplicar el valor de$N$ y presione todos los puntos medios de los intervalos creados previamente.

Si desea tener una solución donde pueda aumentar el número de puntos más gradualmente, puede continuar leyendo ese libro y aprender sobre las secuencias de van der Corput e inversas radicales. Ver secuencias de baja discrepancia en Wikipedia, proporciona todos los detalles.

Actualización: para resolver por$z$, define la suma parcial

$$S_k = \frac1N \sum_{n=1}^k f\Bigl( \frac{2n-1}{2N} \Bigr).$$ Encontrar $k$ tal que $$S_k \le \frac12 S_N < S_{k+1},$$ e interpolar para encontrar $$z_N = \frac{2k-1}{2N} + \frac{S_N/2 - S_k}{N(S_{k+1}-S_k)}.$$ Esta interpolación supone que $f(\cdot)$es continuo Si adicionalmente$f(\cdot)$ es dos veces diferenciable, entonces esta aproximación integrando la expansión de segundo orden para incorporar $S_{k-1}$ y $S_{k+2}$y resolviendo una ecuación cúbica para $z$.