Project: A Statistical Profile of a City's Housing Prices

What you're building

This is the last stop in the Statistics From Zero for Data Work path, and it works a little differently from the eleven articles before it. There is no new statistical idea to introduce here. Instead, you are going to take everything from data types through sampling bias and use it on one dataset, the way a working analyst would: load the data, describe it carefully, question your own summary numbers instead of trusting them blindly, and write up what you found in plain language. By the end you will have produced something you could genuinely hand to someone else - a short, honest statistical profile of a housing market.

The dataset is a set of 500 home sales in a fictional city called Millbrook, spread across five neighborhoods. It is generated with code, not scraped from a real listings site, and that is deliberate. Real housing data needs cleaning, licensing checks and verification before anyone could responsibly publish numbers from it, and this is not a data-cleaning tutorial. Generating the data ourselves has a real teaching advantage too: because we write the rules that produce the prices, we know the ground truth. When a statistic tells us something, we can check whether that something is actually true of how the data was built. That is the fastest way to learn when to trust a summary number and when to question it.

Everything that follows works just as well on a real dataset of home sales, salaries, delivery times, hospital stays, or anything else where one numeric column is the thing people care about and several other columns might explain it. Treat Millbrook as a stand-in for whatever dataset you eventually point this process at.

Here is the route we will take, in order:

• Load the data and classify every column by data type (Part 2)
• Find the centre of the price column with mean, median and mode, overall and by neighborhood (Part 3)
• Measure spread with range, variance, standard deviation and IQR (Part 4)
• Compute quartiles and draw box plots to compare neighborhoods (Part 5)
• Look at the shape of the price distribution and measure its skew (Part 6)
• Check how well the normal distribution's 68-95-99.7 rule fits the data (Part 7)
• Find outliers with the IQR rule and the z-score rule, and decide what to do with them (Part 8)
• Measure correlation between price and the other numeric columns, and spot a misleading one (Part 9)
• Draw a sample from the data and see what sampling bias does to the picture (Part 10)
• Write the findings up as a short statistical profile

Setting up: tools and versions

You need three libraries: pandas and numpy for the calculations, and matplotlib for three plots. This article assumes pandas 2.x, numpy 1.24 or later and matplotlib 3.7 or later, run on Python 3.10 or later. Almost everything shown here has worked the same way in pandas since at least version 1.0, so a somewhat older install will very likely behave the same; the two methods worth double-checking on an old version are Series.quantile and Series.skew, both of which have been part of pandas for a long time.

Run the code blocks in the order they appear, in the same notebook or script, because each one uses variables created by the one before it. If you want your own run to match the reasoning in this article as closely as possible, keep the random seeds exactly as written. Changing a seed will change the specific numbers you see - a different mean, a different set of flagged outliers - but it will not change the general patterns this article talks about, because those patterns come from how the data-generating rules were written, not from the particular random draw.

Build the dataset: five neighborhoods in the fictional city of Millbrook

Millbrook is invented for this project. It has five neighborhoods with different typical price levels, and 500 recorded home sales with the kind of columns a real listings dataset would have: where the home is, how many bedrooms and bathrooms it has, its size, when it was built, what condition it is in, and what it sold for. We will build the price itself out of the other columns plus some random noise, and we will deliberately slip in a handful of very expensive outlier sales, because real housing markets have those too and we need some in our data to practise finding them.

Here is the full column list, with the data type we are aiming to give each one - keep this list nearby, because Step 1 will ask you to justify it from scratch:

• neighborhood - nominal categorical (no natural order)
• condition - ordinal categorical (fair < good < excellent)
• bedrooms - discrete numeric
• bathrooms - discrete numeric
• size_sqft - continuous numeric
• year_built - discrete numeric
• price - continuous numeric, the column we care most about

import numpy as np
import pandas as pd

np.random.seed(42)
n = 500

neighborhoods = ["Riverside", "Oak Hill", "Downtown", "Maple Park", "Westgate"]
neighborhood_weights = [0.25, 0.20, 0.15, 0.25, 0.15]
neighborhood = np.random.choice(neighborhoods, size=n, p=neighborhood_weights)

# typical price per square foot in each neighborhood (made up, for this fictional city)
price_per_sqft = {
    "Riverside": 220, "Oak Hill": 180, "Downtown": 260,
    "Maple Park": 150, "Westgate": 200,
}

condition = np.random.choice(["fair", "good", "excellent"], size=n, p=[0.2, 0.5, 0.3])
condition_multiplier = {"fair": 0.90, "good": 1.00, "excellent": 1.15}

bedrooms = np.random.choice([1, 2, 3, 4, 5], size=n, p=[0.05, 0.25, 0.35, 0.25, 0.10])
bathrooms = np.clip(np.round(bedrooms / 2 + np.random.normal(0, 0.5, n)), 1, 4).astype(int)

base_size = 500 + bedrooms * 280
size_sqft = base_size + np.random.normal(0, 150, n)
size_sqft = np.clip(size_sqft, 400, None).round(0)

year_built = np.random.randint(1950, 2024, size=n)

Walk through this slowly, because every line is one of the data-type decisions from Part 2 made concrete. neighborhood is drawn from a fixed list with probabilities - that makes it nominal categorical: Downtown is not "more" than Riverside, the categories just happen to have different typical prices in our pricing table below. condition is drawn from three ordered labels, which makes it ordinal: excellent really is better than good, which is better than fair, even though we have not yet turned those labels into numbers. bedrooms and bathrooms come from np.random.choice over whole numbers, which is what makes them discrete - a home cannot have 2.5 bedrooms. size_sqft starts as a count of square feet tied to the number of bedrooms, then has normal noise added, which is why it ends up continuous (and we round it only because floor plans are not usually quoted to the millimetre). year_built is an integer year, discrete in the sense that it only takes whole values, though we will treat it more like an ordinary number later since arithmetic on years is meaningful.

base_price = np.array([
    size_sqft[i] * price_per_sqft[neighborhood[i]] * condition_multiplier[condition[i]]
    for i in range(n)
])

# multiplicative noise: keeps price positive and gives it a mild right skew
noise = np.random.lognormal(mean=0, sigma=0.12, size=n)
price = base_price * noise

# a handful of unusually expensive sales - waterfront estates, penthouse conversions,
# the kind of listing that shows up in every real market and pulls the mean around
luxury_idx = np.random.choice(n, size=6, replace=False)
price[luxury_idx] = price[luxury_idx] * np.random.uniform(2.5, 4.0, size=6)

housing = pd.DataFrame({
    "neighborhood": neighborhood,
    "condition": condition,
    "bedrooms": bedrooms,
    "bathrooms": bathrooms,
    "size_sqft": size_sqft,
    "year_built": year_built,
    "price": price.round(0),
})

print(housing.shape)
print(housing.head())

The price formula multiplies size by a neighborhood rate and a condition multiplier, then multiplies the whole thing by random noise drawn from a lognormal distribution instead of adding noise the way we did for size. This matters: additive noise can push a value below zero, which makes no sense for a price, while multiplicative noise centred near 1 only ever stretches or shrinks the base price, keeping it positive and giving the whole column a gentle right-leaning shape even before we do anything else. Then we pick six listings at random out of the 500 and inflate them by two and a half to four times their computed price. Those six are standing in for the genuinely unusual sales that show up in almost any real housing dataset - the mansion, the odd legal dispute, the property with development rights attached - and we are putting them there on purpose so that later sections have real outliers to find instead of having to imagine them.

Run the print statements and you should see housing.shape report (500, 7), and housing.head() show the first five rows with five text/number columns of neighborhood and condition names next to small integers for bedrooms, bathrooms and year, and larger numbers for size and price. The exact neighborhood names and prices in those first five rows depend on the random seed, but the shape and the column layout will always look like this.

Step 1: confirm the data types before you calculate anything

Part 2 of this path made the case that the first question to ask about any column is not "what's the average" but "what kind of thing is this, and which statistics are even meaningful for it". It is worth doing that formally here, because it is the step people skip under time pressure, and skipping it is exactly how a mean gets computed on something that should never have been averaged.

print(housing.dtypes)

pandas will report neighborhood and condition as object (plain text), bedrooms, bathrooms and year_built as integers, and size_sqft and price as floating point numbers. That is pandas describing storage, not meaning - it has no idea that condition has a natural order and neighborhood does not. That distinction lives in your head and in how you choose to encode the columns, so let's make it explicit.

condition_order = ["fair", "good", "excellent"]
housing["condition"] = pd.Categorical(housing["condition"], categories=condition_order, ordered=True)

print(housing["condition"].dtype)
print((housing["condition"] > "fair").sum())

Converting condition to an ordered pandas Categorical does two useful things. First, it documents, inside the dataframe itself, that excellent is better than good is better than fair - anyone who loads this file later and checks the dtype will see that order instead of having to guess it from the label text. Second, it makes comparisons like housing["condition"] > "fair" work correctly, counting every good and excellent home, which a plain text column would refuse to do (or would do wrong, comparing strings alphabetically, where "excellent" sorts before "fair" for entirely the wrong reason).

neighborhood stays as plain text, and that is the right call, not an oversight. It is nominal: there is no meaningful way to say Riverside is "more" than Oak Hill. If you ever find yourself converting a nominal column to 0, 1, 2, 3, 4 so that a function will accept it, stop - a few sections from now, in the correlation step, we will see exactly why that habit produces numbers that look precise and mean nothing.

bedrooms, bathrooms and year_built are discrete numeric columns: you cannot have 2.7 bedrooms, and a home was built in one specific year, not somewhere between two years. size_sqft and price are continuous: in principle they could take any value in a range, and we only see rounded numbers because that is how homes are measured and listed. This distinction is what will tell us, later, that a mode is a genuinely useful summary for bedrooms but close to meaningless for price.

Step 2: find the centre of the price column

A hand-worked warm-up

Before touching the full 500-row dataset, it helps to work through a tiny version of the same question by hand, so you can check the computer's answer against your own arithmetic later. Suppose seven homes sold for these prices, already sorted: 310,000, 325,000, 340,000, 355,000, 362,000, 390,000, and 1,200,000 - that last one is our stand-in for a luxury outlier.

The mean adds everything up and divides by the count. The sum is 3,282,000, and dividing by 7 gives a mean of about 468,857. The median, since there are seven values, is simply the middle one once they are sorted - the fourth value, 355,000. Notice the gap: the mean sits roughly 114,000 above the median, and the entire reason is the single outlier at the end. Remove that one value and the mean of the remaining six drops to 345,333, much closer to their median of about 350,500. This is the same lesson Part 3 built around: the mean listens to every value, including the extreme ones, while the median only cares about position in the sorted order. For a column like price, which almost always has a few unusually expensive outliers, that difference is not a technicality - it changes which number you should lead with when you describe "the typical home".

There is no meaningful mode in that list of seven prices, because no value repeats. That is normal for continuous data recorded to the nearest dollar - the mode is far more useful for a column like bedrooms, where only a handful of distinct values exist and some of them genuinely occur far more often than others.

Now the full dataset

print(housing["price"].mean())
print(housing["price"].median())
print(housing["bedrooms"].mode())

Run that and look at the gap between the first two numbers. By construction, this dataset has six listings inflated to two and a half to four times their computed price, out of 500 - a little over one percent of the rows. That is a small fraction, but because the mean is sensitive to extreme values and those six are extreme by design, you should see the mean sitting noticeably above the median, the same direction as the hand-worked example, just with a smaller gap in percentage terms because the outliers are a smaller share of the data here than they were in our toy example of seven.

The bedrooms mode is a case where you can predict the answer before running the code and then check your prediction. Bedrooms were drawn with probabilities 0.05, 0.25, 0.35, 0.25 and 0.10 for one through five bedrooms - three-bedroom homes were given the single highest probability. With 500 draws, the law of large numbers makes it very likely, though not mathematically certain, that 3 comes out as the most frequent value. If your run shows something else, it is not a bug; it is a reminder that "most likely category" and "guaranteed winner" are not the same thing, especially with a sample size in the hundreds rather than the thousands.

Breaking the centre down by neighborhood

price_by_neighborhood = housing.groupby("neighborhood")["price"].agg(["mean", "median", "std", "count"])
print(price_by_neighborhood.sort_values("median", ascending=False))

This groupby splits the price column into five smaller columns, one per neighborhood, and computes mean, median, standard deviation and count for each. Because we set Downtown's price-per-square-foot rate to 260, the highest of the five, and Maple Park's to 150, the lowest, the gap between the rates (110 dollars per square foot) is large compared with the noise we added (roughly plus or minus 12 percent from the lognormal term, on top of ordinary variation in size). You should see Downtown sitting at the top of this table and Maple Park at the bottom, with Riverside, Westgate and Oak Hill somewhere in between in that order, because that is the order of their price-per-square-foot rates. If one neighborhood's std column looks dramatically larger than the rest, that is worth a second look: it probably means one of the six injected luxury outliers happened to land in that neighborhood by chance, since we chose which rows to inflate at random, independent of neighborhood.

Step 3: measure the spread

A hand-worked warm-up

Take five home sizes in square feet: 1200, 1350, 1400, 1550 and 2100. The range is simply the largest minus the smallest: 2100 minus 1200 is 900 square feet.

For variance, first find the mean: (1200 + 1350 + 1400 + 1550 + 2100) divided by 5 is 1520. Each value's deviation from that mean is -320, -170, -120, 30 and 580 (these always sum to zero, which is a good way to check your arithmetic). Square each deviation - 102400, 28900, 14400, 900 and 336400 - and add them up to get 483000. Divide by 5 (treating these five homes as the entire population we care about) and you get a population variance of 96600, whose square root, the population standard deviation, is about 310.8 square feet. If instead you treat these five as a sample meant to represent a larger population, Part 4's correction divides by 4 instead of 5, giving a sample variance of 120750 and a sample standard deviation of about 347.5. The sample version is always a little larger, because dividing by a smaller number produces a bigger result, and that correction exists precisely because a small sample tends to understate how spread out the full population really is.

For the interquartile range, sort the five values (they already are) and find the 25th and 75th percentiles using the same linear interpolation method pandas uses by default. With five values, the 25th percentile lands exactly on the second value, 1350, and the 75th percentile lands exactly on the fourth value, 1550. The IQR is 1550 minus 1350, which is 200 square feet - a much narrower view of "typical spread" than the range of 900, because the IQR ignores the smallest and largest 25 percent of the data entirely.

Now the full dataset

print("price range:", housing["price"].max() - housing["price"].min())
print("price variance (sample):", housing["price"].var())
print("price std (sample):", housing["price"].std())
print("size range:", housing["size_sqft"].max() - housing["size_sqft"].min())
print("size std (sample):", housing["size_sqft"].std())

pandas' .var() and .std() use the sample formula by default, dividing by n minus 1 rather than n, exactly like the second half of the hand-worked example above. With 500 rows that correction barely matters - dividing by 499 instead of 500 changes almost nothing - but it is worth knowing it is happening, because on a dataset of 5 or 10 rows the difference would be large enough to change your conclusions, as the worked example just showed.

Expect the price standard deviation to look large in absolute terms, and expect it to look large even relative to a typical price, because those six inflated outliers contribute disproportionately to variance - squaring a deviation punishes large deviations far more than small ones, so a handful of very large deviations can dominate the sum even though they are a tiny fraction of the rows. The size_sqft standard deviation, by contrast, should look much more modest relative to its own mean, because size was built from ordinary additive normal noise with no injected outliers. Comparing these two standard deviations side by side is a preview of the shape work in the next two sections: a big standard deviation is not automatically bad news, but it is always a prompt to ask why the data is spread out that much, not just how much.

Step 4: percentiles, quartiles and a box plot by neighborhood

quartiles = housing["price"].quantile([0.25, 0.5, 0.75])
print(quartiles)

Q1 = housing["price"].quantile(0.25)
Q3 = housing["price"].quantile(0.75)
print("IQR:", Q3 - Q1)

quantile(0.25) finds the value below which a quarter of the sorted prices fall, quantile(0.75) finds the value below which three quarters fall, and the difference between them is the interquartile range - the width of the middle 50 percent of the market, ignoring the cheapest quarter and the most expensive quarter. Because the expensive outliers sit well above the 75th percentile, the IQR is a far steadier description of "how spread out typical prices are" than either the full range or the standard deviation, both of which the outliers can stretch a long way.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 5))
housing.boxplot(column="price", by="neighborhood", ax=ax)
ax.set_title("Price distribution by neighborhood")
ax.set_ylabel("Price")
plt.suptitle("")
plt.show()

Recall from Part 5 how to read this: each box spans the IQR for that neighborhood, with a line inside marking the median, whiskers stretching out to the most extreme values that are not flagged as outliers, and individual dots beyond the whiskers for anything the box plot's own rule - usually 1.5 times the IQR beyond the box edges - considers unusually far from the rest of that group. You should see Downtown's box sitting highest on the price axis and Maple Park's sitting lowest, matching the price-per-square-foot ordering from Step 2. You should also see a small number of dots scattered above some of the boxes - these are very likely some or all of the six luxury listings we built into the data, each showing up as an outlier within whichever neighborhood it happened to land in, which is exactly the kind of flag this plot is designed to raise.

Step 5: look at the shape, not just the summary numbers

fig, ax = plt.subplots(figsize=(8, 5))
ax.hist(housing["price"], bins=30, edgecolor="white")
ax.set_xlabel("Price")
ax.set_ylabel("Number of listings")
ax.set_title("Distribution of housing prices in Millbrook")
plt.show()

Expect a tall cluster of bars somewhere in the few-hundred-thousand range, tapering off to the right, with a thin scatter of bars stretching much further right than the main cluster before the histogram runs out of data entirely. That long thin right-hand tail with almost nothing in it is the visual signature of right skew, the same shape Part 6 introduced using income and company size as everyday examples. Housing prices are a classic case: most homes cluster in a fairly ordinary range, and a small number of much more expensive properties stretch the tail out without there being any corresponding cluster of extremely cheap homes to balance it on the left.

print("price skew:", housing["price"].skew())
print("size skew:", housing["size_sqft"].skew())

pandas' built-in .skew() computes a standard numerical measure of skewness: zero for a perfectly symmetric distribution, positive for a right-leaning tail like the one you just saw in the histogram, negative for a left-leaning tail. You should see a clearly positive number for price - this dataset was built to have a long right tail, both from the lognormal noise and from the six injected outliers - and a number close to zero for size_sqft, which was built from ordinary additive normal noise and has no injected outliers of its own. This is a useful pairing to keep in mind: the same dataset can hold one column that behaves close to symmetric and another, right next to it, that does not, and nothing about the column names or data types alone would have told you that in advance. You had to look.

Step 6: check the normal distribution's 68-95-99.7 rule against this data

Part 7 introduced the rule that, for data shaped like a normal distribution, about 68 percent of values fall within one standard deviation of the mean, about 95 percent within two, and about 99.7 percent within three. The rule is a consequence of the normal distribution's specific bell shape - it is not true of distributions in general, and the whole point of this section is to show a column where it clearly does not hold, right next to one where it comes close.

price_mean = housing["price"].mean()
price_std = housing["price"].std()

for k in [1, 2, 3]:
    low = price_mean - k * price_std
    high = price_mean + k * price_std
    share = housing["price"].between(low, high).mean()
    print(f"price within {k} std: {share:.1%}")

Here is the reasoning to apply once you have your own numbers in front of you, rather than a specific percentage to expect, since the exact figures depend on exactly where the random outliers landed. Those six inflated listings do two things to the standard deviation: they make it considerably larger than it would be without them, because squaring a huge deviation contributes a huge amount to the variance sum. A standard deviation that is inflated this way tends to swallow up the ordinary, non-outlier spread of the data almost entirely - so do not be surprised if the "within 1 standard deviation" share you compute comes out well above 68 percent, sometimes a good deal above it. That is not a coincidence and not an error in the code; it is what happens when a few extreme values stretch the ruler (the standard deviation) so much that almost everything else fits comfortably inside one stretched unit of it. The clean 68-95-99.7 progression is a property of symmetric, bell-shaped data, and a handful of one-sided outliers is enough to break that property even though it barely changes the median at all.

size_mean = housing["size_sqft"].mean()
size_std = housing["size_sqft"].std()

for k in [1, 2, 3]:
    low = size_mean - k * size_std
    high = size_mean + k * size_std
    share = housing["size_sqft"].between(low, high).mean()
    print(f"size within {k} std: {share:.1%}")

Run this second version on size_sqft, which has no injected outliers and was built from additive normal noise around a base tied to bedroom count. You should see shares considerably closer to the textbook 68, 95 and 99.7 percent than the price column produced, though not a perfect match - size_sqft is a mixture of five slightly different normal-ish distributions (one per bedroom count) rather than one single clean normal distribution, and it was also clipped at a 400 square foot floor, which trims the distribution's left tail a little. Close-but-not-exact is a perfectly realistic outcome; it is the clear mismatch on the price column that is the real lesson of this section.

It is also worth computing a z-score for one specific listing by hand, the way Part 7 did, because it gives you a concrete sense of what the standard deviation is actually measuring. Suppose the most expensive listing in your run sold for 1,350,000 and your computed mean and standard deviation for price came out to roughly 480,000 and 190,000 (these three numbers are illustrative placeholders - substitute your own output). Its z-score would be (1,350,000 minus 480,000) divided by 190,000, which comes to about 4.58. A value more than four and a half standard deviations from the mean is extreme by any standard, in any roughly bell-shaped dataset - which is exactly why the next section uses a z-score cutoff as one of two ways to flag outliers automatically.

Step 7: find the outliers, then decide what they mean

Part 8 made a point of separating two tasks that are easy to collapse into one: finding unusual values, and deciding what to do about them. Finding them is mechanical. Deciding is judgment, and it depends on facts a spreadsheet cannot supply on its own - is this value a mistake, or a real and important part of the story?

Q1 = housing["price"].quantile(0.25)
Q3 = housing["price"].quantile(0.75)
IQR = Q3 - Q1
lower_fence = Q1 - 1.5 * IQR
upper_fence = Q3 + 1.5 * IQR

iqr_outliers = housing[(housing["price"] < lower_fence) | (housing["price"] > upper_fence)]
print(len(iqr_outliers))
print(iqr_outliers[["neighborhood", "price"]].sort_values("price", ascending=False))

This is the same 1.5-times-IQR rule the box plot drew for you visually in Step 4, written out as code so you can inspect the flagged rows directly. You should see at least the six listings we deliberately inflated, and quite possibly a couple more besides - in a right-skewed dataset like this one, some genuinely large but not artificially manipulated sales can sit just past the fence too, simply because the upper tail is long and gradually thinning rather than cleanly separated from the rest. That is expected, and it is the reason Part 8 insisted on treating a flag as the start of an investigation, not an automatic deletion order: some of these rows are the injected outliers, and some might just be honestly large, expensive homes that belong in the dataset exactly as they are.

housing["price_zscore"] = (housing["price"] - price_mean) / price_std
z_outliers = housing[housing["price_zscore"].abs() > 3]
print(len(z_outliers))

Compare the count from this z-score rule with the count from the IQR rule above. It would not be surprising if the z-score rule flags fewer rows. The reason is almost circular, and worth sitting with: the z-score rule measures distance using the mean and standard deviation, but the same six outliers that we want to catch are also busy inflating that mean and standard deviation in the first place. A standard deviation that has been stretched by the outliers it is trying to detect becomes a less sensitive ruler for detecting them - an outlier has to be even more extreme to register three stretched standard deviations away. The IQR rule avoids this problem because the median and quartiles barely move when a handful of extreme values are added, which is the same robustness property that made the median a better "typical price" than the mean back in Step 2.

In a real project, this is the point where you would stop writing code and start making phone calls, or at least opening the original listing. Is the three-million-dollar row a genuine waterfront estate, in which case it belongs in the dataset and possibly deserves its own line in your eventual report? Or is it a data entry error - a decimal point in the wrong place, a price field that accidentally includes the lot next door - in which case it needs correcting or removing before you say anything else about the market. The statistics in this section can only tell you where to look. They cannot tell you what you will find when you get there.

Step 8: correlation, and a correlation that would mislead you

numeric_cols = ["bedrooms", "bathrooms", "size_sqft", "year_built", "price"]
corr_matrix = housing[numeric_cols].corr()
print(corr_matrix.round(2))

Part 9 defined correlation as a number between -1 and 1 that measures how closely two numeric columns move together in a straight-line sense, with 1 meaning perfectly together, -1 meaning perfectly opposite, and 0 meaning no linear relationship at all. Expect size_sqft and price to show a strong positive correlation, likely somewhere well above 0.7, because price was built directly from size multiplied by a rate - that relationship is baked into the data by construction, so a weak correlation here would actually be a sign something went wrong in the code above. Expect bedrooms and price to show a positive but probably more moderate correlation, and expect year_built and price to come out close to zero, because year_built was never used anywhere in the price formula - any correlation you see there is sampling noise, not a real relationship, and it should shrink further if you regenerate the dataset with a larger n.

Bedrooms predicts price mostly because bedrooms predicts size, and size predicts price - not because adding a bedroom, by itself, adds a fixed amount of value.

That is worth unpacking, because it is the single most useful habit Part 9 tried to build: before treating a correlation as a direct cause-and-effect relationship, ask whether a third variable could be doing the real work. In this dataset we can actually check, because we wrote the generator ourselves. Look back at the code: bedrooms feeds directly into base_size, and base_size feeds directly into price. Bedrooms never appears in the price formula on its own. So the correlation you see between bedrooms and price is entirely mediated through size - bedrooms is a useful stand-in for "probably a bigger home", and it is the size, combined with the neighborhood rate, that actually sets the price. In a real dataset you cannot inspect the generator, but you can still ask the same question: is there a third column, like size or location, that plausibly drives both of the two columns I am correlating?

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(housing["size_sqft"], housing["price"], alpha=0.5)
ax.set_xlabel("Size (sqft)")
ax.set_ylabel("Price")
ax.set_title("Size vs price")
plt.show()

You should see an upward-sloping cloud of points - bigger homes generally cost more, exactly as the formula guarantees - that also fans out wider as size increases, with the spread of prices among the larger homes visibly bigger than the spread among the smaller ones. That fanning-out is a direct consequence of using multiplicative noise: a 10 percent random swing on a 300,000 dollar base home is a smaller dollar amount than the same 10 percent swing on a 900,000 dollar base home. You should also see a few points sitting well above the main cloud, at a given size, well above where the rest of the points at that size level sit - these are your six luxury outliers again, visible here as points whose price is much higher than their size alone would predict, which is a genuinely different signal from "this home is simply bigger than average".

One more caution belongs here, and it connects straight back to Step 1. It would be technically possible to label-encode neighborhood as the numbers 0 through 4 - say, alphabetically - and feed that numeric column into .corr() alongside price. pandas would compute a number without complaining. That number would be close to meaningless, because neighborhood has no real order; an alphabetical encoding imposes an artificial one, and the correlation coefficient would be measuring the strength of a straight-line relationship with a sequence that exists only in how you happened to sort the names. If you want to relate a nominal column like neighborhood to price properly, the groupby comparison from Step 2 - comparing means or medians across groups - is the honest tool, not a correlation coefficient computed on an invented numbering.

Step 9: draw a sample, then draw a biased one

Everything so far has treated all 500 rows of housing as if they were the entire market - the full population of interest. Part 10 made the point that in most real projects you are not that lucky: you have a sample, drawn somehow from a much larger population, and the question of how it was drawn matters as much as what it contains. We can simulate both halves of that lesson here, because we happen to have the full 500-row "population" sitting in memory.

population_mean = housing["price"].mean()
population_median = housing["price"].median()

sample = housing.sample(n=40, random_state=1)
sample_mean = sample["price"].mean()
sample_median = sample["price"].median()

print(f"population mean: {population_mean:.0f}, sample mean: {sample_mean:.0f}")
print(f"population median: {population_median:.0f}, sample median: {sample_median:.0f}")

housing.sample(n=40, random_state=1) picks 40 rows completely at random, with every row equally likely to be chosen - an honest, unbiased sample of the population we already have in full. You should see the sample's mean and median land somewhere in the neighborhood of the population's mean and median, but not exactly on them. That gap is ordinary sampling variability, not a mistake: any random sample of 40 out of 500 will, by chance, include a slightly different mix of cheap and expensive homes than the full set does, and a smaller sample has more room for that kind of chance imbalance than a larger one. Try changing random_state to a few different numbers and rerun - you should see the sample statistics shift around each time, sometimes landing closer to the population values and sometimes further away, which is exactly the lesson: a sample estimate is a reasonable guess, not a guarantee.

biased_sample = housing[housing["neighborhood"] == "Downtown"].sample(n=20, random_state=1)
print("biased sample mean:", biased_sample["price"].mean())
print("population mean:", population_mean)

This second sample is not a random draw from the whole population - it only ever looks at Downtown listings before sampling. Because Downtown carries the highest price-per-square-foot rate of the five neighborhoods, this sample's mean should come out clearly higher than the population mean, and that gap is not sampling noise that would shrink if you drew a bigger sample. It is sampling bias: every row in this sample comes from a slice of the market that is systematically more expensive than the city as a whole, so no amount of additional sampling within Downtown alone would ever bring the estimate back down to the real citywide figure. In a real project, this is exactly the failure mode to watch for when your data comes from a single listings site that happens to cover some neighborhoods better than others, a single real-estate agency that specializes in a particular price range, or a time window that happened to catch a local construction boom in one part of town. The arithmetic on a biased sample is not wrong. The sample itself is answering a narrower question than the one you think you are asking.

Step 10: write the statistical profile

Everything above produced numbers and plots. A statistical profile turns them into a short piece of writing someone without your code could actually use. The discipline worth building here is to let each sentence be earned by a specific number you computed, rather than writing what sounds reasonable and backfilling a statistic to match it afterward.

A solid profile of a dataset like this one tends to cover six things, roughly in this order:

• Scope: how many listings, over what neighborhoods, and - importantly - what is not covered (a point we will come back to below)
• Typical price: the median, with a short note on why the median rather than the mean, if the data is skewed
• Spread: the IQR or a plain-language range for "most homes", separate from the full min-to-max range which outliers can stretch a long way
• Shape and outliers: whether the distribution is skewed, and an honest account of any flagged outliers - investigated, not silently dropped or silently kept
• What relates to price, and what probably causes what: the strongest numeric relationships, with a caveat about mediating variables like size standing in for bedrooms
• Limits of the data: who or what might be under-represented, and what that means for how far the conclusions generalize

Across the 500 recorded sales, the median price was [your median here]; we lead with the median rather than the mean because six sales - roughly one in eighty - were well above the rest of the market and would otherwise pull a single "average price" figure upward in a way that does not reflect what a typical buyer actually pays.

That sentence is written with brackets on purpose: plug in the median you actually computed in Step 2, not a number copied from this article, because your own run of the code - with its own random seed, or with real data in place of Millbrook's - will produce its own figures. The sentence's structure, though, is reusable: state the robust number, then state in one clause why you are leading with it instead of the more commonly reported mean. That one clause is usually the difference between a report that merely states numbers and a report that teaches the reader something about the market.

The scope and limits sections are the ones people skip, and they are often the most important for whoever reads the report next. For Millbrook, an honest scope note might read: these are 500 recorded sales across five neighborhoods; the data does not distinguish between resales and new construction, does not include homes that were listed but never sold, and does not cover rentals at all. For a real dataset, the equivalent note might mention which listings site the data came from, which date range it covers, and which kinds of properties - say, homes sold through private brokers rather than public listings - might be systematically missing. Step 9's biased-sample exercise is the whole reason this note belongs in the report: a number computed correctly from a narrow or skewed slice of the market can still mislead anyone who assumes it describes the whole market.

Common mistakes in a project like this one

• Reporting the mean price without checking the skew first, so a handful of expensive outliers quietly define what counts as "typical"
• Computing a correlation or an average on a nominal categorical column that was label-encoded into numbers for convenience, producing a precise-looking number that has no real meaning
• Treating every value flagged by the IQR rule or a z-score cutoff as an error to delete, instead of a prompt to look closer at what the row actually represents
• Using the z-score rule alone on skewed data without noticing that the same outliers you are hunting for have already inflated the mean and standard deviation you are measuring them against
• Reading a correlation between two columns as proof that one causes the other, without asking whether a third variable - like size standing behind both bedrooms and price - is doing the real work
• Generalizing from a sample that was convenient to collect (one neighborhood, one listings site, one time window) to a claim about the whole city, without naming that limitation anywhere in the write-up
• Quoting a sample statistic as if it were an exact population value, instead of acknowledging the sampling variability that Step 9 demonstrated directly
• Trusting summary numbers without ever plotting the data - different distributions can share very similar means, medians and even standard deviations while looking completely different once you actually look at them, which is exactly why the histogram and box plot earned their own steps in this project rather than being optional extras

Summary and what's next

This project did not introduce one new statistical idea. It took the data type classification from Part 2, the mean, median and mode from Part 3, the range, variance, standard deviation and IQR from Part 4, the percentiles and box plots from Part 5, the skew and shape vocabulary from Part 6, the normal distribution's 68-95-99.7 rule from Part 7, the outlier detection methods from Part 8, the correlation and causation caution from Part 9, and the population, sample and sampling bias distinction from Part 10, and ran all of it, in order, on one dataset with a story behind it. That order was not arbitrary: data types first, because they decide which statistics are even valid; centre and spread next, because they anchor everything else; shape and outliers after that, because they tell you whether to trust the centre and spread you just computed; relationships between columns after that, because they are easy to over-interpret once you are not careful about what causes what; and sampling last, because it is the step that asks whether any of the preceding numbers generalize beyond the rows you happened to have.

That sequence - types, centre, spread, shape, outliers, relationships, sampling, write-up - is reusable on almost any dataset where one numeric column matters to someone: delivery times for an operations team, test scores for a school, response times for a support queue, salaries for an HR review. The code changes very little between those cases; the judgment calls, about what a mean obscures and what a correlation does not prove, stay exactly the same.

This is also the end of the Statistics From Zero for Data Work path. The honest next step, now that you can describe a dataset properly, is inferential statistics: estimating a population value from a sample with a stated margin of uncertainty, and testing whether a difference you see - between two neighborhoods, two time periods, two groups of customers - is likely to be real or could plausibly be explained by chance alone. Those ideas build directly on the sampling and distribution concepts from Parts 7 and 10, and on the descriptive habits practised across this entire project. That territory is not covered by the twelve parts you have just finished, but everything in them is the foundation it stands on, and you now have a complete, working example of that foundation in action - with a dataset and a set of results you built and checked yourself, rather than took on faith.