How to use the pymc.rnormal function in pymc

if c not in area_list:
                    continue

                alpha[c] = mc.rnormal(0., sigma_true[3]**-2.)
                sum_c += alpha[c]
                last_c = c
            if last_c >= 0:
                alpha[last_c] -= sum_c

            alpha[r] = mc.rnormal(0., sigma_true[2]**-2.)
            sum_r += alpha[r]
            last_r = r
        if last_r >= 0:
            alpha[last_r] -= sum_r

        alpha[sr] = mc.rnormal(0., sigma_true[1]**-2.)
        sum_sr += alpha[sr]
        last_sr = sr
    if last_sr >= 0:
        alpha[last_sr] -= sum_sr

    return alpha

def test_age_pattern_model_sim():
    # simulate normal data
    a = pl.arange(0, 100, 5)
    pi_true = .0001 * (a * (100. - a) + 100.)
    sigma_true = .025*pl.ones_like(pi_true)

    p = pl.maximum(0., mc.rnormal(pi_true, 1./sigma_true**2.))

    # create model and priors
    vars = {}

    vars.update(age_pattern.age_pattern('test', ages=pl.arange(101), knots=pl.arange(0,101,5), smoothing=.1))

    vars['pi'] = mc.Lambda('pi', lambda mu=vars['mu_age'], a=a: mu[a])
    vars.update(rate_model.normal_model('test', vars['pi'], 0., p, sigma_true))

    # fit model
    m = mc.MCMC(vars)
    m.sample(2)

vars['mu_age'].trace() * fe_usa_1990 * re_usa_1990)


    ### Prediction case 3: random effect not constant, zero fixed effect coefficients

    # set random seed to make randomness reproducible
    pl.np.random.seed(12345)
    pred = covariate_model.predict_for(d, d.parameters['p'],
                                         'NAHI', 'male', 2005,
                                         'CAN', 'male', 1990,
                                         0., vars, 0., pl.inf)

    # test that the predicted value is as expected
    pl.np.random.seed(12345)
    fe = pl.exp(0.)
    re = pl.exp(mc.rnormal(0., vars['sigma_alpha'][3].trace()**-2))
    assert_almost_equal(pred.mean(0),
                        (vars['mu_age'].trace().T * fe * re).T.mean(0))

def test_log_normal_model_sim(N=16):
    # simulate negative binomial data
    pi_true = 2.
    sigma_true = .1

    n = pl.array(pl.exp(mc.rnormal(10, 1**-2, size=N)), dtype=int)
    p = pl.exp(mc.rnormal(pl.log(pi_true), 1./(sigma_true**2 + 1./n), size=N))

    # create model and priors
    vars = dict(mu_age=mc.Uniform('mu_age', 0., 1000., value=.01),
                sigma=mc.Uniform('sigma', 0., 10000., value=1000.))
    vars['mu_interval'] = mc.Lambda('mu_interval', lambda mu=vars['mu_age']: mu*pl.ones(N))
    vars.update(rate_model.log_normal_model('sim', vars['mu_interval'], vars['sigma'], p, 1./pl.sqrt(n)))

    # fit model
    m = mc.MCMC(vars)
    m.sample(1)

def simulated_age_intervals(data_type, n, a, pi_age_true, sigma_true):
    # choose age intervals to measure
    age_start = pl.array(mc.runiform(0, 100, n), dtype=int)
    age_start.sort()  # sort to make it easy to discard the edges when testing
    age_end = pl.array(mc.runiform(age_start+1, pl.minimum(age_start+10,100)), dtype=int)

    # find truth for the integral across the age intervals
    import scipy.integrate
    pi_interval_true = [scipy.integrate.trapz(pi_age_true[a_0i:(a_1i+1)]) / (a_1i - a_0i) 
                        for a_0i, a_1i in zip(age_start, age_end)]

    # generate covariates that add explained variation
    X = mc.rnormal(0., 1.**2, size=(n,3))
    beta_true = [-.1, .1, .2]
    beta_true = [0, 0, 0]
    Y_true = pl.dot(X, beta_true)

    # calculate the true value of the rate in each interval
    pi_true = pi_interval_true*pl.exp(Y_true)

    # simulate the noisy measurement of the rate in each interval
    p = pl.maximum(0., mc.rnormal(pi_true, 1./sigma_true**2.))

    # store the simulated data in a pandas DataFrame
    data = pandas.DataFrame(dict(value=p, age_start=age_start, age_end=age_end,
                                 x_0=X[:,0], x_1=X[:,1], x_2=X[:,2]))
    data['effective_sample_size'] = pl.maximum(p*(1-p)/sigma_true**2, 1.)

    data['standard_error'] = pl.nan

xmax = .07

replicates = 1000

residuals = [[], []]
coverage = [[], []]

### @export 'neg-binom-sim-study'
pi_true = .025
delta_true = 5.
n_pred = 1.e9

for i in range(replicates):
    print '\nsimulation replicate %d' % i
    ## generate simulated data
    n = pl.array(pl.exp(mc.rnormal(10, 1**-2, size=16)), dtype=int)
    k = pl.array(mc.rnegative_binomial(n*pi_true, delta_true), dtype=float)
    r = k/n


    ## setup negative binomial model
    pi = mc.Uniform('pi', lower=0, upper=1, value=.5)
    delta = mc.Uninformative('delta', value=100.)

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

    @mc.deterministic
    def pred(pi=pi, delta=delta):
        return mc.rnegative_binomial(pi*n_pred, delta) / float(n_pred)

reload(dismod3)

colors = ['#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f0', '#ffff33']

### @export 'initialize'
df = pandas.read_csv('ssas_mx.csv', index_col=None)

ages = pl.arange(101)
knots = [0, 15, 60, 100]
import scipy.interpolate
Y_true = pl.exp(scipy.interpolate.interp1d(knots, pl.log([1.2, .3, .6, 1.5]), kind='linear')(ages))

N = 50
tau = .1**-2
X = pl.array(mc.runiform(pl.arange(0., 100., 100./N), 100./N + pl.arange(0., 100., 100./N), size=N), dtype=int)
Y = mc.rnormal(Y_true[X], tau)

### @export 'initial-rates'
pl.figure(figsize=(17., 11), dpi=72)

dismod3.graphics.plot_data_bars(df, 'talk')
pl.semilogy([0], [.1], '-')

pl.title('All-cause mortality rate\nin 1990 for females\nin sub-Saharan Africa, Southern.', size=55)
pl.ylabel('Rate (Per PY)', size=48)
pl.xlabel('Age (Years)', size=48)

pl.subplots_adjust(.1, .175, .98, .7)
pl.axis([-5, 105, 2.e-4, .8])
pl.xticks(size=30)
pl.yticks(size=30)
pl.grid()

pl.plot(r, n, 'o',
            mew=0, alpha=.25)

    decorate_funnel()
    pl.title(r'$\alpha=%d$, $\beta=%d$' % (alpha, beta))

pl.figure(**book_graphics.half_page_params)
pl.subplots_adjust(wspace=.4)
pl.subplot(2,2,1)
plot_beta_binomial_funnel(20., 980.)

pl.subplot(2,2,2)
plot_beta_binomial_funnel(2000., 98000.)

pl.subplot(2,1,2)
n = pl.exp(mc.rnormal(log_ppc_n, 1**-2, size=16))
n = pl.array(n.round()+1, dtype=int)
k = pl.minimum(n, mc.rnegative_binomial(.02*n, ppc_example_dispersion))
k = pl.array(k, dtype=float)
r = k/n

alpha = mc.Uninformative('alpha', value=1.)
beta = mc.Uninformative('beta', value=99.)
pi = mc.Beta('pi', alpha, beta, value=.01*pl.ones(16))
@mc.potential
def obs(pi=pi):
    return mc.binomial_like(k, n, pi)
@mc.deterministic
def pred(pi=pi, alpha=alpha, beta=beta):
    return [mc.rbinomial(n, pi), mc.rbinomial(n, mc.rbeta(alpha, beta))]

mcmc = mc.MCMC([alpha, beta, pi, obs, pred])

def plot_beta_binomial_funnel(alpha, beta):
    pi_true = alpha/(alpha+beta)
    pi = mc.rbeta(alpha, beta, size=10000)

    n = pl.exp(mc.rnormal(10, 2**-2, size=10000))
    k = mc.rbinomial(pl.array(n, dtype=int), pi)
    r = k/n
    pl.vlines([pi_true], .1*n.min(), 10*n.max(),
              linewidth=2, linestyle='-', color='w', zorder=9)
    pl.vlines([pi_true], .1*n.min(), 10*n.max(),
              linewidth=1, linestyle='--', color='black', zorder=10)
    pl.plot(r, n, 'ko',
            mew=0, alpha=.25)

    pl.semilogy(schiz['r'], schiz['n'], 'ks', mew=1, mec='white', ms=4,
                label='Observed Values')

    pl.xlabel('Rate (Per PY)')
    pl.ylabel('Study Size (PY)')
    pl.axis([-.0001, .0101, 50., 1500000])
    pl.title(r'$\alpha=%d$, $\beta=%d$' % (alpha, beta))

    @mc.deterministic
    def predicted(mu=mu, sigma_explained=sigma_explained, sigma_e=sigma_e):
        return mc.rnormal(mu, 1 / (sigma_explained**2. + sigma_e**2.))