from __future__ import division

from matplotlib import pyplot as plt
import numpy as np
from scipy.stats import gumbel_r

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

fig, ax = plt.subplots(figsize=(8, 6))

n = 100
xs = np.linspace(-5, 5, n)

ax.plot(xs, logistic_probability(xs));

ax.set_xlim(-5, 5);
ax.set_xlabel('$x$');

ax.set_ylabel('$P(Y = 1 | X = x)$');

N = 500

eps0 = gumbel_r.rvs(size=(N, n))
eps1 = gumbel_r.rvs(size=(N, n))

U0 = eps0
U1 = xs + eps1

simulated_ps = (U1 > U0).mean(axis=0)

fig, ax = plt.subplots(figsize=(8, 6))

ax.plot(xs, logistic_probability(xs), c='k', lw=1.5, label='Exact');
ax.plot(xs, simulated_ps, c='r', label='Simulated');

ax.set_xlim(-5, 5);
ax.set_xlabel('$x$');

ax.set_ylabel('$P(Y = 1 | X = x)$');

ax.legend(loc=2);