Estudio: Encontrar las estimaciones de máxima verosimilitud de los parámetros de una función de densidad - actualizado

Question

ACTUALIZADO

Estoy tratando de encontrar el máximo de estimación de la probabilidad de una función de distribución de probabilidad dada a continuación

$$\begin{equation} g(x)=\frac{1}{\Gamma \left( \alpha \right)\gamma^{2\alpha}2^{\alpha-1}}x^{2\alpha-1}\exp\left\{{-\frac{x^2}{2\gamma^{2}}}\right\}I_{{\rm I\!R}^{+}}(x) \end{equation}$$

donde $\alpha >0$ es la forma de parámetros, $\sigma >0$ es el parámetro de escala.

La función de probabilidad está dada entonces por

$$\begin{equation} L(\alpha,\gamma/x)=\prod\limits_{i=1}^{n}\frac{1}{\Gamma \left( \alpha \right)\gamma^{2\alpha}2^{\alpha-1}}x_i^{2\alpha-1}\exp\left\{{-\frac{x_i^2}{2\gamma^{2}}}\right\} \end{equation}$$

Por lo tanto, la probabilidad de la función es entonces $$\begin{equation} L(\alpha,\gamma/x)=\frac{1}{[\Gamma \left( \alpha \right)]^{n}\gamma^{2\alpha n}{2^{n\alpha-n}}} \exp\left\{{-\frac{1}{2\gamma^{2}}\sum\limits_{i=1}^{n}x_{i}^{2}}\right\}\left(\prod\limits_{i=1}^{n}x_{i}\right)^{2\alpha-1} \end{equation}$$

Ahora, la función de verosimilitud logarítmica denotado por $\ell$ es

$$\begin{equation} \begin{aligned} \ell &=\log[L(\alpha,\gamma/x)]\\ &=-n\log(\Gamma \left( \alpha \right))-2\alpha n \log(\gamma)-n(\alpha-1)\log(2)-\frac{1}{2\gamma^{2}}\sum\limits_{i=1}^{n}x_{i}^{2}+(2\alpha-1)\sum\limits_{i=1}^{n}\log(x_{i}) \end{aligned} \end{equation}$$

Las entradas de la puntuación de la función está dada por

$$\begin{equation} \begin{aligned} \frac{\partial \ell}{\partial \alpha}=-n\psi(\alpha)-2n\log(\gamma)-n\log(2)+2\sum\limits_{i=1}^{n}\log(x_{i}) \end{aligned} \end{equation}$$ donde $\psi(\alpha)$ es la función digamma y

$$\begin{equation} \begin{aligned} \frac{\partial \ell}{\partial \gamma}=-\frac{2\alpha n}{\gamma}+\frac{\sum\limits_{i=1}^{n}x_{i}^{2}}{\gamma^{3}} \end{aligned} \end{equation}$$

La configuración de estas dos ecuaciones a cero y resolver simultáneamente los resultados en un máximo de estimaciones de probabilidad (MLE) de los parámetros, $\hat{\alpha}$$\hat{\gamma}$. Sin embargo, las ecuaciones obtenidas mediante la configuración de las anteriores derivadas parciales a cero no son en forma cerrada y los valores de los parámetros de $\alpha$ $\gamma$ debe encontrarse mediante métodos iterativos.

Matriz de información de Fisher se define como $I_{ij}=-E\left\{\frac{\partial^{2} \ell}{\partial \tau_i \partial \tau_j} \log[L(x_i, \vec{\tau})]\ \right\}$ donde$\tau_1=\alpha$$\tau_2=\gamma$. Por lo tanto, la matriz de información para la gamma-distribución de rayleigh es dada por,

$$\begin{equation} I=n \left[ \begin{array}{cc} \psi_{1}(\alpha) & 2/\gamma\\ 2/\gamma & 4\alpha/\gamma^2 \end{array} \right] \end{equation}$$

Estoy tratando de usar Fisher Scoring para encontrar Emv de los parámetros. Aquí está mi código de MATLAB. Yo primero generar 1000 observaciones aleatorias a partir de la gamma-distribución y de ejecutar este código. Mis valores de partida y el resto se da en el código.

clear all;
clc;

%Simulate 1000 sample from Gamma Distribution
n=1000;
alpha=3;
lambda=0.05;
x=gamrnd(alpha,1/lambda,1,n);

figure(1)
histfit(x,8,'gam');

sumlogx=sum(log(x)); sumxsquare=sum(x.^2);

%Initial Values
alpha=mean(x)^2/var(x);
gam=mean(x)/var(x);
theta=[alpha; gam];
S=Inf;

while sum(abs(S) > 10^(-5)) > 0
    S=[-n*psi(theta(1))-2*n*log(theta(2))-n*log(2)+2*sumlogx;...
        (-2*theta(1)*n/theta(2))+(sumxsquare/(theta(2)^3))];
    FIM=n*[psi(1, theta(1)), 2/theta(2);...
        2/theta(2), 4*theta(1)/(theta(2)^2)];
    theta=theta + FIM\S;
end

alpha_hat=theta(1)
gam_hat=theta(2)

fprintf('alpha_hat=%g, gamma_hat=%g \n', theta(1),theta(2))

Pero por algunas razones que no puedo entender, me estoy poniendo "Error en la utilización de psi X debe ser no negativo." error. Mi $\alpha$ de los valores están siendo negativa en la iteración de alguna manera, y no sé cómo solucionarlo!

Yo también soy ejecución de Newton-Raphson cuyo código de MATLAB es la siguiente

clear all;
clc;

%Simulate 100 sample from Gamma Distribution
n=1000;
alpha=3;
lambda=0.05;
x=gamrnd(alpha,1/lambda,1,n);

figure(1)
histfit(x,8,'gam');

sumlogx=sum(log(x)); sumxsquare=sum(x.^2);

%tuning parameters scale=gamma; shape=alpha
itermin=10^-7;
maxiter=10^7;
sc_init=0.000001;
sh_init=0.000001;
converged=[0;0;sc_init;sh_init];

% pdf
pdf=@(x,gam,alpha) 1/(gamma(alpha)*(gam^(2*alpha))*(2^(alpha-1)))*(x^(2*alpha-1))*exp(-(x^2)/(2*(gam^2)));

%score function is the first partial derivative of the log likelihood function
score=@(gam,alpha) -n*psi(alpha)-2*n*log(gam)-n*log(2)+2*sumlogx;

%Hessian function is the negative of the 2nd
hessian=@(gam,alpha) psi(1, alpha);

sc_loop=2; 
scale_hat=zeros(1,maxiter); 
scale_hat(1)=sc_init;

while 1==1
sh_loop=2;
shape_hat=zeros(1,maxiter);
shape_hat(1)=sh_init;

while 1==1
%calculate chat as chat_prev+score(chat_prev)/hessian(chat_prev)
shape_hat(sh_loop)=shape_hat(sh_loop-1)+score(scale_hat(sc_loop-1),shape_hat(sh_loop-1))/hessian(scale_hat(sc_loop-1),shape_hat(sh_loop-1));
%test for a convergence
if abs(shape_hat(sh_loop)-shape_hat(sh_loop-1))<itermin
    break %the process converged to a c value
elseif sh_loop>maxiter
    disp(['max iteration on \alpha achieved:', num2str(maxiter)]);
    return
end
sh_loop=sh_loop+1;
end

scale_hat(sc_loop)=(sum(x.^shape_hat(sh_loop-1))/n)^(1/shape_hat(sh_loop-1));
 %test for a convergence
  if abs(scale_hat(sc_loop)-scale_hat(sc_loop-1))<itermin
        break %the process converged to a gamma value
  end

  converged=[converged,[sc_loop-1;sh_loop-1;scale_hat(sc_loop);shape_hat(sh_loop)]];
  sc_loop=sc_loop+1;
end

%final display
disp(repmat('-',[1,30])),disp(' Iteration Scale Shape'),disp(repmat('-',[1,30]))
disp(num2str(converged','%6.4f')),disp(repmat('-',[1,30]))
disp(['Real values: gamma=', num2str(gam),',alpha=',num2str(alpha)])

Estoy recibiendo el mismo "Error en la utilización de psi, X debe ser no negativo." error.

Podría usted ayudarme? Algo está mal con psi función y no sé. Tal vez debería usar la aproximación, pero no estoy seguro de cuánta de la información que voy a perder!

Answer 1

[Nota: Esta es mi respuesta a la Dec. 19, 2014, versión de la que se trate.]

Si se opera el cambio de variable $y=x^2$ en la densidad de $$f_X(x|\alpha,\beta\sigma)=\frac{1}{\Gamma \left( \alpha \right)\beta^{\alpha}}\exp\left\{{-\frac{x^2}{2\sigma^{2}}\frac{1}{\beta}}\right\}\frac{x^{2\alpha-1}}{2^{\alpha-1}\sigma^{2\alpha}}\mathbb{I}_{{\mathbb{R}}^{+}}(x)$$ the Jacobian is given by $\dfrac{\text{d}y}{\text{d}x}= 2x = 2y^{1/2}$ y por lo tanto $$\begin{align*} f_Y(y|\alpha,\beta,\sigma)&=\frac{1}{\Gamma \left( \alpha \right)\beta^{\alpha}}\exp\left\{{-\frac{y}{2\sigma^{2}}\frac{1}{\beta}}\right\}\frac{y^{\frac{2\alpha-1}{2}}}{2^{\alpha-1}\sigma^{2\alpha}}\frac{1}{2 y^{1/2}}\mathbb{I}_{{\mathbb{R}}^{+}}(y)\\ &=\frac{1}{\Gamma \left( \alpha \right)\beta^{\alpha}}\exp\left\{{-\frac{y}{2\sigma^{2}}\frac{1}{\beta}}\right\}\frac{y^{{\alpha-1}}}{2^{\alpha}\sigma^{2\alpha}}\mathbb{I}_{{\mathbb{R}}^{+}}(y) \end{align*}$$ Esto demuestra que

1. Este es un estándar $\mathcal{G}(\alpha,2\sigma^2\beta)$ modelo, es decir, se observan $$(x_1^2,\ldots,x_n^2)=(y_1,\ldots,y_n)\stackrel{\text{iid}}{\sim}\mathcal{G}(\alpha,\eta);$$
2. el modelo es parametrizada ya que sólo $\eta=2\sigma^2\beta$ puede ser identificado;
3. EM no es necesario encontrar el MLE de $(\alpha,\eta)$, que no está disponible en forma cerrada, pero la solución de $$\hat\eta^{-1}=\bar{y}/\hat{\alpha}n\qquad\log(\hat{\alpha})-\psi(\hat{\alpha})=\log(\bar{y})-\frac{1}{n}\sum_{i=1}^n\log(y_i)$$ where $\psi(\cdot)$ es la di-de la función gamma. Este artículo de Thomas Minka indica rápido aproximaciones a la resolución de la ecuación anterior.

Answer 2

En el caso del algoritmo EM, los valores iniciales puede ser fijado de manera arbitraria ya que las iteraciones se garantiza la convergencia a la máxima:

Hemos visto que tanto la E y la M pasos del algoritmo EM están aumentando el valor de un bien definido obligada en el registro de probabilidad de la función y que el EM ciclo de cambio de los parámetros del modelo de tal manera como para hacer que el registro de probabilidad para aumentar (a menos que ya esté en un máximo, en cuyo caso los parámetros permanecen sin cambios).

Hay varias estrategias para recoger los valores iniciales que mejorará el rendimiento general de los algoritmos:

Tenga en cuenta que el algoritmo EM toma muchas más iteraciones para llegar a (aproximado) la convergencia en comparación con el K-means el algoritmo, y que cada ciclo se requiere significativamente más cálculo. Por tanto, es común para ejecutar el K-means el algoritmo de con el fin de encontrar un adecuado inicialización de un modelo de mezcla de Gaussianas que es posteriormente adaptado con EM. Las matrices de covarianza puede muy bien ser inicializado para la muestra de covarianzas de los clusters encontrados por el K-means el algoritmo, y la mezcla de los coeficientes se pueden establecer para las fracciones de puntos de datos asignado a la respectivos grupos.

Obispo, C. M. (2006). Reconocimiento de patrones y aprendizaje de máquina (Vol. 1, pág. 740). Nueva York: springer.

Usted puede leer sobre esto en el capítulo 9 del libro.

El uso de EM para la estimación de variables latentes en el contexto de Modelos Gaussianos Mezcla de dos componentes ( $\Delta_i = 0$ $\Delta_i = 1$ ), loglikelihood está dada por

(Lo siento, tengo prisa y no puedo escribir esto por mí mismo)

ESLII da algunos consejos sobre cómo seleccionar los valores iniciales:

Una buena manera de construir inicial estimación de $\mu_1$ $\mu_2$ es simplemente elegir dos de las $y_i$ al azar. Tanto en $\sigma_1^2$ $\sigma_2^2$ puede ser igual al total de la varianza de la muestra. La proporción de mezcla $\hat{\pi}$ puede ser puesto en marcha en el valor de 0.5.

Hastie, T., Tibshirani, R.,, Friedman, J. (2008). Los elementos de aprendizaje estadístico: minería de datos, inferencia y predicción. Springer.

Insisto, se puede elegir arbitrariamente los valores iniciales, pero algunos van a converger más rápido.