differentiation of arbitrary matrix expressions

On Saturday, June 29, 2013 09:06:58 AM Stavros Macrakis wrote:
> Thanks for your question.
> 
> Your code does seem to run into a Maxima bug; a simpler version of it is:
> 
> declare([X,Y],nonscalar);
> integrate(f(x),x,0,X);     <<< does not give an error (arguable, but OK)
> integrate(f(x),x,0,X.Y);  <<< gives an error that X.Y is non-real, but X.Y
> can in fact be real/scalar
> 
> There is a simple workaround for this:
> 
> P : tau^(-2) *Phi(Q)^y*(1-Phi(Q))^(1-y)*phi(theta/tau);
> subst(X.beta+theta,Q,P)
Given
phi(x) := exp(-x^2/2)/sqrt(2*%pi);
Phi(x) := integrate(phi(z), z, -inf, x);
I get
P : tau^(-2) *Phi(Q)^y*(1-Phi(Q))^(1-y)*phi(theta/tau);
Maxima encountered a Lisp error:
 Error in PROGN [or a callee]: Bind stack overflow.
Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.

Also, I'm not sure what I'm supposed to do with P; I'll want to do stuff like 
get the derivative with respect to beta.

> 
> Sorry for the inconvenience.
> 
> You can reorganize the result using factor.
> 
> Re: The matrix() function seems designed for matrices of known dimensions
> 
> matrix(...) is a notation for writing explicit matrices in terms of their
> entries, yes.  Not sure what your question is here.
The impliciit question was whether matrix would help with my problem.  I take 
it the answer is no.  I was concerned this was the only way to tell maxima it 
has a matrix; I gather the answer is no because of the nonscalar declarations.

Is there a way to say a matrix has dimension i x j where i and j are unknown 
but fixed values?  I presume it would need that information to decide whether, 
e.g., X.beta was a scalar.
> 
> Re:
> It would also be nice if diff(Phi(x), x) gave me phi(x) instead of the
> fully expanded solution.
> 
> Hmm, not sure what you want exactly here.  You could define
> 
> Phi(x) := integrate('phi(z),z,-inf,x)$
> 
> ... in which case phi(z) will remain unevaluated, and diff(Phi(x),x) will
> get you back phi.  Is that what you have in mind?
This is more of a cosmetic issue; I guess sometimes I want to see the phi 
literally, but others, e.g., taking derivatives, I may want to differentiate 
through it.

Another way of saying that is that it would be nice if the solutions arrived 
at were expressed in terms of phi if possible.
> 
> I think you'll have to tell us more about what sort of result you're trying
> to get for us to help you more.
Ultimately, I have a log likelihood that involves nested sums of terms 
involving matrices and vectors, like X[i, j].beta+theta[i] where X[i, j] is 
itself a row vector.  I want to differentiate with respect to, e.g., beta and 
theta_i.

I did it by hand, but the process was quite error-prone, and I'm probably 
going to need to do it again with even more complex expressions.

Ross

> 
>                 -s
> 
> On Fri, Jun 28, 2013 at 2:42 PM, Ross Boylan <ross at biostat.ucsf.edu> wrote:
> > 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
> > ______________________________**_________________
> > Maxima mailing list
> > Maxima at math.utexas.edu
> > http://www.math.utexas.edu/**mailman/listinfo/maxima<http://www.math.utex
> > as.edu/mailman/listinfo/maxima>