A Monte Carlo-Based Model for Predicting Student Performance

1. Introduction

Predicting student academic performance has emerged as a key research area in Educational Data Mining (EDM) and Learning Analytics, motivated by the need to identify at-risk learners early, deliver personalized interventions, and improve educational outcomes. Accurate predictive models help institutions support timely decision-making, enhance student retention, and improve academic success rates in educational settings [1,2,3].

Student performance prediction typically entails forecasting future academic outcomes—such as grade point average (GPA), course grades, or pass/fail status—using information derived from demographic attributes, academic history, engagement indicators, and learning management system interactions [3]. A wide range of supervised machine learning techniques, including Decision Trees, Support Vector Machines, Random Forests, k-Nearest Neighbors, and Logistic Regression, have been extensively applied to this problem using structured student datasets [4,5,6]. Ensemble learning approaches often achieve higher accuracy by capturing nonlinear relationships among features; however, their performance remains sensitive to data imbalance, feature selection, and noisy or incomplete data.

Recent studies have explored advanced modeling paradigms such as deep neural networks, explainable learning frameworks, heterogeneous ensemble systems, and privacy-preserving federated learning approaches [7,8,9]. While these methods offer improved predictive performance in controlled settings, they typically assume deterministic data representations and produce point estimates without explicitly modeling uncertainty. Moreover, many deep models require large volumes of high-quality data and substantial computational resources, limiting their scalability and interpretability across diverse educational contexts.

These limitations motivate exploration of probabilistic and stochastic modeling frameworks capable of handling uncertainty, variability, and complex data distributions inherent in real-world educational systems. Here, Monte Carlo methods provide a theoretically grounded and computationally scalable alternative.

Monte Carlo algorithms constitute a broad class of randomized computational techniques that rely on stochastic sampling to obtain approximate numerical solutions to analytically intractable problems [10,11]. These methods are particularly effective in high-dimensional settings, where deterministic numerical techniques suffer from exponential computational growth, commonly referred to as the “curse of dimensionality.”

Monte Carlo methods estimate expected values and multidimensional integrals by averaging over randomly generated samples. Consider the estimation of the expectation of a measurable function $f\left(X\right)$:

$$\mathbb{E}\left[f\right(X\left)\right]=\int f\left(x\right)p\left(x\right)dx,$$

where $X$is a random variable with probability density function $p\left(x\right)$. When direct analytical evaluation is infeasible, $N$independent samples $X_1,X_2,\dots ,X_N\sim p\left(x\right)$are generated, and the Monte Carlo estimator is defined as

$${\widehat{\mu }}_N={\displaystyle \frac{1}{N}}\sum _{i=1}^{N}f(X_i).$$

By the Law of Large Numbers (LLN), the estimator converges almost surely to the true expectation as $N\to \mathrm{\infty }$. Furthermore, the Central Limit Theorem (CLT) guarantees asymptotic normality of the estimator, enabling statistical error analysis and confidence interval construction [11,12]. The root-mean-square error (RMSE) of the estimator decreases at a rate

$$\mathrm{RMSE}\left({\widehat{\mu }}_N\right)=\mathcal{O}(1/\sqrt{N}),$$

which is independent of the dimensionality of the underlying problem. This property makes Monte Carlo methods particularly suitable for complex, high-dimensional systems.

A central challenge in Monte Carlo simulation is variance control, as high estimator variance slows convergence. Advanced variance-reduction techniques—such as importance sampling, control variates, and stratified sampling are commonly employed to improve efficiency and numerical stability [11,12]. These techniques have proven effective across domains involving uncertainty, including finance, engineering, and increasingly, predictive analytics.

In predictive modeling, Monte Carlo methods offer a distinct advantage by generating probabilistic forecasts rather than single deterministic predictions. By producing distributions of possible outcomes with associated likelihoods, Monte Carlo-based models provide richer insights into uncertainty and risk, which are essential for educational early-warning systems and academic decision support.

The main contributions of this paper are summarized as follows:

• Novel Monte Carlo Framework: We propose a Monte Carlo–based GPA prediction framework that explicitly incorporates probabilistic sampling, enabling reliable predictions under uncertainty and heterogeneous student data.
• Variance Reduction Techniques: The model integrates advanced variance-reduction strategies, improving convergence speed and stability for high-dimensional educational datasets.
• Comparative Evaluation: We evaluate the proposed framework using real-world datasets from PSAU’s engineering and medical faculties, demonstrating superior performance relative to state-of-the-art machine learning and deep learning models.
• Probabilistic Insight: Unlike deterministic predictors, the Monte Carlo model generates predictive distributions, providing actionable insights into uncertainty and risk in student academic outcomes.
• Scalability and Practical Deployment: The approach is computationally efficient and can be deployed in real educational settings, supporting early warning and targeted intervention systems.

The remainder of this paper is organized as follows. Section II reviews recent studies on student performance prediction. Section III presents the proposed Monte Carlo-based model. Section IV describes the experimental setup and evaluation metrics. Section V discusses the results and comparative analysis. Section VI concludes the paper and outlines recommendations for future research.

2. Previous Studies

While most research on predicting student performance traditionally relies on machine learning algorithms such as Support Vector Machines (SVM), decision trees, and linear regression , recent studies continue to corroborate the effectiveness of these approaches in educational contexts. For example, a recent systematic review highlights the extensive use of classical ML models (e.g., SVM, decision trees, logistic regression) for forecasting academic success and identifying at-risk students, demonstrating that such models often achieve high classification accuracy and enable early interventions based on behavioral and institutional attributes [13]. Empirical studies further show the effectiveness of these techniques: classification algorithms, including Random Forests, SVM, and nearest neighbors, have been successfully applied to predict student exam performance, helping institutions identify students who are likely to succeed or fail [14]. Additional research confirms that both traditional ML and deep learning models (e.g., neural networks) provide consistent performance prediction when evaluated on real academic datasets [15].

In contrast, although Monte Carlo methods have been extensively applied across a wide range of domains, including physics (as examples [16,17]), finance (as examples [18,19]), engineering (as examples [20,21]), healthcare (as examples [22,23]), and risk analysis (as examples [24,25,26]). Their adoption in the field of educational data mining remains limited. In particular, relatively few studies have explicitly employed Monte Carlo methods for predicting students’ performance outcomes, as in [27,28,29,30,31].

In [27], the authors develope a Monte Carlo simulation predictive model that allows prospective students to estimate the number of applicants to the Department of Hadith Sciences for the upcoming academic year. The research technique employed in this study uses a Monte Carlo simulation to apply this research. The number of students in the Department of Hadith Sciences who enrolled from 2019 to 2023 serves as the data used to predict the number of new student registrations. They develop a simulation system to forecast the number of new applicants to the Hadith Sciences program using the classical simulation methodology.

In [28], the authors develop a generic computer-based tool that uses Monte Carlo stochastic simulation and historical data to analyze and predict student performance within a given curriculum. They model subject completion as a sequence of probabilistic events; their model computes critical metrics defined by Spain's National Quality Evaluation and Accreditation Agency (ANECA), specifically the graduation rate (students finishing on time or with one extra year) and the efficiency rate (the ratio of required to actual credits taken). They design an experiment to test how changes in regulations, such as the maximum number of credits a student can enroll in, affect overall success before implementing new curricula. Unlike previous studies focused on individual students, this simulator provides a "workshop" to optimize entire academic programs and identify which features most influence indicators, as demonstrated by its specific application to the Bachelor in Engineering Technologies at the Technical University of Madrid

In [29], the authors examine factors that may affect students' motivation to learn scientific computing in the undergraduate systems engineering program at the University of Córdoba (Colombia). They conducted a survey of 117 students and then applied multidimensional linear regression to identify a function that captures the regular patterns associated with motivation and its influencing factors. They used the Monte Carlo method in exploring the multidimensional space of independent variables to estimate the probability that a student attains one of ten motivation levels, ranging from completely demotivated to highly motivated to learn scientific computing.

In [30], the authors examine the divergent mechanisms by which teachers' support for students affects students' performance on low- and high-stakes assessments. They used Bayesian structural equation modeling (BSEM) with the Markov Chain Monte Carlo (MCMC) algorithm to evaluate the structural relationships among teacher support, student behavior, and [formative and summative] academic performance. We also assessed moderated relationships, focusing on students' gender and interactions with teacher support. They used Bayesian structural equation modeling with the Markov Chain Monte Carlo algorithm and Ghanaian data to capture the cumulative nature of student learning and performance assessment. They evaluated student performance during daily schooling (formative assessment) and at the end of their academic program (summative assessment).

In [31], the authors study factors that can affect student performance in online learning by developing an empirical model. They tested the reliability of their model linguist reviews using construct reliability coefficients and confirmatory factor analysis (CFA), and a one-sample Kolmogorov-Smirnov test with the Monte Carlo method.

3. The Model

Figure 1 presents the proposed Monte Carlo–based framework for predicting student performance. The model follows a structured pipeline that begins with data acquisition from the real educational system and proceeds through probabilistic modeling, stochastic simulation, and performance evaluation. Each stage of the model is designed to capture the uncertainty inherent in student behavior and academic outcomes.

3.1. Data Collection from the Real System

The first stage of the proposed model involves gathering raw data from the real educational environment. This dataset includes historical records related to students’ academic performance, such as assessment scores, course completion status, and other relevant academic indicators. The collected data serve as the empirical foundation upon which the probabilistic model is constructed.

3.2. Data Preprocessing

Before probabilistic modeling, the raw data undergoes preprocessing to ensure consistency, completeness, and suitability for analysis. This step includes handling missing values, removing outliers, and transforming variables where necessary. The preprocessed dataset is then divided into two subsets: a training dataset, used to estimate the underlying probability distributions, and a testing dataset, reserved exclusively for evaluating model performance.

The training data are used to construct the probability model and perform Monte Carlo simulations, while the testing data remain unseen during training to provide an unbiased evaluation of predictive accuracy. This separation ensures that model validation reflects real-world predictive performance.

3.3. Probability of Occurrence Estimation

Using the training dataset, the probability of occurrence for each observed outcome is estimated empirically. Let $X$denote a discrete random variable representing a student performance metric. The probability mass function is computed as:

$$p\left(x\right)={\displaystyle \frac{\mathrm{frequency}\mathrm{of}x}{\mathrm{total}\mathrm{number}\mathrm{of}\mathrm{observations}}}$$

This empirical probability estimation captures the observed distribution of student outcomes without imposing parametric assumptions.

3.4. Cumulative Distribution Function Construction

Based on the estimated probabilities, the cumulative distribution function (CDF) is constructed. The CDF is defined as:

$$F_X\left(x\right)=P(X\le x)$$

The CDF plays a central role in the Monte Carlo simulation process, as it enables the transformation of uniformly distributed random numbers into samples that follow the empirical distribution of student performance [32]

3.5. Random Number Generation

Random number generation constitutes a core component of the proposed Monte Carlo–based simulation framework. At this stage, pseudo-random numbers uniformly distributed over the interval $\left[0,1\right]$are generated and used as stochastic inputs for sampling from the empirical cumulative distribution function (CDF) of student performance.

In computational simulations, true randomness is typically approximated using pseudo-random number generators (PRNGs), which produce deterministic sequences of numbers that exhibit statistical properties consistent with independent realizations from a target distribution. Most modern PRNGs are designed to generate sequences that closely approximate the continuous uniform distribution over $\left[0,1\right]$, making them suitable for Monte Carlo simulation and probabilistic modeling tasks [33,34]

Let $U$be a random variable uniformly distributed on the unit interval, i.e. :

$$U\sim \mathcal{U}\left(\mathrm{0,1}\right).$$

Accordingly, a sequence of independent pseudo-random numbers $\left\{u_1,u_2,\dots ,u_N\right\}$is generated, where each $u_i\in \left[\mathrm{0,1}\right]$represents a stochastic realization. These values serve as inputs for inverse CDF sampling from the empirical distribution of student performance outcomes.

From an algorithmic perspective, a wide range of well-established PRNGs exist for generating uniform random numbers, including Linear Congruential Generators (LCGs), Mersenne Twister, Xorshift, and Permuted Congruential Generators (PCG). These algorithms differ in terms of period length, equidistribution properties, and computational efficiency. However, contemporary scientific programming environments and software platforms (e.g., MATLAB, Python, R, Java) implement high-quality PRNGs by default, which have been extensively validated for Monte Carlo applications.

Since the proposed model relies exclusively on uniformly distributed pseudo-random numbers and does not require specialized stochastic properties beyond uniformity and independence, the internal PRNG provided by the software environment is sufficient. Therefore, the design of a custom random number generation algorithm is outside the scope of this work. The focus instead is placed on how the generated uniform samples are transformed via the empirical CDF to model uncertainty in student performance predictions.

3.6. Monte Carlo Prediction via Distribution Transformation

The generated random numbers are mapped to predicted student performance values using inverse transformation and empirical distribution techniques. By repeatedly sampling from the CDF, the model produces a large number of simulated outcomes, thereby estimating the probabilistic behavior of student performance under uncertainty. This Monte Carlo process allows the model to generate predictions that reflect variability and uncertainty rather than single-point estimates.

To transform the uniformly distributed random variables into samples following the empirical distribution of student performance, the inverse cumulative distribution function (inverse CDF) method is employed [34,35] as follows:

Let $F_X\left(x\right)$denote the cumulative distribution function of the discrete random variable $X$, constructed from the training data as:

$$F_X\left(x\right)=\sum _{t\le x}p\left(t\right).$$

The inverse CDF, also known as the quantile function, is defined as:

$$F_X^{-1}\left(u\right)=\mathrm{i}\mathrm{n}\mathrm{f}\{x\in \mathbb{R}:F_X(x)\ge u\},u\in \left[\mathrm{0,1}\right].$$

A simulated student performance value $\widehat{x}$is then obtained by:

$$\widehat{x}=F_X^{-1}\left(u\right),u\sim \mathcal{U}\left(\mathrm{0,1}\right).$$

By repeatedly applying this transformation for a large number of samples $N$, the model generates a set of simulated outcomes:

$$\{{\widehat{x}}_1,{\widehat{x}}_2,\dots ,{\widehat{x}}_N\}.$$

These simulated values approximate the empirical distribution of student performance and enable probabilistic prediction under uncertainty.

3.7. Model Evaluation and Result Analysis

The model outcomes are then compared with the testing dataset to assess predictive performance through a designed and implemented experiment. Evaluation metrics such as prediction error, distributional similarity, or classification accuracy (depending on the application context) are computed. This comparison quantifies the model’s ability to generalize beyond the training data. The final output of this stage represents the evaluation results, which indicate the effectiveness of the proposed Monte Carlo–based prediction model.

The expected predicted performance can be estimated as:

$$\mathbb{E}\left[X\right]\approx {\displaystyle \frac{1}{N}}\sum _{i=1}^N{\widehat{x}}_i,$$

while the variability of predicted outcomes is captured by the sample variance:

$$\mathrm{V}\mathrm{a}\mathrm{r}\left(X\right)\approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N({\widehat{x}}_i-\bar{x}{)}^2,\bar{x}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\widehat{x}}_i.$$

This formulation allows the proposed model to provide not only point predictions but also uncertainty-aware estimates of student performance.

Beyond experiment-specific diagnostics, a broader evaluation is conducted through a benchmark to assess the overall predictive capability of the proposed modeling framework relative to deterministic alternatives.

4. Experimental Setup and Evaluation

The proposed model was implemented in Python. We used Python because it is a widely adopted, high-level, interpreted programming language that supports scientific computing and data analysis. Owing to its extensive ecosystem of numerical and statistical libraries, it makes it an effective and flexible environment for implementing Monte Carlo simulations and probabilistic modeling techniques[36,37].

This section details the experimental design and implementation of the proposed model. Specifically, it describes the dataset characteristics, the preprocessing and data transformation procedures, and the key implementation details, while ensuring consistency with the model framework introduced in Section 3.

4.1. Dataset Description

The dataset used in this study consists of real academic records collected from the Primary Year program at Prince Sattam bin Abdulaziz University (PSAU) over two consecutive academic years: 2015/2016 ($Y_1$) and 2016/2017 ($Y_2$). In the primary year, there are two major disciplines: one for students who decide to study engineering and science, and the other for students who decide to study medicine and medical science. The total number of students in the academic year $Y_1$was 1224 students (675 in the engineering discipline and 537 in the medical science discipline), while the total number of students in the academic year $Y_2$was 959 students (554 in the engineering discipline and 405 in the medical science).

To evaluate the predictive capability of the proposed model, the dataset was temporally partitioned. Records from the $Y_1$academic year were used to construct the empirical distributions required by the model (training phase), while data from the $Y_2$ academic year were reserved for testing and validation. This split reflects a realistic deployment scenario in which historical academic data are used to predict future student outcomes. In addition, we will train and test our model on 3 categories of datasets: individual-discipline datasets and the entire dataset.

In the PSAU Primary Year system, students complete two academic semesters, after which their Grade Point Average (GPA) is computed. This GPA determines eligibility and placement into specific colleges, making it a critical variable for both students and academic administrators. Consequently, accurate modeling of GPA-related outcomes is essential for informed decision-making, academic planning, and early intervention strategies.

The dataset comprises all student data during this year, including the final GPA attributes, making it well-suited for assessing the flexibility of the proposed probabilistic model in handling heterogeneous data types.

Table 1 presents a sample of part of the original dataset after de-identifying the data to protect student privacy.

4.2. Data Preprocessing

Data preprocessing was performed to align the raw academic records with the requirements of the proposed Monte Carlo simulation framework. At this stage, student records were organized and categorized according to academic performance levels.

Specifically, students were grouped into six discrete grade categories based on GPA ranges defined by the Primary Year academic regulations. Each GPA value was mapped to a corresponding grade label. This discretization step serves two purposes:

• It reduces continuous GPA values into interpretable performance classes, and
• It facilitates the construction of empirical probability distributions required for Monte Carlo sampling.

To demonstrate the model’s capability to process both numerical and textual representations, grade information was retained in textual form while preserving numerical GPA values for computation. Frequency counts for each grade category were then calculated.

The preprocessing algorithm, illustrated in Figure 2, is responsible for transforming raw academic data into a structured form suitable for probabilistic modeling. The process begins by reading the cleaned dataset from an Excel file and extracting the Grade Point Average (GPA) values into a one-dimensional array.

Let the set of GPA values be defined as:

$$G=\{g_1,g_2,\dots ,g_N\}$$

where $N$denotes the total number of students.

The GPA values are optionally sorted to support descriptive analysis and distribution inspection.

Each GPA value is then classified into a discrete grade category based on predefined institutional thresholds using an iterative conditional structure. The classification rule is formally defined as:

$$f\left(g\right)=\left\{\begin{array}{cc}A,& g\ge 4.5\\ B,& 3.75\le g<4.5\\ C,& 2.75\le g<3.75\\ D,& 2.0\le g<2.75\\ F,& 0<g<2.0\\ W,& g=0\end{array}\right.$$

The resulting preprocessed dataset, which includes grade labels and their associated frequencies, is summarized in Table 2 for the results for the first year 2015/2016 ($Y_1$) (the training dataset) and in Table 3 for the results for the first year 2015/2016 ($Y_1$) (the training dataset).

4.3. Data Transformation and Probabilistic Modeling

Following preprocessing, the dataset was transformed into a probabilistic representation compatible with the Monte Carlo model process described in Section 3.

For each grade category, the probability of occurrence was computed by dividing its frequency by the total number of students. Based on these probabilities, the cumulative distribution function (CDF) was constructed. The CDF provides a monotonic mapping between uniformly distributed random numbers and grade categories, enabling inverse transform sampling.

Additionally, random-digit intervals were assigned to each grade according to cumulative probability ranges. These intervals define the mapping rule by which a generated uniform random number is associated with a specific grade outcome. This transformation step is essential for translating stochastic realizations into predicted academic performance categories.

The processing algorithm, shown in Figure 3, converts the grade frequencies obtained during preprocessing into a probabilistic representation. Let $f_i$represent the frequency of grade category $i$, where:

$$i\in \{A,B,C,D,F,O\}$$

The probability of occurrence for each grade is computed as:

$$P_i={\displaystyle \frac{f_i}{N}}$$

These probabilities form a discrete probability distribution:

$$P=\{P_A,P_B,P_C,P_D,P_F,P_O\}$$

To enable Monte Carlo predicting mode, the cumulative distribution function (CDF) is computed as:

$$CD{F}_k=\sum _{i=1}^k{P}_i$$

This cumulative probability array defines the random digit assignment intervals, mapping uniformly distributed random values to discrete grade outcomes

The fully processed probabilistic dataset is presented in Table 4 and Table 5, and Table 6, which include grade frequencies, probabilities, cumulative probabilities, and random digit assignments.

4.4. Predicting the Students’ Performance

The prediction stage, illustrated in Figure 3, applies the probabilistic model constructed in the previous steps. Pseudo-random numbers uniformly distributed over the interval [0,1] were generated using Python’s built-in random number generator.

Each generated random value was compared against the predefined random-digit intervals obtained from the empirical cumulative distribution function (CDF). Based on this comparison, the value was mapped to its corresponding grade category, yielding a single predicted GPA. This process was repeated iteratively so that the total number of generated GPA values equaled the number of students in the test dataset, ensuring a one-to-one correspondence between predicted and observed student records during model evaluation.

In accordance with the prediction algorithm illustrated in Figure 3, the model was evaluated using three distinct test datasets representing different student populations: engineering students, medical students, and the entire student cohort. This experimental design enables an assessment of the model’s predictive behavior across academically diverse groups. Repeated application of this algorithm yielded a set of predicted GPAs and an overall simulated distribution of student performance outcomes. Table 7 presents representative samples of the predicted grades generated through this procedure, and Table 8 the final model predicted grade compared with the actual results.

4.5. Experimental Evaluation

This section evaluates the outcomes of the experiment by examining the agreement between the predicted and observed grade distributions produced by the proposed model. The evaluation is organized into two complementary components. First, global agreement is quantified using distribution-based error metrics and goodness-of-fit measures that capture the overall similarity between predicted and empirical distributions. Second, residual diagnostics are performed to investigate category-level deviations and identify localized prediction behavior.

4.5.1. Distributional Evaluation and Error Metrics

Agreement between predicted and observed distributions is quantified using count-based errors (RMSE and MAE) and distribution-level similarity measures. The Chi-square goodness-of-fit statistic tests whether predicted counts are consistent with observed counts. Distribution distances include KL divergence (Actual || Predicted), Jensen–Shannon divergence (JS), total variation distance (TV), and Bhattacharyya distance. Lower error/distance indicates closer agreement; for Chi-square, a larger p-value indicates a closer fit. Table 9 presents the distribution-level evaluation metrics (Monte Carlo vs. actual), while Figure 4 and Figure 5, and 6 visualize predicted and observed grade distributions for each cohort. Figure 7 and Figure 8, and 9 overlay plots that highlight distribution-shape agreement, and Figure 10 provides direct category-level comparisons in a bar chart.

Figure 4. Engineering Students: Predicted vs. actual grade distribution (overlay).

Figure 4. Engineering Students: Predicted vs. actual grade distribution (overlay).

Figure 5. Medical Students: Predicted vs. actual grade distribution (overlay).

Figure 5. Medical Students: Predicted vs. actual grade distribution (overlay).

Figure 6. Entire Students: Predicted vs. actual grade distribution (overlay).

Figure 6. Entire Students: Predicted vs. actual grade distribution (overlay).

Figure 7. Engineering Students: Predicted vs. actual grade distribution (bar chart).

Figure 7. Engineering Students: Predicted vs. actual grade distribution (bar chart).

Figure 8. Medical Students: Predicted vs. actual grade distribution (bar chart).

Figure 8. Medical Students: Predicted vs. actual grade distribution (bar chart).

Figure 9. Entire Students: Predicted vs. actual grade distribution (bar chart).

Figure 9. Entire Students: Predicted vs. actual grade distribution (bar chart).

Figure 10. Summary of key distribution metrics across cohorts.

Figure 10. Summary of key distribution metrics across cohorts.

4.5.2. Residual Analysis and Diagnostic Assessment

To ensure that the model has captured the underlying data patterns and that the remaining error is merely random noise, residuals were computed per grade category as:

$$e_i=y_i-{\widehat{y}}_i$$

Where:

$e_i$ is residual (The remaining error is merely random noise

$y_i$ is the actual value and

Figure 11. Engineering Students: Residual plot (Actual − Predicted).

Figure 11. Engineering Students: Residual plot (Actual − Predicted).

Figure 12. Engineering Students: Residual plot (Actual − Predicted).

Figure 12. Engineering Students: Residual plot (Actual − Predicted).

Figure 13. Entire Students: Residual plot (Actual − Predicted).

Figure 13. Entire Students: Residual plot (Actual − Predicted).

Table 10 presents the analysis results. Positive residuals indicate underprediction by the model, while negative residuals indicate overprediction. Residual plots provide a diagnostic view of where the model deviates from the observed distribution (see Figures 11,12, and 13).

Residual analysis is a critical diagnostic method in statistics and machine learning used to evaluate the appropriateness of a model by analyzing the residuals—the differences between observed values and values predicted by the model. It checks key assumptions like linearity, constant variance (homoscedasticity), and independence of errors. A random pattern in residual plots indicates a good fit, while structured patterns suggest model deficiencies [38].

Normality of residuals was assessed using the Shapiro–Wilk test [39]. Since all p-values exceed 0.05, the null hypothesis of normality cannot be rejected, implying that deviations are consistent with stochastic variation rather than structural model misspecification. Bias and balance evaluation based on residual sign symmetry and zero mean indicates calibrated probabilistic interval construction and stable stochastic sampling behavior.

Residual magnitudes are concentrated primarily in extreme and transition grade categories, which is expected in Monte Carlo sampling due to the higher variability associated with low-probability outcomes. No systematic clustering or directional drift was observed. These findings indicate that the empirical probability representation captures the underlying distribution adequately.

Overall, diagnostic assessment confirms statistical neutrality, absence of bias, and residual patterns consistent with well-behaved stochastic simulation models, reinforcing confidence in the proposed probabilistic forecasting framework.

4.5.3. Results Discussion

The Monte Carlo model captures the overall distributional structure across cohorts while exhibiting cohort-dependent deviations in specific grade categories. Engineering students show the largest absolute deviations in the A and W categories, which typically have lower frequency and higher relative variability. Medical students show pronounced deviations in the B, C, and D categories, reflecting distributional drift between academic years. For the entire cohort, aggregation amplifies mid-range category deviations (C and D), consistent with heterogeneous subgroup behavior. Residual diagnostics identify the dominant error contributors and can guide refinement of empirical probabilities and interval assignment.

The evaluation results demonstrate that the proposed model effectively captures the probabilistic structure of student performance outcomes. By relying on empirical distributions and stochastic sampling, the model provides a realistic representation of uncertainty that is often overlooked by deterministic prediction approaches.

These findings suggest that Monte Carlo–based modeling constitutes a viable and interpretable alternative for academic performance prediction, particularly in scenarios where uncertainty, variability, and limited data availability play a critical role.

5. Model Evaluation

To evaluate the proposed Monte Carlo simulation framework, a linear regression baseline will be adopted rather than more complex machine-learning algorithms because:

• The available datasets consist of aggregated grade-frequency counts rather than individual-level feature vectors; applying instance-level machine-learning classifiers to aggregated data risks introducing ecological inference bias, whereby relationships observed at the group level may not hold at the individual level [40].
• Regression at the aggregate level remains appropriate when the target of inference is group-level behaviour rather than individual prediction, as parameters estimated from aggregated observations may still provide valid insights about aggregate outcomes when interpreted consistently with the data level [41]. This ensures methodological alignment between the Monte Carlo model output (distribution predictions) and the comparator model output.
• Simple linear models are widely accepted as interpretable and computationally efficient baselines in predictive modelling, providing transparent parameterization and diagnostic insight before introducing additional complexity. Empirical evidence demonstrates that simpler linear approaches may rival or outperform more complex machine-learning techniques in certain applied settings while offering clearer interpretability and reduced computational cost. Moreover, rigorous evaluation practice emphasises the importance of selecting appropriate baseline models and metrics to avoid misleading conclusions regarding model generalizability and performance [42].

Taken together, these considerations justify the use of linear regression as a principled baseline for evaluating the Monte Carlo model. The regression benchmark provides an interpretable deterministic reference operating in the same aggregated output space, allowing performance improvements attributable to stochastic simulation to be assessed without introducing methodological inconsistencies or ecological inference bias.

5.1. Implementing the Linear Regression Prediction Model

In this section, a linear regression predictor will be implemented as a deterministic benchmark for the proposed Monte Carlo simulation model. Both approaches are evaluated in the same output space (cohort-level grade-frequency distributions). Using Table 2 (section 4) as the training dataset (Y₁) and Table 3 (section 4) as the test dataset (Y₂), the benchmark learns a mapping from Y₁ grade counts to Y₂ grade counts and predicts the expected distribution in the subsequent year. Figure 14 shows the algorithm used for implementing the Linear Regression Prediction Model.

For each cohort, an ordinary least squares (OLS) linear regression model is fitted at the grade-category level as:

$${\widehat{Y}}_2=\alpha Y_1+\beta$$

where $Y_1$ denotes the observed grade frequency in Year 1, ${\widehat{Y}}_2$ is the predicted grade frequency in Year 2, α is the slope coefficient, and β is the intercept. Parameters are estimated by minimizing the squared error over grade categories.

5.2. Results, Evolution, and Comparisons

Table 11 presents training (Y₁), test (Y₂), and model prediction results (Regression and Monte Carlo), and Figure 15 and Figure 17, and 18 visualizing these results

Figure 15. Engineering Students: Linear Regression vs Monte Carlo.

Figure 15. Engineering Students: Linear Regression vs Monte Carlo.

Figure 16. Medical Students: Linear Regression vs Monte Carlo.

Figure 16. Medical Students: Linear Regression vs Monte Carlo.

Figure 17. Entire Students: Linear Regression vs Monte Carlo.

Figure 17. Entire Students: Linear Regression vs Monte Carlo.

Prediction accuracy is measured using RMSE and MAE over grade counts. Goodness-of-fit is evaluated with R² and the Chi-square statistic (χ²) by comparing predicted and observed grade-frequency distributions. Lower RMSE/MAE and χ² values show closer agreement; higher R² indicates better explained variance. Chi-square p-values offer a statistical measure of distributional consistency. Table 12 prsents the performance results and Table 13 present linear regression parameters, while Figure 18 shows the error summary (RMSE and MAE) for Regression vs Monte Carlo across cohorts.

Regarding residual comparison, Table 14 presents the residuals by grade category (Actual − Predicted) for both Linear Regression and Monte Carlo models, while Figure 19 and Figure 20, and Figure 21 visualize this comparison.

5.3. Discussion

Across cohorts, the linear regression benchmark provides a simple deterministic mapping from Year-1 to Year-2 grade counts. However, regression tends to smooth distribution changes and can misrepresent category migration when year-to-year dynamics are nonlinear or stochastic. The Monte Carlo model, by contrast, is designed to preserve distributional structure through probabilistic sampling and thus can better capture shifts in the tails (e.g., A and W) and mid-range categories depending on cohort behavior.

Engineering cohort: Regression and Monte Carlo show comparable magnitude errors, with differences concentrated in extreme categories. Medical cohort: Monte Carlo better reflects distribution drift, particularly in categories with sharper year-to-year changes. Entire cohort: aggregation amplifies deviations; regression may fit the global trend but often underestimates variability, whereas Monte Carlo maintains distributional fidelity. Overall, the benchmark comparison supports the use of probabilistic simulation when cohort performance evolves stochastically and when preserving distribution shape is the primary modeling objective.

However, the experiment demonstrates that the proposed Monte Carlo–driven predictive framework achieves performance comparable to established statistical baselines while providing a more accurate representation of distributional uncertainty. This validates the methodological viability of integrating stochastic simulation into predictive modeling pipelines. The study, therefore, contributes an alternative modeling perspective that broadens the analytical toolkit for predicting academic performance and supports the future development of hybrid simulation–learning architectures.

6. Limitations

This study presents several limitations that frame the interpretation of the proposed modeling framework. The approach operates on aggregated grade-frequency data, which supports distribution-level forecasting but restricts incorporation of individual-level explanatory variables that could enhance predictive precision. Additionally, the Monte Carlo sampling process relies on empirical distributions derived from historical observations, assuming temporal stability; structural changes in academic environments may therefore require periodic recalibration to maintain model validity. The comparative evaluation is limited to a linear regression baseline due to data aggregation constraints, preventing broader benchmarking against feature-driven machine learning classifiers. Finally, as a stochastic simulation method, the model produces probabilistic estimates subject to sampling variability, implying that outputs should be interpreted as uncertainty-aware guidance rather than deterministic predictions. Collectively, these considerations define the operational scope of the proposed framework and highlight directions for future methodological enhancement.

7. Conclusions and Recommendations

This study developed and evaluated a Monte Carlo–based probabilistic framework for predicting cohort-level academic performance distributions and benchmarked its performance against a linear regression baseline. Monte Carlo captures uncertainty through stochastic sampling from empirical distributions, a methodology widely used for modelling variability and quantifying statistical outcomes in complex systems [11,43]. Experimental results demonstrated that the proposed model reproduces structural characteristics of observed grade distributions while allowing localized deviations consistent with stochastic variability. Such findings align with research showing that probabilistic simulation provides advantages over deterministic modelling when uncertainty and system variability must be explicitly represented [44].

The comparison with linear regression highlights complementary roles between deterministic and probabilistic approaches. Linear models provide interpretable mappings and are commonly employed as baseline benchmarks in predictive analysis [45,46]. However, deterministic scaling cannot fully represent nonlinear migration across grade categories. The Monte Carlo approach, by preserving distributional behaviour through sampling, demonstrated improved capability in capturing structural variability. Accordingly, this model should not be viewed as an alternative replacement to machine learning algorithms, but rather as a scientific complementary addition that enhances modelling capability in uncertainty-driven environments. Integration of stochastic simulation with data-driven learning frameworks has been identified as a promising research direction bridging statistical simulation and modern ML systems [47].

Based on these findings, several recommendations emerge. First, probabilistic simulation should be adopted when forecasting cohort-level outcomes where variability representation is critical. Second, linear regression should be retained as a transparent baseline benchmark for validating modelling improvements. Third, hybrid architectures that combine stochastic simulation with machine learning predictors should be explored to utilize the strengths of both paradigms. Finally, enriching datasets with behavioural or academic covariates may further enhance predictive accuracy and model generalizability [48].

In conclusion, the proposed Monte Carlo model constitutes a meaningful scientific contribution to predictive modelling research by extending beyond deterministic estimation toward uncertainty-aware forecasting. Its integration alongside traditional machine learning methods provides a reliable and interpretable pathway for modelling educational performance dynamics and other complex stochastic systems.