2026-10-17
In brevity, Monte Carlo Simulation is the usage of historical data to predict the probabilities of future data across multiple, intersecting variables.
It is especially useful in fields like:
Basically, any data that has a lot of unpredictable variables.
To show the usefulness of Monte Carlo Simulation, we’ll be using a modified version of the iris dataset from the “datasets” package from RStudio. To simplify things, we’ll be reducing the dataset to three random variables.
Pretend that a flower’s value is determined solely by the dimensions of its petal and the width of its sepal. Specifically, the product of these three dimensions determine the flower’s value.
Assuming we continue using the same batch of seeds from the dataset of flowers that were grown, how do we determine the approximate profit margin of each batch of 150 flowers?
Btw, you can scroll vertically on any slide with text overflow.
Firstly, we’ll start by trimming down the original iris dataset into the random variables:
We’ll use the dlpyr library to make things simple. NA values (there should be none) are omitted.
Code:
suppressPackageStartupMessages(library(dplyr))
data("iris")
df = tibble(PW = iris$Petal.Width,
PL = iris$Petal.Length,
SW = iris$Sepal.Width) %>% na.omit()
We’ll also add in a new column that displays the product of the previous three columns. This’ll be a visual indicator for now, but it gives us a sense of each flower’s value.
It’s a basic formula, something like this: \[ (PetalWidth)(PetalLength)(SepalWidth) = FlowerPrice \]
Code:
df = mutate(df,
FP = round((PW * PL * SW), 2))
# df will be used for further processing due to shorter variable names.
df_present = rename(df,
"Petal Width" = PW,
"Petal Length" = PL,
"Sepal Width" = SW,
"Flower Price" = FP)
The following shows a few columns from the new modified iris dataset:
# A tibble: 6 × 4
`Petal Width` `Petal Length` `Sepal Width` `Flower Price`
<dbl> <dbl> <dbl> <dbl>
1 0.2 1.4 3.5 0.98
2 0.2 1.4 3 0.84
3 0.2 1.3 3.2 0.83
4 0.2 1.5 3.1 0.93
5 0.2 1.4 3.6 1.01
6 0.4 1.7 3.9 2.65
Not bad, but we can see that the flower price can range from below a dollar to over two dollars.
With the preprocessing finished, let’s first visualize the three dimensions for each flower.
A colorscale is used to show the relative price of each flower; lighter colors are higher priced.
We’ll use the Plotly library to achieve this.
Code:
suppressPackageStartupMessages(library(plotly))
xax = list(
title = "Petal Width",
titlefont = list(family = "Noto Sans")
)
yax = list(
title = "Petal Length",
titlefont = list(family = "Noto Sans")
)
zax = list(
title = "Sepal Width",
titlefont = list(family = "Noto Sans")
)
# Referenced plotly.com's documentation on 3d scatter plots
fig = plot_ly(df, x = ~PW, y = ~PL, z = ~SW,
type = "scatter3d",
mode = "markers",
marker = list(color = ~FP,
colorscale = "Viridis",
showscale = TRUE)) %>%
add_markers() %>%
layout(scene = list(xaxis = xax,
yaxis = yax,
zaxis = zax),
annotations = list(
x = 1,
y = 1,
text = "Relative Flower Price",
xref = "paper",
yref = "paper",
showarrow = FALSE),
showlegend = FALSE)
fig
[Backup Image of 3D Plotly] 
The 3D plot shows a more comprehensive overview of the flower dataset, especially with the range of flower prices going from below ten dollars to over 50 dollars. That is an especially significant range.
However, with the major spread of flower prices, it is still difficult to determine how much profit a batch of 150 flowers will yield.
As such, we’ll use Monte Carlo simulation to get a probable range of profit margins.
Before we start, we need to obtain the mean and standard deviation (sample) of each random variable.
To simplify things, we’ll use two vectors: one for mean and another for standard deviation.
In case you forgot, here are the formulas for mean and standard deviation:
Mean: \[ \bar{x} = \frac{\sum^n_{i=1}x_i}{n} \]
Standard Deviation: \[ s = \sqrt{\frac{\sum(x_i-\bar{x})^2}{n-1}} \]
Code for Prepping Monte Carlo:
means = c(mean(df$PW), mean(df$PL), mean(df$SW)) stddev = c(sd(df$PW), sd(df$PL), sd(df$SW)) means
[1] 1.199333 3.758000 3.057333
stddev
[1] 0.7622377 1.7652982 0.4358663
This code calculated the mean and standard deviation of each flower dimension, both in the order of Petal Width, Petal Length, and Sepal Width.
Decent numbers for parameters.
Let’s simulate 100 seed batches, that’s 15,000 flowers.
For each random variable, we’ll use the mean and add/subtract by a certain value depending on a RNG function and the standard deviation of the variable.
This will obtain a random, probable value for each random variable for all 15000 flowers.
Code:
set.seed(777) # Seed for consistency; jackpot number
n = 15000
PW_trial = rnorm(n, means[1], stddev[1])
PL_trial = rnorm(n, means[2], stddev[2])
SW_trial = rnorm(n, means[3], stddev[3])
df_mt = tibble(PW_MT = PW_trial,
PL_MT = PL_trial,
SW_MT = SW_trial)
We’ll use ggPlot2 to see how the trials turned out.
Code:
suppressPackageStartupMessages(library(ggplot2))
mt_P = ggplot(data = df_mt,
aes(
x = PW_MT,
y = PL_MT
)
) +
geom_point() +
theme_bw() +
labs(x = "Petal Widths",
y = "Petal Lengths",
title = "Monte Carlo Simulation of Petal Area")
That is definitely a pretty normal distribution, however, there are unrealistic negative values. Negative values will be trimmed out using this code:
df_mt[df_mt < 0] = NA df_mt = na.omit(df_mt)
We’ll use the same mutate function to obtain the flower prices for the monte carlo simulation.
df_mt = mutate(df_mt,
FP_mt = round((PW_MT * PL_MT * SW_MT), 2))
Using the Monte Carlo simulation, we can see the comprehensive overview of flower prices.
Code:
mt_Prices = ggplot(df_mt,
aes(
x = FP_mt
)
) +
geom_histogram(binwidth = 5,
fill = "lemonchiffon",
aes(
y = after_stat(count / sum(count))
)
) +
labs(x = "Flower Prices ($)",
y = "Chance of Price (%)",
title = "Estimated Flower Prices Using Monte Carlo") +
scale_y_continuous()
This is quite the telling plot. A major bulk of flowers will be worth under $25, meaning that the unit cost of growing each flower should be under $25 if we hope to profit off these flowers.
To actually make things useful for a business owner, it would be better if the histogram displayed revenue per batch of seeds instead of individual flower prices.
df_mt = mutate(df_mt,
Profits = FP_mt * 150
)
Code:
ProfitCall = ggplot(df_mt,
aes(
x = Profits
)
) +
geom_histogram(binwidth = 1024,
fill = "gold",
) +
labs(x = "Revenue Per Batch ($)",
y = "Number of Batches",
title = "Estimated Batch Revenues Using Monte Carlo") +
scale_y_continuous()
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0 1002 1936 2282 3195 13268
Looking at the histogram and summary output, the results are quite disappointing, well, depending on how you look at it.
Just for fun, here’s the 3D plot from earlier but with Monte Carlo Simulation.
Code:
library(plotly)
xax = list(
title = "Petal Width",
titlefont = list(family = "Noto Sans")
)
yax = list(
title = "Petal Length",
titlefont = list(family = "Noto Sans")
)
zax = list(
title = "Sepal Width",
titlefont = list(family = "Noto Sans")
)
# Referenced plotly.com's documentation on 3d scatter plots
fig = plot_ly(df_mt, x = ~PW_MT, y = ~PL_MT, z = ~SW_MT,
type = "scatter3d",
mode = "markers",
marker = list(color = ~Profits,
colorscale = "Viridis",
showscale = TRUE)) %>%
add_markers() %>%
layout(scene = list(xaxis = xax,
yaxis = yax,
zaxis = zax),
annotations = list(
x = 1,
y = 1,
text = "Batch Revenue ($)",
xref = "paper",
yref = "paper",
showarrow = FALSE),
showlegend = FALSE)
fig
[Backup Image of 3D Plotly] 
It’s a pretty nice change from the first 3D plot, now showing that the highest possible batch price stems from two moderately high Petal & Sepal Widths but with a very high Petal Length. However, it’s important to notice that obtaining the max value for any dimension will not yield a high batch revenue if either of the two other dimensions are low.
I hope this little showcase showed how useful Monte Carlo can be.
I mostly used this article to define Monte Carlo, but I’m sure you can find better sources if you decide to delve into this rabbit hole