combineIndepDstns handles multivatriate dstns
HARKutilities.combineIndepDstns previously could not handle multivariate
distributions as inputs.  Now it can, which will simplify construction
of income process in agg shocks models.
@@ -739,14 +739,16 @@ def combineIndepDstns(*distributions):
'''
Given n lists (or tuples) whose elements represent n independent, discrete
probability spaces (probabilities and values), construct a joint pmf over
- all combinations of these independent points.
+ all combinations of these independent points. Can take multivariate discrete
+ distributions as inputs.
Parameters
----------
distributions : [np.array]
Arbitrary number of distributions (pmfs). Each pmf is a list or tuple.
- For each pmf, the first vector is probabilities and the second is values.
- For each pmf, this should be true: len(X_pmf[0]) = len(X_pmf[1])
+ For each pmf, the first vector is probabilities and all subsequent vectors
+ are values. For each pmf, this should be true:
+ len(X_pmf[0]) == len(X_pmf[j]) for j in range(1,len(distributions))
Returns
-------
@@ -759,7 +761,7 @@ def combineIndepDstns(*distributions):
Discrete points for the joint discrete probability mass function.
Written by Nathan Palmer
- Latest update: 31 August 2015 by David Low
+ Latest update: 5 July August 2017 by Matthew N White
'''
# Very quick and incomplete parameter check:
for dist in distributions:
@@ -767,8 +769,10 @@ def combineIndepDstns(*distributions):
# Get information on the distributions
dist_lengths = ()
+ dist_dims = ()
for dist in distributions:
dist_lengths += (len(dist[0]),)
+ dist_dims += (len(dist)-1,)
number_of_distributions = len(distributions)
# Initialize lists we will use
@@ -788,23 +792,22 @@ def combineIndepDstns(*distributions):
# Now we are ready to tile.
# We don't use the np.meshgrid commands, because they do not
# easily support non-symmetric grids.
- Xmesh = np.tile(dist[1].reshape(dist_newshape),dist_tiles)
- Pmesh = np.tile(dist[0].reshape(dist_newshape),dist_tiles)
- # Now flatten the tiled arrays.
- flatX = Xmesh.ravel()
- flatP = Pmesh.ravel()
+ # First deal with probabilities
+ Pmesh = np.tile(dist[0].reshape(dist_newshape),dist_tiles) # Tiling
+ flatP = Pmesh.ravel() # Flatten the tiled arrays
+ P_temp += [flatP,] #Add the flattened arrays to the output lists
- # Add the flattened arrays to the output lists.
+ # Then loop through each value variable
+ for n in range(1,dist_dims[dd]+1):
+ Xmesh = np.tile(dist[n].reshape(dist_newshape),dist_tiles)
+ flatX = Xmesh.ravel()
X_out += [flatX,]
- P_temp += [flatP,]
# We're done getting the flattened X_out arrays we wanted.
# However, we have a bunch of flattened P_temp arrays, and just want one
# probability array. So get the probability array, P_out, here.
- P_out = np.ones_like(X_out[0])
- for pp in P_temp:
- P_out *= pp
+ P_out = np.prod(np.array(P_temp),axis=0)
assert np.isclose(np.sum(P_out),1),'Probabilities do not sum to 1!'
return [P_out,] + X_out
