Una cosa a tener en cuenta es que asintóticamente, los MLE van a ser normales (bajo algunas condiciones). Suponiendo que esas condiciones se apliquen, entonces si N/M es realmente grande y las observaciones se han dado aleatoriamente a las computadoras (para que tengan funciones de probabilidad igualmente inciertas) esas probabilidades deberían ser casi cuadráticas cerca de sus máximos (con una segunda derivada casi constante a través de los subconjuntos), y los estimadores deberían ser casi normales con varianza constante, y en esas circunstancias podría argumentarse que tiene sentido simplemente promediarlos.
[Tienes razón en preocuparte por sesgo, por supuesto. ... si el estimador tiene cierto sesgo en n=N/M, eso no se eliminará promediando. Sin embargo, es posible que pueda aprovechar técnicas de bootstrap para reducir el sesgo.]
De todas formas, incluso si esas cosas no son ciertas (no podemos asumir la distribución asintótica, por ejemplo), no significa que no podamos hacer algo acerca de MLE de manera distribuida, al menos en algunas situaciones particulares. Aquí hay un ejemplo de algo más que se podría hacer (pero ve los comentarios al final).
Si las probabilidades son continuas y suaves (/no muy onduladas), podemos evaluar toda la función de log-probabilidad en algún intervalo (sobre una cuadrícula bastante fina) y pasar esa información, sumando las probabilidades al final.
Por ejemplo, digamos que N es 10^10 y M es 1000. Pasamos 10^7 observaciones a cada computadora y producimos la log-probabilidad en una cuadrícula fina de valores (digamos algo del orden de 10^3 a 10^4 valores), combinando todas las log-probabilidades de alguna manera adecuada (con el objetivo de mantener la precisión y no consumir demasiado tiempo, como sumar en pares en paralelo, lo que llevará $\sim\log_2(N/M)$ ciclos de suma al final). Esto nos permite lidiar con casos como la multimodalidad, pero necesitamos tener alguna idea sensata de dónde colocar nuestras cuadrículas para comenzar, lo que puede requerir algún preprocesamiento (o algún conocimiento de dominio muy específico). Entonces, cuando identificamos la ubicación del máximo general en el último paso, es posible que podamos suavizar nuestra función en la cuadrícula discreta tratándola nuevamente como cuadrática en el vecindario inmediato del óptimo, y así obtener una estimación considerablemente mejor que la precisión de la cuadrícula.
Esta idea implica pasar de vuelta (y posteriormente sumar) vectores para la probabilidad en lugar del argmax solamente.
Este tipo de enfoque se puede adaptar a algunas otras circunstancias distintas a las suaves y continuas, pero ese es un caso muy común con el que lidiar.
Este enfoque no debería sufrir el problema de retener el sesgo en $n=N/M$ que obtendrías al simplemente promediar los MLE.
Aquí hay un ejemplo ilustrativo sobre el parámetro de forma de una densidad gamma:
![enter image description here]()
El negro representa el log-likelihood completo hecho en la cuadrícula, el verde y el azul son los log-likelihoods para dos mitades del conjunto de datos, y el rojo es su suma. Si graficas el rojo sobre el negro, coinciden (como era de esperarse).
Al ajustar una cuadrática en el pico sugiere que el máximo está en aproximadamente 0.4603 (el máximo de la cuadrícula es 0.46, el verdadero parámetro era 0.5).
En la práctica, podrías hacerlo en dos pasos ... sobre una cuadrícula gruesa primero, luego sobre una cuadrícula mucho más fina cerca del máximo.
Este enfoque solo sería razonable para un problema con muy pocos parámetros. (No se mencionó la dimensión en la pregunta, pero en algunas situaciones la dimensión es solo 1 o 2 o 3). No recomendaría esto cuando hay más de un número muy pequeño de dimensiones.
