differentiation of arbitrary matrix expressions

I'm interested in working with things I know are matrices and vectors, 
but whose exact dimensions I don't want to specify.  I want to 
differentiate these things.  I have expressions like (using pseudo 
tex-ish notation) 1/tau^2 prod_i prod_j Phi(X_{ij} beta+theta_i)^y_{ij} 
phi(theta_i/tau)
where X_{ij} is itself a matrix (a row vector for the j'th partner of 
respondent i).
Phi is the cumulative standard normal distribution and phi is its density.

I tried to start simple, without the ij.  So now
phi(x) := exp(-x^2/2)/sqrt(2*%pi);
Phi(x) := integrate(phi(z), z, -inf, x);
declare(X, nonscalar)
declare(beta, nonscalar)
P(beta, theta, tau, y) := tau^(-2) *Phi(X . beta+theta)^y*(1-Phi(X . 
beta + theta))^(1-y)*phi(theta/tau);
log(P(beta, theta, tau, y))
#0: Phi(x=X . beta+theta)
#1: P(beta=beta,theta=theta,tau=tau,y=y)
  -- an error. To debug this try: debugmode(true);

Phi(X . beta + theta);
#0: Phi(x=X . beta+theta)
  -- an error. To debug this try: debugmode(true);

Phi(X*beta); --works

The matrix() function seems designed for matrices of known dimensions

It would also be nice if diff(Phi(x), x) gave me phi(x) instead of the 
fully expanded solution.

BTW,
declare(X[i,j], nonscalar);  --error

Thanks for any help.  cc's appreciated, as I'm having mail "issues".

Ross Boylan