Hvordan kan vi lære noe av et datasett?
Ikke alt dekkes i forelesninger. Noe er bare i boka
La oss anta at vi har lest inn prisdata fra nordpool og lagt dem i en csv-fil. Så begynner vi med prisen på strømmen.
import pandas as pd
df_pris = pd.read_csv("data/nordpool-dec-24.csv")
df_forbruk = pd.read_csv("data/meteringvalues-dec-2024.csv", sep=";", decimal=",")
nettopris = sum(df_forbruk['KWH 60 Forbruk']*df_pris["value"])/1000
med_moms = nettopris*1.25
print(f"Spotpris på strømmen inkl. moms er {med_moms:.2f} kr")Spotpris på strømmen inkl. moms er 2899.55 krFørst må vi ha værdata. Jeg antar vi har lagret det i en fil weather_blindern.json.
import json
with open("data/weather_blindern.json") as ifile:
json_content = json.load(ifile)
my_list = []
for data in json_content["data"]:
my_list += data["observations"]
df_temp = pd.DataFrame(my_list)
df_temp = df_temp[df_temp["timeOffset"] == "PT6H"]
print(df_temp.head()) elementId value unit \
1 mean(air_temperature P1D) 7.2 degC
3 mean(air_temperature P1D) 7.8 degC
5 mean(air_temperature P1D) 4.2 degC
7 mean(air_temperature P1D) -3.5 degC
9 mean(air_temperature P1D) -3.5 degC
level timeOffset \
1 {'levelType': 'height_above_ground', 'unit': '... PT6H
3 {'levelType': 'height_above_ground', 'unit': '... PT6H
5 {'levelType': 'height_above_ground', 'unit': '... PT6H
7 {'levelType': 'height_above_ground', 'unit': '... PT6H
9 {'levelType': 'height_above_ground', 'unit': '... PT6H
timeResolution timeSeriesId performanceCategory exposureCategory \
1 P1D 0 C 2
3 P1D 0 C 2
5 P1D 0 C 2
7 P1D 0 C 2
9 P1D 0 C 2
qualityCode
1 NaN
3 NaN
5 NaN
7 NaN
9 NaN from datetime import datetime
def date_func(row):
date_object = datetime.strptime(row["Fra"], "%d.%m.%Y %H:%M")
return date_object.date()
df_forbruk["date"] = df_forbruk.apply(date_func, axis=1)
forbruk_per_dag = df_forbruk.groupby("date").sum()
x = df_temp["value"]
y = forbruk_per_dag["KWH 60 Forbruk"]/var/folders/qn/3_cqp_vx25v4w6yrx68654q80000gp/T/ipykernel_59245/910005881.py:7: FutureWarning:
The default value of numeric_only in DataFrameGroupBy.sum is deprecated. In a future version, numeric_only will default to False. Either specify numeric_only or select only columns which should be valid for the function.
import matplotlib.pyplot as plt
plt.scatter(x, y)
plt.xlabel("Utetemperatur [$^\circ$C]")
plt.ylabel("Strømforbruk [kWh]")Text(0, 0.5, 'Strømforbruk [kWh]')
import numpy as np
dag_varmepumpe_alternativer = np.arange(12, 21)
for i, dag_varmepumpe in enumerate(dag_varmepumpe_alternativer):
plt.subplot(3,3, i+1)
forbruk_per_dag["is_varmepumpe"] = [0] * dag_varmepumpe+[1] * (31-dag_varmepumpe)
x = df_temp["value"]
y = forbruk_per_dag["KWH 60 Forbruk"]
plt.scatter(x, y, c=forbruk_per_dag["is_varmepumpe"])
plt.title(f"{dag_varmepumpe} desember")
plt.tight_layout()
Så stirrer vi på grafene og ser etter om det er noe valg som gir oss to klare grupper med datapunkter.
Fra strøm i naturfag på vgs/ungdomsskolen vet vi at \[U = RI \implies I = \frac{1}{R}U\]
Kan plotte:
Strøm \(I\) er en lineær funksjon av spenning \(U\), med invers motstand som proporsjonalitetskonstant.
Vi går på laben og finner en Ohmsk motstand, og måler, da forventer vi å finne noe slikt:
Siden stigningen er 3.33…, så er motstanden \(0.3 \mathrm{A}\).
Dermed kan vi plotte en linje sammen med datapunktene
Men er det slik måleserier ser ut i virkeligheten?
Vi får nok heller noe som ser slik ut:
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
R = 0.3
U = np.random.random(10) * 5
I = 1.0 / R * U + np.random.normal(0, 0.5, 10) # Current calculation
plt.scatter(U, I)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
R = 0.3
U = np.random.random(10) * 5
I = 1.0 / R * U + np.random.normal(0, 0.5, 10) # Current calculation
p = np.polyfit(U, I, 10)
x = np.linspace(np.min(U), np.max(U), 100)
plt.scatter(U, I)
plt.plot(x, np.polyval(p, x))/var/folders/qn/3_cqp_vx25v4w6yrx68654q80000gp/T/ipykernel_59245/659525506.py:8: RankWarning:
Polyfit may be poorly conditioned
Nå må vi forholde oss til verden slik den er. Det er alltid en viss usikkerhet vi ikke klarer å forklare uansett hvor god modell vi har. \[y=f(x) + \varepsilon\]
Der \(\varepsilon\) representeter en usikkerhet som vi ikke klarer å plukke opp med vår modell av fenomenet.
Vi ønsker å finne parametre \(\beta_0\) og \(\beta_1\) i funksjonen \[f(x) = \beta_1 x + \beta_0\]
Slik at når vi sier at \(\hat{y} = f(x)\), så minimerer vi
\[\sum_{i=0}^N (y_i - \hat{y}_i)^2\]
Dette heter minste kvadraters lineær regresjon.
import numpy as np
p = np.polyfit(U, I, 1)
print(p)[3.34088808 0.17706321]plt.scatter(U, I)
I_fit = np.polyval(p, U)
plt.plot(U, I_fit)
Skriv opp for deg selv, hva som er sammenhengen mellom den rette linja og datapunktene her om vi har gjort lineær regresjon.
def f(x):
return 2*x + 3
N = 10
slopes = []
for i in range(N):
x = np.linspace(0, 10, 10)
y = f(x) + np.random.normal(0, 1.5, 10)
p = np.polyfit(x, y, 1)
slopes.append(p[0])
y_fit = np.polyval(p, x)
plt.plot(x, y_fit, linewidth=0.75, color="k")
plt.plot(x, f(x), color="r")
plt.plot(x, y_fit, linewidth=0.75, color="k", label="Samples with n=10")
plt.plot(x, f(x), color="r", label="Population mean")
plt.legend()
a) Lag deg en affin funsjon (f. eks \(f(x) = 2x + 3\)). Trekk 10 tilfeldige tall, beregn \(f(x) + \varepsilon\) der du setter inn et tilfeldig normalfordelt tall som \(\varepsilon\) med np.random.normal(...). Verifiser med et plott at ting fungerer som de skal
b) Gjør en lineærtilpasning med np.polyfit(x, y, 1). Plott lineærtilpasningen sammen med de litt spredte punktene dine.
c) Gjenta mange ganger (f. eks 10000 i en for-løkke), og lagre stigningstallet som estimeres for hver enkelt linje. Plott et histogram over stigningstallene, med plt.hist.
def f(x):
return 2*x + 3
x = np.linspace(0, 10, 10)
y = f(x) + np.random.normal(0, 1.5, 10)
p = np.polyfit(x, y, 1)
y_hat = np.polyval(p, x)
plt.plot(x, y, "o")
plt.plot(x, y_hat)
_ = plt.figure()<Figure size 672x480 with 0 Axes>N = 10000
slopes = []
for i in range(N):
x = np.linspace(0, 10, 10)
y = f(x) + np.random.normal(0, 1.5, 10)
p = np.polyfit(x, y, 1)
slopes.append(p[0])
_ = plt.hist(slopes, bins=30)
N = 10000
slopes = []
for i in range(N):
x = np.linspace(0, 10, 100)
y = f(x) + np.random.normal(0, 1.5, 100)
p = np.polyfit(x, y, 1)
slopes.append(p[0])
_ = plt.hist(slopes, 30)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.savefig("figures/est_beta1.png")
Hva skjedde med fordelingen av stigningstall?
a) Hva skjedde med fordelingen av stigningstall ved lineær regresjon da vi økte fra 10 til 100 punkter?
b) Hvordan skal vi tolke dette?
Fra boka har vi følgende formel for forventet feil i stigningtallet: \[\mathrm{SE}(\hat{\beta}_1) = \sqrt{\frac{\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2}}\]
der \(\sigma\) er standardavviket til den tilfeldige variablen \(\varepsilon\) i \(y = f(x) + \varepsilon\)
Vi trakk tilfeldige tall med \(\sigma=1.5\) og hhv. 10 og 100 jevnt fordelte tall fra 0 til 10. Altså:
x = np.linspace(0, 10, 10)
se = np.sqrt(1.5**2/np.sum((x-np.mean(x))**2))
print(f"Standard feil for n=10: {se}")
x = np.linspace(0, 10, 100)
se = np.sqrt(1.5**2/np.sum((x-np.mean(x))**2))
print(f"Standard feil for n=100: {se}")Standard feil for n=10: 0.1486301082920587
Standard feil for n=100: 0.051444481273168134def est_dist(x, mu, std):
return 1/np.sqrt(2*np.pi*std**2)*np.exp(-(x-mu)**2/(2*std**2))
for n, color in zip([10, 100], ["red", "blue"]):
x = np.linspace(0, 10, n)
se = np.sqrt(1.5**2/np.sum((x-np.mean(x))**2))
est_vals = np.linspace(1.5, 2.5, 100)
plt.plot(est_vals, est_dist(2, est_vals, se), "--", color=color, label=f"n={n}")
plt.legend()
def est_dist(x, mu, std):
return 1/np.sqrt(2*np.pi*std**2)*np.exp(-(x-mu)**2/(2*std**2))
N = 10000
for n, color in zip([10, 100], ["red", "blue"]):
slopes = []
for i in range(N):
x = np.linspace(0, 10, n)
y = f(x) + np.random.normal(0, 1.5, n)
p = np.polyfit(x, y, 1)
slopes.append(p[0])
_ = plt.hist(slopes, 30, color=color, density=True)
se = np.sqrt(1.5**2/np.sum((x-np.mean(x))**2))
est_vals = np.linspace(1.5, 2.5, 50)
plt.plot(est_vals, est_dist(2, est_vals, se), "--", color=color, label=f"n={n}")
plt.legend()
Men om vi har målt noe, så har vi jo bare én måleserie. Hva kan vi gjøre da?
np.random.seed(16)
x = np.linspace(0, 10, 100)
y = f(x) + np.random.normal(0, 1, 100)
slopes = []
for i in range(1000):
selection = np.random.randint(0, len(x), 50)
x_c = x[selection]
y_c = y[selection]
p = np.polyfit(x_c, y_c, 1)
slopes.append(p[0])
_ = plt.hist(slopes, bins=30)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.savefig("figures/est_bootstrap.png")
Se nøye på figurene, hva skjer her?