For the following, we will assume normality of residuals, i.i.d data and so on. These are not trivial assumptions and should always be tested, but for our discussion, let's just say we know it's nice, Gaussian data.
Question 1: No, you cannot estimate such probabilities directly from a model as there is no such output. However, let's keep in mind that when estimating a Y for specific X1 and X2 inputs on a linear regression, we are creating a conditions distribution. Something to help visualize this:
So, what we're trying to do is this. We have a linear regression, which gives us conditional distributions of Y dependent on values of X. Now, given this conditional distribution, you want to figure out the probability that Y falls within a certain range. This is pretty straightforward now, because you have the distribution and therefore the rough PDF (probability density function), so you can figure out whatever probabilities you'd like. Of course, be careful and take it with a grain of salt - you are estimating uncertainty with many assumptions, finite data etc.
In R, lm() objects have various parameters/outputs. I'm sure you're familiar with using predict() on an lm() object to generate a prediction on new data using an existing linear regression model. That single prediction corresponds to the mean of the conditional distribution at those X values. I hope this makes sense. The thing you're missing is the standard deviation of the (we've assumed) Gaussian distribution of data points around the mean of that conditional distribution. This "standard deviation" we're looking for is the "Residual Standard Error", and can be found by using the summary() function on your lm() object.
So you've assumed normality, and you have a mean and standard deviation for the distribution of points at a particular value of X. From here, just use these values in R or in some quick, free online tool to calculate out the PDF of your normal distribution, and then apply whatever values you're interested in seeing the probability for. For example, you determine that at X1=3 and X2=5, the point prediction of the lm() (the mean) is 5.6, with Residual Standard Error of 2. You can now easily figure out, given this distribution, what the probability is of a value between 3 and 6, or whatever you're interested in.
Question 2:
As far as I know, the way to get "uncertainty" metrics from XGBoost is less straightforward and more theoretically debated. If you're doing classification, it's easy because XGBoost outputs a class probability. But for regression, honestly I'm not sure. You might consider training multiple models and then make a mini-distribution of their predictions to get a sense of the uncertainty/disagreement.
If you use RandomForest, it's a bit more straightforward because in RF, the trees are independent, so you could (instead of generate entire new models), look at the distribution of predictions from individual trees to get a sense of the uncertainty in the prediction.
I say "uncertainty" because nice, clean statistical assumptions are likely to start breaking down, and you probably won't get super clean normal distributions in the outputs. So to get a "probability estimate" of a certain value, you might have to resort to just looking at the outputs of the trees and saying "okay, 95% of the trees estimated values over my X threshold, so I will say it's a 95% probability of X threshold being met". Maybe someone can add some comments or responses here, because I'm pretty certain that the logic I just described is theoretically/mathematically incorrect, in the sense that you probably can't make those kinds of assumptions based on the outputs of RF. I think from a practical point of view, I wouldn't bother trying to make specific probability predictions for regression from RF. Maybe just "the prediction is very uncertain here, so it's likely to be a much greater range of values, and the range is probably reasonably reflected by the range of outputs from individual trees".
I hope this helps!
