code: def predictive_probability_multistate(M_c, X_L_list, X_D_list, Y, Q):
      
      logprobs = [float(predictive_probability(M_c, X_L, X_D, Y, Q))
          for X_L, X_D in zip(X_L_list, X_D_list)]
      return logmeanexp(logprobs)
docs: Returns the predictive probability, averaged over each sample.