Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Análisis bayesiano para un modelo de regresión logística

Este ejemplo muestra cómo hacer inferencias bayesianas para un modelo de regresión logística usando.slicesample

Las inferencias estadísticas generalmente se basan en la estimación de máxima verosimilitud (MLE). MLE elige los parámetros que maximizan la probabilidad de los datos, y es intuitivamente atractiva. En MLE, se supone que los parámetros son desconocidos pero fijos, y se estiman con cierta confianza. En las estadísticas bayesianas, la incertidumbre sobre los parámetros desconocidos se cuantifica utilizando la probabilidad para que los parámetros desconocidos se consideren variables aleatorias.

Inferencia bayesiana

La inferencia bayesiana es el proceso de análisis de modelos estadísticos con la incorporación de conocimientos previos sobre los parámetros del modelo o modelo. La raíz de tal inferencia es el teorema de Bayes:

Por ejemplo, supongamos que tenemos observaciones normales

donde se conoce Sigma y la distribución previa de Theta es

En esta fórmula MU y tau, a veces conocido como hiperparámetros, también se conocen. Si observamos muestras de, podemos obtener la distribución posterior para Theta comonX

El siguiente gráfico muestra el anterior, probabilidad, y posterior para Theta.

rng(0,'twister');  n = 20; sigma = 50; x = normrnd(10,sigma,n,1); mu = 30; tau = 20; theta = linspace(-40, 100, 500); y1 = normpdf(mean(x),theta,sigma/sqrt(n)); y2 = normpdf(theta,mu,tau); postMean = tau^2*mean(x)/(tau^2+sigma^2/n) + sigma^2*mu/n/(tau^2+sigma^2/n); postSD = sqrt(tau^2*sigma^2/n/(tau^2+sigma^2/n)); y3 = normpdf(theta, postMean,postSD); plot(theta,y1,'-', theta,y2,'--', theta,y3,'-.') legend('Likelihood','Prior','Posterior') xlabel('\theta') 

Datos del experimento del coche

En algunos problemas simples como el ejemplo de inferencia media anterior normal, es fácil averiguar la distribución posterior en un formulario cerrado. Pero en general los problemas que implican los Priores no conjugados, las distribuciones posteriores son difíciles o imposibles de calcular analíticamente. Consideraremos la regresión logística como un ejemplo. Este ejemplo implica un experimento para ayudar a modelar la proporción de automóviles de varios pesos que fallan una prueba de kilometraje. Los datos incluyen observaciones de peso, número de coches probados y número fallido. Trabajaremos con una versión transformada de los pesos para reducir la correlación en nuestras estimaciones de los parámetros de regresión.

% A set of car weights weight = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'; weight = (weight-2800)/1000;     % recenter and rescale % The number of cars tested at each weight total = [48 42 31 34 31 21 23 23 21 16 17 21]'; % The number of cars that have poor mpg performances at each weight poor = [1 2 0 3 8 8 14 17 19 15 17 21]'; 

Modelo de regresión logística

La regresión logística, un caso especial de un modelo lineal generalizado, es adecuada para estos datos, ya que la variable de respuesta es binomial. El modelo de regresión logística se puede escribir como:

donde X es la matriz de diseño y b es el vector que contiene los parámetros del modelo. En MATLAB®, podemos escribir esta ecuación como:

logitp = @(b,x) exp(b(1)+b(2).*x)./(1+exp(b(1)+b(2).*x)); 

Si tiene algún conocimiento previo o algunos antecedentes no informativos están disponibles, puede especificar las distribuciones de probabilidad previas para los parámetros del modelo. Por ejemplo, en este ejemplo, usamos los antecedentes normales para la interceptación y la pendiente, es decir,b1b2

prior1 = @(b1) normpdf(b1,0,20);    % prior for intercept prior2 = @(b2) normpdf(b2,0,20);    % prior for slope 

Por el teorema de Bayes, la distribución posterior conjunta de los parámetros del modelo es proporcional al producto de la probabilidad y los Priores.

post = @(b) prod(binopdf(poor,total,logitp(b,weight))) ...  % likelihood             * prior1(b(1)) * prior2(b(2));                  % priors 

Tenga en cuenta que la constante de normalización para el posterior en este modelo es analíticamente inintractable. Sin embargo, incluso sin conocer la constante de normalización, puede visualizar la distribución posterior, si conoce el rango aproximado de los parámetros del modelo.

b1 = linspace(-2.5, -1, 50); b2 = linspace(3, 5.5, 50); simpost = zeros(50,50); for i = 1:length(b1)     for j = 1:length(b2)         simpost(i,j) = post([b1(i), b2(j)]);     end; end; mesh(b2,b1,simpost) xlabel('Slope') ylabel('Intercept') zlabel('Posterior density') view(-110,30) 

Este posterior es alargado a lo largo de una diagonal en el espacio de parámetros, indicando que, después de mirar los datos, creemos que los parámetros están correlacionados. Esto es interesante, ya que antes de recopilar cualquier dato asumimos que eran independientes. La correlación proviene de la combinación de nuestra distribución previa con la función de probabilidad.

Muestreo de rebanada

Los métodos de Montecarlo se utilizan a menudo en el análisis de datos bayesianos para resumir la distribución posterior. La idea es que, incluso si no se puede calcular la distribución posterior analíticamente, se puede generar una muestra aleatoria a partir de la distribución y utilizar estos valores aleatorios para estimar la distribución posterior o estadísticas derivadas como la media posterior, mediana, desviación estándar, etc. El muestreo de rebanada es un algoritmo diseñado para muestrear a partir de una distribución con una función de densidad arbitraria, conocida sólo hasta una constante de proporcionalidad-exactamente lo que se necesita para el muestreo de una distribución posterior complicada cuya normalización constante es desconocida. El algoritmo no genera muestras independientes, sino más bien una secuencia de Markovian cuya distribución estacionaria es la distribución de destino. Por lo tanto, el muestreador de rebanada es un algoritmo de cadena de Markov Monte Carlo (MCMC). Sin embargo, difiere de otros algoritmos MCMC conocidos, ya que solo se especifica la necesidad de escalado posterior, no se necesitan distribuciones marginales ni propuestas.

Este ejemplo muestra cómo utilizar el muestreador de sectores como parte de un análisis Bayesiano del modelo de regresión logística de prueba de kilometraje, incluida la generación de una muestra aleatoria a partir de la distribución posterior para los parámetros del modelo, el análisis de la salida del muestreador y hacer inferencias sobre los parámetros del modelo. El primer paso es generar una muestra aleatoria.

initial = [1 1]; nsamples = 1000; trace = slicesample(initial,nsamples,'pdf',post,'width',[20 2]); 

Análisis de salida de muestreador

Después de obtener una muestra aleatoria del muestreador de sectores, es importante investigar cuestiones como la convergencia y la mezcla, para determinar si la muestra puede tratarse razonablemente como un conjunto de realizaciones aleatorias de la distribución posterior de destino. Mirar las parcelas de traza marginal es la forma más sencilla de examinar la salida.

subplot(2,1,1) plot(trace(:,1)) ylabel('Intercept'); subplot(2,1,2) plot(trace(:,2)) ylabel('Slope'); xlabel('Sample Number'); 

De estas parcelas se desprende que los efectos de los valores iniciales de los parámetros toman un tiempo para desaparecer (quizás 50 o más muestras) antes de que el proceso empiece a parecer estacionario.

También es útil en la comprobación de la convergencia para utilizar una ventana móvil para calcular estadísticas como la media de la muestra, la mediana o la desviación estándar para la muestra. Esto produce un trazado más suave que los seguimientos de muestra sin procesar y puede facilitar la identificación y comprensión de cualquier no estacionariedad.

movavg = filter( (1/50)*ones(50,1), 1, trace); subplot(2,1,1) plot(movavg(:,1)) xlabel('Number of samples') ylabel('Means of the intercept'); subplot(2,1,2) plot(movavg(:,2)) xlabel('Number of samples') ylabel('Means of the slope'); 

Dado que se trata de promedios móviles sobre una ventana de 50 iteraciones, los primeros 50 valores no son comparables al resto de la gráfica. Sin embargo, el resto de cada parcela parece confirmar que el parámetro posterior significa haber convergido a estacionariedad después de 100 o así iteraciones. También es evidente que los dos parámetros están correlacionados entre sí, de acuerdo con la trama anterior de la densidad posterior.

Dado que el período de asentamiento representa muestras que no se pueden tratar razonablemente como realizaciones aleatorias de la distribución de destino, probablemente es aconsejable no usar los primeros 50 o así valores al principio de la salida del muestreador de sectores. Podría simplemente eliminar esas filas de la salida, sin embargo, también es posible especificar un período de "burn-in". Esto es conveniente cuando ya se conoce una longitud de quemado adecuada, tal vez de las corridas anteriores.

trace = slicesample(initial,nsamples,'pdf',post, ...                                      'width',[20 2],'burnin',50); subplot(2,1,1) plot(trace(:,1)) ylabel('Intercept'); subplot(2,1,2) plot(trace(:,2)) ylabel('Slope'); 

Estas parcelas no parecen mostrar ninguna estacionariedad, lo que indica que el período de combustión ha hecho su trabajo.

Sin embargo, hay un segundo aspecto de las parcelas de traza que también se deben explorar. Mientras que el rastro para la intercepción parece el ruido de alta frecuencia, el rastro para la pendiente parece tener un componente de la frecuencia más baja, indicando allí la autocorrelación entre los valores en las iteraciones adyacentes. Todavía podríamos calcular la media de esta muestra autocorrelacionada, pero a menudo es conveniente reducir los requisitos de almacenamiento eliminando la redundancia en la muestra. Si esto eliminara la autocorrelación, también nos permitiría tratarlo como una muestra de valores independientes. Por ejemplo, puede diluir la muestra manteniendo sólo cada valor 10.

trace = slicesample(initial,nsamples,'pdf',post,'width',[20 2], ...                                                 'burnin',50,'thin',10); 

Para comprobar el efecto de este adelgazamiento, es útil estimar las funciones de autocorrelación de muestra de las trazas y utilizarlas para comprobar si las muestras se mezclan rápidamente.

F    =  fft(detrend(trace,'constant')); F    =  F .* conj(F); ACF  =  ifft(F); ACF  =  ACF(1:21,:);                          % Retain lags up to 20. ACF  =  real([ACF(1:21,1) ./ ACF(1,1) ...              ACF(1:21,2) ./ ACF(1,2)]);       % Normalize. bounds  =  sqrt(1/nsamples) * [2 ; -2];       % 95% CI for iid normal  labs = {'Sample ACF for intercept','Sample ACF for slope' }; for i = 1:2     subplot(2,1,i)     lineHandles  =  stem(0:20, ACF(:,i) , 'filled' , 'r-o');     lineHandles.MarkerSize = 4;     grid('on')     xlabel('Lag')     ylabel(labs{i})     hold on     plot([0.5 0.5 ; 20 20] , [bounds([1 1]) bounds([2 2])] , '-b');     plot([0 20] , [0 0] , '-k');     hold off     a  =  axis;     axis([a(1:3) 1]); end 

Los valores de autocorrelación en el primer desfase son significativos para el parámetro de intercepción, y más aún para el parámetro Slope. Podríamos repetir el muestreo usando un parámetro de simplificación más grande para reducir aún más la correlación. Sin embargo, para los fines de este ejemplo, seguiremos usando el ejemplo actual.

Inferencia para los parámetros del modelo

Como era de esperar, un histograma de la muestra imita la trama de la densidad posterior.

subplot(1,1,1) hist3(trace,[25,25]); xlabel('Intercept') ylabel('Slope') zlabel('Posterior density') view(-110,30) 

Puede utilizar un histograma o una estimación de densidad de suavizado del kernel para resumir las propiedades de distribución marginal de las muestras posteriores.

subplot(2,1,1) hist(trace(:,1)) xlabel('Intercept'); subplot(2,1,2) ksdensity(trace(:,2)) xlabel('Slope'); 

También puede calcular estadísticas descriptivas como la media posterior o los percentiles de las muestras aleatorias. Para determinar si el tamaño de la muestra es lo suficientemente grande para lograr una precisión deseada, es útil supervisar la estadística deseada de las trazas como una función del número de muestras.

csum = cumsum(trace); subplot(2,1,1) plot(csum(:,1)'./(1:nsamples)) xlabel('Number of samples') ylabel('Means of the intercept'); subplot(2,1,2) plot(csum(:,2)'./(1:nsamples)) xlabel('Number of samples') ylabel('Means of the slope'); 

En este caso, parece que el tamaño de la muestra de 1000 es más que suficiente para dar buena precisión para la estimación media posterior.

bHat = mean(trace) 
 bHat =     -1.6931    4.2569  

Resumen

El™ de herramientas de estadísticas y aprendizaje automático ofrece una variedad de funciones que le permiten especificar las probabilidades y antecedentes fácilmente. Se pueden combinar para derivar una distribución posterior. La función le permite realizar análisis bayesianos en MATLAB utilizando la simulación de cadena Monte Carlo de Markov.slicesample Se puede utilizar incluso en problemas con distribuciones posteriores que son difíciles de muestrear utilizando generadores de números aleatorios estándar.