Ah, then I suspect that it may be the list that is slowing things up here. In general, it would be more efficient to keep everything as tensors. So if you have a set of N columns in your data (each representing a var), then you can define N parameters like this (instead of the loop you are currently using):
coefs = pm.Normal('coefs', mu=0, sigma=`, shape=N)
Or you can use the new dims/coords feature, which would allow you to “name” each parameter as you like. Either way, your sum should be something like this:
sums = pm.math.dot(coefs, data_columns.T)