spearman - statistics with python pdf

Trazar intervalos de confianza para la estimación de máxima verosimilitud (3)

Intento escribir código para generar intervalos de confianza para la cantidad de libros diferentes en una biblioteca (y también para producir una trama informativa).

Mi primo está en la escuela primaria y cada semana recibe un libro de su maestro. Luego lo lee y lo devuelve a tiempo para obtener otro la próxima semana. Después de un tiempo comenzamos a notar que estaba obteniendo libros que había leído antes y esto se volvió gradualmente más común con el tiempo.

Supongamos que el número verdadero de libros en la biblioteca es N y el maestro elige uno al azar (con reemplazo) para darle cada semana. Si en la semana t el número de ocasiones en las que ha recibido un libro que ha leído es x, entonces puedo generar un cálculo de probabilidad máxima para el número de libros en la biblioteca, siguiendo https://math.stackexchange.com/questions/ 615464 / how-many-books-are-in-a-library .

Ejemplo: considere una biblioteca con cinco libros A, B, C, D y E. Si recibe libros [A, B, A, B, B, D] en siete semanas sucesivas, entonces el valor de x (el número de duplicados) será [0, 0, 1, 1, 2, 3, 3] después de cada una de esas semanas, lo que significa que después de siete semanas, habrá recibido un libro que ya ha leído en tres ocasiones.

Para visualizar la función de verosimilitud (asumiendo que he entendido cuál es la correcta) he escrito el siguiente código, que creo traza la función de verosimilitud. El máximo es alrededor de 135, que es de hecho la estimación de máxima verosimilitud de acuerdo con el enlace MSE anterior.

from __future__ import division import random import matplotlib.pyplot as plt import numpy as np #N is the true number of books. t is the number of weeks.unk is the true number of repeats found t = 30 unk = 3 def numberrepeats(N, t): return t - len(set([random.randint(0,N) for i in xrange(t)])) iters = 1000 ydata = [] for N in xrange(10,500): sampledunk = [numberrepeats(N,t) for i in xrange(iters)].count(unk) ydata.append(sampledunk/iters) print "MLE is", np.argmax(ydata) xdata = range(10, 500) print len(xdata), len(ydata) plt.plot(xdata,ydata) plt.show()

La salida se ve como

Mis preguntas son estas:

¿Hay una manera fácil de obtener un intervalo de confianza del 95% y trazarlo en el diagrama?
¿Cómo se puede superponer una curva suavizada sobre la trama?
¿Hay alguna forma mejor en que se debería haber escrito mi código? No es muy elegante y también es bastante lento.

Encontrar el intervalo de confianza del 95% significa encontrar el rango del eje x de modo que el 95% de las veces la estimación de máxima verosimilitud empírica que obtenemos mediante muestreo (que debería ser teóricamente 135 en este ejemplo) caerá dentro de ella. La respuesta que @mbatchkarov ha dado no hace esto correctamente actualmente.

Ahora hay una respuesta matemática en https://math.stackexchange.com/questions/656101/how-to-find-a-confidence-interval-for-a-maximum-likelihood-estimate .

Aquí hay una respuesta a su primera pregunta y un puntero a una solución para el segundo:

plot(xdata,ydata) # calculate the cumulative distribution function cdf = np.cumsum(ydata)/sum(ydata) # get the left and right boundary of the interval that contains 95% of the probability mass right=argmax(cdf>0.975) left=argmax(cdf>0.025) # indicate confidence interval with vertical lines vlines(xdata[left], 0, ydata[left]) vlines(xdata[right], 0, ydata[right]) # hatch confidence interval fill_between(xdata[left:right], ydata[left:right], facecolor=''blue'', alpha=0.5)

Esto produce la siguiente figura:

Trataré de responder la pregunta 3 cuando tenga más tiempo :)

La manera simple (numérica) de obtener un intervalo de confianza es simplemente ejecutar el script muchas veces, y ver cuánto varía su estimación. Puede usar esa desviación estándar para calcular el intervalo de confianza.

En aras del tiempo, otra opción es ejecutar una serie de ensayos con cada valor de N (utilicé 2000), y luego usar submuestreo aleatorio de esos ensayos para obtener una estimación de la desviación estándar del estimador. Básicamente, esto implica seleccionar un subconjunto de las pruebas, generar su curva de probabilidad usando ese subconjunto, y luego encontrar el máximo de esa curva para obtener su estimador. Hace esto en muchos subconjuntos y esto le da un conjunto de estimadores, que puede usar para encontrar un intervalo de confianza en su estimador. Mi guion completo es el siguiente:

import numpy as np t = 30 k = 3 def trial(N): return t - len(np.unique(np.random.randint(0, N, size=t))) def trials(N, n_trials): return np.asarray([trial(N) for i in xrange(n_trials)]) n_trials = 2000 Ns = np.arange(1, 501) results = np.asarray([trials(N, n_trials=n_trials) for N in Ns]) def likelihood(results): L = (results == 3).mean(-1) # boxcar filtering n = 10 L = np.convolve(L, np.ones(n) / float(n), mode=''same'') return L def max_likelihood_estimate(Ns, results): i = np.argmax(likelihood(results)) return Ns[i] def max_likelihood(Ns, results): # calculate mean from all trials mean = max_likelihood_estimate(Ns, results) # randomly subsample results to estimate std n_samples = 100 sample_frac = 0.25 estimates = np.zeros(n_samples) for i in xrange(n_samples): mask = np.random.uniform(size=results.shape[1]) < sample_frac estimates[i] = max_likelihood_estimate(Ns, results[:,mask]) std = estimates.std() sterr = std * np.sqrt(sample_frac) # is this mathematically sound? ci = (mean - 1.96*sterr, mean + 1.96*sterr) return mean, std, sterr, ci mean, std, sterr, ci = max_likelihood(Ns, results) print "Max likelihood estimate: ", mean print "Max likelihood 95% ci: ", ci

Hay dos inconvenientes para este método. Una es que, dado que toma muchas submuestras del mismo conjunto de pruebas, sus estimaciones no son independientes. Para limitar el efecto de esto, solo usé el 25% de los resultados para cada subconjunto. Otro inconveniente es que cada submuestra es solo una fracción de sus datos, por lo que las estimaciones derivadas de estos subconjuntos tendrán más variaciones que las estimaciones derivadas de ejecutar el script completo muchas veces. Para dar cuenta de esto, calculé el error estándar como la desviación estándar dividida por la raíz cuadrada de 4, ya que tenía cuatro veces más datos en mi conjunto de datos completo que en una de las submuestras. Sin embargo, no estoy lo suficientemente familiarizado con la teoría de Monte Carlo para saber si esto es matemáticamente sólido. Ejecutar mi script varias veces parecía indicar que mis resultados eran razonables.

Por último, utilicé un filtro de caja en las curvas de probabilidad para suavizar un poco. Idealmente, esto debería mejorar los resultados, pero incluso con el filtrado todavía había una considerable cantidad de variabilidad en los resultados. Al calcular el valor para el estimador general, no estaba seguro de si sería mejor calcular una curva de probabilidad de todos los resultados y usar el máximo de eso (esto es lo que terminé haciendo), o usar la media de todos los estimadores de subconjuntos. Usar la media de los estimadores de subconjuntos podría ayudar a cancelar parte de la rugosidad en las curvas que quedan después del filtrado, pero no estoy seguro de esto.

Parece que estás bien en la primera parte, así que abordaré tu segundo y tercer punto.

Hay muchas formas de ajustar curvas suaves, con scipy.interpolate y splines, o con scipy.optimize.curve_fit . Personalmente, prefiero curve_fit , porque puedes suministrar tu propia función y dejar que se ajuste a tus parámetros.

Alternativamente, si no desea aprender una función paramétrica, puede hacer un suavizado de ventana móvil simple con numpy.convolve .

En cuanto a la calidad del código: no estás aprovechando la velocidad de numpy, porque estás haciendo cosas en python puro. Escribiría tu código (existente) de esta manera:

from __future__ import division import numpy as np import matplotlib.pyplot as plt # N is the true number of books. # t is the number of weeks. # unk is the true number of repeats found t = 30 unk = 3 def numberrepeats(N, t, iters): rand = np.random.randint(0, N, size=(t, iters)) return t - np.array([len(set(r)) for r in rand]) iters = 1000 ydata = np.empty(500-10) for N in xrange(10,500): sampledunk = np.count_nonzero(numberrepeats(N,t,iters) == unk) ydata[N-10] = sampledunk/iters print "MLE is", np.argmax(ydata) xdata = range(10, 500) print len(xdata), len(ydata) plt.plot(xdata,ydata) plt.show()

Probablemente sea posible optimizar esto aún más, pero este cambio hace que el tiempo de ejecución de su código sea de ~ 30 segundos a ~ 2 segundos en mi máquina.