¿Por qué la regresión logística no aprende las ponderaciones correctas para una variable Y definida como logit de X y ponderaciones?

Question

Estoy intentando ver si la implementación de la regresión logística en mi es correcta. Para ello primero empecé a generar algunos datos aleatorios X y especificar algunos pesos predefinidos.

Cuando entreno el modelo utilizando una implementación estándar de regresión logística, el modelo no aprende las ponderaciones correctas especificadas previamente. ¿Estoy haciendo algo mal en mi aplicación?

Este es mi código:

import numpy as np
import statsmodels.api as sm

X = sm.add_constant(np.random.randn(1000,3)*10)
w = np.array([0.25, 0.1, 0.3, -1.5])

def logit(x):
    return 1. / (1. + np.exp(x))

y = logit(np.dot(X,w) + np.random.randn(X.shape[0])*0.0001) # y ~ logit(X.w + Random Noise)
model = sm.Logit(y>0.5, X) # Label y as 1 if y > 0.5 else 0
res = model.fit()
print res.summary2()

Este es el resultado que obtengo de mi modelo:

Warning: Maximum number of iterations has been exceeded.
         Current function value: 0.000000
         Iterations: 35
                          Results: Logit
==================================================================
Model:              Logit            Pseudo R-squared: 1.000      
Dependent Variable: y                AIC:              8.0000     
Date:               2016-02-21 18:08 BIC:              27.6310    
No. Observations:   1000             Log-Likelihood:   -6.5515e-06
Df Model:           3                LL-Null:          -692.70    
Df Residuals:       996              LLR p-value:      4.3502e-300
Converged:          0.0000           Scale:            1.0000     
No. Iterations:     35.0000                                       
------------------------------------------------------------------
           Coef.    Std.Err.    z    P>|z|     [0.025     0.975]  
------------------------------------------------------------------
const     -28.4569 1024.9463 -0.0278 0.9779  -2037.3148  1980.4009
x1        -13.8000  411.2444 -0.0336 0.9732   -819.8241   792.2242
x2        -40.6588 1204.5892 -0.0338 0.9731  -2401.6104  2320.2927
x3        201.4672 5953.2122  0.0338 0.9730 -11466.6144 11869.5488
==================================================================

¿Qué estoy haciendo mal? ¿No debería el modelo aprender los coeficientes como las ponderaciones que especifiqué anteriormente?

Answer 1

No estás simulando un experimento Bernoulli. Además, tu y se puede predecir perfectamente y debería haber una advertencia o excepción de predicción perfecta.

Si cambio a un experimento Bernoulli, y aumento el tamaño de la muestra para que la estimación se acerque más a los valores verdaderos, y añado una semilla para la replicabilidad, entonces obtengo una buena coincidencia.

Creo que tu función Logit tiene el signo contrario al que se utiliza en los modelos estadísticos, es decir, utiliza ceros como caso de "éxito".

import numpy as np
import statsmodels.api as sm

np.random.seed(12345)
X = sm.add_constant(np.random.randn(20000,3)*10)
w = np.array([0.25, 0.1, 0.3, -1.5])

def logit(x):
    return 1. / (1. + np.exp(x))

y = logit(np.dot(X,w))
from scipy import stats
#yo = stats.bernoulli.rvs(y, size=y.shape[0])   # opposite sign
yo = 1 - stats.bernoulli.rvs(y, size=y.shape[0])
model = sm.Logit(yo, X)
res = model.fit()
print res.summary()

que produce

Optimization terminated successfully.
         Current function value: 0.084227
         Iterations 11
                         Results: Logit
================================================================
Model:              Logit            Pseudo R-squared: 0.878    
Dependent Variable: y                AIC:              3377.0925
Date:               2016-02-21 21:23 BIC:              3408.7065
No. Observations:   20000            Log-Likelihood:   -1684.5  
Df Model:           3                LL-Null:          -13858.  
Df Residuals:       19996            LLR p-value:      0.0000   
Converged:          1.0000           Scale:            1.0000   
No. Iterations:     11.0000                                     
------------------------------------------------------------------
         Coef.    Std.Err.      z       P>|z|     [0.025    0.975]
------------------------------------------------------------------
const    0.2499     0.0446     5.6100   0.0000    0.1626    0.3373
x1       0.1007     0.0050    20.2428   0.0000    0.0909    0.1104
x2       0.2964     0.0085    34.6648   0.0000    0.2796    0.3131
x3      -1.5137     0.0373   -40.5711   0.0000   -1.5868   -1.4405
================================================================

adición ampliar comentarios

El proceso subyacente de generación de datos de Logit y muchas otras funciones de enlace para experimentos Bernoulli también puede representarse mediante un latente y una ecuación de observación que limita lo que se puede observar de la variable latente.

Como señala dsaxton en el comentario, la distribución de la variable latente es la distribución logística.

Creo que el código Python se explica por sí mismo:

np.random.seed(987654)

linpred = np.dot(X,w)
r = stats.logistic.rvs(size=X.shape[0])    
y = linpred + r  # latent variable
yo = y >= 0  # observed binary variable

model = sm.Logit(yo, X)
res = model.fit()
print res.summary2()

Los parámetros estimados son (he ajustado un poco la semilla para obtener un "buen" sorteo aleatorio).

>>> print(res.params)
[ 0.26791279  0.09486941  0.30443008 -1.49688738]

La formulación de variable latente se utiliza a menudo en econometría porque se generaliza fácilmente, por ejemplo, a modelos truncados como Tobit o Logit o Probit ordenados.

Answer 2

Primero una cosa terminológica menor, pero la función que has definido se llama en realidad sigmoide y no logit, que es la inversa de la función sigmoidea. En segundo lugar te falta el signo menos y en su lugar debería ser $1 / (e^{- w^T x} + 1)$ .

Lo más importante es que estás añadiendo un término de error al sigmoide cuando no existe tal cosa en la regresión logística. La forma en que debe generar $Y_i$ es por muestreo a partir de una distribución Bernoulli con $1 / (e^{- w^T x_i} + 1)$ como probabilidad de éxito. Si hace esto, debería obtener las estimaciones correctas en su simulación.

Answer 3

Creo que el problema está en tu definición de "logit" (que en realidad es la función logística, la inversa de logit): debería ser $exp(-x)$ no $exp(x)$ .