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.

Adaptación de distribuciones univariadas personalizadas

En este ejemplo se muestra cómo utilizar la función estadísticas y herramientas de aprendizaje automático™ para adecuar las distribuciones personalizadas a los datos univariados.mle

Utilizando, puede calcular las estimaciones de parámetros de máxima verosimilitud y estimar su precisión, para muchos tipos de distribuciones más allá de aquellas para las que el cuadro de herramientas proporciona funciones de ajuste específicas.mle

Para ello, debe definir la distribución utilizando una o varias funciones M. En los casos más simples, puede escribir código para calcular la función de densidad de probabilidad (PDF) para la distribución que desea ajustar y hará la mayor parte del trabajo restante por usted.mle Este ejemplo cubre esos casos. En los problemas con los datos censurados, también debe escribir código para calcular la función de distribución acumulativa (CDF) o la función de supervivencia (SF). En algunos otros problemas, puede ser ventajoso definir la función de log-verosimilitud (LLF) sí mismo. La segunda parte de este ejemplo, abarca ambos casos.Adaptación de distribuciones univariadas personalizadas, parte 2

Adaptación de distribuciones personalizadas: Un ejemplo de Poisson truncado cero

Los datos de recuento a menudo se modelan mediante una distribución de Poisson, y puede usar la función estadísticas y Toolbox de aprendizaje automático para ajustarse a un modelo de Poisson.poissfit Sin embargo, en algunas situaciones, los recuentos que son cero no se registran en los datos, por lo que el ajuste de una distribución de Poisson no es sencillo debido a esos "ceros faltantes". Este ejemplo mostrará cómo ajustar una distribución de Poisson a datos truncados cero, utilizando la función.mle

En este ejemplo, usaremos datos simulados de una distribución de Poisson truncada cero. Primero, generamos algunos datos aleatorios de Poisson.

rng(18,'twister'); n = 75; lambda = 1.75; x = poissrnd(lambda,n,1); 

A continuación, eliminamos todos los ceros de los datos para simular el truncamiento.

x = x(x > 0); length(x) 
 ans =      65  

Este es un histograma de estos datos simulados. Tenga en cuenta que los datos se ven razonablemente como una distribución de Poisson, excepto que no hay ceros. Los ajustaremos con una distribución que es idéntica a un Poisson en los enteros positivos, pero que no tiene ninguna probabilidad en cero. De esta manera, podemos estimar el parámetro de Poisson lambda mientras se hace una contabilización de los "ceros faltantes".

hist(x,[0:1:max(x)+1]); 

El primer paso es definir la distribución de Poisson truncada cero por su función de probabilidad (PF). Crearemos una función para calcular la probabilidad para cada punto en x, dado un valor para la lambda de parámetro medio de la distribución de Poisson. El PF para un Poisson truncado cero es sólo el PF de Poisson habitual, renormalizado de modo que suma a uno. Con cero truncamiento, la renormalización es sólo 1-PR {0}. La forma más sencilla de crear una función para el PF es utilizar una función anónima.

pf_truncpoiss = @(x,lambda) poisspdf(x,lambda) ./ (1-poisscdf(0,lambda)); 

Para simplificar, hemos asumido que todos los valores x dados a esta función serán enteros positivos, sin comprobaciones. La comprobación de errores, o una distribución más complicada, probablemente llevaría más de una sola línea de código, lo que sugiere que la función debe definirse en un archivo independiente.

El siguiente paso es proporcionar una primera conjetura aproximada razonable para el parámetro lambda. En este caso, sólo usaremos la media de la muestra.

start = mean(x) 
 start =      2.2154  

Proporcionamos con los datos, y con la función anónima, utilizando el parámetro ' pdf '.mle (El Poisson es discreto, así que esto es realmente una función de probabilidad, no un PDF.) Dado que el parámetro medio de la distribución de Poisson debe ser positivo, también se especifica un límite inferior para lambda. Devuelve la estimación de máxima verosimilitud de lambda y, opcionalmente, los intervalos de confianza aproximados del 95% para los parámetros.mle

[lambdaHat,lambdaCI] = mle(x, 'pdf',pf_truncpoiss, 'start',start, 'lower',0) 
 lambdaHat =      1.8760   lambdaCI =      1.4990     2.2530  

Observe que la estimación del parámetro es menor que la media de la muestra. Así es como debería ser, porque la estimación de máxima verosimilitud representa los ceros que faltan que no están presentes en los datos.

También podemos calcular una estimación de error estándar para lambda, utilizando la aproximación de varianza de muestra grande devuelta por.mlecov

avar = mlecov(lambdaHat, x, 'pdf',pf_truncpoiss); stderr = sqrt(avar) 
 stderr =      0.1923  

Suministro de valores adicionales a la función de distribución: Un ejemplo normal truncado

A veces también ocurre que los datos continuos se truncan. Por ejemplo, es posible que las observaciones mayores que algún valor fijo no se registren debido a limitaciones en la forma en que se recopilan los datos. En este ejemplo se muestra cómo ajustar una distribución normal a los datos truncados, utilizando la función.mle

Para este ejemplo, simulamos datos de una distribución normal truncada. En primer lugar, generamos algunos datos normales aleatorios.

n = 75; mu = 1; sigma = 3; x = normrnd(mu,sigma,n,1); 

A continuación, eliminamos cualquier observación que caiga más allá del punto de truncamiento, xTrunc. A lo largo de este ejemplo, asumiremos que xTrunc es conocido y no necesita ser estimado.

xTrunc = 4; x = x(x < xTrunc); length(x) 
 ans =      64  

Este es un histograma de estos datos simulados. Los ajustaremos con una distribución que es idéntica a una normal para x < xTrunc, pero que tiene una probabilidad cero por encima de xTrunc. De esta manera, podemos estimar los parámetros normales MU y Sigma mientras que la contabilización de la "cola perdida".

hist(x,(-10:.5:4)); 

Como en el ejemplo anterior, definiremos la distribución normal truncada por su PDF, y crearemos una función para calcular la densidad de probabilidad para cada punto en x, valores dados para los parámetros MU y Sigma. Con el punto de truncamiento fijo y conocido, el PDF para una normal truncada es sólo el pdf normal habitual, truncado, y luego ponderaciones renormalizadas para que se integre a uno. La renormalización es sólo el CDF evaluado en xTrunc. Para simplificar, supondremos que todos los valores x son menores que xTrunc, sin comprobación. Usaremos una función anónima para definir el PDF.

pdf_truncnorm = @(x,mu,sigma) normpdf(x,mu,sigma) ./ normcdf(xTrunc,mu,sigma); 

El punto de truncamiento, xTrunc, no se está calculando, por lo que no está entre los parámetros de distribución en la lista de argumentos de entrada de la función PDF. xTrunc tampoco forma parte del argumento de entrada de vector de datos. Con una función anónima, podemos simplemente hacer referencia a la variable xTrunc que ya se ha definido en el espacio de trabajo, y no hay necesidad de preocuparse de pasarlo como un argumento adicional.

También necesitamos proporcionar una conjetura inicial aproximada para las estimaciones de parámetros. En este caso, debido a que el truncamiento no es demasiado extremo, la media de la muestra y la desviación estándar probablemente funcionarán bien.

start = [mean(x),std(x)] 
 start =      0.2735    2.2660  

Proporcionamos con los datos, y con la función anónima, utilizando el parámetro ' pdf '.mle Dado que Sigma debe ser positivo, también especificamos límites de parámetros más bajos. Devuelve las estimaciones de máxima verosimilitud de MU y Sigma como un solo Vector, así como una matriz de intervalos de confianza aproximados del 95% para los dos parámetros.mle

[paramEsts,paramCIs] = mle(x, 'pdf',pdf_truncnorm, 'start',start, 'lower',[-Inf 0]) 
 paramEsts =      0.9911    2.7800   paramCIs =     -0.1605    1.9680     2.1427    3.5921  

Observe que las estimaciones de MU y Sigma son bastante más grandes que la media de la muestra y la desviación estándar. Esto se debe a que el ajuste del modelo ha representado la cola superior "ausente" de la distribución.

Podemos calcular una matriz de covarianza aproximada para las estimaciones de parámetros utilizando.mlecov La aproximación es generalmente razonable en muestras grandes, y los errores estándar de las estimaciones pueden ser aproximados por las raíces cuadradas de los elementos diagonales.

acov = mlecov(paramEsts, x, 'pdf',pdf_truncnorm) stderr = sqrt(diag(acov)) 
 acov =      0.3452    0.1660     0.1660    0.1717   stderr =      0.5876     0.4143  

Ajuste de una distribución más complicada: Una mezcla de dos normales

Algunos datasets exhiben bimodalidad, o incluso multimodalidad, y el ajuste de una distribución estándar a estos datos normalmente no es apropiado. Sin embargo, una mezcla de distribuciones unimodales simples a menudo puede modelar estos datos muy bien. De hecho, incluso puede ser posible dar una interpretación a la fuente de cada componente de la mezcla, basándose en el conocimiento específico de la aplicación.

En este ejemplo, ajustaremos una mezcla de dos distribuciones normales a algunos datos simulados. Esta mezcla podría describirse con la siguiente definición constructiva para generar un valor aleatorio:

  First, flip a biased coin.  If it lands heads, pick a value at random   from a normal distribution with mean mu_1 and standard deviation   sigma_1. If the coin lands tails, pick a value at random from a normal   distribution with mean mu_2 and standard deviation sigma_2.

Para este ejemplo, vamos a generar datos a partir de una mezcla de distribuciones t de Student en lugar de usar el mismo modelo que estamos encajando. Este es el tipo de cosas que puede hacer en una simulación Monte-Carlo para probar cuán robusto es un método de ajuste a las desviaciones de las suposiciones del modelo que se ajustan. Aquí, sin embargo, encajaremos solo un conjunto de datos simulado.

x = [trnd(20,1,50) trnd(4,1,100)+3]; hist(x,-2.25:.5:7.25); 

Como en los ejemplos anteriores, definiremos el modelo para ajustarlos creando una función que calcule la densidad de la probabilidad. El PDF para una mezcla de dos normales es sólo una suma ponderada de los PDF de los dos componentes normales, ponderados por la probabilidad de la mezcla. Este PDF es lo suficientemente simple para crear usando una función anónima. La función toma seis entradas: un vector de datos en el que evaluar el PDF y los cinco parámetros de la distribución. Cada componente tiene parámetros para su media y desviación estándar; la probabilidad de la mezcla hace un total de cinco.

pdf_normmixture = @(x,p,mu1,mu2,sigma1,sigma2) ...                          p*normpdf(x,mu1,sigma1) + (1-p)*normpdf(x,mu2,sigma2); 

También necesitaremos una suposición inicial para los parámetros. Cuantos más parámetros tenga un modelo, más importa un punto de partida razonable. Para este ejemplo, comenzaremos con una mezcla igual (p = 0,5) de normales, centrada en los dos cuartes de los datos, con desviaciones estándar iguales. El valor inicial de la desviación estándar proviene de la fórmula para la varianza de una mezcla en términos de la media y la varianza de cada componente.

pStart = .5; muStart = quantile(x,[.25 .75]) sigmaStart = sqrt(var(x) - .25*diff(muStart).^2) start = [pStart muStart sigmaStart sigmaStart]; 
 muStart =      0.5970    3.2456   sigmaStart =      1.2153  

Por último, necesitamos especificar los límites de cero y uno para la probabilidad de mezcla, y los límites inferiores de cero para las desviaciones estándar. Los elementos restantes de los vectores de límites se establecen en + inf o-inf, para indicar que no hay restricciones.

lb = [0 -Inf -Inf 0 0]; ub = [1 Inf Inf Inf Inf]; paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, ...                           'lower',lb, 'upper',ub) 
Warning: Maximum likelihood estimation did not converge.  Iteration limit exceeded.   paramEsts =      0.3523    0.0257    3.0489    1.0546    1.0629  

El valor predeterminado para las distribuciones personalizadas es 200 iteraciones.

statset('mlecustom') 
 ans =     struct with fields:            Display: 'off'       MaxFunEvals: 400           MaxIter: 200            TolBnd: 1.0000e-06            TolFun: 1.0000e-06        TolTypeFun: []              TolX: 1.0000e-06          TolTypeX: []           GradObj: 'off'          Jacobian: []         DerivStep: 6.0555e-06       FunValCheck: 'on'            Robust: []      RobustWgtFun: []            WgtFun: []              Tune: []       UseParallel: []     UseSubstreams: []           Streams: {}         OutputFcn: []  

Invalide ese valor predeterminado, utilizando una estructura de opciones creada con la función.statset Aumente también el límite de evaluación de funciones.

options = statset('MaxIter',300, 'MaxFunEvals',600); paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, ...                           'lower',lb, 'upper',ub, 'options',options) 
 paramEsts =      0.3523    0.0257    3.0489    1.0546    1.0629  

Parece que las iteraciones finales a la convergencia sólo importaba en los últimos dígitos del resultado. No obstante, siempre es una buena idea asegurarse de que se ha alcanzado la convergencia.

Por último, podemos trazar la densidad ajustada contra un histograma de probabilidad de los datos sin procesar, para comprobar el ajuste visualmente.

bins = -2.5:.5:7.5; h = bar(bins,histc(x,bins)/(length(x)*.5),'histc'); h.FaceColor = [.9 .9 .9]; xgrid = linspace(1.1*min(x),1.1*max(x),200); pdfgrid = pdf_normmixture(xgrid,paramEsts(1),paramEsts(2),paramEsts(3),paramEsts(4),paramEsts(5)); hold on plot(xgrid,pdfgrid,'-') hold off xlabel('x') ylabel('Probability Density') 

Uso de funciones anidadas: Un ejemplo normal con precisiones desiguales

A veces sucede cuando se recopilan datos de que cada observación se realizó con una precisión o confiabilidad diferentes. Por ejemplo, si varios experimentadores hacen un número de mediciones independientes de la misma cantidad, pero cada uno notifica sólo el promedio de sus mediciones, la confiabilidad de cada punto de datos reportado dependerá del número de observaciones crudas que fueron en él. Si los datos sin procesar originales no están disponibles, una estimación de su distribución debe tener en cuenta el hecho de que los datos que están disponibles, los promedios, cada uno tienen diferentes desviaciones. Este modelo tiene realmente una solución explícita para las estimaciones de parámetros de máxima verosimilitud. Sin embargo, a los efectos de la ilustración, utilizaremos para estimar los parámetros.mle

Supongamos que tenemos 10 puntos de datos, donde cada uno es en realidad el promedio de cualquier lugar de 1 a 8 observaciones. Esas observaciones originales no están disponibles, pero sabemos cuántos había para cada uno de nuestros puntos de datos. Necesitamos estimar la media y la desviación estándar de los datos sin procesar.

x = [0.25 -1.24 1.38 1.39 -1.43 2.79 3.52 0.92 1.44 1.26]; m = [   8     2    1    3     8    4    2    5    2    4]; 

La varianza de cada punto de datos es inversamente proporcional al número de observaciones que entraron en él, por lo que utilizaremos 1/m para ponderar la varianza de cada punto de datos en un ajuste de máxima verosimilitud.

w = 1 ./ m 
 w =    Columns 1 through 7      0.1250    0.5000    1.0000    0.3333    0.1250    0.2500    0.5000    Columns 8 through 10      0.2000    0.5000    0.2500  

En el modelo que estamos encajando aquí, podríamos definir la distribución por su PDF, pero el uso de un archivo PDF de registro es algo más natural, porque el PDF normal es de la forma

  c .* exp(-0.5 .* z.^2),

y tendría que tomar el registro del PDF de todos modos, para calcular la log-verosimilitud.mle Así que en su lugar, crearemos una función que calcula el PDF de registro directamente.

La función log PDF tiene que computar el registro de la densidad de probabilidad para cada punto en x, valores dados para MU y Sigma. También tendrá que tener en cuenta las diferentes ponderaciones de desviación. A diferencia de los ejemplos anteriores, esta función de distribución es un poco más complicada que una sola línea, y se escribe más claramente como una función independiente en su propio archivo. Dado que la función de PDF de registro necesita los recuentos de observación como datos adicionales, la forma más directa de lograr este ajuste es usar funciones anidadas.

Hemos creado un archivo independiente para una función llamada wgtnormfit.m. Esta función contiene la inicialización de datos, una función anidada para el PDF de registro en el modelo normal ponderado y una llamada a la función para ajustarse realmente al modelo.mle Debido a que Sigma debe ser positivo, debemos especificar límites de parámetros más bajos. La llamada a devuelve las estimaciones de máxima verosimilitud de MU y Sigma en un solo vector.mle

type wgtnormfit.m 
 function paramEsts = wgtnormfit %WGTNORMFIT Fitting example for a weighted normal distribution.  %   Copyright 1984-2012 The MathWorks, Inc.   x = [0.25 -1.24 1.38 1.39 -1.43 2.79 3.52 0.92 1.44 1.26]'; m = [   8     2    1    3     8    4    2    5    2    4]';      function logy = logpdf_wn(x,mu,sigma)         v = sigma.^2 ./ m;         logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);     end  paramEsts = mle(x, 'logpdf',@logpdf_wn, 'start',[mean(x),std(x)], 'lower',[-Inf,0]);  end 

En, pasamos un manejador a la función anidada, utilizando el parámetro ' logpdf '.wgtnormfit.mmlelogpdf_wn Esa función anidada hace referencia a los recuentos de observaciones, m, en el cálculo del PDF de registro ponderado. Dado que el vector m se define en su función principal, tiene acceso a él, y no hay necesidad de preocuparse de pasar m en como un argumento de entrada explícito.logpdf_wn

Tenemos que proporcionar una primera suposición aproximada para las estimaciones de parámetros. En este caso, la media de la muestra no ponderada y la desviación estándar deben estar bien, y eso es lo que usa.wgtnormfit.m

start = [mean(x),std(x)] 
 start =      1.0280    1.5490  

Para ajustar el modelo, ejecutamos la función de ajuste.

paramEsts = wgtnormfit 
 paramEsts =      0.6244    2.8823  

Observe que la estimación de MU es inferior a dos tercios de la media de la muestra. Así es como debería ser: la estimación se ve influenciada más por los puntos de datos más fiables, es decir, los que se basaron en el mayor número de observaciones crudas. En este conjunto de datos, esos puntos tienden a extraer la estimación de la media de la muestra no ponderada.

Uso de una transformación de parámetros: El ejemplo normal (continuación)

En la estimación de máxima verosimilitud, los intervalos de confianza para los parámetros se calculan normalmente utilizando una aproximación normal de gran tamaño para la distribución de los estimadores. Esto es a menudo una suposición razonable, pero con tamaños de muestra pequeños, a veces es ventajoso para mejorar esa aproximación normal mediante la transformación de uno o más parámetros. En este ejemplo, tenemos un parámetro Location y un parámetro scale. Los parámetros de escala a menudo se transforman en su registro, y lo haremos aquí con Sigma. En primer lugar, crearemos una nueva función de registro PDF y luego calcularemos de nuevo las estimaciones usando esa parametrización.

La nueva función de registro PDF se crea como una función anidada dentro de la función wgtnormfit2.m. Como en el primer ajuste, este archivo contiene la inicialización de datos, una función anidada para el PDF de registro en el modelo normal ponderado y una llamada a la función para ajustarse realmente al modelo.mle Dado que Sigma puede ser cualquier valor positivo, log (Sigma) es ilimitado, y ya no necesitamos especificar límites inferiores o superiores. Además, la llamada a en este caso devuelve las estimaciones de parámetros y los intervalos de confianza.mle

type wgtnormfit2.m 
 function [paramEsts,paramCIs] = wgtnormfit2 %WGTNORMFIT2 Fitting example for a weighted normal distribution (log(sigma) parameterization).  %   Copyright 1984-2012 The MathWorks, Inc.   x = [0.25 -1.24 1.38 1.39 -1.43 2.79 3.52 0.92 1.44 1.26]'; m = [   8     2    1    3     8    4    2    5    2    4]';      function logy = logpdf_wn2(x,mu,logsigma)         v = exp(logsigma).^2 ./ m;         logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);     end  [paramEsts,paramCIs] = mle(x, 'logpdf',@logpdf_wn2, 'start',[mean(x),log(std(x))]);  end 

Tenga en cuenta que utiliza el mismo punto de partida, transformado a la nueva parametrización, es decir, el registro de la desviación estándar de la muestra.wgtnormfit2.m

start = [mean(x),log(std(x))] 
 start =      1.0280    0.4376  
[paramEsts,paramCIs] = wgtnormfit2 
 paramEsts =      0.6244    1.0586   paramCIs =     -0.2802    0.6203     1.5290    1.4969  

Dado que la parametrización utiliza log (Sigma), tenemos que volver a la escala original para obtener un intervalo de estimación y confianza para Sigma. Observe que las estimaciones para MU y Sigma son las mismas que en el primer ajuste, porque las estimaciones de máxima verosimilitud son invariables a la parametrización.

muHat = paramEsts(1) sigmaHat = exp(paramEsts(2)) 
 muHat =      0.6244   sigmaHat =      2.8823  
muCI = paramCIs(:,1) sigmaCI = exp(paramCIs(:,2)) 
 muCI =     -0.2802     1.5290   sigmaCI =      1.8596     4.4677