def likefunc(self,n,s,dels,delb):
     """log-likelihood function for individual parameter points in the model.
     Contains the two nuisance parameters dels and delb, which
     parameterise the systematic errors. Marginalise these out to be
     Bayesian, or profile them to be pseudo-frequentist (they still
     have priors).
     The parameter 's' (signal mean) should then be the only free 
     parameter left. 
     Args:
     n - observed number of events
     i - which signal region we are currently looking at
     dels - systematic error parameter for signal
     delb - systematic error parameter for background
     s - expected number of events due to signal
     ssys - estimated gaussian uncertainty on expected number of signal events (effectively a prior)
     b - expected number of events due to background
     bsys - estimated gaussian uncertainty on expected number of background events (effectively a prior)
     bstat - estimated "statistical" gaussian uncertainty on expected number of background events (also effectively a prior)
     K - signal efficiency scaling factor
     """
     #bsystot = np.sqrt(self.bsys**2 + self.bstat**2)                             # assume priors are independent
     siglike = logpoissonlike(n,self.sK*s*(1+dels*self.ssys)+self.b*(1+delb*self.bsystot))  # poisson signal + background log likelihood
     #Need to change the scaling of the prior to match the simulated data.
     #Makes no difference to inferences.
     Pdels = pymc.normal_like(dels,0,1) #+ 0.5*np.log(2*np.pi)            #standard normal gaussian log prior on dels
     Pdelb = pymc.normal_like(delb,0,1) #+ 0.5*np.log(2*np.pi)            #standard normal gaussian log prior on delb
     
     if siglike + Pdels + Pdelb < -1e200:
         print dels, delb
         print siglike,Pdels,Pdelb, self.sK*s*(1+dels*self.ssys)+self.b*(1+delb*self.bsystot)
         raise
     
     return siglike + Pdels + Pdelb

def X(value=X_true, K=K, A=A, mu = mu_x_init, tau = tau_x_init):
    """Autoregression"""

    # Initial data
    logp=normal_like(value[:K], mu, tau)

    # Difference equation
    for i in xrange(K,T):
        logp += normal_like(value[i], sum(A[:K]*value[i-K:i]), 1.)

    return logp

 def X_obs(pi=pi, sigma=sigma, value=X):
     logp = mc.normal_like(pl.array(value).ravel(), 
                           (pl.ones([N,J*T])*pl.array(pi).ravel()).ravel(), 
                           (pl.ones([N,J*T])*pl.array(sigma).ravel()).ravel()**-2)
     return logp
     
     logp = pl.zeros(N)
     for n in range(N):
         logp[n] = mc.normal_like(pl.array(value[n]).ravel(),
                                  pl.array(pi+beta).ravel(),
                                  pl.array(sigma).ravel()**-2)
     return mc.flib.logsum(logp - pl.log(N))

def plot_ratio_analysis(data_samples=(100,), dataset_samples=(100,), datasets=100):
    x, y = np.meshgrid(data_samples, dataset_samples)
    z = np.empty(x.shape, dtype=np.float)

    for i, data_sample in enumerate(data_samples):
        for j, dataset_sample in enumerate(dataset_samples):
            data = np.random.randn(x[j, i])
            errors = []
            sl_sum = 0
            pt_sum = 0
            for rep in range(1, 200):
                # Chose two random mu pts
                mu1 = (np.random.rand()-.5) * 3
                mu2 = (np.random.rand()-.5) * 3

                # Evaluate true likelihood
                pt1 = pm.normal_like(data, mu=mu1, tau=1)
                pt2 = pm.normal_like(data, mu=mu2, tau=1)

                ptr = pt1 / pt2
                pt_sum += pt1
                pt_sum += pt2

                #print ptr

                # Evaluate synth likelihood
                ps1 = synth_likelihood(data, mu1, 1, dataset_samples=y[j, i], samples=datasets)
                ps2 = synth_likelihood(data, mu2, 1, dataset_samples=y[j, i], samples=datasets)

                sl_sum += ps1
                sl_sum += ps2

                pts = ps1 / ps2
                #print pts

                errors.append((pts - ptr)**2)
            print pt_sum
            print sl_sum
            z[j, i] = np.mean(errors)
            print x[j, i], y[j,i], z[j, i]

    print x
    print y
    print z
    cont = plt.contourf(x, y, z)

    plt.colorbar(cont)
    plt.xlabel('Number of samples per dataset')
    plt.ylabel('Size of input data.')

def logdoublenormal(x,mean,sigmaP,sigmaM):
    #mean is measured value
    #x is computed theory value
    #sigmaP and sigmaM are distances from mean to upper and lower 1 sigma (68%)
    #confidence limits.
    if x==None: return -1e300
    if x>=mean:
        tauP = 1./sigmaP**2
        #need to remove the normalisation factor so we get the same normalisation
        #for each half of the likelihood.
        loglike = pymc.normal_like(x,mean,tauP) - pymc.normal_like(mean,mean,tauP)
    if x<mean:
        tauM = 1./sigmaM**2
        loglike = pymc.normal_like(x,mean,tauM) - pymc.normal_like(mean,mean,tauM)
    return loglike

 def x(N=N, mu=moo, tau=tau, n=n, value=np.log(data)):
    k = N-n
    dev = (value[0]-mu)*np.sqrt(tau)   
    out = gammaln(N+1) - gammaln(k) + (k-1)*np.log(pm.utils.normcdf(dev)) + pm.normal_like(value, mu, tau)
    if np.isnan(out):
        raise ValueError
    return out

 def obs(f=rate_stoch,
         age_indices=age_indices,
         age_weights=age_weights,
         value=d_val,
         tau=1./(d_se)**2):
     f_i = dismod3.utils.rate_for_range(f, age_indices, age_weights)
     return mc.normal_like(value, f_i, tau)

def get_likelihood_M0(map_M0, x, pwr, sigma, tau, obstype):
    A0 = get_variables_M0(map_M0)
    A0 = curve_fit_M0(x, pwr, A0, sigma)
    if obstype == '.logiobs':
        return pymc.normal_like(pwr, get_spectrum_M0(x, A0), tau)
    else:
        return pymc.lognormal_like(pwr, get_spectrum_M0(x, A0), tau)

def get_likelihood_M2(map_M2, x, pwr, sigma, tau, obstype):
    A2 = get_variables_M2(map_M2)
    A2 = curve_fit_M2(x, pwr, A2, sigma)
    if obstype == '.logiobs':
        return pymc.normal_like(pwr, get_spectrum_M2(x, A2), tau)
    else:
        return pymc.lognormal_like(pwr, get_spectrum_M2(x, A2), tau)

def get_likelihood_M1(map_M1, x, pwr, sigma, tau, obstype):
    A1 = get_variables_M1(map_M1)
    A1 = curve_fit_M1(x, pwr, A1, sigma)
    if obstype == '.logiobs':
        return pymc.normal_like(pwr, get_spectrum_M1(x, A1), tau)
    else:
        return pymc.lognormal_like(pwr, get_spectrum_M1(x, A1), tau)

    def covariate_constraint(mu=vars['mu_age'], alpha=vars['alpha'], beta=vars['beta'],
                             U_all=U_all,
                             X_sex_max=X_sex_max,
                             X_sex_min=X_sex_min,
                             lower=np.log(model.parameters[name]['level_bounds']['lower']),
                             upper=np.log(model.parameters[name]['level_bounds']['upper'])):
        log_mu_max = np.log(mu.max())
        log_mu_min = np.log(mu.min())

        alpha = np.array([float(x) for x in alpha])
        if len(alpha) > 0:
            for U_i in U_all:
                log_mu_max += max(0, alpha[U_i].max())
                log_mu_min += min(0, alpha[U_i].min())

        # this estimate is too crude, and is causing problems
        #if len(beta) > 0:
        #    log_mu_max += np.sum(np.maximum(X_max*beta, X_min*beta))
        #    log_mu_min += np.sum(np.minimum(X_max*beta, X_min*beta))

        # but leaving out the sex effect results in strange problems, too
        log_mu_max += X_sex_max*float(beta[sex_index])
        log_mu_min += X_sex_min*float(beta[sex_index])

        lower_violation = min(0., log_mu_min - lower)
        upper_violation = max(0., log_mu_max - upper)
        return mc.normal_like([lower_violation, upper_violation], 0., 1.e-6**-2)

def multi_normal_like(values, vec_mu, tau):
    """logp for multi normal"""
    logp = 0
    for i in range(len(vec_mu)):
        logp += pm.normal_like(values[i,:], vec_mu[i], tau)

    return logp

    def mixture(value=1., gamma=gamma, pi=[0.2, 0.8], mu=[-2., 3.],
                sigma=[0.01, 0.01]):
        """
        The log probability of a mixture of normal densities.

        :param value:       The point of evaluation.
        :type value :       float
        :param gamma:       The parameter characterizing the SMC one-parameter
                            family.
        :type gamma :       float
        :param pi   :       The weights of the components.
        :type pi    :       1D :class:`numpy.ndarray`
        :param mu   :       The mean of each component.
        :type mu    :       1D :class:`numpy.ndarray`
        :param sigma:       The standard deviation of each component.
        :type sigma :       1D :class:`numpy.ndarray`
        """
        # Make sure everything is a numpy array
        pi = np.array(pi)
        mu = np.array(mu)
        sigma = np.array(sigma)
        # The number of components in the mixture
        n = pi.shape[0]
        # pymc.normal_like requires the precision not the variance:
        tau = np.sqrt(1. / sigma ** 2)
        # The following looks a little bit awkward because of the need for
        # numerical stability:
        p = np.log(pi)
        p += np.array([pymc.normal_like(value, mu[i], tau[i])
                       for i in range(n)])
        p = math.fsum(np.exp(p))
        # logp should never be negative, but it can be zero...
        if p <= 0.:
            return -np.inf
        return gamma * math.log(p)

 def obs(f=vars['rate_stoch'],
         age_indices=age_indices,
         age_weights=age_weights,
         value=np.log(dm.value_per_1(d)),
         tau=se**-2, data=d):
     f_i = rate_for_range(f, age_indices, age_weights)
     return mc.normal_like(value, np.log(f_i), tau)

        def smooth_gamma(gamma=flat_gamma, knots=knots, tau=smoothing**-2):
            # the following is to include a "noise floor" so that level value
            # zero prior does not exert undue influence on age pattern
            # smoothing
            gamma = gamma.clip(pl.log(pl.exp(gamma).mean()/10.), pl.inf)  # only include smoothing on values within 10x of mean

            return mc.normal_like(pl.sqrt(pl.sum(pl.diff(gamma)**2 / pl.diff(knots))), 0, tau)

def plot_erroranalysis(data_samples=(10,), dataset_samples=(10,), datasets=200):
    x = dataset_samples
    y = np.empty(x.shape, dtype=np.float)

    for data_sample in data_samples:
        data = np.random.randn(data_sample)

        for i, dataset_sample in enumerate(dataset_samples):
            errors = []
            sl_sum = 0
            pt_sum = 0
            for rep in range(1, 400):
                # Chose two random mu pts
                mu1 = 0
                mu2 = (np.random.rand()-.5)

                # Evaluate true likelihood
                pt1 = pm.normal_like(data, mu=mu1, tau=1**-2)
                pt2 = pm.normal_like(data, mu=mu2, tau=1**-2)

                ptr = pt1 / pt2
                pt_sum += pt1
                pt_sum += pt2

            #print ptr

                # Evaluate synth likelihood
                ps1 = synth_likelihood(data, mu1, 1, dataset_samples=x[i], samples=datasets)
                ps2 = synth_likelihood(data, mu2, 1, dataset_samples=x[i], samples=datasets)

                sl_sum += ps1
                sl_sum += ps2

                pts = ps1 / ps2
            #print pts

                errors.append((pts - ptr)**2)
            print pt_sum
            print sl_sum
            y[i] = np.mean(errors)

        plt.plot(x, y, label='%i' % data_sample)

    plt.xlabel('Number of samples per dataset')
    plt.ylabel('MSE')
    plt.legend()

    def r_like(b1=beta1_summ, n=obs_summ['n']):
        """Likelihood for correlation coefficients of summarized data"""
        # Convert slope to r
        rho = b1*stdev_phe/stdev_iq

        # Fisher transformation to allow for normality assumption
        eps = np.arctan(rho) - np.arctan(obs_summ['correlation'])
        # Difference should be mean-zero
        return normal_like(eps, mu=np.zeros(len(n)), tau=n-3)

        def obs(value=logit_val, logit_se=logit_se,
                X=covariates(d),
                alpha=alpha, beta=beta, gamma=gamma, sigma=sigma,
                age_indices=age_indices,
                age_weights=age_weights):

            # calculate study-specific rate function
            mu = predict_logit_rate(X, alpha, beta, gamma)
            mu_i = rate_for_range(mu, age_indices, age_weights)
            
            tau_i = 1. / (sigma**2 + logit_se**2)
            logp = mc.normal_like(x=value, mu=mu_i, tau=tau_i)
            return logp

      def model(params=params, vars=vars, paramnames=paramnames, filters=filters,
            value=1.0):
         # Set the parameters in the model
         for i,param in enumerate(paramnames):
            if debug:
               print "setting ",param, " to ",params[i]
            self.model.parameters[param] = params[i]

         logp = 0
         numpts = 0
         for i,f in enumerate(filters):
            mod,err,mask = self.model(f, self.sn.data[f].MJD)
            m = mask*self.sn.data[f].mask
            if not np.sometrue(m):
               continue
            numpts += np.sum(m)
            tau = np.power(vars[i] + np.power(self.sn.data[f].e_mag,2),-1)
            logp += pymc.normal_like(self.sn.data[f].mag[m],mod[m],tau[m])
         #if numpts < len(paramnames):
         #   return -np.inf
         return logp

 def deriv_sign_rate(f=rate,
                     age_indices=age_indices,
                     tau=1.e14,
                     deriv=deriv, sign=sign):
     df = pl.diff(f[age_indices], deriv)
     return mc.normal_like(pl.absolute(df) * (sign * df < 0), 0., tau)