# clt.py - hoe goed voorspelt het gemiddelde van n dagen? (model voor de zonnestroom in april)
import numpy as np

rng = np.random.default_rng(1)
EMAX, A, B = 7.6, 2.2, 1.035               # heldere dag 7,6 kWh; vorm van de verdeling (scheef)
mu = EMAX * A / (A + B)                   # verwachting: 5,17 kWh (PVGIS voor april)
sigma = EMAX * np.sqrt(A * B / ((A + B) ** 2 * (A + B + 1)))
dagen = EMAX * rng.beta(A, B, size=(20000, 30))   # 20 000 keer een maand van 30 dagen

print("mu = %.2f kWh, sigma = %.2f kWh" % (mu, sigma))
print(" n   gemiddelde   s van de gemiddelden   sigma/sqrt(n)   P(gem < 4 kWh)")
for n in (1, 2, 7, 30):
    gem = dagen[:, :n].mean(axis=1)       # het gemiddelde van de eerste n dagen van elke maand
    print("%2d   %8.3f   %14.3f   %18.3f   %11.4f"
          % (n, gem.mean(), gem.std(ddof=1), sigma / np.sqrt(n), np.mean(gem < 4.0)))