Symbolic Conditional Expectation in Mathematica

2019-04-15 04:20发布

I want to implement the conditional expectation operator. I will use the capital epsilon E to denote the operator. I expect at least the following input (underscore means subscript)

E_2[a]
E_2[x_1]
E_2[x_1 + y_5]
E_1[3 a + b - 4 + 2 x_]
E_6[x_5 x_7]
E_t[x_t]
E_t[3 x_{t - 1} x_{t + 2}]

to produce the following output

a
x_1
E_2[y_5] + x_1
-4 + 3 a + b + 2 E_2[x_5]
E_6[x_7] x_5
x_t
3 E_t[x_{t + 2}] x_{t - 1}

The examples above are not the only input/output pairs I need to produce, but rather serve as a test and an illustration for the syntax I prefer.

I've got this far. ce means Conditional Expectation, its third component is whether the "expectation propagation" is finalized or not (otherwise the infinite recursion occurs in the product rule), mv stands for Measurable Variable.

Notation[Subscript[E, t_][y_]  ==> ce[y_, t_, False]];
Notation[Subscript[E, t_][y_] <==  ce[y_, t_, _]];

Notation[Subscript[x_, t_] <==> mv[x_, t_]];

(* Atomic Elements and Measurable Variables *)
ce[x_, t_, _] := x /; (AtomQ[x] || Head[x] === mv && 0 <= t - x[[2]]);

(* Distribution over Addition *)
ce[x_ + y__, t_, s_] := ce[x, t, s] + Plus @@ (ce[#, t, s] & /@ {y});

(* Distribution over Product *)
ce[x__Times, t_, False] := Module[{v, m, n},

   (* All Variables in the Product *)
   v = List @@ x;

   (* Measurable Among Them *)
   m = Select[v, AtomQ[#] || Head[#] === mv && 0 <= t - #[[2]] &];

   (* The Rest is not Measurable *)
   n = Complement[v, m];

   Times[Times @@ m, ce[Times @@ n, t, True]]

];       

1条回答
2楼-- · 2019-04-15 05:09

I think I can get you close to what you want; I'm not going to do it all, though, as it can be tricky, but I'll point you in the right direction.

First of all, using subscripts to denote different variables is tricky in Mathematica as it interprets E0 as Subscript[E,0] and both E and Subscript are reserved. (As Sjoerd said, E = 2.718....) To get Mathematica to recognize <anything><something> as a distinct symbol, you need to load the Notations package via <<Notations`. Then using the Notations Palette, Symbolize Subscript[E,0]. (As a word of caution, don't try to do that without using the palette to set up the code correctly, otherwise it may not work.)

Once all of your variables are symbolized, as needed, you need to set up the appropriate transformation rules. The first two are simplest, enter

E_0[a] = a
E_0[x_0] = x_0

Rule 3 and 4:

E_0[x_Plus]:=Distribute[E_0[x]]
E_0[x_Times]:=Distribute[E_0[x], Times]

Those were the easy ones, the next three require a different kind of association, neither Set nor SetDelayed will work here as the outer symbol being evaluated is Dt, and you can't associate new rules with it as it is Protected. However, there are two methods of associating such expressions with an internal symbol: UpSet (^=) (or UpSetDelayed (^:=)) or TagSet (/:). I prefer to use TagSet as it is more explicit, but either should work.

Rule 5 and 6:

E_0 /: Dt[ E_0[ x_ ], y_ ] := E_0[ Dt[x,y] ]

This will also get you close to rule 7, but adding this alongside rules 3 and 4 causes a Recursion Limit error as it bounces back and forth trying to figure out how to evaluate it. Instead, replace rule 3 and 4 with

E_0[x_ + y__]:= E_0[x] + Plus@@( E_0 /@ {y} )
E_0[x_ y__ ] := E_0[x] Times@@( E_0 /@ {y} )

Which puts definite limits on the recursion. As far as rule 7 is concerned, you get this

E_0[D[x_1[t_1,q_0], t_1]] E_0[Dt[t_1, y_0]] 
+ E_0[D[x_1[t_1,q_0], q_0]] E_0[Dt[q_0,y]]

which is a consequence of the Dt rule and rule 4. To get E_0 not to distribute over D and Dt is left as an exercise.


Edit: I'd like to make a few comments on the solution code you've provided. First, clever use of a Boolean to stop the recursion, and it works well with your Notation. I'd suggest several changes to you product distribution, though. First, I'd use x__Times instead of the condition (/; Head[x] == Times) as it is easier to read, and I believe (but haven't tested) it may be faster, i.e. less overhead to process it. Second, replace you use of Table with List @@ x, where @@, called Apply, replaces Times with List, and it is again easier to read and write. For your definition of n, consider using Complement; I don't know if it is faster, but I tend to prefer set theoretic constructs for this type of thing. Lastly, unless you need a variable to be reevaluated whenever it is used, do not use SetDelayed (:=), use Set (=). By using :=, m is evaluated twice, and v is evaluated 3 times!

Pros and Cons: The chief reasons to do this is ease of use and readability. By defining your own objects and how they behave, you give yourself a lot of flexibility and simplify your code. That alone makes it worth it. However, I've had difficulty doing this in the past, and such a setup can be persnickety, and I'd recommend thorough testing. Secondly, by adding these extra layers, you may slow down your code, so be careful if this is to be used in mission critical applications. Additionally, you have to include Notation every time you use it, and the palette will become annoying at some point. Although, the palette can be dealt with by setting AutoLoadNotationPalette = False prior to loading the Notation package.

查看更多
登录 后发表回答