Я пытаюсь реализовать байесовский PCA, используя библиотеку PyMC для python. Но я застрял там, где определяю координаты более низкого измерения...
Модель
x = Wz + e
где x — вектор наблюдения, W — матрица преобразования, а z — вектор координат меньшего размера.
Сначала я определяю распределение для матрицы преобразования W. Каждый столбец взят из нормального распределения (для простоты нулевое среднее значение и тождественная ковариация).
def W_logp(value):
logLikes = np.array([multivariate_normal.logpdf(value[:,i], mean=np.zeros(dimX), cov=1) for i in range(0, dimZ)])
return logLikes.sum()
def W_random():
W = np.zeros([dimX, dimZ])
for i in range(0, dimZ):
W[:,i] = multivariate_normal.rvs(mean=np.zeros(dimX), cov=1)
return W
w0 = np.random.randn(dimX, dimZ)
W = pymc.Stochastic(
logp = W_logp,
doc = 'Transformation',
name = 'W',
parents = {},
random = W_random,
trace = True,
value = w0,
dtype = float,
rseed = 116.,
observed = False,
cache_depth = 2,
plot = False,
verbose = 0)
Затем я хочу определить распределение для z, которое снова является многомерным нормальным (нулевое среднее и тождественная ковариация). Однако мне нужно нарисовать z для каждого наблюдения отдельно, а W общее для всех. Итак, я попытался
z = pymc.MvNormal('z', np.zeros(dimZ), np.eye(dimZ), size=N)
Однако pymc.MvNormal не имеет параметра size. Так что выдает ошибку. Следующим шагом будет
m = Data.mean(axis=0) + np.dot(W, z)
obs = pymc.MvNormal('Obs', m, C, value=Data, observed=True)
Я не приводил спецификацию для C выше, так как сейчас это не имеет значения. Есть идеи как реализовать?
Спасибо
РЕДАКТИРОВАТЬ
После ответа Криса Фоннесбека я изменил свой код следующим образом
numD, dimX = Data.shape
dimZ = 3
mm = Data.mean(axis=0)
tau = pymc.Gamma('tau', alpha=10, beta=2)
tauW = pymc.Gamma('tauW', alpha=20, beta=2, size=dimZ)
@pymc.deterministic(dtype=float)
def C(tau=tau):
return (tau)*np.eye(dimX)
@pymc.deterministic(dtype=float)
def CW(tau=tauW):
return np.diag(tau)
W = [pymc.MvNormal('W%i'%i, np.zeros(dimZ), CW) for i in range(dimX)]
z = [pymc.MvNormal('z%i'%i, np.zeros(dimZ), np.eye(dimZ)) for i in range(numD)]
mu = [pymc.Lambda('mu%i'%i, lambda W=W, z=z: mm + np.dot(np.array(W), np.array(z[i]))) for i in range(numD)]
obs = [pymc.MvNormal('Obs%i'%i, mu[i], C, value=Data[i,:], observed=True) for i in range(numD)]
model = pymc.Model([tau, tauW] + obs + W + z)
mcmc = pymc.MCMC(model)
Но на этот раз он пытается выделить огромное количество памяти (более 8 ГБ) при запуске pymc.MCMC(model)
с numD=45
и dimX=504
. Даже когда я пробую это только с numD=1
(то есть создаю только 1 z
, mu
и obs
), он делает то же самое. Есть идеи, почему?