Modelling Multivariate Normal RV, with each row in data is sum of these RVs at different number of times

Questions
1

Hi,
I have simulated data using the following code

# Actual mean and covariance of multivariate normal
μ_actual = np.array([20.0, 4.0])
sigmas_actual = np.array([10, 1.5])
Rho_actual = np.matrix([[1.0, 0.4], [0.4, 1.0]])
cov_actual = np.diag(sigmas_actual) * Rho_actual * np.diag(sigmas_actual)
# simulating data
num_days = [30, 66, 15, 90, 30, 25, 32, 36, 53, 67, 23]
readings_sim = np.zeros((len(num_days), 2))
for i, days in enumerate(num_days):
    ew = np.random.multivariate_normal(μ_actual, cov_actual, size=days).sum(0)
    readings_sim[i] = ew

This is how the data looks like
image

if we sum these Normal RVs to n times. The new RV will be a normal distribution with
new mean = nmean and
new cov= nactual cov

keeping that in mind I have coded the following model

with pm.Model() as model1:
  mu1 = pm.Normal('mu1',0, 10)
  mu2 = pm.Normal('mu2', 0, 10)
  mu = tt.stack([mu1, mu2])
  sd_dist = pm.HalfNormal.dist(10)
  chol, corr, std = pm.LKJCholeskyCov("chol", 
                                      eta=1, n=2, 
                                      sd_dist=sd_dist, 
                                      compute_corr=True
                                      )
  cov_mat = pm.Deterministic("cov",chol.dot(chol.T))
  for i,n in enumerate(num_days):
    obs = pm.MvNormal("obs_"+str(i), mu = n* mu, cov= n * cov_mat, observed=readings_sim[i])
 
  trace = pm.sample(tune=500,target_accept=0.9)

The way I have handled the observed value does not look pythonic at all. This data is actually a simulation of water and electricity consumption data. later I need to run this at a customer level. therefore looking for a better solution.
thanks in advance