diff --git a/examples/prompteval/demo_prompteval_dataloader.py b/examples/prompteval/demo_prompteval_dataloader.py
new file mode 100644
index 0000000..8ac1f58
--- /dev/null
+++ b/examples/prompteval/demo_prompteval_dataloader.py
@@ -0,0 +1,41 @@
+"""Demo: Loading and inspecting the PromptEvaldataset using the nearest_neighbors package.
+
+Run this script from the root of the repository with:
+    python examples/demo_prompteval_dataloader.py
+"""
+
+from nearest_neighbors.datasets.dataloader_factory import NNData
+
+print("====== PromptEval Data Demo ======")
+
+# Show dataset help info
+print("\nDataset Help Information:")
+NNData.help("prompteval")
+
+# Load with default settings
+print("\n====== Example 1: Default Parameters ======")
+loader = NNData.create("prompteval", seed=42)
+data, mask = loader.process_data_scalar()
+
+print("Observed Ratings Matrix (with NaNs):")
+print(data)
+print("\nMask Matrix (True = rating present, False = missing):")
+print(mask)
+
+# Load with custom sampling
+print("\n====== Example 2: Custom Parameters ======")
+loader_custom = NNData.create("movielens", seed=42, sample_users=100, sample_movies=50)
+data_custom, mask_custom = loader_custom.process_data_scalar()
+
+print("Custom-Sized Ratings Matrix:")
+print(data_custom)
+print("\nCustom Mask Matrix:")
+print(mask_custom)
+
+# View saved internal state
+print("\n====== Full State Dictionary Keys ======")
+state = loader.get_full_state_as_dict(include_metadata=True)
+print("Top-level keys:")
+print(list(state.keys()))
+print("\nCustom parameters used:")
+print(state.get("custom_params", {}))
diff --git a/src/nearest_neighbors/datasets/prompteval/__init__.py b/src/nearest_neighbors/datasets/prompteval/__init__.py
new file mode 100644
index 0000000..342ca37
--- /dev/null
+++ b/src/nearest_neighbors/datasets/prompteval/__init__.py
@@ -0,0 +1 @@
+from .loader import PromptEvalDataLoader  # noqa: F401
diff --git a/src/nearest_neighbors/datasets/prompteval/loader.py b/src/nearest_neighbors/datasets/prompteval/loader.py
new file mode 100644
index 0000000..30cc019
--- /dev/null
+++ b/src/nearest_neighbors/datasets/prompteval/loader.py
@@ -0,0 +1,255 @@
+"""Dataset loader for the PromptEval (MMLU) dataset.
+
+Source1: https://huggingface.co/datasets/PromptEval
+Source2: https://github.com/kyuseongchoi5/EfficientEval_BayesOpt
+Paper Reference for transformation implementation:
+    Felipe Maia Polo et al.
+    "Efficient multi-prompt evaluation of LLMs."
+    Neurips, 2024.
+    https://arxiv.org/pdf/2405.17202
+"""
+
+from nearest_neighbors.datasets.dataloader_base import NNDataLoader
+from nearest_neighbors.datasets.dataloader_factory import register_dataset
+import numpy as np
+import pandas as pd
+import pickle
+from typing import Any
+import logging
+from joblib import Memory
+
+
+memory = Memory(".joblib_cache", verbose=2)
+logger = logging.getLogger(__name__)
+
+params = {
+    "tasks": (
+        list[str],
+        None,
+        "List of tasks to evaluate on. By default, returns all.",
+    ),
+    "models": (
+        list[str],
+        None,
+        "List of models to evaluate by. By default, returns all.",
+    ),
+    "seed": (int, None, "Random seed for reproducibility"),
+    "propensity": (float, None, "Proportion of data to keep"),
+}
+
+
+@register_dataset("prompteval", params)
+class PromptEvalDataLoader(NNDataLoader):
+    """Data from the PromptEval study formatted into a matrix or tensor.
+    To initialize with default settings, use: NNData.create("prompteval").
+
+    """
+
+    urls = {
+        "full_data": "https://github.com/kyuseongchoi5/EfficientEval_BayesOpt/tree/main/data/MMLU/data_all.pkl"
+    }
+
+    def __init__(
+        self,
+        tasks: list[str] | None = None,
+        models: list[str] | None = None,
+        seed: int | None = None,
+        propensity: float = 1.0,  # Default to 1.0 (keeping all data)
+        **kwargs: Any,
+    ):
+        """Initializes the PromptEval data loader.
+
+        Args:
+        ----
+            tasks: benchmark tasks to evaluate on. Default: None (use all tasks).
+            models: models that are evaluated on for each tasks. Default: None (use all models).
+            seed: Random seed for reproducibility. Default: None
+            propensity: Proportion of data to keep. Default: 1.0
+            kwargs: Additional keyword arguments.
+
+        """
+        super().__init__(
+            **kwargs,
+        )
+        self.tasks = tasks
+        self.models = models
+        self.propensity = propensity
+        if seed is not None:
+            np.random.seed(
+                seed=seed
+            )  # instantiate random seed if provided but do it only once here
+
+    def process_data_scalar(self) -> tuple[np.ndarray, np.ndarray]:
+        """Processes the data into scalar setting. This implementation is applicable when generating (template * example) matrix, while fixing model and task.
+
+        Returns
+        -------
+            data: 2d processed data matrix of floats (in this case, each entry is boolean as the metric is correctness)
+            mask: Mask for processed data
+
+        """
+        if not self.tasks or not self.models:
+            raise ValueError("Tasks and models must be specified")
+
+        model = self.models[0]  # Use the first model
+        task = self.tasks[0]  # Use the first task
+        propensity = self.propensity
+
+        assert len(self.models) == 1, "Only one model is supported in scalar mode"
+        assert len(self.tasks) == 1, "Only one task is supported in scalar mode"
+
+        full_data = self._load_data()
+
+        df = pd.DataFrame(full_data[0][model])
+        df_subject = df[df["subject"] == task]
+        num_examples = len(
+            df_subject["correctness"]
+        )  # decide the column dimension of the data matrix
+
+        temp_example = np.zeros(
+            (len(full_data), num_examples)
+        )  # num_templates * num_examples
+
+        for j in range(len(full_data)):
+            # j : per prompt template
+            df = pd.DataFrame(full_data[j][model])
+            print(f"Loading {model} for {j}th prompt template with {task} subject")
+            df_subject = df[df["subject"] == task]
+            temp_example[j, :] = df_subject["correctness"]
+
+        temp_example_df = pd.DataFrame(temp_example)
+
+        # Create a mask of the original data according to the propensity
+        original_mask = (
+            temp_example_df.notna()
+        )  # This just tells us what's naturally present
+
+        n_rows = temp_example_df.shape[0]
+        n_cols = temp_example_df.shape[1]
+        n_rows_keep = int(n_rows * propensity)
+        n_cols_keep = int(n_cols * propensity)
+
+        # Randomly select which rows/columns to keep
+        rows_keep_indices = np.random.choice(n_rows, n_rows_keep, replace=False)
+        cols_keep_indices = np.random.choice(n_cols, n_cols_keep, replace=False)
+
+        # Create a new mask starting with all False
+        propensity_mask = pd.DataFrame(
+            False, index=original_mask.index, columns=original_mask.columns
+        )
+
+        # Set the randomly selected rows/columns to True
+        for i in rows_keep_indices:
+            for j in cols_keep_indices:
+                propensity_mask.iloc[i, j] = True
+
+        data = temp_example_df.to_numpy()
+        mask = propensity_mask.to_numpy()
+        self.data = data
+        self.mask = mask
+        return data, mask
+
+    def process_data_distribution(self) -> tuple[np.ndarray, np.ndarray]:
+        """Process the data into distributional setting.
+
+        Returns
+        -------
+            data: task * model * template matrix of floats (in this case, each entry average of correctness across examples)
+            mask: Mask for processed data
+
+        """
+        if not self.tasks or not self.models:
+            raise ValueError("Tasks and models must be specified")
+
+        models = self.models
+        tasks = self.tasks
+        propensity = self.propensity
+
+        full_data = self._load_data()
+
+        task_model_temp = np.zeros((len(tasks), len(models), len(full_data)))
+
+        for j in range(len(full_data)):
+            # j : per prompt template
+            for k, model in enumerate(models):
+                # model : per model
+                df = pd.DataFrame(full_data[j][model])
+                for l, subject in enumerate(tasks):
+                    print(
+                        f"Loading {model} for {j}th prompt template with {subject} subject"
+                    )
+                    df_subject = df[df["subject"] == subject]
+                    task_model_temp[l, k, j] = np.mean(df_subject["correctness"])
+
+        task_model_temp_df = pd.DataFrame(task_model_temp)
+
+        # Create a mask of the original data according to the propensity
+        original_mask = (
+            task_model_temp_df.notna()
+        )  # This just tells us what's naturally present
+
+        n_rows = task_model_temp_df.shape[0]
+        n_cols = task_model_temp_df.shape[1]
+        n_rows_keep = int(n_rows * propensity)
+        n_cols_keep = int(n_cols * propensity)
+
+        # Randomly select which rows/columns to keep
+        rows_keep_indices = np.random.choice(n_rows, n_rows_keep, replace=False)
+        cols_keep_indices = np.random.choice(n_cols, n_cols_keep, replace=False)
+
+        # Create a new mask starting with all False
+        propensity_mask = pd.DataFrame(
+            False, index=original_mask.index, columns=original_mask.columns
+        )
+
+        # Set the randomly selected rows/columns to True
+        for i in rows_keep_indices:
+            for j in cols_keep_indices:
+                propensity_mask.iloc[i, j] = True
+
+        data = task_model_temp_df.to_numpy()
+        mask = propensity_mask.to_numpy()
+        self.data = data
+        self.mask = mask
+        return data, mask
+
+    def get_full_state_as_dict(self, include_metadata: bool = False) -> dict:
+        """Returns the full state as a dictionary. For HeartSteps, this includes the data, masking matrix, and the custom parameters (if include_metadata == True
+
+        If the data and mask are None, then the data has not been processed yet. Call process_data_scalar() or process_data_distribution() to process the data first.
+
+        Args:
+            include_metadata (bool): Whether to include metadata in the dictionary. Default: False. The metadata for HeartSteps is currently empty.
+
+        """
+        full_state = {
+            "data": self.data,
+            "mask": self.mask,
+        }
+        return full_state
+
+    @classmethod
+    @memory.cache
+    def _load_data(cls) -> Any:
+        """Load the MMLU full dataset.
+
+        Returns:
+            full_data: Any Python object stored in the pickle file
+
+        """
+        logger.info("Retrieving MMLU full dataset from url...")
+        full_data_path = cls.urls["full_data"]
+
+        # GitHub URLs can't be directly downloaded; you would need to use raw content
+        # or download from releases. Here we're assuming the file is downloaded locally.
+        try:
+            # Using standard pickle module instead of pandas
+            with open(full_data_path, "rb") as f:
+                full_data = pickle.load(f)
+        except Exception as e:
+            logger.error(f"Error loading pickle file: {e}")
+            raise ValueError(
+                f"Could not load data from {full_data_path}. Make sure the file exists."
+            )
+
+        return full_data
diff --git a/src/nearest_neighbors/datasets/synthetic_data/loader.py b/src/nearest_neighbors/datasets/synthetic_data/loader.py
index a4a539b..cef3265 100644
--- a/src/nearest_neighbors/datasets/synthetic_data/loader.py
+++ b/src/nearest_neighbors/datasets/synthetic_data/loader.py
@@ -1,6 +1,7 @@
 from nearest_neighbors.datasets.dataloader_base import NNDataLoader
 from nearest_neighbors.datasets.dataloader_factory import register_dataset
 import numpy as np
+import math as math
 from typing import Any
 
 params = {
@@ -186,8 +187,90 @@ def _make_mcar(self) -> None:
         self.availability_mask = A
 
     def _make_mnar(self) -> None:
-        raise NotImplementedError("MNAR yet to be implemented")
-        # TODO: Aashish/Caleb/Kyuesong/Tatha: come back to this later
+        """Makes values missing not at random (MNAR), specifically staggered adoption (non-positivity & confounded)"""
+        U = self.row_latent
+
+        N = self.num_rows
+        T = self.num_cols
+
+        missing_mask = np.zeros((N, T))
+        pre_Masking = np.zeros((N, T))
+
+        # Divide units into 3 groups
+        g1_inds = np.arange(0, N // 3)
+        g2_inds = np.arange(N // 3, 2 * N // 3)
+        g3_inds = np.arange(2 * N // 3, N)
+
+        gamma_1 = [2, 0.7, 1, 0.7]
+        gamma_2 = [2, 0.2, 1, 0.2]
+
+        # TODO: make beta a parameter (currently hardcoded)
+        # first group adopts at the first 30% of the time period
+        # second group adopts at the first 70% of the time period
+        beta = [0.3, 0.7]
+
+        T1_lower = math.floor(T ** beta[0])
+        T2_lower = math.floor(T ** beta[1])
+
+        for i in range(N):
+            if i in g1_inds:
+                pre_Masking[i, :] = np.concatenate(
+                    (np.ones(T1_lower), np.zeros(T - T1_lower))
+                )
+                for t in range(T - T1_lower):
+                    pre_Masking[i, (t + T1_lower)] = (
+                        np.random.binomial(  # each units' adoption time probability is affected by their neighbors
+                            1,
+                            self._expit(
+                                gamma_1[0]
+                                + (0.99**t) * gamma_1[1] * U[i - 1]
+                                + gamma_1[2] * U[i]
+                                + (0.99**t) * gamma_1[3] * U[i + 1]
+                            ),
+                            1,
+                        )
+                    )
+                pre_A = pre_Masking[i, :]
+                if len([i for i in range(len(pre_A)) if pre_A[i] == 0]) == 0:
+                    missing_mask[i, :] = pre_A
+                elif len([i for i in range(len(pre_A)) if pre_A[i] == 0]) > 0:
+                    adopt_time = min([i for i in range(len(pre_A)) if pre_A[i] == 0])
+                    missing_mask[i, :] = np.concatenate(
+                        (np.ones(adopt_time), np.zeros(T - adopt_time))
+                    )
+            elif i in g2_inds:
+                pre_Masking[i, :] = np.concatenate(
+                    (np.ones(T2_lower), np.zeros(T - T2_lower))
+                )
+                for t in range(T - T2_lower):
+                    pre_Masking[i, (t + T2_lower)] = np.random.binomial(
+                        1,
+                        self._expit(
+                            gamma_2[0]
+                            + (1.01**t) * gamma_2[1] * U[i - 1]
+                            + gamma_2[2] * U[i]
+                            + (1.01**t) * gamma_2[3] * U[i + 1]
+                        ),
+                        1,
+                    )
+                pre_A = pre_Masking[i, :]
+                if len([i for i in range(len(pre_A)) if pre_A[i] == 0]) == 0:
+                    missing_mask[i, :] = pre_A
+                elif len([i for i in range(len(pre_A)) if pre_A[i] == 0]) > 0:
+                    adopt_time = min([i for i in range(len(pre_A)) if pre_A[i] == 0])
+                    missing_mask[i, :] = np.concatenate(
+                        (np.ones(adopt_time), np.zeros(T - adopt_time))
+                    )
+            elif i in g3_inds:
+                missing_mask[i, :] = np.ones(T)
+
+        data_obs = self.data_noisy.copy()
+        data_obs[missing_mask.astype(bool)] = np.nan
+        A = ~missing_mask.astype(
+            bool
+        )  # A = NOT M, i.e. A_ij = 1 if Y_ij is observed, 0 if missing
+        self.data_obs = data_obs
+        self.availability_mask = A
 
     def get_full_state_as_dict(self, include_metadata: bool = False) -> dict:
         """Returns the full state of this object as a dictionary"""
diff --git a/src/nearest_neighbors/fit_method.py b/src/nearest_neighbors/fit_method.py
new file mode 100644
index 0000000..0246c3d
--- /dev/null
+++ b/src/nearest_neighbors/fit_method.py
@@ -0,0 +1,149 @@
+from .nnimputer import FitMethod
+from .nnimputer import NearestNeighborImputer
+from .data_types import DistributionKernelMMD
+import numpy.typing as npt
+import numpy as np
+
+
+class DirectOptimization(FitMethod):
+    """Non-cross-validation fit method. Analytically optimizes the squared MMD error."""
+
+    def __init__(
+        self, row: int, column: int, kernel: str, eta_cand: npt.NDArray, delta: float
+    ):
+        """Initialize the fit method with additional parameters.
+
+        Args:
+            kernel (str): Kernel to use for the MMD
+            eta_cand (npt.NDArray): Candidate distance thresholds
+            delta (float): Significance level
+            row (int): target row index
+            column (int): target column index
+
+        """
+        supported_kernels = ["exponential"]
+
+        if kernel not in supported_kernels:
+            raise ValueError(
+                f"Kernel {kernel} is not supported. Supported kernels are {supported_kernels}"
+            )
+
+        self.kernel = kernel
+        self.eta_cand = eta_cand
+        self.delta = delta
+        self.row = row
+        self.column = column
+
+    def fit(
+        self,
+        data_array: npt.NDArray,
+        mask_array: npt.NDArray,
+        imputer: NearestNeighborImputer,
+    ) -> float:
+        """Analytically optimizes the squared MMD error.
+
+        Args:
+            data_array (npt.NDArray): Data array
+            mask_array (npt.NDArray): Mask array
+            imputer (NearestNeighborImputer): Imputer
+
+        Returns:
+            float: Best distance threshold
+
+        """
+        # Initialize sup_kern outside conditional
+        sup_kern = 1  # Default value
+        if self.kernel == "exponential":
+            sup_kern = 1  # TODO: need to change for general kernels
+
+        delta = self.delta
+        eta_cand = self.eta_cand
+        row = self.row
+        column = self.column
+
+        n_rows, n_cols = data_array.shape[0], data_array.shape[1]
+        n = data_array[0, 0].shape[0]  # number of samples per distribution
+        data_type = DistributionKernelMMD(self.kernel)
+
+        row_distances = np.zeros(n_rows)
+        for i in range(n_rows):
+            # Get columns observed in both row i and row
+            overlap_columns = np.logical_and(mask_array[row], mask_array[i])
+
+            if not np.any(overlap_columns):
+                row_distances[i] = np.inf
+                continue
+
+            # Calculate distance between rows
+            for j in range(n_cols):
+                if (
+                    not overlap_columns[j] or j == column
+                ):  # Skip missing values and the target column
+                    continue
+                row_distances[i] += data_type.distance(
+                    data_array[row, j], data_array[i, j]
+                )
+            row_distances[i] /= np.sum(overlap_columns)
+
+        perf = []
+
+        for eta in eta_cand:
+            neighborhood = np.where(
+                (row_distances < eta) * (mask_array[:, column]) == 1
+            )[0]  # Set of neighbors: (i) within eta distance (ii) observed
+
+            if (
+                sum(np.isin(neighborhood, row)) == 1
+            ):  # Pretending as if (row, column) entry is missing
+                neighborhood = np.delete(neighborhood, np.where(neighborhood == row)[0])
+
+            if (
+                len(neighborhood) == 0
+            ):  # Default (null) output when there is zero neighbor
+                perf.append(10**5)  # Avoid selecting such eta without neighbors
+            else:
+                overlap = []
+                for neighbor in neighborhood:
+                    overlap.append(np.sum(mask_array[row, :] * mask_array[neighbor, :]))
+
+                Bias = (
+                    8
+                    * np.exp(1 / np.exp(1))
+                    * sup_kern
+                    * np.log(2 * n_rows / delta)
+                    / (np.sqrt(2 * np.log(2) * np.min(overlap)))
+                )
+                Variance = 4 * sup_kern * (np.log(n) + 1.5) / (n * len(neighborhood))
+
+                perf.append(eta + Bias + Variance)
+
+        if not perf:  # Handle case when perf list is empty
+            return float(
+                "inf"
+            )  # Return infinity as a default value when no valid threshold is found
+
+        eta_star = eta_cand[np.argmin(perf)]
+        return eta_star
+
+
+class CrossValidation(FitMethod):
+    """Cross-validation fit method. Uses cross-validation to find the best distance threshold."""
+
+    def fit(
+        self,
+        data_array: npt.NDArray,
+        mask_array: npt.NDArray,
+        imputer: NearestNeighborImputer,
+    ) -> float:
+        """Uses cross-validation to find the best distance threshold.
+
+        Args:
+            data_array (npt.NDArray): Data array
+            mask_array (npt.NDArray): Mask array
+            imputer (NearestNeighborImputer): Imputer
+
+        Returns:
+            float: Best distance threshold
+
+        """
+        return 0.0  # TODO: Implement cross-validation
diff --git a/src/nearest_neighbors/nnimputer.py b/src/nearest_neighbors/nnimputer.py
index 5976930..6258829 100644
--- a/src/nearest_neighbors/nnimputer.py
+++ b/src/nearest_neighbors/nnimputer.py
@@ -128,6 +128,8 @@ def fit(
         """Find the best distance threshold for the given data.
 
         Args:
+            row (int): Row index
+            column (int): Column index
             data_array (npt.NDArray): Data matrix
             mask_array (npt.NDArray): Mask matrix
             imputer (NearestNeighborImputer): Imputer object
diff --git a/src/nearest_neighbors/simulations/mcar.py b/src/nearest_neighbors/simulations/mcar.py
index 4d0c7d1..62fda56 100644
--- a/src/nearest_neighbors/simulations/mcar.py
+++ b/src/nearest_neighbors/simulations/mcar.py
@@ -139,3 +139,54 @@ def expit(x: np.ndarray) -> np.ndarray:
     # Data[Masking == 0] = Y0[Masking == 0]
     Data: np.ndarray = np.array(Y)
     return Data, Theta, Masking
+
+
+def gendata_dist_mcar(
+    N: int, T: int, n: int, d: int, p: float, seed: int
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    """Generates multivariate Gaussian data of multiple measurements with latent dimension r = 2.
+
+    Args:
+        N (int): Number of users.
+        T (int): Number of time periods.
+        n (int): Number of samples per distribution.
+        d (int): Dimension of the data.
+        p (float): Probability of an entry being observed.
+        seed (int): Random seed for reproducibility.
+
+    Returns:
+        Data (np.ndarray): Generated data matrix of shape (N, T, n, d).
+        Masking (np.ndarray): Masking matrix indicating observed entries of shape (N, T).
+        True Mean (np.ndarray): True mean of the data of shape (N, T, d).
+        True Covariance (np.ndarray): True covariance of the data of shape (N, T, d, d).
+
+    """
+    np.random.seed(seed=seed)
+
+    ## Data Matrix (N * T * n * d)
+    Data = np.zeros((N, T, n, d))
+    true_Mean = np.zeros((N, T, d))
+    true_Cov = np.zeros((N, T, d, d))
+
+    u_1 = np.random.uniform(-1, 1, N)
+    u_2 = np.random.uniform(0.2, 1, N)
+
+    v_1 = np.random.uniform(-2, 2, T)
+    v_2 = np.random.uniform(0.5, 2, T)
+
+    even_ones = np.repeat([0, 1], int(d / 2))
+    odd_ones = np.repeat([1, 0], int(d / 2))
+
+    for i in range(N):
+        for t in range(T):
+            m_it = u_1[i] * v_1[t] * (even_ones - odd_ones)
+            c_it = np.diag(u_2[i] * v_2[t] * (0.5 * even_ones + odd_ones))
+            true_Mean[i, t, :] = m_it
+            true_Cov[i, t, :, :] = c_it
+            dat_mat = np.random.multivariate_normal(m_it, c_it, size=n)
+            Data[i, t, :, :] = dat_mat
+
+    Masking = np.zeros((N, T))
+    Masking = np.reshape(np.random.binomial(1, p, (N * T)), (N, T))
+
+    return Data, Masking, true_Mean, true_Cov