Notebook on nbviewer (original) (raw)

y = df.Salary

yMean = y.mean() ySD = y.std() agammaShRa = gammaShRaFromModeSD(ySD/2 , 2*ySD )

x1 = df.Pos.cat.codes.values Nx1Lvl = len(df.Pos.cat.categories)

x2 = df.Org.cat.codes.values Nx2Lvl = len(df.Org.cat.categories)

cellSDs = df.groupby(['Org','Pos']).std().dropna() medianCellSD = cellSDs.median().squeeze() sdCellSD = cellSDs.std().squeeze() sgammaShRa = gammaShRaFromModeSD(medianCellSD, 2*sdCellSD)

with pm.Model() as model2:

#a0 = pm.Normal('a0', yMean, tau=1/(ySD*5)**2)
a0_tilde = pm.Normal('a0_tilde', mu=0, sd=1)
a0 = pm.Deterministic('a0', yMean + ySD*5*a0_tilde)
    
a1SD = pm.Gamma('a1SD', agammaShRa[0], agammaShRa[1])
#a1 = pm.Normal('a1', 0.0, tau=1/a1SD**2, shape=Nx1Lvl)
a1_tilde = pm.Normal('a1_tilde', mu=0, sd=1, shape=Nx1Lvl)
a1 = pm.Deterministic('a1', 0.0 + a1SD*a1_tilde)

a2SD = pm.Gamma('a2SD', agammaShRa[0], agammaShRa[1])
#a2 = pm.Normal('a2', 0.0, tau=1/a2SD**2, shape=Nx2Lvl)
a2_tilde = pm.Normal('a2_tilde', mu=0, sd=1, shape=Nx2Lvl)
a2 = pm.Deterministic('a2', 0.0 + a2SD*a2_tilde)
    
a1a2SD = pm.Gamma('a1a2SD', agammaShRa[0], agammaShRa[1])
#a1a2 = pm.Normal('a1a2', 0.0, 1/a1a2SD**2, shape=(Nx1Lvl, Nx2Lvl))
a1a2_tilde = pm.Normal('a1a2_tilde', mu=0, sd=1, shape=(Nx1Lvl, Nx2Lvl))
a1a2 = pm.Deterministic('a1a2', 0.0 + a1a2SD*a1a2_tilde)

mu = a0 + a1[x1] + a2[x2] + a1a2[x1, x2]
         
sigmaMode = pm.Gamma('sigmaMode', sgammaShRa[0], sgammaShRa[1]) 
sigmaSD = pm.Gamma('sigmaSD', sgammaShRa[0], sgammaShRa[1]) 
sigmaRa = pm.Deterministic('sigmaRa', (sigmaMode + tt.sqrt(sigmaMode**2 + 4*sigmaSD**2)) / (2*sigmaSD**2))
sigmaSh = pm.Deterministic('sigmaSh', (sigmaMode * sigmaRa))
sigma = pm.Gamma('sigma', sigmaSh, sigmaRa, shape=(Nx1Lvl, Nx2Lvl))
ySigma = pm.Deterministic('ySigma', tt.maximum(sigma, medianCellSD/1000))
    
nu = pm.Exponential('nu', 1/30.)
    
like = pm.StudentT('like', nu, mu, lam=1/(ySigma[x1, x2])**2, observed=y)