Esta página es para la versión anterior. La página correspondiente en inglés ha sido eliminada en la versión actual.

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

En este ejemplo se muestra cómo realizar inferencias bayesianas para un modelo de regresión logística mediante .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 son intuitivamente atractivos. 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 como variables aleatorias.

Inferencia bayesiana

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

$$P(\mathrm{parameters|data}) = \frac
 {P(\mathrm{data|parameters}) \times P(\mathrm{parameters})}
 {P(\mathrm{data})}
 \propto \mathrm {likelihood} \times \mathrm{prior}$$

Por ejemplo, supongamos que tenemos observaciones normales

$$X|\theta \sim N(\theta, \sigma^2)$$

donde sigma es conocida y la distribución previa para theta es

$$\theta \sim N(\mu, \tau^2)$$

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

$$\theta|X \sim N\left(\frac{\tau^2}{\sigma^2/n + \tau^2} \bar X +
 \frac{\sigma^2/n}{{\sigma^2/n + \tau^2}} \mu,
 \frac{(\sigma^2/n)\times \tau^2}{\sigma^2/n +
 \tau^2}\right)$$

El siguiente gráfico muestra el anterior, la probabilidad y el 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 de experimentos de coche

En algunos problemas simples, como el ejemplo anterior de inferencia media normal, es fácil averiguar la distribución posterior en una forma cerrada. Pero en general los problemas que implican antecedentes 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 coches de varios pesos que no superan una prueba de kilometraje. Los datos incluyen observaciones de peso, número de coches probados y número fallado. 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:

$$P(\mathrm{failure}) = \frac{e^{Xb}}{1+e^{Xb}}$$

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 anteriores para los parámetros del modelo. Por ejemplo, en este ejemplo, usamos priores 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 

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

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 la posterior en este modelo es analíticamente intratable. 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 se alarga a lo largo de una diagonal en el espacio de parámetros, lo que indica que, después de mirar los datos, creemos que los parámetros están correlacionados. Esto es interesante, ya que antes de recopilar los datos 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 rebanadas

Los métodos de Monte Carlo se utilizan a menudo en el análisis de datos bayesianos para resumir la distribución posterior. La idea es que, incluso si no puede calcular la distribución posterior analíticamente, puede generar una muestra aleatoria a partir de la distribución y utilizar estos valores aleatorios para estimar la distribución posterior o las estadísticas derivadas, como la media posterior, la mediana, desviación estándar, etc. El muestreo de cortes es un algoritmo diseñado para muestrear desde 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 constante de normalización es desconocido. El algoritmo no genera muestras independientes, sino una secuencia de Markovian cuya distribución estacionaria es la distribución de destino. Por lo tanto, el sampler de rebanadas es un algoritmo Markov Chain Monte Carlo (MCMC). Sin embargo, difiere de otros algoritmos MCMC conocidos porque sólo se especifica la necesidad posterior escalada -- no se necesitan propuestas ni distribuciones marginales.

En este ejemplo se muestra cómo utilizar el muestreador de cortes 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 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 la salida del muestreador

Después de obtener una muestra aleatoria del muestreador de rebanadas, es importante investigar problemas 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 objetivo. Mirar las gráficas de trazado marginales 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 gráficas se desprende que los efectos de los valores iniciales del parámetro tardan un tiempo en desaparecer (quizás cincuenta o más muestras) antes de que el proceso comience a verse estacionario.

También es útil para comprobar la convergencia 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 una gráfica más suave que las trazas de muestra sin procesar, y puede facilitar la identificación y comprensión de cualquier no estacionario.

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 medias 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 las medias posteriores del parámetro han convergido a la estacionariedad después de más de 100 iteraciones. También es evidente que los dos parámetros están correlacionados entre sí, de acuerdo con la gráfica anterior de la densidad posterior.

Dado que el período de liquidación representa muestras que no se pueden tratar razonablemente como realizaciones aleatorias de la distribución de destino, es probable que sea aconsejable no utilizar los primeros 50 valores más o menos al principio de la salida del muestreador de cortes. Simplemente podría eliminar esas filas de la salida, sin embargo, también es posible especificar un período de "quemado". Esto es conveniente cuando ya se conoce una longitud de quemado adecuada, tal vez de ejecuciones 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 trazas no parecen mostrar ninguna no estacionalidad, lo que indica que el período de quemado ha hecho su trabajo.

Sin embargo, hay un segundo aspecto de las trazas que también deben explorarse. Mientras que el seguimiento de la interceptación parece ruido de alta frecuencia, el seguimiento de la pendiente parece tener un componente de frecuencia más baja, lo que indica que hay autocorrelación entre los valores en iteraciones adyacentes. Todavía podríamos calcular la media de este ejemplo autocorrelacionados, pero a menudo es conveniente reducir los requisitos de almacenamiento eliminando la redundancia en el ejemplo. Si esto elimina la autocorrelación, también nos permitiría tratar esto como una muestra de valores independientes. Por ejemplo, puede adelgazar la muestra manteniendo solo cada 10o valor.

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 muestras a partir 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 desajuste son significativos para el parámetro de interceptación, y aún más para el parámetro slope. Podríamos repetir el muestreo usando un parámetro de adelgazamiento 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 gráfica 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 núcleo 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 como para lograr la precisión deseada, es útil supervisar la estadística deseada de las trazas en 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 una buena precisión para la estimación media posterior.

bHat = mean(trace) 
 bHat =     -1.6931    4.2569  

Resumen

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