# How to implement a formal expectation operator over an unknown distribution

GROUPS:
 Sune Jespersen 1 Vote I need to do some simplification of an expression involving averages over a stochastic variable (in order to verify a long analytical calculation). The easiest way to do that, I figured, were if I could implement an operator which would basically be short-hand for the averaging procedure, with all the appropriate properties. Then of course this operator would be present in the final expression, which is fine, and would enable me to compare easily with my own calculations. So assuming I use x for the stochastic variable, I tried defining av usingav[y_ + z_] := av[y] + av[z]  av[c_ y_] := c av[y] /; FreeQ[c, x] av[c_] := c /; FreeQ[c, x]  Then when I writeD[av[x y], y]I get av[x]which is fine, but when I writeD[av[Exp[-x y]], y]I get-E^(-x y) y  instead of -y av[Exp[-x y]] as I want; i.e., the av is removed somehow.I tried using UpValue for teaching Mathematica that it could interchange differentiation and av, but apparently that is not the problem. I might be going about this entirely the wrong way, but I'd be grateful for any input. Note the builtin Expectation function does not accomplish it either.
4 years ago
8 Replies
 I think I got it:In[1]:= av /: D[av[f___], x_] := av[D[f, x]]In[2]:= av[y_ + z_] := av[y] + av[z]In[3]:= av[c_ y_] := c av[y] /; FreeQ[c, x]In[4]:= av[c_] := c /; FreeQ[c, x]In[5]:= D[av[x y], x]Out[5]= yIn[6]:= D[av[Exp[-x y]], x]Out[6]= -y av[E^(-x y)]Sune
4 years ago
 Well, turns out there's still a problem:In[9]:= D[Log[av[Exp[-b x]]], b]Out[9]= -((E^(-b x) x)/av[E^(-b x)])instead of-(av[(E^(-b x) x)]/av[E^(-b x)])What's going on?Sune
4 years ago
 Chip Hurst 2 Votes I have a guess as to what's going on.Your definition of derivative for av disobeys the chain rule. I think Mathematica is first applying the chain rule, then your rule.We can 'hack' our way around this (kind of ugly) Unprotect[D]; ClearAttributes[D, Protected]  (* avoid Mathematica's chain rule *) D[f_[HoldPattern[av][g_]], x_] := Module[{u},   av[D[g, x]] (D[f[u], u] /. u -> av[g]) ]  In[1]:= D[Log[av[Exp[-b x]]], b]Out[1]= -(av[E^(-b x) x]/av[E^(-b x)])
4 years ago
 Hi chip.Thanks for your reply. I think you're right that it has to do with the chain rule, but I don't exactly understand it. When applying the chain rule, it still needs to apply the derivative to the av operator, and then my rule would apply, wouldn't it?I guess your hack could be formulated as a upvalue also, which if so probably would be better practice.@Peter: I didn't want to introduce a new operator, because I would like to take advantage of functions such as Series, which call D.S
4 years ago
 When applying the chain rule, it still needs to apply the derivative to the av operator, and then my rule would apply, wouldn't it?When applying the chain rule, Mathematica would only need to evaluate D[av(u), u], which is 1 (df/dx = df/du du/dx).Now that I look at it again, I think both answers (yours and Mathematica's) are correct. In fact, can't we use both results to conclude av(x) is the identity function?Your way (special D rule): D[Log[av[f(x)]], x] == D[av[f(x)], x] / av[f(x)] == av[D[f(x), x]] / av[f(x)] == av[f'(x)] / av[f(x)].Mathematica (chain rule): D[Log[av[f(x)]], x] == D[av[f(x)], x] / av[f(x)] == (D[av(u), u] /. u -> f(x)) D[f(x), x] / av[f(x)] == 1 * f'(x) / av[f(x)].This means av[f'(x)] == f'(x), and letting f(x) == x^2/2 gives us av(x) == x?