
Re: defining averages over unknown PDF
Posted:
Nov 6, 2013 12:26 AM


On Mon, Nov 04, 2013 at 11:16:59PM 0500, Sune wrote: > Dear all. > > I want to do some symbolic manipulations of an expression involving averages over a stochastic variable with an unknown density. Therefore, I figured I could define a function av, which would correspond to the average over this unknown parameter density function. > I did as follows: > av[y_ + z_, x_] := av[y, x] + av[z, x]? > av[c_ y_, x_] := c av[y, x] /; FreeQ[c, x] > av[c_, x_] := c /; FreeQ[c, x] > > So these are basic properties of the average over the distribution of X. Some things work okay, for example > In[52]:= av[Exp[x y], x]? > Out[52]= av[E^(x y), x] > and > In[79]:= D[av[x y, x], x]? > Out[79]= y > and > In[80]:= D[av[x y, x], y]? > Out[80]= av[x, x]. > > However, the most vital part for my calculations does not work: > In[81]:= D[av[Exp[x y], x], y]? > Out[81]= E^(x y) x > > It should have been av[Exp[x y] x,x]. > > Any clues to what I'm doing wrong? I'm thinking that I need to specify some rules for differentiation, but I don't know how. But then I'm wondering how come it got the other expressions for differentiation right.
Ahh, the subtle treacheries of partial differentiation. Note that by your definition,
In[71]:= av[Exp[x y] + h, x]  av[Exp[x y], x]
Out[71]= h
So that
In[72]:= Limit[(av[Exp[x y] + h, x]  av[Exp[x y], x])/h, h > 0]
Out[72]= 1
So both your "correct" and "incorrect" answers are consistent with the chain rule and and the above computation of partial derivatives. So why is D computing the partial derivative in such a stupid way? Well, it isn't, at least not directly. D correctly computes the partial derivative as
f'[x] * Derivative[1, 0][av][f[x], x] + Derivative[0, 1][av][f[x],x]
But now Derivative helpfully tries compute these partials using pure functions, and then your definitions kick in, giving 1 and 0 for the partials. In particular, your third definitions means av[#1,#2]& === #1, and you're doomed.
So you want to abort the automatic differentiation rules with your own custom rule, which you can do with the following syntax:
av /: D[av[f_, x_], y_] /; x =!= y := av[D[f, y], x]
In[65]:= D[av[Exp[x y], x], y]
Out[65]= av[E^(x y) x, x]
 Itai Seggev Mathematica Algorithms R&D 2173980700

