def create_expanded_design_for_mixing(design,
                                      draw_list,
                                      mixing_pos,
                                      rows_to_mixers):
    """
    Parameters
    ----------
    design : 2D ndarray.
        All elements should be ints, floats, or longs. Each row corresponds to
        an available alternative for a given individual. There should be one
        column per index coefficient being estimated.
    draw_list : list of 2D ndarrays.
        All numpy arrays should have the same number of columns (`num_draws`)
        and the same number of rows (`num_mixers`). All elements of the numpy
        arrays should be ints, floats, or longs. Should have as many elements
        as there are lements in `mixing_pos`.
    mixing_pos : list of ints.
        Each element should denote a column in design whose associated index
        coefficient is being treated as a random variable.
    rows_to_mixers : 2D scipy sparse array.
        All elements should be zeros and ones. Will map the rows of the design
        matrix to the particular units that the mixing is being performed over.
        Note that in the case of panel data, this matrix will be different from
        `rows_to_obs`.

    Returns
    -------
    design_3d : 3D numpy array.
        Each slice of the third dimension will contain a copy of the design
        matrix corresponding to a given draw of the random variables being
        mixed over.
    """
    if len(mixing_pos) != len(draw_list):
        msg = "mixing_pos == {}".format(mixing_pos)
        msg_2 = "len(draw_list) == {}".format(len(draw_list))
        raise ValueError(msg + "\n" + msg_2)

    # Determine the number of draws being used. Note the next line assumes an
    # equal number of draws from each random coefficient's mixing distribution.
    num_draws = draw_list[0].shape[1]
    orig_num_vars = design.shape[1]

    # Initialize the expanded design matrix that replicates the columns of the
    # variables that are being mixed over.
    arrays_for_mixing = design[:, mixing_pos]
    expanded_design = np.concatenate((design, arrays_for_mixing),
                                     axis=1).copy()
    design_3d = np.repeat(expanded_design[:, None, :],
                          repeats=num_draws,
                          axis=1)

    # Multiply the columns that are being mixed over by their appropriate
    # draws from the normal distribution
    for pos, idx in enumerate(mixing_pos):
        rel_draws = draw_list[pos]
        # Note that rel_long_draws will be a dense, 2D numpy array of shape
        # (num_rows, num_draws).
        rel_long_draws = rows_to_mixers.dot(rel_draws)
        # Create the actual column in design 3d that should be used.
        # It should be the multiplication of the draws random variable and the
        # independent variable associated with the param that is being mixed.
        # NOTE THE IMPLICIT ASSUMPTION THAT ONLY INDEX COEFFICIENTS ARE MIXED.
        # Also, the final axis is selected on because the final axis sepecifies
        # the particular variable being multiplied by the draws. We select with
        # orig_num_vars + pos since the variables being mixed over were added,
        # in order so we simply need to start at the first position after all
        # the original variables (i.e. at orig_num_vars) and iterate.
        design_3d[:, :, orig_num_vars + pos] *= rel_long_draws

    return design_3d
Parameters
    ----------
    design : 2D ndarray.
        All elements should be ints, floats, or longs. Each row corresponds to
        an available alternative for a given individual. There should be one
        column per index coefficient being estimated.
    draw_list : list of 2D ndarrays.
        All numpy arrays should have the same number of columns (`num_draws`)
        and the same number of rows (`num_mixers`). All elements of the numpy
        arrays should be ints, floats, or longs. Should have as many elements
        as there are lements in `mixing_pos`.
    mixing_pos : list of ints.
        Each element should denote a column in design whose associated index
        coefficient is being treated as a random variable.
    rows_to_mixers : 2D scipy sparse array.
        All elements should be zeros and ones. Will map the rows of the design
        matrix to the particular units that the mixing is being performed over.
        Note that in the case of panel data, this matrix will be different from
        `rows_to_obs`.

    Returns
    -------
    design_3d : 3D numpy array.
        Each slice of the third dimension will contain a copy of the design
        matrix corresponding to a given draw of the random variables being
        mixed over.