How to find the contributing probability of individual features in a machine learning model? In some machine learning models there is a predict_proba function for the overall class probability based on all of the features in the model. Is there a way to find the individual probability contribution of each feature? For example, if there is a 85% chance of raining based on three features: cloud coverage, humidity, and temperature. Is there a way to determine that dense cloud coverage by itself means there is a 60% chance of rain, high humidity means there is a 70% chance of rain, and high temperature means there is a 10% chance of rain, but these three features combined produce an 85% chance of rain? Looking for a general direction or materials to read. Thanks in advance.

You can look at what is known as the "marginal effect" of variables.  This doesn't quite give you what you want, because all marginal effects do not sum to the probability, but it is a way to say "how much does the probability change on average when this variable changes".  This is a technique that can apply to models more generally, though you typically see them used with GLMs like logistic regression because people want to do inference on the quantity.
Working with your example, let's say one of our variables is a binary indicator for "Is Cloudy Outside".  To compute the marginal effect of cloudiness, simple compute the average difference in predicted probabilities when all observations have the variable "Is Cloudy = 1" and when all observations have "Is Cloudy=0".
Mathematically, Let $f$ be a model which produces estimated probabilities $\hat{p}$ from a given set of data $\{ (\mathbf{x}_1, w_1, y_1), \dots, (\mathbf{x}_n,  y_n) \}$. Here, $w$ is the variable for which we want to compute the marginal effect, $\mathbf{x}$ are other variables in the model, and $y$ is the label.
The marginal effect would then be
$$  \dfrac{1}{n} \sum f(\mathbf{x}_i, w_i=1) - f(\mathbf{x}_i, w_i=0) \>.$$
This is easy to do in code.  Here is an example in R.  Let's first make a problem
set.seed(0)
N <- 10000
x <- rbinom(N, ,m 0.5)
w <- rbinom(N, ,m 0.5)
y <- rbinom(N, 1, 0.25 + 0.05*x + 0.075*w)

From the code, we already know the marginal effect of $w$, its 0.075. So on average, changing $w$ from 0 to 1 will yield a 7.5% increase in probability.  Let's fit a model to these data using fit <- glm(y~x+w, family = binomial()) and compute the marginal effect
p_w_1 <- predict(fit, newdata = list(x=x, w=rep(1,N)), type = 'response')
p_w_0 <- predict(fit, newdata = list(x=x, w=rep(0,N)), type = 'response')
effect = mean(p_w_1 - p_w_0)
effect
>[1] 0.06910256


Our estimated marginal effect of the variable $w$ is ~7%. This approach lines up with what libraries would give us
library(marginaleffects)
library(tidyverse)
marginaleffects(fit) %>% 
  summary()
 Term  Effect Std. Error z value   Pr(>|z|)   2.5 %  97.5 %
1    x 0.04832   0.009159   5.276 1.3218e-07 0.03037 0.06627
2    w 0.06895   0.009130   7.552 4.2956e-14 0.05105 0.08684

Model type:  glm 
Prediction type:  response

You could apply these approaches to any model you like to get an estimated marginal effect.
