Week 5: The first regressors

Goal of today: a script that trains two regressors (linear regression and random forest) and shows that they forecast tomorrow's max temperature better than any single weather model:

Eger / models: 974 complete days
  ecmwf                 0.82 °C
  gfs                   1.56 °C
  icon                  1.28 °C
  average of sources    0.93 °C
  linear                0.68 °C  <-- best
  random_forest         0.82 °C

How to work: do the steps in order. After each step, run the code and compare with the expected output. Go on only when yours looks the same. (Your numbers will differ a little, because you have more days.)


What you need to know

Machine learning means: instead of writing the rule yourself, you show the computer many examples and it finds the rule. Our rule: "three forecasts → the real temperature".

Regression or classification? If the answer is a number (a temperature, a price), it is regression. If the answer is a category (rain / no rain, spam / not spam), it is classification. We do regression. A model that does regression is called a regressor.

Features and target. Each complete day is one example. Here is one day (one row) of the history set; the whole table has about 970 such rows:

features (inputs, X) target (answer, y)
ecmwf gfs icon measured
11.1 11.4 12.4 12.2

The model learns from many rows like this. Later we give it only the features (tomorrow's forecasts) and it gives back its own forecast for the target.

fit and predict. In scikit-learn every model has the same two methods: - model.fit(X, y): learn from examples (features + answers). - model.predict(X): forecast answers for new features.

So trying another model is usually one extra line.

Training and testing. A student who learned the answers of last year's exam by heart is not necessarily good at the subject. A model is the same: on the days it learned from, it can look perfect just by remembering them. So we keep some days aside (the test days). The model never sees them while learning, and we measure the error only on them.

Overfitting is exactly that: the model learned the training days by heart, including their random noise, and it does badly on new days. A model with a small training error and a big test error is overfitting.

Split by date, never randomly. We train on the older 80% of the days and test on the newest 20%. This is how the model will be used in real life: learn from the past, forecast the future. A random split would put Monday into training and Tuesday into testing. Weather on neighbouring days is similar, so the test would be too easy and the error too optimistic.

MAE (mean absolute error) from week 4: the average of |forecast − measured|, in °C. Lower is better. Example: forecasts 20, 22, 25 and measured 21, 22, 23 → errors 1, 0, 2 → MAE = 1.0 °C.

Baseline. Before believing a model is good, compare it with simple solutions that need no learning: each source alone, and the average of the three sources. A model that cannot beat these is useless.

Linear regression learns one weight for each input and a constant:

forecast = w1·ecmwf + w2·gfs + w3·icon + b

It chooses w1, w2, w3, b so that the error on the training days is as small as possible. This way it learns two useful things automatically: which source to trust more (bigger weight) and systematic errors (if the models are always a bit too warm, the constant b corrects it).

Decision tree and random forest. A decision tree is a list of yes/no questions:

is ecmwf > 21.3?
 ├─ yes: is icon > 25.0?  ─ yes: 26.1   no: 23.4
 └─ no:  is gfs  > 12.0?  ─ yes: 15.2   no:  8.7

Each "leaf" at the end gives the average of the training days that ended up there. A random forest builds many such trees (we use 200), each on a random part of the data, and averages their answers. It can learn curved, non-linear rules. But its answers are always averages of training days, so it can never predict a value higher than the hottest training day. We will see this in step 7.


Step 1: Install scikit-learn

scikit-learn (imported as sklearn) is the most used Python library for classic machine learning. Add it to requirements.txt:

requirements.txt

requests
tzdata
pandas
flask
scikit-learn

Install it (with the venv active):

pip install -r requirements.txt

Try it:

python -c "import sklearn; print(sklearn.__version__)"

You should see a version number, e.g. 1.7.2.


Step 2: ml.py: the data set

Create a new file ml.py with this content:

"""Machine learning: learn how to combine several forecasts into a better one.

Run it directly to compare all models:   python ml.py
"""
from config import CITIES, SETS
from db import load_table
def build_dataset(city, sources):
    """Complete days only: every forecast AND the measured value must be present."""
    return load_table(city, sources).dropna()

load_table() (week 3) gives one row per day. Some rows are not complete: tomorrow has no measured value yet, and a source may be missing on some day. dropna() drops every row that has an empty value (NaN). A model cannot learn from a row without an answer.

Try it:

python
>>> from ml import build_dataset
>>> table = build_dataset("eger", ["ecmwf", "gfs", "icon"])
>>> table.head()
            ecmwf   gfs  icon  measured
day
2024-02-05   11.1  11.4  12.4      12.2
2024-02-06   12.4  13.1  13.0      14.7
2024-02-07   11.0  12.2  13.4      12.6
2024-02-08   12.0  11.3  12.4      11.9
2024-02-09   11.4  10.8  13.3      11.6
>>> len(table)
974
>>> exit()

head() shows the first 5 rows. 974 complete days: that is our data set. The first three columns are the features, measured is the target.


Step 3: Train and test parts

Add this function to the end of ml.py:

def split_by_date(table, test_share=0.2):
    """Older days for training, the newest days for testing. Never shuffle time series!"""
    cut = int(len(table) * (1 - test_share))
    return table.iloc[:cut], table.iloc[cut:]

Try it:

python
>>> from ml import build_dataset, split_by_date
>>> table = build_dataset("eger", ["ecmwf", "gfs", "icon"])
>>> train, test = split_by_date(table)
>>> len(train), train.index.min(), train.index.max()
(779, '2024-02-05', '2026-03-24')
>>> len(test), test.index.min(), test.index.max()
(195, '2026-03-25', '2026-10-05')
>>> exit()

779 days for learning, 195 days for testing, and every test day is after every training day.


Step 4: Measuring the error

Add import numpy as np at the top of ml.py, so the beginning looks like this:

"""Machine learning: learn how to combine several forecasts into a better one.

Run it directly to compare all models:   python ml.py
"""
import numpy as np

from config import CITIES, SETS
from db import load_table

Then add this function to the end:

def mae(predicted, measured):
    """Mean absolute error: on average, how many degrees we are wrong."""
    return float(np.mean(np.abs(np.asarray(predicted) - np.asarray(measured))))

numpy is the library for fast calculations on lists of numbers (pandas is built on it). np.asarray() turns a list or a pandas column into a numpy array, then -, np.abs() and np.mean() work on all elements at once. No loop needed.

Try it with the example from What you need to know:

python
>>> from ml import mae
>>> mae([20, 22, 25], [21, 22, 23])
1.0
>>> exit()

Step 5: Baselines: how good are the forecasts alone?

Add this function to the end of ml.py:

def baseline_scores(test, sources):
    """How good are the forecasts without any machine learning?"""
    scores = {source: mae(test[source], test["measured"]) for source in sources}
    scores["average of sources"] = mae(test[sources].mean(axis=1), test["measured"])
    return scores

Try it:

python
>>> from ml import build_dataset, split_by_date, baseline_scores
>>> sources = ["ecmwf", "gfs", "icon"]
>>> train, test = split_by_date(build_dataset("eger", sources))
>>> baseline_scores(test, sources)
{'ecmwf': 0.8179487179487177, 'gfs': 1.556923076923077, 'icon': 1.281025641025641, 'average of sources': 0.9324786324786325}
>>> exit()

ECMWF is the best source in Eger: 0.82 °C average error. Note that the simple average (0.93) is worse than ECMWF alone: averaging blindly mixes in the bad sources too. Our models must beat 0.82.


Step 6: The first model: linear regression

Before writing it into ml.py, try a model by hand, to see each step:

python
>>> from sklearn.linear_model import LinearRegression
>>> from ml import build_dataset, split_by_date, mae
>>> sources = ["ecmwf", "gfs", "icon"]
>>> train, test = split_by_date(build_dataset("eger", sources))
>>> model = LinearRegression()
>>> model.fit(train[sources], train["measured"])
LinearRegression()
>>> model.coef_
array([0.46460461, 0.0556095 , 0.47428087])
>>> model.intercept_
np.float64(-0.44572442289740266)

What happened: - LinearRegression() creates an empty, untrained model. - fit(X, y): X is the three forecast columns of the training days, y is the measured column. The model learns. - coef_ are the learned weights, intercept_ is the constant. So the model learned this rule for Eger:

forecast = 0.46·ecmwf + 0.06·gfs + 0.47·icon − 0.45

GFS gets almost no weight: the model found that GFS adds little that the other two do not already know.

Now let it forecast the test days, which it has never seen:

>>> predicted = model.predict(test[sources])
>>> predicted[:3]
array([15.33821823, 13.66369479, 11.79994104])
>>> test["measured"].head(3)
day
2026-03-25    15.4
2026-03-26    15.3
2026-03-27    11.2
Name: measured, dtype: float64
>>> mae(predicted, test["measured"])
0.6831933850344467
>>> exit()

0.68 °C, better than the best single source (0.82). The model learned, from the past, how to combine the forecasts.


Step 7: The second model: random forest

Try it the same way:

python
>>> from sklearn.ensemble import RandomForestRegressor
>>> from ml import build_dataset, split_by_date, mae
>>> sources = ["ecmwf", "gfs", "icon"]
>>> train, test = split_by_date(build_dataset("eger", sources))
>>> forest = RandomForestRegressor(n_estimators=200, min_samples_leaf=5, random_state=0)
>>> forest.fit(train[sources], train["measured"])
RandomForestRegressor(min_samples_leaf=5, n_estimators=200, random_state=0)
>>> mae(forest.predict(test[sources]), test["measured"])
0.8226063969605073

0.82 °C: no better than ECMWF alone. Why? Look at the hottest days:

>>> train["measured"].max(), test["measured"].max()
(np.float64(34.9), np.float64(38.6))
>>> forest.predict(test[sources]).max()
np.float64(33.8)
>>> exit()

The hottest training day was 34.9 °C, but the test period (summer 2026) had days up to 38.6 °C. The forest never predicted more than 33.8 °C: its answers are averages of training days, so it cannot go beyond what it has seen. On the 18 test days above 33 °C, its error is 2.06 °C, while linear regression's is 0.76 °C. Linear regression just continues the straight line, and here the relationship really is almost a straight line.


Step 8: All models in one function

Now put what you did by hand into ml.py. Add the two imports at the top, and the constant MIN_DAYS after the imports. The beginning of ml.py should look like this:

"""Machine learning: learn how to combine several forecasts into a better one.

Run it directly to compare all models:   python ml.py
"""
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression

from config import CITIES, SETS
from db import load_table

MIN_DAYS = 30  # do not train on fewer complete days than this

Then add these two functions to the end:

def make_models():
    linear = LinearRegression()
    forest = RandomForestRegressor(n_estimators=200, min_samples_leaf=5, random_state=0)
    return {"linear": linear, "random_forest": forest}


def evaluate(city, data_set):
    """Train every model on the older days, measure the error on the newest days."""
    sources = SETS[data_set]
    table = build_dataset(city, sources)
    if len(table) < MIN_DAYS:
        return None, len(table)
    train, test = split_by_date(table)
    scores = baseline_scores(test, sources)
    for name, model in make_models().items():
        model.fit(train[sources], train["measured"])
        scores[name] = mae(model.predict(test[sources]), test["measured"])
    return scores, len(table)

Try it:

python
>>> from ml import evaluate
>>> evaluate("eger", "models")
({'ecmwf': 0.8179487179487177, 'gfs': 1.556923076923077, 'icon': 1.281025641025641, 'average of sources': 0.9324786324786325, 'linear': 0.6831933850344467, 'random_forest': 0.8226063969605073}, 974)
>>> evaluate("eger", "providers")
(None, 0)
>>> exit()

The same numbers as in steps 5–7. The providers set (your live data) has no complete days yet. That is normal. In the next weeks it grows; by week 8 it will have enough.


Step 9: A readable report for all cities

Add these to the end of ml.py:

def print_scores(scores):
    best = min(scores, key=scores.get)
    for name, value in scores.items():
        mark = "  <-- best" if name == best else ""
        print(f"  {name:20s} {value:5.2f} °C{mark}")
if __name__ == "__main__":
    for city in CITIES:
        for data_set in SETS:
            scores, days = evaluate(city, data_set)
            print(f"\n{CITIES[city]['name']} / {data_set}: {days} complete days")
            if scores is None:
                print(f"  not enough data yet (need {MIN_DAYS})")
            else:
                print_scores(scores)

Run it:

python ml.py

You should see:

Budapest / models: 974 complete days
  ecmwf                 0.64 °C
  gfs                   1.18 °C
  icon                  0.70 °C
  average of sources    0.67 °C
  linear                0.58 °C  <-- best
  random_forest         0.73 °C

Budapest / providers: 0 complete days
  not enough data yet (need 30)

Eger / models: 974 complete days
  ecmwf                 0.82 °C
  gfs                   1.56 °C
  icon                  1.28 °C
  average of sources    0.93 °C
  linear                0.68 °C  <-- best
  random_forest         0.82 °C

Eger / providers: 0 complete days
  not enough data yet (need 30)

In both cities, linear regression beats the best weather model: Budapest 0.64 → 0.58 °C, Eger 0.82 → 0.68 °C.


Your complete ml.py

Compare with yours. If something does not work, copy this one.

ml.py

"""Machine learning: learn how to combine several forecasts into a better one.

Run it directly to compare all models:   python ml.py
"""
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression

from config import CITIES, SETS
from db import load_table

MIN_DAYS = 30  # do not train on fewer complete days than this


def build_dataset(city, sources):
    """Complete days only: every forecast AND the measured value must be present."""
    return load_table(city, sources).dropna()


def split_by_date(table, test_share=0.2):
    """Older days for training, the newest days for testing. Never shuffle time series!"""
    cut = int(len(table) * (1 - test_share))
    return table.iloc[:cut], table.iloc[cut:]


def mae(predicted, measured):
    """Mean absolute error: on average, how many degrees we are wrong."""
    return float(np.mean(np.abs(np.asarray(predicted) - np.asarray(measured))))


def baseline_scores(test, sources):
    """How good are the forecasts without any machine learning?"""
    scores = {source: mae(test[source], test["measured"]) for source in sources}
    scores["average of sources"] = mae(test[sources].mean(axis=1), test["measured"])
    return scores


def make_models():
    linear = LinearRegression()
    forest = RandomForestRegressor(n_estimators=200, min_samples_leaf=5, random_state=0)
    return {"linear": linear, "random_forest": forest}


def evaluate(city, data_set):
    """Train every model on the older days, measure the error on the newest days."""
    sources = SETS[data_set]
    table = build_dataset(city, sources)
    if len(table) < MIN_DAYS:
        return None, len(table)
    train, test = split_by_date(table)
    scores = baseline_scores(test, sources)
    for name, model in make_models().items():
        model.fit(train[sources], train["measured"])
        scores[name] = mae(model.predict(test[sources]), test["measured"])
    return scores, len(table)


def print_scores(scores):
    best = min(scores, key=scores.get)
    for name, value in scores.items():
        mark = "  <-- best" if name == best else ""
        print(f"  {name:20s} {value:5.2f} °C{mark}")


if __name__ == "__main__":
    for city in CITIES:
        for data_set in SETS:
            scores, days = evaluate(city, data_set)
            print(f"\n{CITIES[city]['name']} / {data_set}: {days} complete days")
            if scores is None:
                print(f"  not enough data yet (need {MIN_DAYS})")
            else:
                print_scores(scores)

Check

Extra

  1. In Eger, ICON gets as much weight as ECMWF (step 6), although its own error is much bigger (1.28 vs 0.82 °C). How is that possible? (Hint: compute (test["icon"] - test["measured"]).mean(). A source that is always about 1 °C too warm is wrong, but still very useful, because its error is predictable and the constant corrects it.)
  2. Add more regressors to make_models() and compare. Each one is a single line: from sklearn.linear_model import Ridge from sklearn.neighbors import KNeighborsRegressor from sklearn.ensemble import GradientBoostingRegressor ... "ridge": Ridge(alpha=1.0), "knn": KNeighborsRegressor(n_neighbors=10), "boosting": GradientBoostingRegressor(random_state=0),
  3. Try a random split (train_test_split from sklearn.model_selection with shuffle=True) and compare the errors with the date split. Which one is more honest?

If something goes wrong

Problem Reason
ModuleNotFoundError: No module named 'sklearn' Install it (step 1), with the venv active. The package is called scikit-learn, the import sklearn.
not enough data yet for models The history set is missing. Run load_history.py (week 3).
ValueError: Input contains NaN You removed dropna() from build_dataset().
ImportError: cannot import name 'evaluate' in the interactive Python You changed ml.py: exit() and start python again.
The numbers differ a little from this page Normal: you have more days than when this page was written.

← Week 4Week 6 →