Factor Analysis with PyMC, due to John Salvatier, https://groups.google.com/forum/?fromgroups#!searchin/pymc/meanReturns/pymc/I6sot0JbVo4/nuI1BpJt4LkJ

fa.py

from numpy import *
import pymc
from scipy import stats
from scipy.stats import distributions as d

#parameters about the da
dimensions = 5
observations = 30
shape = (dimensions, observations)


# first generate some fake data
#generate a factor (just one)
trueFactorMagnitudes = d.norm(  loc = 0, scale = 1).rvs(observations)
# make up a factor loading
trueFactorLoadings = d.norm( loc = 1, scale = .2).rvs (5)
#make up some error scales
trueErrorSds = d.gamma.rvs(5, scale = .05, size = dimensions)

#make up the actual data
data = trueFactorMagnitudes[newaxis, :] * trueFactorLoadings[:, newaxis] + d.norm(loc = 0, scale = trueErrorSds[:, newaxis]).rvs(shape)


factors = pymc.Normal("factormagnitudes",mu = zeros(observations), tau = ones(observations), value = mean(data, axis = 0))
loadings = pymc.Normal("factorloadings",mu = zeros(dimensions), tau = ones(dimensions)*.001,  value = array([1,1,1,1,1]))
returnSDs = pymc.Gamma("residualsds", alpha = ones(dimensions) * 5 , beta = ones(dimensions) * .05, value = array([1,1,1,1,1]))

@pymc.deterministic
def returnPrecisions ( stdev = returnSDs):
    precisions = (ones(shape) * (stdev**-2)[:, newaxis]).ravel()

    return precisions


@pymc.deterministic
def meanReturns (factors = factors, loadings = loadings):

    means = factors[newaxis, :] * loadings[:,newaxis]

    return means.ravel()

returns = pymc.Normal ("returns", mu = meanReturns, tau = returnPrecisions, observed = True, value = data)

if __name__=="__main__":
    import fa
    from pymc import MCMC
    M = MCMC(fa)