Monte Carlo methods are computational algorithms that rely on repeated random sampling to solve problems that might be deterministic.

They usually follow a particular pattern:

1. Define a domain of possible inputs.
2. Generate inputs randomly from a probability distribution over the domain.
3. Perform a deterministic computation of the outputs.
4. Aggregate the results.

We know that the probability of getting an even number on a die is:

\[P(Even) = \frac{Number\ of\ even\ outcomes}{Total\ number\ of\ outcomes} = \frac{3}{6} = 0.5\]

Monte Carlo Approach:

1. Simulate dice rolls N times.
2. Count the number of even rolls X.
3. Estimate probability: \[ \hat{P} = \frac{X}{N} \]

For our experiment, we will simulate a random person’s characteristics like age, sex, resting blood pressure, cholesterol level, etc. and determine how likely it is for them to have heart disease. We will he using the “heart_cleveland_upload.csv” database.

First, we will fit a logistic regression model using our dataset to see how each variable actually relates to heart disease.

model = glm(condition ~ age + sex + cp + trestbps + chol + fbs +
               restecg + thalach + exang + oldpeak + slope + ca + thal,
             data = heart, family = binomial)

We will use the following R code to generate a plot for the coefficients we got. Longer bars to the right increase the likelihood of heart disease, while longer bares to the left decrease it.

coefs <- data.frame(variable = names(coef(model)), estimate = coef(model))

ggplot(coefs, aes(x = variable, y = estimate)) + geom_col() + 
  coord_flip() + labs(x = "Variable", y = "Coefficient Estimate")

Now we will use the following code to generate a simulation for the probability of 20,000 people having heart disease using the Monte Carlo Method.

sim = 20000
diseased = 0

for(i in 1:sim){
  randAge = sample(heart$age, 1, replace = TRUE)
  randSex = sample(heart$sex, 1, replace = TRUE)
  randCp = sample(heart$cp, 1, replace = TRUE)
  randTrest = sample(heart$trestbps, 1, replace = TRUE)
  randChol = sample(heart$chol, 1, replace = TRUE)
  randFbs = sample(heart$fbs, 1, replace = TRUE)
  randEcg = sample(heart$restecg, 1, replace = TRUE)
  randThalach = sample(heart$thalach, 1, replace = TRUE)
  randEx = sample(heart$exang, 1, replace = TRUE)
  randOld = sample(heart$oldpeak, 1, replace = TRUE)
  randSlope = sample(heart$slope, 1, replace = TRUE)
  randCa = sample(heart$ca, 1, replace = TRUE)
  randThal = sample(heart$thal, 1, replace = TRUE)
  
  stats = coefs$estimate[1] + coefs$estimate[2] * randAge + 
          coefs$estimate[3] * randSex + coefs$estimate[4] * randCp +
          coefs$estimate[5] * randTrest + coefs$estimate[6] * randChol +
          coefs$estimate[7] * randFbs + coefs$estimate[8] * randEcg +
          coefs$estimate[9] * randThalach + coefs$estimate[10] * randEx +
          coefs$estimate[11] * randOld + coefs$estimate[12] * randSlope +
          coefs$estimate[13] * randCa + coefs$estimate[14] * randThal
  
  prob = 1 / (1 + exp(-stats))
  
  if(prob > 0.5){
    diseased = diseased + 1
  }
}
probHeart = diseased/sim
print(probHeart)

## [1] 0.462

flush.console()

Now we will use the formula \[ P = \frac{Number\ of\ people\ with\ heart\ conditions}{Total\ number\ of\ people}\] to find out the actual probability that someone has heart disease.

## [1] 0.4612795

As we can see, this value 0.46128 is quite close to the one we obtained by running a Monte Carlo Simulation.

We have demonstrated that Monte Carlo methods can accurately estimate real-world probabilities. Our simulation of 20,000 random patients yielded a heart disease probability of approximately 0.46, which closely matches the actual rate of 0.46128 in our dataset. This validates both our Monte Carlo approach and the Law of Large Numbers.

Here are some more interesting plots related to our data. On the left, it shows how age differs between people with and without heart disease. On the right, we can see how the 3 key variables – Age, Cholesterol, and Blood Pressure – relate to heart disease: