How to create a bivariate normal distribution I have a data set like this:
df

Income     Education_in_years
40,000     10
50,000      9
70,000     12
30,000      5
100,000    20

I would like to create a bivariate distribution from this and try to guess probability of income given eduction in years.
I can build the linear model as follows:
lin <- lm(Income~Eduction_in_years, data=df)

I could come up with a formula like this:
Income = a*education_in_years + e

What I would like to do is, given the model, create a bivariate normal, and run a 'what-if' analysis to determine income given eduction_in_years.
Can somebody walk me through how to create a bivariate normal distribution from this dataset and determine the income level given the education?
 A: Simple linear regression is based on a bivariate normal model. However, in practice, we assume the independent variables are fixed/known and we are finding the conditional distribution of Y at each X.  This is what you are looking for, so just apply simple regression and you will get your answer. However, if you would like to assume the X and Y are random and follow a joint normal distribution, then you need to do some more work.
Have you taken a look a the Bivariate Normal density? It specifically has parameters for the marginal means, variances, and the correlation between the two variables. You can estimate each of these quantities from your current data to get an estimated joint normal distribution. For a more sophisticated approach, you can create a prediction region for the next observation of X and Y and see what prediction probability the formula assigns to different . This is not as simple as just "plugging in" your estimates, but it will account for the error in using estimated parameters, which yhou are not doing by just plugging in. This paper has a formula for multivariate prediction regions, see equation (2.2)
