def _create_intermediate_nodes(d, c, f_name, h_name, verbose=False):
    """Returns 'intermediate' nodes; i.e., false-alarm and hit probabilites.

    """
    f = pm.Deterministic(
        f_name,
        (1 + T.erf((-d/2 - c)/T.sqrt(2))) / 2
    )
    h = pm.Deterministic(
        h_name,
        (1 + T.erf((d/2 - c)/T.sqrt(2))) / 2
    )

    return f, h

    def __init__(self, random_state=None, mu=0.0, sigma=1.0):
        super(Normal, self).__init__(mu=mu,
                                     sigma=sigma,
                                     random_state=random_state,
                                     optimizer=None)

        # pdf
        self.pdf_ = (
            (1. / np.sqrt(2. * np.pi)) / self.sigma *
            T.exp(-(self.X - self.mu) ** 2 / (2. * self.sigma ** 2))).ravel()
        self.make_(self.pdf_, "pdf")

        # -log pdf
        self.nnlf_ = bound(
            T.log(self.sigma) + T.log(np.sqrt(2. * np.pi)) +
            (self.X - self.mu) ** 2 / (2. * self.sigma ** 2),
            np.inf,
            self.sigma > 0.).ravel()
        self.make_(self.nnlf_, "nnlf")

        # cdf
        self.cdf_ = 0.5 * (1. + T.erf((self.X - self.mu) /
                                      (self.sigma * np.sqrt(2.)))).ravel()
        self.make_(self.cdf_, "cdf")

        # ppf
        self.ppf_ = (self.mu +
                     np.sqrt(2.) * self.sigma * T.erfinv(2. * self.p - 1.))
        self.make_(self.ppf_, "ppf", args=[self.p])

 def __init__(self, x, mu, sigma, *args, **kwargs):
     super(Normal, self).__init__(*args, **kwargs)
     self._logp = bound(-(x - mu)**2 / (2 * sigma**2) + T.log(1 / T.sqrt(sigma**2 * 2 * np.pi)), sigma > 0)
     self._cdf = 0.5 * (1 + T.erf((x - mu)/(sigma*T.sqrt(2))))
     self._add_expr('x', x)
     self._add_expr('mu', mu)
     self._add_expr('sigma', sigma)

 def get_output_for(self, input, **kwargs):
     eps = 1e-7
     x_axis = theano.shared(np.arange(0, 600, dtype='float32')).dimshuffle('x',0)
     # This needs to be clipped to avoid NaN's!
     sigma = T.exp(T.clip(input[:,1].dimshuffle(0,'x'), -10, 10))
     #theano_printer.print_me_this("sigma", sigma)
     x = (x_axis - input[:,0].dimshuffle(0,'x')) / (sigma * np.sqrt(2).astype('float32'))
     return (T.erf(x) + 1)/2

def theano_mu_sigma_erf(mu_erf, sigma_erf, eps=1e-7):
    x_axis = theano.shared(np.arange(0, 600, dtype='float32')).dimshuffle('x',0)
    if sigma_erf.ndim==0:
        sigma_erf = T.clip(sigma_erf.dimshuffle('x','x'), eps, 1)
    elif sigma_erf.ndim==1:
        sigma_erf = T.clip(sigma_erf.dimshuffle(0,'x'), eps, 1)
    x = (x_axis - mu_erf.dimshuffle(0,'x')) / (sigma_erf * np.sqrt(2).astype('float32'))
    return (T.erf(x) + 1)/2

 def forward_theano(self, x):
     abs_x = tt.abs_(x)
     y = tt.switch(abs_x < self.c, tt.erf(x / 2.**0.5),
                   (((self.beta**2 - 4 * self.alpha *
                     (self.gamma - abs_x))**0.5
                    - self.beta) /
                    (2 * self.alpha)) * tt.sgn(x))
     return y

def lognormal_cdf_math(x, mu, sigma, eps=1e-12):
    # wikipedia claims cdf is
    # .5 + .5 erf( log(x) - mu / sqrt(2 sigma^2))
    #
    # the maximum is used to move negative values and 0 up to a point
    # where they do not cause nan or inf, but also don't contribute much
    # to the cdf.
    return 0.5 + 0.5 * tensor.erf((tensor.log(tensor.maximum(x, eps)) - mu) / tensor.sqrt(2 * sigma ** 2))

def cdf(x, location=0, scale=1):
    location = T.cast(location, theano.config.floatX)
    scale = T.cast(scale, theano.config.floatX)

    div = T.sqrt(2 * scale ** 2 + epsilon)
    div = T.cast(div, theano.config.floatX)

    erf_arg = (x - location) / div
    return .5 * (1 + T.erf(erf_arg + epsilon))

    def get_output_for(self, input, **kwargs):
        mu = input[0]
        sigma = input[1]

        x_range = T.arange(0, self.max_support).dimshuffle('x', 0)
        mu = T.repeat(mu, self.max_support, axis=1)
        sigma = T.repeat(sigma, self.max_support, axis=1)
        x = (x_range - mu) / (sigma * T.sqrt(2.) + 1e-16)
        cdf = (T.erf(x) + 1.) / 2.
        return cdf

 def logp(self, value):
     tau = self.tau
     sd = self.sd
     mu = self.mu
     alpha = self.alpha
     return bound(
         tt.log(1 +
         tt.erf(((value - mu) * tt.sqrt(tau) * alpha) / tt.sqrt(2)))
         + (-tau * (value - mu)**2
         + tt.log(tau / np.pi / 2.)) / 2.,
         tau > 0, sd > 0)

    def get_output_for(self, input, **kwargs):
        mu = input[0]
        sigma = input[1]
        w = input[2]
        if self.log:
            sigma = T.exp(sigma)

        x_range = T.arange(0, 600).dimshuffle('x', 0, 'x')
        mu = mu.dimshuffle(0, 'x', 1)
        sigma = sigma.dimshuffle(0, 'x', 1)
        x = (x_range - mu) / (sigma * T.sqrt(2.) + 1e-16)
        cdf = (T.erf(x) + 1.) / 2.  # (bs, 600, n_mix)
        cdf = T.sum(cdf * w.dimshuffle(0, 'x', 1), axis=-1)
        return cdf

    def get_output_for(self, input, **kwargs):
        if input.ndim > 3:
            # input: (batch, time, axis, verti, horiz)
            # needs: (batch, time, pixels)
            input = input.flatten(ndim=3)

        eps=1e-7
        clipped_input = T.clip(input, eps, 1-eps)
        mu = T.sum(clipped_input, axis=2).dimshuffle(0,1,'x')

        sigma = T.sqrt(T.sum(clipped_input * (1-clipped_input), axis=2).dimshuffle(0,1,'x') + eps)
        x_axis = theano.shared(np.arange(0, 600, dtype='float32')).dimshuffle('x','x',0)
        x = (x_axis - mu) / sigma
        return (T.erf(x) + 1)/2

def cdf(sample, location=0, scale=1):
    """Return a theano expression representing the values of the cumulative
    density function of a Gaussian distribution.

    Parameters
    ----------

    sample : Theano variable
        Array of shape ``(n,)`` where ``n`` is the number of samples.

    location : Theano variable
        Scalar representing the mean of the distribution.

    scale : Theano variable
        Scalar representing the standard deviation of the distribution.

    Returns
    -------

    l : Theano variable
        Array of shape ``(n,)`` where each entry represents the cumulative
        density of the corresponding sample.


    Examples
    --------

    >>> import theano
    >>> import theano.tensor as T
    >>> import numpy as np
    >>> from breze.learn.utils import theano_floatx
    >>> sample, mean, std = T.vector(), T.scalar(), T.scalar()
    >>> c = cdf(sample, mean, std)
    >>> f_c = theano.function([sample, mean, std], c)

    >>> X, = theano_floatx(np.array([-1, 0, 1]))
    >>> cs = f_c(X, 0.1, 1.2)
    >>> np.allclose(cs, [0.17965868, 0.46679324, 0.77337265])
    True
    """
    location = T.cast(location, theano.config.floatX)
    scale = T.cast(scale, theano.config.floatX)

    div = T.sqrt(2 * scale ** 2 + epsilon)
    div = T.cast(div, theano.config.floatX)

    erf_arg = (sample - location) / div
    return .5 * (1 + T.erf(erf_arg + epsilon))

    def s_expectation_lt_thresh(self, x, thresh):
        """
        return \int_{-inf}^{thresh} (thresh-y)*p(y|x) dy
        
        p(y | x) = gaussian with center mu(x) and variance sigma(x)**2

        """

        mu = self.s_mean(x)
        sigma = tensor.sqrt(
                tensor.maximum(self.s_variance(x),
                    self.min_variance))        
        a = 0.5 * (mu - thresh) 
        delta = (thresh - mu) / (sqrt(2) * sigma)
        sbar = sigma / sqrt(2 * pi)
        rval = sbar * tensor.exp(-delta ** 2) - a * (1 + tensor.erf(delta))
        rval = tensor.maximum(rval, 1e-7)
        
        if rval.dtype != self.dtype:
            raise TypeError('rval dtype', rval.dtype)
        return rval

    def __init__(self, mu=0.0, sigma=1.0):
        """Constructor.

        Parameters
        ----------
        * `mu` [float]:
            The distribution mean.

        * `sigma` [float]:
            The distribution standard deviation.
        """
        super(Normal, self).__init__(mu=mu, sigma=sigma)

        # pdf
        self.pdf_ = (
            (1. / np.sqrt(2. * np.pi)) / self.sigma *
            T.exp(-(self.X - self.mu) ** 2 / (2. * self.sigma ** 2))).ravel()
        self._make(self.pdf_, "pdf")

        # -log pdf
        self.nll_ = bound(
            T.log(self.sigma) + T.log(np.sqrt(2. * np.pi)) +
            (self.X - self.mu) ** 2 / (2. * self.sigma ** 2),
            np.inf,
            self.sigma > 0.).ravel()
        self._make(self.nll_, "nll")

        # cdf
        self.cdf_ = 0.5 * (1. + T.erf((self.X - self.mu) /
                                      (self.sigma * np.sqrt(2.)))).ravel()
        self._make(self.cdf_, "cdf")

        # ppf
        self.ppf_ = (self.mu +
                     np.sqrt(2.) * self.sigma * T.erfinv(2. * self.p - 1.))
        self._make(self.ppf_, "ppf", args=[self.p])

def cdf(z):
  """Cumulative distribution function via erf (Error function)"""
  return (numpy.float32(1) + T.erf(z)) / numpy.float32(2)

def std_cdf(x):
    """
    Calculates the standard normal cumulative distribution function.
    """
    return 0.5 + 0.5 * tt.erf(x / tt.sqrt(2.))

def normcdf(X, nu = 0, sigma=1):
    return 0.5 * (1 + TT.erf( (X-nu) / (sigma * 2**0.5)))

 def Phi(x):
     erfarg = (x - avg) / (std * SQRT2)
     rval = 0.5 * (1. + T.erf(erfarg))
     return rval.astype(dtype)

def cdf(x, miu=0.0, variance=1.0):
    return 1.0 / 2 * (1.0 + T.erf((x - miu) / T.sqrt(2 * variance)))