Overview
A real estate agency wants a formula for what a house is worth. They hand over 506 homes around Boston, each described by eleven numbers — the crime rate on the street, the rooms under the roof, how far it is to work, what the neighbours pay in tax — and the median value the market actually put on them.
Fitting a linear model to that takes four lines. The rest of the project is the part that matters: working out whether the numbers that come out mean anything. I split the data 90/10 before looking at it, fitted on 455 homes and kept 51 sealed until the end. The model explains 71.8% of the variation and misses by about $3,700 on homes it had never seen, against a spread of $9,200 in the prices themselves — roughly five times better than answering "the average" to every question.
Then the bootstrap: refit the model a thousand times on resampled data and watch which coefficients hold still. Two of them do not. The proportion of industrial land and the share of pre-war housing both wander across zero, which is a statistician's way of saying the data cannot tell you whether they matter at all. A regression will always hand you eleven numbers — it will not volunteer that two of them are noise.
What I built
Describing the data before modelling it
- Histogram, box plot and violin plot of the target variable, and a scatter of price against the average number of rooms — the one relationship visible without a model.
- The finding that shapes everything afterwards: every home worth more than $50,000 was recorded as exactly $50,000, so the top of the distribution is a ceiling in the data rather than a feature of the market.
The model, and how much of it to trust
- An OLS regression on the 455 training homes, with the intercept and the coefficient on room count read back in plain language: about $3,690 per additional room, holding the other ten variables still.
- A 1,000-iteration bootstrap for the confidence intervals, testing each coefficient by whether zero falls inside it. Industrial land lands at (−0.084, +0.142) and fails; room count lands at (+2.21, +5.44) and holds.
- Evaluation on the held-out tenth: MSE 13.61, RMSE $3,689, against a baseline of 73.71 for predicting the mean every time.
The fitted model, its coefficients and the fifty-one held-out homes are all in the block above — it runs in the browser, on the real numbers.
The model, running
Coefficients from project3.ipynb
This really is the fitted regression, not a demonstration of one: the eleven β values below were printed by statsmodels on the training split and copied over unchanged. Price a house with them, read what the fit believes, then watch it miss.
Estimated median value
$23,236
Exactly the median house of the dataset.
What is moving the price
Every slider is sitting on the median of the dataset, so there is nothing to push the price either way. Move one.
Boston, 1978, 506 homes. A teaching dataset with a well-known history — including a variable on the racial composition of each town, which is excluded here, as it was in the assignment.
Inside the notebook
Read-only · outputs from project3.ipynb
Five cells from the assignment: the split and the shape of the prices, the one relationship you can see by eye, the regression itself, a thousand bootstrap refits deciding which coefficients survive, and the test set having the last word.
data = pd.read_csv("housing_data.csv")
# The split has to be reproducible: seed it with my student number
student_number = 4914
np.random.seed(student_number)
train_indices = np.random.choice(len(data), int(0.9 * len(data)), replace=False)
test_indices = [i for i in range(len(data)) if i not in train_indices]
train_data, test_data = data.iloc[train_indices], data.iloc[test_indices]
…
plt.hist(data['medv'], bins=20, edgecolor='black', color='yellow')
plt.title('Distribution of Housing Prices (medv)')
plt.show()plt.scatter(data['rm'], data['medv'], color='blue', alpha=0.5)
plt.title('Scatter Plot of Housing Prices vs. Average Number of Rooms')
plt.xlabel('Average Number of Rooms (rm)')
plt.ylabel('Housing Prices (medv) in thousands of USD')
plt.show()X_train = sm.add_constant(train_data[['crim','zn','indus','nox','rm','age',
'dis','rad','tax','ptratio','lstat']])
y_train = train_data['medv']
model = sm.OLS(y_train, X_train).fit()
print(model.summary())n_iterations = 1000
bootstrap_coeffs = []
for _ in range(n_iterations):
sample_indices = np.random.choice(train_indices, size=len(train_indices), replace=True)
resampled_model = sm.OLS(y_train.loc[sample_indices], X_train.loc[sample_indices]).fit()
bootstrap_coeffs.append(resampled_model.params)
bootstrap_coeffs = np.array(bootstrap_coeffs)
# Is 0 inside the interval? Then the variable is not telling us anything.
for i, name in [(3, "indus"), (5, "rm")]:
lo = np.percentile(bootstrap_coeffs[:, i], 2.5)
hi = np.percentile(bootstrap_coeffs[:, i], 97.5)
print(f"95% CI for {name}: ({lo:.4f}, {hi:.4f})")Y_pred = model.predict(sm.add_constant(test_data[cols]))
mse = np.mean((y_test - Y_pred) ** 2)
print(f"Mean Squared Error (MSE): {mse}")
print(f"Root Mean Squared Error (RMSE): {np.sqrt(mse)}")
print(f"Standard Deviation of MEDV: {data['medv'].std()}")
# The honest baseline: predict the mean, every time
baseline_mse = np.mean((y_test - y_test.mean()) ** 2)
print(f"Baseline MSE: {baseline_mse}")