How to use the pymc.potential function in pymc

    @mc.potential
    def obs(pi=pi, delta=delta):
        return mc.negative_binomial_like(r*n, pi*n, delta)

    @mc.potential(name='mu_age_derivative_potential_%s'%name)
    def mu_age_derivative_potential(mu_age=mu_age,
                                    increasing_a0=pl.clip(parameters['increasing']['age_start']-ages[0], 0, len(ages)),
                                    increasing_a1=pl.clip(parameters['increasing']['age_end']-ages[0], 0, len(ages)),
                                    decreasing_a0=pl.clip(parameters['decreasing']['age_start']-ages[0], 0, len(ages)),
                                    decreasing_a1=pl.clip(parameters['decreasing']['age_end']-ages[0], 0, len(ages))):
        mu_prime = pl.diff(mu_age)
        inc_violation = mu_prime[increasing_a0:increasing_a1].clip(-pl.inf, 0.).sum()
        dec_violation = mu_prime[decreasing_a0:decreasing_a1].clip(0., pl.inf).sum()
        return -1.e12 * (inc_violation**2 + dec_violation**2)

        @mc.potential(name='deriv_sign_{%d,%d,%d,%d}^%s' % (deriv, sign, age_start, age_end, rate))
        def deriv_sign_rate(f=rate,
                            age_indices=age_indices,
                            tau=10000.,
                            deriv=deriv, sign=sign):
            df = np.diff(f[age_indices], deriv)
            return -tau * np.dot(df**2, (sign * df < 0))
        return [deriv_sign_rate]

        @mc.potential(name='deriv_sign_{%d,%d,%d,%d}^%s' % (deriv, sign, age_start, age_end, str(rate)))
        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)
        return [deriv_sign_rate]

        @mc.potential(name='hierarchical_potential_%s'%stoch_key)
        def hier_potential(r1=vars[stoch_key]['rate_stoch'], r2=world_rate,
                           c1=vars[stoch_key]['conf'], c2=world_confidence):
            return mc.normal_like(np.diff(r1) - np.diff(r2), 0., c1 + c2)
        vars[stoch_key]['h_potential'] = hier_potential

            @mc.potential(name='zero-%d-%d^%d' % (age_start, age_end, rf.id))
            def zero_rate(f=rf.vars['Erf_%d'%rf.id], age_start=age_start, age_end=age_end, tau=1./(1e-4)**2):
                return mc.normal_like(f[range(age_start, age_end+1)], 0.0, tau)
            rf.vars['prior'] += [zero_rate]

            @mc.potential(name='unimodal_{%d,%d}^%s' % (age_start, age_end, rate))
            def unimodal_rate(f=rate, age_indices=age_indices, tau=1000.):
                df = np.diff(f[age_indices])
                sign_changes = pl.find((df[:-1] > NEARLY_ZERO) & (df[1:] < -NEARLY_ZERO))
                sign = np.ones(len(age_indices)-2)
                if len(sign_changes) > 0:
                    change_age = sign_changes[len(sign_changes)/2]
                    sign[change_age:] = -1.
                return -tau*np.dot(np.abs(df[:-1]), (sign * df[:-1] < 0))
            priors += [unimodal_rate]

            @mc.potential(name='unimodal_{%d,%d}^%s' % (age_start, age_end, str(rate)))
            def unimodal_rate(f=rate, age_indices=age_indices, tau=1.e5):
                df = pl.diff(f[age_indices])
                sign_changes = pl.find((df[:-1] > NEARLY_ZERO) & (df[1:] < -NEARLY_ZERO))
                sign = pl.ones(len(age_indices)-2)
                if len(sign_changes) > 0:
                    change_age = sign_changes[len(sign_changes)/2]
                    sign[change_age:] = -1.
                return -tau*pl.dot(pl.absolute(df[:-1]), (sign * df[:-1] < 0))
            priors += [unimodal_rate]

        @mc.potential(name='age_coeffs_potential_%s' % key)
        def gamma_potential(gamma=gamma, mu_gamma=mu_gamma, tau_gamma=1./sigma_gamma[param_mesh]**2, param_mesh=param_mesh):
            return mc.normal_like(gamma[param_mesh], mu_gamma[param_mesh], tau_gamma)

# TODO: Is there a security issue with passing the output of str(x) to exec?
    function_namespace.update([(str(param), param) for param in params])
    param_kwarg_names = ','.join([str(param) + '=' + str(param) for param in params])
    param_arg_names = ','.join([str(param) for param in params])

    function_namespace.update({'dataset_error_funcs': dataset_error_funcs})
    # Now we have to do some metaprogramming to get the variable names to bind properly
    # This code doesn't yet allow custom distributions for the error
    error_func_code = """def error({0}):
    result = zeros_like(dataset_names, dtype='float')
    for idx in range(len(dataset_names)):
        result[idx] = divide(square(dataset_error_funcs[idx]({1})).mean(), dataset_est_variances[idx])
    return -result.sum()""".format(param_kwarg_names, param_arg_names)
    error_func_code = compile(error_func_code, '', 'exec')
    exec(error_func_code, function_namespace)
    error = pymc.potential(function_namespace['error'])
    mod = pymc.Model([function_namespace[str(param)] for param in params] + [error])
    return mod, datasets