Simplify expectation taking code

[sbenthall]

The code that gets the expected values of distributions in the solver can be simplified by making generic methods on the distribution classes.

https://github.com/econ-ark/HARK/blob/master/HARK/ConsumptionSaving/ConsAggShockModel.py#L633-L667

[sbenthall]

duplicate #621

[sbenthall]

Mathematica syntax is:
Expectation[f[x],x [Distributed] LogNormal[mu,sigma]]

[sbenthall]

This code is calculating the expectation of:

• a mathematical function: \psi ^ (1 - \rho) * vPfunc_Next ? (a function over all the values in the distribution)
• over the Income distribution.... IncomeDstn

[sbenthall]

Try to mirror the Sympy syntax:

https://docs.sympy.org/latest/modules/stats.html?highlight=expectation#sympy.stats.Expectation

[sbenthall]

https://gist.github.com/mnwhite/a77387cf2633f0a2404342d39fdb5284

[mnwhite]

I'm going to take a stab at this tonight or tomorrow morning. It's easier than I thought.

…

On Thu, Jun 18, 2020 at 1:06 PM Sebastian Benthall ***@***.***> wrote: https://gist.github.com/mnwhite/a77387cf2633f0a2404342d39fdb5284 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#625 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADKRAFPCBZHMOC5NUIQLNQDRXJCPZANCNFSM4MJX676A> .

[llorracc]

By "this" you mean working out sympy syntax for some expression(s) in your gist code?

[mnwhite]

I'm making a non-sympy expectations function. Here's the docstring, which won't format correctly here:

`
def calcExpectation(func,values,dstn):
'''
Calculate the expectation of a stochastic function at an array of values.

Parameters
----------
func : function
    The function to be evaluated, which can take M+N identically shaped arrays
    as arguments and returns 1 array as an output (of the same shape).  The
    first M arguments are non-stochastic, representing the inputs passed in
    the argument values.  The latter N arguments are stochastic, where N is
    the dimensionality of dstn.
values : np.array or [np.array]
    One or more arrays of input values for func, representing the non-stochastic
    arguments.  If the array(s) are 1D, they are interpreted as grids over
    each of the M non-stochastic arguments of the function; the expectation
    will be evaluated at the tensor product set, so the output will have shape
    (values[0].size,values[1].size,...,values[M].size).  Otherwise, the arrays
    in values must all have the same shape, and the expectation is computed
    at f(values[0],values[1],...,values[M],*dstn).
dstn : DiscreteDistribution
    The distribution over which the function is to be evaluated.  It should
    have dimension N, representing the last N arguments of func.
    
Returns
-------
f_exp : np.array
    The expectation of the function at the queried values.
'''
pass

`

[mnwhite]

See the calcExpectations branch for a completely untested draft of the function. It's what I was able to write in an hour before needing to jump on baby duty. Will test it on nonsense functions to make sure it works, then will pastebin an edited version of the ToddlerCSmodel.py thingy where it's implemented.

…

On Thu, Jun 18, 2020 at 5:07 PM Christopher Llorracc Carroll < ***@***.***> wrote: By "this" you mean working out sympy syntax for some expression(s) in your gist code? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#625 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADKRAFOT4HIAKIWXAEMER3LRXJ6ZHANCNFSM4MJX676A> .

[sbenthall]

Right on. I'll be glad to review this PR when it's ready.

In my view, it makes sense to make the problem easier by breaking it down.

Refactoring the solution code into methods/functions/objections that create an abstraction layer of mathematical operations over the lower level implementation (which for now will be matrix manipulations but later may be symbolic) is very difficult for me because I don't have the math in my head. It makes tons of sense for @mnwhite to do that, if he has time.

Massaging these functions into a smooth, Mathematica/Sympy inspired syntax is much easier (for me) once that first step is completed.

My only quibble with what Matt's proposed is a style issue that it pervasive across the whole library and perhaps can't be helped. But I should bring it up nonetheless. According to PEP 8, the recommendation is:

• Class names should normally use the CapWords convention.
• Function names should be lowercase, with words separated by underscores as necessary to improve readability. Variable names follow the same convention as function names. So lower_case_with_underscores

mixedCase is allowed for "backwards compatibility", which was a nice allowance given the variety of styles used in Python. But PEP 8 was published in 2001. So it's a little silly for a project launched almost 15 years later to be using the deprecated convention.

I realize that this is a sticky problem but I'd volunteer to do the mass conversion before the 1.0 release.

[mnwhite]

The capitalization style is intentional; to my recollection, Nate Palmer told us this was the standard for functions when we began the project. The lack of underscores in variable names has a different origin.

…

On Fri, Jun 19, 2020 at 10:00 AM Sebastian Benthall < ***@***.***> wrote: Right on. I'll be glad to review this PR when it's ready. In my view, it makes sense to make the problem easier by breaking it down. Refactoring the solution code into methods/functions/objections that create an abstraction layer of mathematical operations over the lower level implementation (which for now will be matrix manipulations but later may be symbolic) is very difficult *for me* because I don't have the math in my head. It makes tons of sense for @mnwhite <https://github.com/mnwhite> to do that, if he has time. Massaging these functions into a smooth, Mathematica/Sympy inspired syntax is much easier (for me) once that first step is completed. My only quibble with what Matt's proposed is a style issue that it pervasive across the whole library and perhaps can't be helped. But I should bring it up nonetheless. According to PEP 8 <https://www.python.org/dev/peps/pep-0008/>, the recommendation is: - Class names should normally use the CapWords convention. - Function names should be lowercase, with words separated by underscores as necessary to improve readability. Variable names follow the same convention as function names. So lower_case_with_underscores mixedCase is allowed for "backwards compatibility", which was a nice allowance given the variety of styles used in Python. But PEP 8 was published in 2001. So it's a little silly for a project launched almost 15 years later to be using the deprecated convention. I realize that this is a sticky problem but I'd volunteer to do the mass conversion before the 1.0 release. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#625 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADKRAFK22MJDWI2XHNAQIL3RXNVQFANCNFSM4MJX676A> .

[sbenthall]

@llorracc today in the meeting has pointed out that it would be nice if the object representing a random variable (for example, the distribution objects) were used both in simulation (when sampling) and in solution (by taking expectations).

[sbenthall]

#720 #884 now merged; expectation taking code is in.

Next step on this is refactoring to use calcExpectations in the solution code throughout the library

[sbenthall]

I think the current implementation of calcExpectation is quite complicated and can be further simplified using this method that comes with numpy

https://numpy.org/doc/stable/reference/generated/numpy.apply_along_axis.html

[sbenthall]

One aspect of the calcExpectation method is the broadcasting of the function over some values that are not part of the probability distribution.

This introduces a lot of complexity as the current method is trying to cover a couple different kinds of broadcasting without relying on the underlying numpy machinery.

One small issue resulting from this is that the documentation is ambiguous as to whether the 'default' value of these values is an empty list or None.

https://github.com/econ-ark/HARK/blob/master/HARK/distribution.py#L1020-L1032

Correspondingly, I've found a small inconvenience coming from whether the output of calcExpectation is an numpy array or a scalar value under conditions where there are no reported values.

Consider:

dd_0_1_20 = distribution.Normal().approx(20)
ce1 = distribution.calcExpectation(dd_0_1_20)

In this code ce1 is a numpy array containing one element, which is very close to zero. However, under these conditions, we might expect the output to be a scalar. Note that the output of numpy.dot is a scalar if both inputs are 1-D arrays.

https://numpy.org/doc/stable/reference/generated/numpy.dot.html

I think it will save the user a lot of headache if we think through how to implement distributions and expectation taking in terms of numpy array broadcasting.

[sbenthall]

The current calcExpectation code handles the case where non-stochastic values are handed in as a set of 1D arrays that are to be interpreted as grids over multiple dimensions. The code intends to handle this by first creating a completed M-dimensional grid which is the tensor product over each dimensional array.

The condition for using this tensor-product mechanism, as opposed to the alternative simpler case where the function simply broadcasts over the values array, is unclear. I don't think the current code/documentation successfully distinguished between
these two cases, as in principle one could broadcast over several 1D arrays.

https://github.com/econ-ark/HARK/blob/master/HARK/distribution.py#L1047-L1048

Instead of combining these two complex operations into a single function, it would be better to break them down into simpler steps.

Creating the tensor product of several grids can be its own function and live elsewhere. That way the expectation-taking code can be focused on doing one thing well.

I'll make this change in a PR.

[sbenthall]

I'm not sure, but I believe that the intention of the grid 'tensor product' work can be done with this built-in numpy function for the Kronecker product:

https://numpy.org/doc/stable/reference/generated/numpy.kron.html#numpy.kron

Or possibly it can be done using this numpy.tensordot function

https://numpy.org/doc/stable/reference/generated/numpy.tensordot.html?highlight=tensor%20product

Because this is a very simple operation, and one that could easily be done prior to the input to calcExpectation, I won't try to preserve the code for the tensor product case in the next PR. Likely the current way this is being done, with np.tile, is inefficient and it is certainly less clear what is going on to the reader.

As calcExpectation has not yet been worked into the code base at large, this reduction in the method's functionality will have no effect on downstream use. In the upcoming, PR, I'll focus on trimming the functionality down and providing thorough tests.

[sbenthall]

Ok, I've tested it and numpy.apply_along_axis automatically broadcasts over additional values args if they are arrays.

https://numpy.org/doc/stable/reference/generated/numpy.apply_along_axis.html#numpy.apply_along_axis

That makes the implementation of calcExpecation very easy.

In general, HARK should be refactored to use the core libraries effectively as much as possible.

[sbenthall]

It's actually quite wonderful. args can be a mix of constants and same-sized arrays and it will act very reasonably.

[Mv77]

@sbenthall had mentioned he'd be keen to hear feedback on the new calcExpectation from other users.

I have been incorporating it into my new portfolio model today. It works great: intuitive to use, surprisingly faster, and makes for clearer code, pushing HARK towards explicit definitions of transition equations.

[sbenthall]

oh fantastic. thanks!

[sbenthall]

Maybe the right numpy tool for making the grids over several scalar function inputs is np.meshgrid

https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html

[Mv77]

I think that's it but I'm not sure what it would be needed for in the current implementation given that the multivariate distribution already has that grid for shocks.

Maybe it could be used for args, so that they aren't required to be of the same shape and one can just pass the 1-d grids instead of n-d meshes. But this depends on whether you want the user or the function to take care of args. It's not obvious what the best practice would be there, since sometimes the meshes are useful beyond taking expectations for, e.g., evaluating value functions and so the user might have them laying around already.

[Mv77]

Still, I believe there are various places in HARK where multiple lines of np.tile(...) could be replaced with single lines of np.meshgrid.

[sbenthall]

There should be a separate issue for this.

…

On Wed, Jan 27, 2021, 11:48 AM Mateo Velásquez-Giraldo < ***@***.***> wrote: Still, I believe there are various places in HARK where multiple lines of np.tile(...) could be replaced with single lines of np.meshgrid. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#625 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAAQZEEUSNHEFLEBW64SBXLS4A7U3ANCNFSM4MJX676A> .

[sbenthall]

Closing this issue. The new calcExpectations function is used in several places in the library. There are further code simplifications along these lines to be made; see #934 for the np.meshgrid issue.