from manim import *
from matplotlib import pyplot as plt
plt.rcParams['figure.figsize'] = [12, 9]
import scipy.optimize as opt
import numpy as np
import math
# jag skapar en lista "guesses" där parametrarna till alla anrop
# av vår modellfunktion kommer sparas - egna anrop såsom anrop via
# "curve_fit"
guesses = []
# här defineras modellfunktionen, en 3e-grads polynom
# "guesses" deklareras som global, dvs ändringar i listan görs globalt
def func(x, a, b, c, d):
global guesses
guesses.append([a, b, c, d])
return a*x*x*x + b*x*x + c*x + d
# lite förberedelser, en lista med x-koordinater, 100 punkter mellan 0 och 10
xdata = np.linspace(0, 10, 100)
# de ursprungliga parametrarna till min 3e-grads funktion som utgångspunkt, dvs
# de SANNA värden som vi ska försöka återskapa ur den brusiga mätdatan senare.
# Jag använde DESMOS för att hitta en funktion som såg rimligt ut
# f(x) = 0.04 x^3 - 0.4 x^2 + 0.2 x + 4.8
parameters = [.04, -0.4, 0.2, 4.8]
# sedan gäller det att beräkna de SANNA y-värden i listan "y"
y = func(xdata, *parameters)
# här skapas normalfördelade slumptal med en standardavvikelse på 0.5
# ( np.random_normal() har en standardavvikelse på 1.0 )
y_noise = .5 * np.random.normal(size=xdata.size)
# och här fuskar vi fram våra simulerade mätpunkter med de slumpmässiga mätfelen
ydata = y + y_noise
# vi ger curve_fit() en startpunkt på f(x) = x^3 + x^2 + x + 1
# därifrån ska sedan de parametrarna bestämmas, som passar bäst till våra mätpunkter
startparams = [1, 1, 1, 1]
# här görs kurvanpassningen: curve_fit() kommer anropa vår modellfunktion ett
# antal gånger på jakt efter paramterar som minimerar summan av de kvadratiska avstånden
# mellan vår mätpunkter och modellfunktionen
optimizedParameters, pcov = opt.curve_fit(f=func, xdata=xdata, ydata=ydata, p0=[1,2,1,1], method="trf", verbose=2)
#=========================================================================
# OCH NU BÖRJAR ANIMERINGEN MED MANIM (efter 3blue1brown, Grant Sanderson)
#
# listan "guesses[]" innehåller nu en serie med alla parametrar som har använts
# i funktionsanrop till vår modellfunktion.
#
# lista[0] är vårt anrop med de SANNA parametrarna
#
# till varje steg i anpassningen så utförde curve_fit() faktiskt 5 anrop av
# vår modellfunktion med bara små ändringar i parametrarna för att internt
# kunna beräkna de partiella derivatorna av summan av avståndskvadraten
#
class curvefit_x3(Scene):
def construct(self):
# Axes() skapar representationen av ett koordinatsystem
ax = Axes(
x_range=[-0.5,np.ceil(max(xdata)),1],
y_range=[-0.5,np.ceil(max(ydata)),1],
x_length=8,
tips= False,
)
ax.add_coordinates()
y_label = ax.get_y_axis_label("y", edge=LEFT, direction=LEFT, buff=0.4)
x_label = ax.get_x_axis_label("x")
# self.play() är den typiska animeringen där objekt ritas 3blue1brown-style
# här ska nu alltså ritas upp koordinatsystemet och texten längs axlarna
self.play(Create(ax),Create(x_label),Create(y_label))
# Dot() kan användas för att rita en punkt i grafen
# VGroup() skapar en sammanhängande lista av objekt, i det här
# fallet en lista med alla våra mätpunkter dvs punkterna som
# skapades med hjälp av de slumpmässiga felen
# ax.c2p(x,y) behövs så att punkterna hamnar på rätt position i koordinatsystemet
datapts = VGroup()
for i in range(len(xdata)):
datapts.add(Dot([ax.c2p(xdata[i],ydata[i])], radius=.05, color=RED))
# nu ritar vi in punkterna i grafen
self.play(Create(datapts))
# self.wait() fryser animeringen för den angivna tiden (i sekunder)
self.wait(1)
# DecimalNumber skapar ett objekt på scenen som visar ett numeriskt värde
# "sum" är namnet på objektet genom vilket vi kan ändra objektets egenskaper
# såsom värdet som ska visas
#
# MathTex() skapar ett text-objekt med LaTeX-formatering av texten
sum = DecimalNumber(0, num_decimal_places=2,group_with_commas=False)
eq = MathTex(r"\sum_{i=0}^n \left(y_i - f(x_i)\right)^2=")
# sum.next_to() placerar talet "sum" till höger om texten "eq"
sum.next_to(eq, RIGHT)
# jag grupperar ihop "eq" och "sum" till en ny grupp "sqsum" så
# att jag kan förflytta dessa två objekt tillsammans på scenen
sqsum = VGroup(eq,sum)
sqsum.to_edge(UP)
self.play(Create(sqsum))
# jag gör en kopia av vår lista över parametrar till vår
# modellfunktionen - detta på grund av att vi nu kommer
# själva anropa modellfunktionen upprepade gånger som varje gång skulle lägga
# till nya element i listan som vi inte vill ha med i animeringen
local_guesses = guesses.copy()
params = startparams
# lambda-funktioner är fortfarande lite suspekta för mig i Python, men
# här används en sådan för att rita upp den första kurvan som
# är vår modellfunktion med startvärden som parametrar som gavs
# som startvärden till curve_fit()
graph = ax.plot(
lambda x: func(x, *params),
color = BLUE,
)
self.play(Create(graph))
# och längst nere på skärmen skrivs det ut den faktiska ekvationen med ytterligare ett MathTex-objekt
eq = MathTex(r"f_(x)= {:.3f}\,x^3 {:+.3f}\,x^2 {:+.3f}\,x {:+.3f}".format(*params), color=BLUE)
eq.to_edge(DOWN, buff=0.1)
self.play(Create(eq))
# nu skapas det en VGroup() som kommer innehålla de gula linjerna
# som ritades mellan varje röd märpunkt och den motsvarande platsen längs
# den blåa funktionsgrafen
# För varje x-koordinat xdata[i] räknas det fram det motsvarande funktionsvärdet fi
# "start" och "end" är de motsvarande start- och slutpunkterna av de gula linjerna
# och värdet i "sum" ökas med kvadraten av linjens längd
#
# detta är alltså den första sådan ritning som då gäller för startvärdet hos parametrarna
fit = VGroup()
for i in range(len(xdata)):
fi = func(xdata[i],*params)
start = ax.c2p(xdata[i],ydata[i])
end = ax.c2p(xdata[i],fi)
fit.add(Line(start=start, end=end, color=YELLOW))
self.add(fit)
self.wait(.2)
sum.set_value(sum.get_value() + (ydata[i]-fi)**2)
# nu stegar vi oss fram genom hela listan av paramterar som modellfunktionen
# anropades med. Eftersom det var 5 anrop för varje steg i kurvanpassningen
# gör vi en slice där vi bara tar vart femte anrop och vi utelämnar det första anropet
for params in local_guesses[1::5]:
# samma procedur som ovan
graph1 = ax.plot(
lambda x: func(x, *params),
color = BLUE,
)
newsum = 0
fit1 = VGroup()
for i in range(len(xdata)):
fi = func(xdata[i],*params)
start = ax.c2p(xdata[i],ydata[i])
end = ax.c2p(xdata[i],fi)
newsum += (ydata[i]-fi)**2
fit1.add(Line(start=start, end=end, color=YELLOW))
eq1 = MathTex(r"f_(x)= {:.3f}\,x^3 {:+.3f}\,x^2 {:+.3f}\,x {:+.3f}".format(*params), color=BLUE)
eq1.to_edge(DOWN, buff=.1)
self.play(sum.animate.set_value(newsum),Transform(graph,graph1),Transform(fit,fit1),Transform(eq,eq1))
self.wait(3)
# färdigt, vi har gått genom hela sekvensen av kurvanpassningen och ritar en blå ram runt
# den slutgiltiga funktionen
srect = SurroundingRectangle(eq, color=BLUE)
self.play(Create(srect))
self.wait(2)
# sista steget är att rita den SANNA funktionen, dvs den som vi ursprungligen började med...
graph2 = ax.plot(
lambda x: func(x, *parameters),
color = RED,
)
self.play(Create(graph2))
self.wait(5)
Comments