# This source code is public domain

import numpy
import matplotlib.pyplot as plt

#Alter von(yr); bis(yr); Groesse(m); Gewicht(kg)
stat="""0; 1; 0.67; 7.7
1; 2; 0.83; 11.9
2; 3; 0.93; 14.4
3; 4; 1.01; 16.6
4; 5; 1.08; 18.7
5; 6; 1.15; 21.0
6; 7; 1.22; 23.9
7; 8; 1.27; 27.0
8; 9; 1.33; 30.4
9; 10; 1.39; 34.1
10; 11; 1.44; 37.8
11; 12; 1.50; 42.6
12; 13; 1.56; 47.5
13; 14; 1.63; 53.5
14; 15; 1.70; 59.9
15; 16; 1.75; 65.8
16; 17; 1.78; 69.8
17; 18; 1.80; 72.9
18; 20; 1.81; 75.5
20; 25; 1.81; 78.0
25; 30; 1.80; 80.7
30; 35; 1.80; 83.4
35; 40; 1.80; 84.8
40; 45; 1.80; 85.1
45; 50; 1.79; 85.7
50; 55; 1.78; 85.8
55; 60; 1.77; 85.9
60; 65; 1.76; 85.2
65; 70; 1.76; 84.4
70; 75; 1.74; 83.3
75; 80; 1.73; 79.0""".split('\n') #Daten von 2009

gewicht = [list(map(float,i.split(';'))) for i in stat]
x = [(i[1]+i[0])*0.5 for i in gewicht[18:-1]]
y = [i[3] for i in gewicht[18:-1]]

N = 3
A = numpy.zeros((N,N))
for i in range(N):
    for j in range(N):
        A[i,j] = sum(xi**(i+j) for xi in x)
b = numpy.zeros((N))
for i in range(N):
    b[i] = sum(xi**(i)*yi for xi,yi in zip(x,y))
c = numpy.linalg.solve(A, b)
xneu = numpy.sort(numpy.linspace(17, 74.5, num=50).tolist() + x)
yneu = numpy.sum([c[i]*xneu**i for i in range(len(c))],axis=0)

fig = plt.figure(figsize=(4.2, 3.2))
y1, = plt.plot(x,y,'o',label='Messpunkte')
y2, = plt.plot(xneu,yneu,'r-',label='Modellfunktion')
plt.xlabel('x (Alter in Jahren)')
plt.ylabel('y (Gewicht in kg)')
order = y1,y2
plt.legend(order,[p.get_label() for p in order],frameon=True, loc='lower right')
plt.grid(True, alpha=0.7)
plt.tight_layout()
plt.savefig('MDKQ3.svg')