Teoría de la probabilidad para la representación físicamente precisa

Introduccion

En el renderizado, a menudo se usa el cálculo de integrales definidas multidimensionales: por ejemplo, para determinar la visibilidad de las fuentes de luz espacial (área de luz), la luminosidad que llega a la región de píxeles, la luminosidad que llega durante un período de tiempo y la irradiación que atraviesa el hemisferio de un punto de superficie. El cálculo de estas integrales generalmente se lleva a cabo utilizando la integración de Monte Carlo, en la cual la integral se reemplaza por la expectativa de un experimento estocástico.

En este artículo hablaré en detalle sobre el proceso básico de integración de Monte Carlo, así como sobre varias técnicas para reducir la variación de la técnica. Esto se hará desde un punto de vista práctico: se supone que el lector no está muy familiarizado con la teoría de la probabilidad, pero que aún desea desarrollar algoritmos de representación efectivos y correctos.

Integrales definidas

Una integral definida es una integral de la forma $$$\ int_ {a} ^ {b} f (x) dx$ donde $$$[a, b]$ Es un segmento (o región), $$$x$ - escalar, y $$$f (x)$ - una función que se puede calcular para cualquier punto del segmento. Como está escrito en Wikipedia , cierta integral es un área con un signo en un plano $$$x$ limitado por horario $$$f$ eje $$$x$ y líneas verticales $$$x = a$ y $$$x = b$ ( Figura 1a ).

Este concepto se extiende lógicamente a un mayor número de dimensiones: para una determinada integral doble, el área con un signo se convierte en un volumen con un signo ( Figura 1b ), y en general para ciertas integrales múltiples se convierte en un volumen multidimensional con un signo .


Figura 1: ejemplos de ciertas integrales.

En algunos casos, el área se puede determinar analíticamente , por ejemplo, para $$$f (x) = 2$ : en el segmento $$$[a, b]$ el área será igual $$$2 (b-a)$ . En otros casos, una solución analítica es imposible, por ejemplo, cuando necesitamos averiguar el volumen de la parte del iceberg sobre el agua ( Figura 1c ). En tales casos $$$f (x)$ a menudo se puede determinar por muestreo .

Integración numérica

Podemos calcular aproximadamente el área de integrales complejas utilizando la integración numérica . Un ejemplo es la suma de Riemann . Esta cantidad se calcula dividiendo el área en formas regulares (por ejemplo, rectángulos), que juntas forman un área similar a un área real. La suma de Riemann se define de la siguiente manera:

$$

$$\ tag {1} S = \ sum_ {i = 1} ^ {n} f (x_i) \ Delta x_i$$

$$$n$ Es el número de subintervalos y $$$\ Delta x_i = \ frac {b-a} {n}$ - el ancho de un subintervalo. Para cada intervalo $$$i$ nosotros tomamos muestras $$$f$ en un punto $$$x_i$ dentro del subintervalo (en la Figura 2, este punto se encuentra al comienzo del subintervalo).


Figura 2: suma de Riemann.

Vale la pena señalar que con el aumento $$$n$ la suma de Riemann converge al valor real de la integral:

$$

$$\ tag {2} \ int_ {a} ^ {b} f (x) dx = \ lim_ {|| \ Delta x || \ to 0} \ sum_ {i = 1} ^ {n} f (x_i) \ Delta x_i$$

La suma de Riemann también se puede utilizar para grandes dimensiones ( Figura 3 ). Sin embargo, aquí nos enfrentamos a un problema: para una función con dos parámetros, el número de subintervalos debería ser mucho mayor si queremos lograr una resolución comparable a la utilizada en el caso bidimensional. Este fenómeno se llama la maldición de las dimensiones , y en las dimensiones superiores se exacerba.


Figura 3: suma de Riemann para una integral doble.

Ahora evaluaremos la precisión de la suma de Riemann para la siguiente función (elegimos intencionalmente una función compleja):

$$

$$\ tag {3} f (x) = \ left | \ sin \ left (\ frac {1} {2} x + \ frac {\ pi} {2} \ right) \ tan \ frac {x} {27} + \ sin \ left (\ frac {3} {5} x ^ 2 \ right) + \ frac {4} {x + \ pi + 1} -1 \ right |$$

Gráfico de función en un segmento $$$[- 2.5,2.5]$ se muestra a continuación. Como referencia, calculamos una cierta integral en Wolfram Alpha $$$\ int _ {- 2.5} ^ {2.5} f (x)$ área de obtención $$$3.12970$ . El gráfico de la derecha muestra la precisión de la integración numérica utilizando la suma de Riemann para aumentar $$$n$ .


Figura 4: Gráfico de funciones y precisión de suma de Riemann. Incluso con pequeños $$$n$ Obtenemos un resultado bastante preciso.

Para tener una idea de la precisión, damos los números: para $$$n = 50$ el error es $$$~ 2 \ veces10 ^ {- 3}$ . En $$$n = 100$ el error es $$$~ 3 \ veces10 ^ {- 4}$ . El siguiente orden de magnitud se obtiene con $$$n = 200$ .

Para obtener más información sobre las cantidades de Riemann, consulte los siguientes recursos:

• Academia Kahn: https://www.khanacademy.org/math/ap-calculus-ab/ab-accumulation-riemann-sums/ab-riemann-sums/v/simple-riemann-approximation-using-rectangles
• Wikipedia: https://en.wikipedia.org/wiki/Riemann_sum

Montecarlo (1)

Al renderizar, casi ninguna (¿y tal vez ninguna?) Las integrales son individuales . Esto significa que nos encontraremos rápidamente con la maldición de las dimensiones. Además, el muestreo de una función a intervalos iguales está sujeto a un muestreo y una distorsión insuficientes : podemos omitir valores importantes de la función u obtener interferencia mutua involuntaria entre la función muestreada y el patrón de muestreo ( Figura 5 ).


Figura 5: las distorsiones conducen a la pérdida de partes de la función muestreada (rojo) y, en este caso, a una interpretación completamente incorrecta de la función.

Estos problemas se resuelven utilizando una técnica llamada integración de Monte Carlo . Similar a la suma de Riemann, también usa el muestreo de funciones en un conjunto de puntos, pero a diferencia del patrón de suma determinista de Riemann, usamos un ingrediente fundamentalmente no determinista : números aleatorios.

La integración de Monte Carlo se basa en la siguiente observación: la integral puede ser reemplazada por la expectativa de un experimento estocástico:

$$

$$\ tag {4} \ int_ {a} ^ {b} f (x) dx = (ba) E \ left [f (X) \ right] \ approx \ frac {ba} {n} \ sum_ {i = 1} ^ {n} f (X)$$

En otras palabras, muestreamos la función $$$n$ veces en puntos aleatorios dentro de un segmento (denotado por una mayúscula $$$X$ ), promedie las muestras y multiplique por el ancho del segmento (para una función unidimensional). Como en el caso de la suma de Riemann, cuando $$$n$ hasta el infinito, el valor promedio de las muestras converge a la expectativa, es decir, al valor verdadero de la integral.

Un poco de teoría de la probabilidad.

Es importante comprender cada uno de los conceptos utilizados aquí. Comencemos con la espera : este es el valor esperado para una sola muestra. Tenga en cuenta que este no es necesariamente un valor posible , que puede parecer contradictorio. Por ejemplo, cuando tiramos el dado, la expectativa es igual a $$$3.5$ - el promedio de todos los resultados posibles: $$$(1 + 2 + 3 + 4 + 5 + 6) /6=21/6=3.5$ .

El segundo concepto son los números aleatorios . Esto puede parecer obvio, pero para la integración de Monte Carlo necesitamos números aleatorios distribuidos uniformemente, es decir. cada valor debe tener la misma probabilidad de generación. Hablaremos más sobre esto más tarde.

El tercer concepto es la desviación y la varianza asociada a ella. Incluso cuando tomamos una pequeña cantidad de números, el valor promedio esperado, así como la expectativa de cada muestra individual, deberían ser los mismos. Sin embargo, cuando calculamos la ecuación 4, rara vez obtenemos ese valor. La desviación es la diferencia entre la expectativa y el resultado del experimento: $$$X-E (X)$ .

En la práctica, esta desviación tiene una distribución interesante:


Este es un gráfico de la distribución normal , o distribución gaussiana : muestra que no todas las desviaciones son igualmente probables. De hecho, aproximadamente el 68,2% de las muestras están en el rango $$$-1 \ sigma..1 \ sigma$ donde $$$\ sigma$ (sigma) es la desviación estándar . La desviación estándar se puede describir de dos maneras:

• La desviación estándar es una medida de la variabilidad de los datos.
• 95% de los puntos de datos están dentro $$$2 \ sigma$ del promedio

Hay dos métodos para determinar la desviación estándar:

1. Desviacion estandar $$$\ sigma = \ sqrt {\ frac {1} {n} \ sum_ {i = 1} ^ {n} \ left (X_i-E \ left [X \ right] \ right) ^ 2}$ : se puede calcular si hay una distribución de probabilidad discreta y se conoce la expectativa $$$E [X]$ . Esto es cierto para cubos en los que $$$X = {1,2,3,4,5,6}$ y $$$E [X] = 3.5$ . Sustituyendo los números, obtenemos $$$\ sigma = 1.71$ .
2. Además, la desviación estándar de las muestras se puede calcular como $$$\ sigma = \ sqrt {\ frac {1} {n-1} \ sum_ {i = 1} ^ {n} \ left (X_i-X \ right) ^ 2}$ . Lea más sobre esto en Wikipedia .

Verificación: ¿es correcto? Si $$$\ sigma = 1.71$ , declaramos que el 68.2% de las muestras están dentro de 1.71 de 3.5. Sabemos que $$${2,3,4,5}$ cumplir este criterio, y $$$1$ y $$$6$ no Cuatro de cada seis son 66.7%. Si nuestro cubo pudiera producir algún valor en el intervalo $$$[1..6]$ , entonces obtendríamos exactamente 68.2%.

En lugar de la desviación estándar, el concepto asociado de varianza , que se define como $$$Var \ left [X \ right] = \ sigma ^ 2$ . Como se usa el cuadrado, la varianza siempre es positiva, lo que ayuda en los cálculos.

Montecarlo (2)

Arriba, calculamos aproximadamente la ecuación 3 usando la suma de Riemann. Ahora repetimos este experimento con la integración de Monte Carlo. Recuerde que la integración de Monte Carlo se define de la siguiente manera:

$$

$$\ tag {5} \ int_ {a} ^ {b} f (x) dx = (ba) E \ left [f (X) \ right] \ approx \ frac {ba} {n} \ sum_ {i = 1} ^ {n} f (X)$$

Traduzcamos esto al código C:

double sum = 0; for( int i = 0; i < n; i++ ) sum += f( Rand( 5 ) - 2.5 ); sum = (sum * 5.0) / (double)n; 

Resultado para valores de $$$n = 2$ antes $$$n = 200$ se muestra en el cuadro a continuación. De esto se puede suponer que la integración de Monte Carlo se manifiesta mucho peor que la suma de Riemann. Un examen más detallado del error dice que con $$$n = 200$ El error promedio de la suma de Riemann es $$$0.0002$ y Monte Carlo $$$0.13$ .


Figura 6: error de Monte Carlo en 2..200 muestras.

En dimensiones superiores, esta diferencia se reduce, pero no se elimina por completo. La ecuación que se muestra a continuación es una versión ampliada de la utilizada anteriormente, que recibe dos parámetros:

$$

$$f (x, y) = \ left | \ sin \ left (\ frac {1} {2} x + \ frac {\ pi} {2} \ right) \ tan \ frac {x} {27} + \ sin \ left (\ frac {1} {6} x ^ 2 \ right) + \ frac {4} {x + \ pi + 1} -1 \ right | \ left | \ sin \ left (1.1y \ right) \ cos \ left (2.3x \ right) \ right | (6)$$

Figura 7: Gráfico de la ecuación anterior.

En el campo de la definición $$$x∈ [-2.5,2.5], y∈ [-2.5,2.5]$ volumen limitado por esta función y plano $$$xy$ es igual $$$6.8685$ . En $$$n = 400$ (20 × 20 muestras) el error de la suma de Riemann es $$$0.043$ . Con el mismo número de muestras, el error promedio de integración de Monte Carlo es $$$0.33$ . Esto es mejor que el resultado anterior, pero la diferencia sigue siendo significativa. Para comprender este problema, estudiaremos la conocida técnica de reducción de dispersión de integración Monte Carlo llamada "estratificación".


Figura 8: impacto de la estratificación; a) muestras con mala distribución; b) muestras con distribución uniforme.

La estratificación aumenta la uniformidad de los números aleatorios. En la Figura 8a , se usan ocho números aleatorios para muestrear la función. Dado que cada número se selecciona al azar, a menudo se distribuyen de manera desigual sobre el dominio de definición. La Figura 8b muestra el efecto de la estratificación: el área de definición se divide en ocho estratos, y se selecciona una posición aleatoria en cada estrato, lo que mejora la uniformidad.

El efecto sobre la varianza es bastante obvio. La figura 9a muestra un gráfico de los resultados con y sin estratificación. La figura 9b muestra el error de valor aproximado. En $$$n = 10$ el error promedio para 8 estratos es $$$0.05$ ; para 20 estratos $$$0.07$ , y para 200 estratos disminuye a $$$0.002$ . Según estos resultados, parece que vale la pena usar una gran cantidad de estratos. Sin embargo, la estratificación tiene desventajas que aumentan con el aumento del número de estratos. En primer lugar, el número de muestras siempre debe ser un múltiplo del número de estratos; en segundo lugar, como en la suma de Riemann, la estratificación sufre la maldición de las dimensiones.


Figura 9: estratificación y varianza: a) un valor aproximado para el número de muestras de n = 2 a n = 200; b) desviación.

Muestra de importancia

En las secciones anteriores, muestreamos las ecuaciones de manera uniforme. La extensión de la función integradora de Monte Carlo nos permite cambiar la situación:

$$

$$\ tag {7} \ int_ {a} ^ {b} f (x) dx = (ba) E \ left [f (X) \ right] \ approx \ frac {ba} {n} \ sum_ {i = 1} ^ {n} \ frac {f (X)} {p (X)}$$

Aqui $$$p (X)$ Es una función de densidad de probabilidad (pdf) : determina la probabilidad relativa de que una variable aleatoria tome un cierto valor.

Para una variable aleatoria uniforme en el intervalo $$$0..1$ , pdf es 1 ( Figura 10 a), y esto significa que cada valor tiene la misma probabilidad de elección. Si integramos esta función sobre $$$[0,0.5]$ entonces obtenemos la probabilidad en $$$0.5$ de que $$$X <\ frac {1} {2}$ . Para $$$X> \ frac {1} {2}$ Obviamente tenemos la misma probabilidad.


Figura 10: Distribuciones de probabilidad. a) pdf constante en el que cada muestra tiene la misma probabilidad de elección; b) pdf, donde las muestras por debajo de 0.5 tienen una mayor probabilidad de selección.

La figura 10b muestra otro pdf. En este caso, la probabilidad de generar un número es menor $$$\ frac {1} {2}$ igual al 70%. Esto se puede implementar utilizando el siguiente fragmento de código:

 float SamplePdf() { if (Rand() < 0.7f) return Rand( 0.5f ); else return Rand( 0.5f ) + 0.5f; } 

Este pdf se define de la siguiente manera:

$$

$$\ tag {8} p (x) = \ left \ {\ begin {matrix} 1.4, si x <\ frac {1} {2} \\ 0.6, de lo contrario \ end {matrix} \ right.$$

Los numeros $$$1.4$ y $$$0.6$ reflejar la necesidad de que la probabilidad $$$x <\ frac {1} {2}$ fue igual al 70%. Al integrar pdf por $$$[0 .. \ frac {1} {2}]$ da $$$1.4 \ veces \ frac {1} {2}$ y $$$0.6 \ veces \ frac {1} {2}$ es igual $$$0.3$ . Esto ilustra un requisito importante para todos los archivos PDF en general: el resultado de la integración de PDF debería ser 1. Otro requisito es que $$$p (x)$ no puede ser cero si $$$f (x)$ distinto de cero: significaría que las partes $$$f$ tener una probabilidad de muestreo cero, lo que obviamente afecta el valor.

Algunos consejos para entender el concepto de pdf:

• Un valor de pdf no describe la probabilidad: por lo tanto, el pdf local puede ser mayor que 1 (por ejemplo, como en el pdf que se acaba de examinar).
• Sin embargo, la integral sobre el dominio de definición de pdf es una probabilidad, lo que significa que la integración de pdf da 1.

Un valor puede interpretarse como la posibilidad relativa de la aparición de un valor específico.

Vale la pena considerar que la distribución normal es una función de distribución de probabilidad: nos da la probabilidad de que alguna variable aleatoria esté en un cierto intervalo. En el caso de una distribución normal, esta variable aleatoria es una desviación de la media. Como cualquier pdf decente, el resultado de integrar la distribución normal es 1.

Por lo tanto, la ecuación 7 nos permite realizar un muestreo no uniforme. Lo compensa dividiendo cada muestra por la probabilidad relativa de su elección. La importancia de esto se muestra en la Figura 11a . El gráfico de funciones muestra un intervalo significativo en el que su valor es $$$0$ . El muestreo en esta área no tiene sentido: no se agrega nada a la suma, simplemente dividimos por un número mayor. Recuerde el iceberg de la Figura 1c : no tiene sentido muestrear la altura en un área grande alrededor del iceberg.


Figura 11: pdf para una función con valores cero.

En la Figura 11b se muestra un pdf que utiliza este conocimiento de la función. Tenga en cuenta que este pdf es en realidad cero para el rango de valores. Esto no lo convierte en un pdf incorrecto: en algunos lugares, la función es cero. Podemos extender esta idea más allá de cero. Las muestras se gastan mejor en aquellos lugares en los que la función tiene valores significativos. De hecho, el pdf ideal es proporcional a la función muestreada . En la Figura 12a se muestra un muy buen pdf para nuestra función. Un pdf aún mejor se muestra en la Figura 12b . En ambos casos, no debemos olvidar normalizarlo para que la integral sea igual a 1.


Figura 12: PDF mejorado para la función en la Figura 11.

El pdf en la Figura 12 plantea dos tareas para nosotros:

1. cómo crear tal pdf;
2. ¿Cómo probar tal pdf?

La respuesta a ambas preguntas es la misma: no necesitamos hacer esto. En muchos casos, la función que queremos integrar es desconocida, y la única forma de determinar los lugares donde es importante es muestrearla, y para esto necesitamos pdf; Situación clásica de "pollo y huevos".

Sin embargo, en otros casos, tenemos una idea aproximada de dónde la función puede dar valores más altos o valores cero. En tales casos, un pdf muy aproximado a menudo es mejor que ningún pdf.

Además, podemos tener la oportunidad de crear pdf sobre la marcha. Un par de muestras dan una idea de la forma de la función, y en base a ella apuntamos las muestras posteriores a aquellos lugares donde esperamos valores altos, que usamos para mejorar el pdf, y así sucesivamente.

En el siguiente artículo, aplicaremos estos conceptos a la implementación de renderizado. Un desafío serio es construir pdf. Examinamos varios casos en los que los archivos PDF ayudan en el muestreo.