code: def get_expected_states(qs, B, policy):
      """
      Compute the expected states under a policy, also known as the posterior predictive density over states
  
      Parameters
      ----------
      qs: ``numpy.ndarray`` of dtype object
          Marginal posterior beliefs over hidden states at a given timepoint.
      B: ``numpy.ndarray`` of dtype object
          Dynamics likelihood mapping or 'transition model', mapping from hidden states at ``t`` to hidden states at ``t+1``, given some control state ``u``.
          Each element ``B[f]`` of this object array stores a 3-D tensor for hidden state factor ``f``, whose entries ``B[f][s, v, u]`` store the probability
          of hidden state level ``s`` at the current time, given hidden state level ``v`` and action ``u`` at the previous time.
      policy: 2D ``numpy.ndarray``
          Array that stores actions entailed by a policy over time. Shape is ``(num_timesteps, num_factors)`` where ``num_timesteps`` is the temporal
          depth of the policy and ``num_factors`` is the number of control factors.
  
      Returns
      -------
      qs_pi: ``list`` of ``numpy.ndarray`` of dtype object
          Predictive posterior beliefs over hidden states expected under the policy, where ``qs_pi[t]`` stores the beliefs about
          hidden states expected under the policy at time ``t``
      """
      n_steps = policy.shape[0]
      n_factors = policy.shape[1]
  
      # initialise posterior predictive density as a list of beliefs over time, including current posterior beliefs about hidden states as the first element
      qs_pi = [qs] + [
          utils.obj_array(n_factors) for t in range(n_steps)
      ]
  
      # get expected states over time
      for t in range(n_steps):
          for control_factor, action in enumerate(policy[t, :]):
              qs_pi[t + 1][control_factor] = B[control_factor][
                  :, :, int(action)
              ].dot(qs_pi[t][control_factor])
  
      return qs_pi[1:]
  
imports: import itertools
  import numpy
  from pymdp.maths import softmax, softmax_obj_arr, spm_dot, spm_wnorm, spm_MDP_G, spm_log_single, spm_log_obj_array
  from pymdp import utils
  import copy
label_desc: Use image node_id 3 for calling a global function with example usage: get_expected_states(qs, B, policy) and returns: qs_pi