Wie werden künstliche Daten für die logistische Regression simuliert?

45

Ich weiß, dass ich etwas in meinem Verständnis der logistischen Regression vermisse und würde mich über jede Hilfe sehr freuen.

Nach meinem Verständnis geht die logistische Regression davon aus, dass die Wahrscheinlichkeit eines 1-Ergebnisses bei den Eingaben eine lineare Kombination der Eingaben ist, die durch eine inverse logistische Funktion geleitet wird. Dies wird im folgenden R-Code veranschaulicht:

#create data:
x1 = rnorm(1000)           # some continuous variables 
x2 = rnorm(1000)
z = 1 + 2*x1 + 3*x2        # linear combination with a bias
pr = 1/(1+exp(-z))         # pass through an inv-logit function
y = pr > 0.5               # take as '1' if probability > 0.5

#now feed it to glm:
df = data.frame(y=y,x1=x1,x2=x2)
glm =glm( y~x1+x2,data=df,family="binomial")

und ich bekomme folgende Fehlermeldung:

Warnmeldungen: 1: glm.fit: Algorithmus konvergierte nicht 2: glm.fit: angepasste Wahrscheinlichkeiten mit 0 oder 1 sind aufgetreten

Ich habe jetzt einige Zeit mit R gearbeitet; genug, um zu wissen, dass wahrscheinlich ich die Schuld bin .. was passiert hier?

zorbar
quelle
2
Die Art und Weise, wie Sie Ihre Daten simulieren, kommt mir komisch vor. Wenn Sie eine Alternative für einen Standard suchen, können Sie hier nachschauen
ocram
@ocram: du hast recht; Dies ist eine doppelte Frage!
user603
2
Ich habe eine fehlerhafte Simulation durchgeführt, wie @ Stéphane Laurent erklärte. Das Problem war jedoch die perfekte Trennung in der logistischen Regression , ein Problem, mit dem ich nicht vertraut war und über das ich ziemlich überrascht war.
Zorbar
@zorbar: es war in meiner antwort auf deine frage (jetzt gelöscht).
user603
1
@ user603: Ich habe wahrscheinlich deine Antwort verpasst;
Trotzdem

Antworten:

63

yich1pr(ich)

> set.seed(666)
> x1 = rnorm(1000)           # some continuous variables 
> x2 = rnorm(1000)
> z = 1 + 2*x1 + 3*x2        # linear combination with a bias
> pr = 1/(1+exp(-z))         # pass through an inv-logit function
> y = rbinom(1000,1,pr)      # bernoulli response variable
> 
> #now feed it to glm:
> df = data.frame(y=y,x1=x1,x2=x2)
> glm( y~x1+x2,data=df,family="binomial")

Call:  glm(formula = y ~ x1 + x2, family = "binomial", data = df)

Coefficients:
(Intercept)           x1           x2  
     0.9915       2.2731       3.1853  

Degrees of Freedom: 999 Total (i.e. Null);  997 Residual
Null Deviance:      1355 
Residual Deviance: 582.9        AIC: 588.9 
Stéphane Laurent
quelle
Du hast recht - ich habe diesen Schritt verpasst. Vielen dank für Deine Hilfe!
Zorbar
1
Ich hatte eine Frage zur Art und Weise, wie Sie Daten simulieren. Wenn wir Daten für die lineare Regression simulieren, simulieren wir auch ein Rauschen (\ epsilon). Ich verstehe, dass die logistische Funktion eine Funktion der Erwartung ist, die selbst den Lärm auslöscht. Ist das der Grund, warum du kein Geräusch in deinem z hast?
Sam
5
mittlere Antwort+Lärm
@ StéphaneLaurent, genau. Ich verstehe völlig. Vielen Dank für Ihre Antwort.
Sam
2

LogisticRegression eignet sich zur Anpassung, wenn Wahrscheinlichkeiten oder Proportionen als Ziele angegeben werden und nicht nur 0/1 Ergebnisse.

import numpy as np
import pandas as pd
def logistic(x, b, noise=None):
    L = x.T.dot(b)
    if noise is not None:
        L = L+noise
    return 1/(1+np.exp(-L))

x = np.arange(-10., 10, 0.05)
bias = np.ones(len(x))
X = np.vstack([x,bias]) # Add intercept
B =  [1., 1.] # Sigmoid params for X

# True mean
p = logistic(X, B)
# Noisy mean
pnoisy = logistic(X, B, noise=np.random.normal(loc=0., scale=1., size=len(x)))
# dichotomize pnoisy -- sample 0/1 with probability pnoisy
dichot = np.random.binomial(1., pnoisy)

pd.Series(p, index=x).plot(style='-')
pd.Series(pnoisy, index=x).plot(style='.')
pd.Series(dichot, index=x).plot(style='.')

Hier haben wir drei mögliche Ziele für die logistische Regression. Dies pist das wahre / Ziel-Verhältnis / die Wahrscheinlichkeit, pnoisydas / die p ist, wobei normales Rauschen in die logarithmische Wahrscheinlichkeitsskala aufgenommen wird, und dichotdas / die als Parameter für das binomische PDF behandelt und daraus abgetastet wird. Sie sollten alle 3 testen. Ich habe festgestellt, dass einige Open-Source-LR-Implementierungen nicht passen p.

Abhängig von Ihrer Anwendung bevorzugen Sie möglicherweise pnoisy.

In der Praxis sollten Sie auch überlegen, wie wahrscheinlich das Rauschen in Ihrer Zielanwendung geformt wird, und versuchen, dies zu emulieren.

user48956
quelle