Fehlerrechnung mit Mathematica

Autor: Julien Kluge

Vorwort

Hey Leute,

dieses Dokument ist gedacht für Leute, die mit Mathematica anfangen wollen für das Praktikum Fehlerrechnung zu lernen. 
Es setzt eigentlich nur voraus, dass ihr eine kleines bisschen Mathematica schonmal gesehen habt und wisst wie man Variablen setzt,
Befehle ausführt und damit rechnet. Ansonsten erkläre ich hier alles von Grund auf.

Warum Fehlerrechnung mit Mathematica in diesem Notebook hier lernen?
Alle Beispiele könnt ihr selbst ausführen, verändern, daran herumprobieren, usw.!

Anders als in der PDF “Kurzeinführung in die Fehlerrechnung für das Physik-Praktikum”, werde ich auch hier die Beispiele weiter
ausführen und gegebenenfalls tiefer in die Materie gehen. Ihr könnt das überspringen, wenn ihr das nicht genauer wissen wollt.
Allerdings ist dieses Notebook auch kein Ersatz. Ich schreibe nicht alles über die Formeln, sondern es geht hier mehr darum,
wie man das in Mathematica einsetzt.

Zu guter Letzt möchte ich noch kurz ein paar Worte darüber verlieren wieso Mathematica und nicht Python, Julia, Excel, etc. ...
Mathematica ist als Sprache explizit auf symbolische Auswertung ausgelegt. Die Zeitersparniss die ihr dadurch bekommt, ganz
zu schweigen von den vielen Möglichkeiten, kann keine andere Sprache bieten die ich bisher kenne. Ich möchte dazu ein kleines
Beispiel aus meiner Forschung erwähnen (wärend ich das hier schreibe, bin ich im zweiten Jahr meines Doktors):
Stellt euch vor, ihr habt ein System an drei Differentialgleichungen, die die Temperatur eines Apparates an drei verschiedenen
Stellen beschreiben. Die DGL’s sind nichtlinear, gekoppelt, haben keine analytische Lösung und enthalten zu allen Übel noch einen
sogenannten Verzögerungsterm (sowas wie ((∂f(t - (Subscript[t, 0])))/(∂t))). Jetzt sind in diesen DGL’s Konstanten, die mein System beschreiben, die ich gerne hätte. 
Mit Unsicherheiten, wenn es geht. Also was muss ich tun:
- Messen meines Systems um drei Datenströme zu bekommen für meine drei DGL’s. Leider kann ich aber nur an zwei der drei Orte messen. Also müssen wir damit klar kommen.
- Ich muss meine DGL’s lösen. Das geht aber nicht analytisch also muss ich das numerisch machen. Dafür muss ich aber Werte für
  meine Konstanten angeben, die ich ja finden will.
- Die Daten müssen also an die DGL’s quasi multidimensional gefittet werden. Um das zu tun muss in jedem “Fit-Schritt” das gesamte DGL System gelöst und
  interpoliert werden, um an den Datenpunkten ausgerechnet zu werden. Damit optimiert man die Parameter durch maximieren der Log-Likelihood des Systems.
- Dabei wird der Einfluss der zweiten DGL nur an die Kopplung der ersten und dritten optimiert da wir hier keine Datenpunkte haben.
- Um Die Unsicherheiten zu bestimmen muss das gelöste und optimierte System merhmals abgeleitet werden, um die Hesse
  Matrix zu bestimmen. In diese Matrix kann man dann die Daten einsetzen. Die Matrix muss anschließend invertiert werden.
- Hat man all das geschafft müssen natürlich auch schöne Plots und Ergebnisse her.
Klingt jetzt erstmal alles super böse. Aber es reicht wohl zu sagen, dass man all diese Schritte, in 5 Zeilen Mathematica code hinsetzen
kann. fünf Zeilen
Dass Mathematica mir das einfach macht und all diese Sachen für mich in fünf Zeilen hinsetzt, ist der Grund dass 80% meiner Arbeit im Doktor
über Mathematica geht (10% geht noch an C, C++, Cuda-C für high-performance Simulationen und 8% Python an Steuerprogramme
für Experimente (ja mir ist klar dass da noch 2% fehlen...über die reden wir nicht)).
Ist natürlich nicht so cool dass Mathematica nicht Open-Source ist, aber ganz ehrlich,
ich würde es nicht missen wollen.

- 02.06.2021
PS: Ich hoffe wenn ihr das hier lest, ist COVID-19 ein Ding der Vergangenheit.

Voraussetzungen

Ihr braucht mindestens Mathematica version 12.3 oder höher.
Ich habe dieses Notebook mit der Version:

$Version

geschrieben. Ebenfalls benutze ich für alle Plots das Scientific-Theme. Das können wir also gleich mal global einstellen:

$PlotTheme="Scientific";

Die Einzelnen Kapitel kann man aufklappen. Einfach auf den Pfeil links daneben drücken.

Mittelwerte & Pythagoreische Addition

Arithmetischer Mittelwert (Mean)

Der arithmetische Mittelwert  der werte  mit  ist gegeben mit:

Wir definieren uns jetzt also einen Datensatz mit allen Werten :

data={1.2,3.4,5.2,3.2,2.5,3.3}

Die Anzahl der Datenpunkte  kann man über den Length[] Befehl finden:

Length[data]

Also ist . Die Summe einer Liste kann man in Mathematica mit Total[] finden:

Total[data]

Wenn man nun beides verbindet bekommt man den Mittelwert mit:

Total[data]/Length[data]

Mathematica kennt allerdings einen fertigen Befehl um das für euch zu machen. Mean[]

Mean[data]

Gewichteter Mittelwert

Manueller gewichteter Mittelwert

Man stelle sich folgende Situation vor: man hat die Werte  und möchte aus denen ein
Mittelwert bilden. Jetzt weiß man aber, dass einige dieser Werte viel genauer und/oder
weniger unsicher sind als andere. Da hätte man es doch lieber, diese Werte stärker in den
Mittelwert eingehen zu lassen. Dafür existiert der gewichtete Mittelwert.
Starten wir damit unsere Daten und Unsicherheiten in eine Liste zu speichern:

data={1.2,3.4,5.2,3.2,2.5,3.3} (*Daten*)
uncert={0.5,0.2,0.8,0.1,0.4,0.2} (*Unsicherheiten*)

Um nun zu gewichten, muss klar sein wonach gewichtet wird. Üblicherweise benutzt man
die instrumentelle Gewichtung. Diese ist das inverse der statistischen Varianz. Die Varianz
ist das Quadrat der Unsicherheit. Also die Gewichte  wobei  die Unsicherheiten
und  demnach die Varianzen zu den Werten  sind.

varnce=uncert^2 (*Varianz*)
gewicht=1/varnce (*Gewichtungen*)

Der gewichtete Mittelwert  und der damit assoziierte Fehler  sind gegeben mit

Und so rechnen wir also beides mit Total[] und Sqrt[] (Sqrt: Quadratwurzel) aus :

Total[gewicht*data]/Total[gewicht]
1/Sqrt[Total[gewicht]]

Unser Ergebnis ist also  und  oder zusammen geschrieben: .

Für die Zukunft macht es Sinn, sich selbst eine Funktion zu schreiben, die das automatisch berechnen kann:

WeightedMean[data_,uncert_]:={Total[(1/uncert^2)*data]/Total[1/uncert^2],1/Sqrt[Total[1/uncert^2]]}

Damit können wir einfach rechnen:

WeightedMean[data,uncert]

Gewichteter Mittelwert mit MeanAround[]

Dieses Kapitel macht erst Sinn, wenn ihr mit Around[] umgehen könnt.
Dieser Befehl wird erst später in diesen Notebook erklärt. Kommt wieder zurück wenn ihr den Befehl könnt.

Wir starten mit den selben Daten wie im manuell berechneten Fall aber formen diese in Around-Objekte:

rawdata={1.2,3.4,5.2,3.2,2.5,3.3}; (*Daten*)
uncert={0.5,0.2,0.8,0.1,0.4,0.2}; (*Unsicherheiten*)
dataAndUncert=Transpose[{rawdata,uncert}] (*Wir benutzen Transpose[] um die Daten mit den Unsicherheiten in eine Matrix zu schreiben. Also Transpose[{{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),...},{(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3]),...}}]={{(Subscript[x, 1]),(Subscript[u, 1])},{(Subscript[x, 2]),(Subscript[u, 2])},{(Subscript[x, 3]),(Subscript[u, 3])},...}*)
data=Around@@@dataAndUncert (*Beides in Around gießen. Die @@@ sind ein Operator. Der wirkt folgendermaßen: Around@@@{{a,b},{x,y}} -> {Around[a,b],Around[x,y]}*)

Hat man erstmal die Daten in der Around-Form, kann man einfach MeanAround[] benutzen. Wenn MeanAround sieht, dass es Around Objekte bekommt, macht es automatisch eine instrumenteller Gewichtung:

MeanAround[data]

Und Fertig.

Pythagoreische Addition

Manuelle phytagoreische Addition

Wenn man beispielsweise mehrere Unsicherheiten  mit  für einen Unsicherheits-Wert  zusammenführen will, muss das geometrisch geschehen.
Dafür benutzt man die phytagoreische Addition. Sie ist folgendermaßen definiert:

Der Grund dafür wird nachher im Distributionstheorie-Kapitel mehrmals aufgeführt.
Wir können also einfach die bekannten Befehle benutzen. 
Beginnend damit Unsicherheiten  in eine Liste zu speichern:

uncert={0.5,0.3,0.2,0.8,0.4} (*Unsicherheiten*)

Dann wird einfach berechnet

Sqrt[Total[uncert^2]]

Phythagoreische Addition mit Around[]

Dieses Kapitel macht erst Sinn, wenn ihr mit Around[] umgehen könnt.
Dieser Befehl wird erst später in diesen Notebook erklärt. Kommt wieder zurück wenn ihr den Befehl könnt.

Wir können uns dem Fakt bedienen dass die Fehlerfortpflanzung (siehe das Kapitel dazu) einer Addition
und die phytagoreische Addition identisch sind:

Around[x1,u1]+Around[x2,u2]+Around[x3,u3]+Around[x4,u4]

Wenn man nun also  setzt, kann man zu beliebigen Werten Unsicherheiten phytagoreisch dazu addieren.
Im folgenden Beispiel addieren wir zu  zwei Unsicherheiten  und  phytagoreisch dazu:

Wert=Around[5.687,0.04]
finalerWert=Wert+Around[0,0.02]+Around[0,0.05]

Das stimmt mit der manuellen Berechnung überein:

finalerWert["Error"]
Sqrt[0.04^2+0.02^2+0.05^2]

Standardabweichung

Ist ein Messwert mit einem zufälligen Fehler belegt, so streuen diese um einen gemeinsamen Punkt (ohne Systematische Abweichungen ist dieser Punkt der Wahrheitswert).
Dabei nimmt diese Streuung die Form einer Gaußschen Normalverteilung (auch Glockenkurve genannt) an.
Mathematica kennt diese Verteilung unter dem befehl NormalDistribution[] und kann die Verteilungsfunktion mit PDF[] aufrufen. Dabei müssen zwei Parameter angegeben werden.
Der Mittelwert  um den die Punkte streuen und eine Breite wie weit die Punkte um den Mittelwert streuen, die sogenannte Standardabweichung :

mittelwert=5;
standardabweichung=2;
distFormel=PDF[NormalDistribution[mittelwert,standardabweichung],x]
Plot[distFormel,{x,mittelwert-5*standardabweichung,mittelwert+5*standardabweichung}]

Die Form einer Glocke erklärt sich leicht, wenn man annimmt, dass Streuungen mit steigendem Abstand vom Mittelwert  zunehmend unwahrscheinlicher werden.
Die Wahrscheinlichkeit selbst irgendwo in dieser Kurve zu landen ist gegeben durch die Fläche unter der Kurve (Also das Integral).
Es gilt also das die Gesamtfläche 100% sein muss:

Integrate[distFormel,{x,-∞,∞}]

Man sieht leicht, dass diese symmetrisch um den Mittelwert  ist. Damit können nun symmetrische Intervalle festgelegt werden in denen eine bekannte
Wahrscheinlichkeit existiert, dass eine Messung dort landet. Diese Intervalle werden mit einem vielfachen  der Standardabweichung  um den Mittelwert 
festgelegt. Die Wahrscheinlichkeit in diesen Intervallen zu landen ergibt sich mit der, in diesem Interval, eingeschlossenen Fläche.

percentages=Integrate[PDF[NormalDistribution[μ,σ],x],{x,μ-n*σ,μ+n*σ},Assumptions->σ∈Reals]

Grid[{{"Interval ±σ","statistische Sicherheit"}}~Join~Transpose[{
Range[1,5],
PercentForm/@(percentages/.n->Range[1.0`5,5.0])
}],Frame->All]

Wir rechnen in der Physik für gewöhnlich mit einer statistischen Unsicherheit von einer Standartabweichung (ein Sigma, ).
Hat man eine Messreihe gemacht, kann die empirische Standartabweichung  davon berechnet werden mit:

Wir berechnen also nun mal beispielhaft an  gemessenen Erdbeschleunigungen:

gs={9.81279,9.83616,9.76557,9.78496,9.84230,9.75096,9.72767,9.84866,9.74724,9.76847};
mittelwertGS=Mean[gs];
anzahl=Length[gs];
stdAbweichung=Sqrt[1/(anzahl-1)*Total[(gs-mittelwertGS)^2]]

Mathematica bietet aber auch gleich einen fertigen Befehl dafür, StandardDeviation[]:

StandardDeviation[gs]

Extra: Der n-1 Faktor:

Dieses Unterkapitel ist extra und nur für die Neugierigen gedacht.

Habt ihr euch mal gefragt wieso in der empirischen Standardabweichung der inverse Faktor  ist und nicht nur  ?
Der Grund dafür liegt darin, dass ihr den Mittelwert  eigentlich auch nicht kennt sondern nur abschätzen könnt. Tatsächlich aber,
ist die Wahl des arithmetischen Mittelwerts hier nicht förderlich, da dieser die Standardabweichung IMMER minimiert.

Wir können das einfach nachprüfen, indem wir einfach dummy daten in die Definition der Standardabweichung mit einem
unbekannten Vorfaktor einsetzen und berechnen bei welchen Mittelwert mu das ganze extremal wird (Ableiten, null setzen und lösen).

data={x1,x2,x3};
std=Sqrt[Vorfaktor*Total[(data-mu)^2]];(*Standardabweichung*)
Solve[D[std,mu]==0,mu]

Und voilà, wir sehen der Arithmetische Mittelwert macht die Standardabweichung extremal (minimal in diesem Fall).
Ahh das ist schlecht. Um dem Abhilfe zu schaffen, kann man tatsächlich den Faktor von  auf  ändern und das gibt
einen besseren Schätzer. Wir können das auch nachprüfen.

Sei eine Normalverteilung mit  und  gegeben.

mean=100;
std=1;
dist=NormalDistribution[mean,std];
data=RandomVariate[dist,5]; (*Diese Funktion verteilt zufällig 5 Datenpunkte anhand der gegebenen Verteilung*)

Also vergleichen wir jetzt zwei Fälle:
unbekannter Mittelwert, n-Faktor
unbekannter Mittelwert, (n-1)-Faktor

Sqrt[1/(Length[data])*Total[(data-Mean[data])^2]]
Sqrt[1/(Length[data]-1)*Total[(data-Mean[data])^2]]

Wir sehen ganz klar, mit Faktor  unterschätzen wir die Standardabweichung, mit nicht. (Natürlich ist das Experiment hier dem Zufall unterworfen. Aber im Mittel wird sich dieses Ergebnis bestätigen)

Es gibt tatsächlich viel interessante Literatur wenn es um diesen  Faktor geht. Oft wird auch von Freiheitsgraden gesprochen. Dadurch dass man den Mittelwert auf einen Wert festmacht, verliert man
einen Freiheitsgrad von den  vielen, die man von den  Datenpunkten hat. Das wird uns beim Fitten auch nochmal entgegenkommen.

Der Vollständigkeit halber: hier nochmal 1millionen mal getestet und gemittelt:

Mean@Table[{Sqrt[1/(Length[#])*Total[(#-Mean[#])^2]],Sqrt[1/(Length[#]-1)*Total[(#-Mean[#])^2]]}&@RandomVariate[NormalDistribution[100,1],5],{i,1,1000000}]

Interessanterweise unterschätzen beide die tatsächliche Standardabweichung. Das hat damit zu tun, dass die empirische Standardabweichung ein minimierer unter den gegebenen Daten ist.
Das kann man mit einem sogenannten Student T-Wert noch korrigieren. Dazu im nächsten Kapitel was. Mit diesem, wird man aber sehen, dass die (n-1) Gewichtung stimmt.

Vertrauensbereich

Mithilfe der Standardabweichung lässt sich, in einer bestimmten statistischen Wahrscheinlichkeit,
die zufällige Messunsicherheit abschätzen. Diese Abschätzung nennt sich Vertrauensbereich  und
ist folgendermaßen definiert:

Der Vorfaktor , Student-T Faktor genannt, ist nicht immer nötig. Im Folgenden die Kapitel zur Unterscheidung:

Fall  manuell

Für viele Datenpunkte lässt sich annehmen, dass  ist und man kann stumpf rechnen:

data={9.81279,9.83616,9.76557,9.78496,9.84230,9.75096,9.72767,9.84866,9.74724,9.76847};
StandardDeviation[data]/Sqrt[Length[data]]

Fall  mit MeanAround[]

Dieses Kapitel macht erst Sinn, wenn ihr mit Around[] umgehen könnt.
Dieser Befehl wird erst später in diesen Notebook erklärt. Kommt wieder zurück wenn ihr den Befehl könnt.

Wir starten mit den selben Daten wie im manuell berechneten Fall. Wir können einfach MeanAround darauf
anwenden und bekommen ein Around-Objekt zurück mit abgschätzten Mittelwert und Vertrauensbereich als Unsicherheit:

data={9.81279,9.83616,9.76557,9.78496,9.84230,9.75096,9.72767,9.84866,9.74724,9.76847};
MeanAround[data]

Fall

Der Student-Faktor (benannt nach William Sealy Gosset der unter dem Pseudonym Student gearbeitet hat) kann auf komplizierter Weise,
von der Student T-Verteilung berechnet werden. Diese Verteilung besitzt nur die Anzahl der Messwerte n (eigentlich Freiheitsgrade: ν = n − 1)
als Parameter und ist um null symmetrisch. Er wird nötig da der empirische Vertrauensbreich (eigentlich Standardabweichung) ebenfalls wieder
ein Extremal der gegebenen Daten und nicht der darunterliegenden Verteilung ist.
In den Einführungspraktika und Grundpraktikum I & II ist es erlaubt die Näherung  für n ≥ 6 zu machen. Will man diese Näherung nicht
eingehen, kann man den tatsächlich relativ leicht in Mathematica berechnen. Die Formel dazu ist relativ kompliziert und nicht näherer
Gegestand dieses Scripts. Es reicht das andere sich die Arbeit gemacht haben das zu implementieren. Man kann die Funktionzur Berechnung
aus dem Wolfram-Function-Repository direkt in das Notebook laden:

StudentTValue=ResourceFunction["StudentTValue"]
Information[StudentTValue]

Wir können diese Funktion nun nutzen. Nehmen wir mal an, wir haben einen Wert dreimal gemessen.
Nun berechnen wir seinen Student-T wert mit dem Argument der Freiheitsgrade :

data={3.2,2.8,3.3};
tWert=StudentTValue[Length[data]-1]

Diese Formel sieht schonmal böse aus. Aber wir können Mathematica einfach nach einen numerische Wert fragen:

N[tWert]

Aha. Unser Vertrauensbreich wäre mit 1 also in etwa 32.1% zu klein gewählt. Vergleichen wir also

StandardDeviation[data]/Sqrt[Length[data]] (*t=1*)
tWert*StandardDeviation[data]/Sqrt[Length[data]] (*t=1.32128*)

Naja das ist doch schonmal ein Unterschied.

Fall  mit MeanAround[]

Dieses Kapitel macht erst Sinn, wenn ihr mit Around[] umgehen könnt.
Dieser Befehl wird erst später in diesen Notebook erklärt. Kommt wieder zurück wenn ihr den Befehl könnt.

Ihr solltet vorher auch das  manuell Kapitel gelesen haben.
Die gegebene Funktion kann uns direkt den StudentT Wert dazu multiplizieren.
Wir schreiben also erst die alten Daten auf und rufen die Funktion ab:

data={9.81279,9.83616,9.76557,9.78496,9.84230,9.75096,9.72767,9.84866,9.74724,9.76847};
value=MeanAround[data]
StudentTValue=ResourceFunction["StudentTValue"]

Und nun wenden wir an:

StudentTValue[Length[data]-1,value]

Die Funktion kann sogar die Rolle von MeanAround übernehmen:

StudentTValue[data]

Kochrezept Messwertberechnung & Around[]

Manuell

Hat man nun sukzessiv die erste Messreihe für einen Wert durchgeführt, gibt es ein Kochrezept wie wir diesen sowie seinen Fehler ordentlich berechnen können.
Folgende Daten sind gegeben:

data={9.81279,9.83616,9.76557,9.78496,9.84230,9.75096,9.72767,9.84866,9.74724,9.76847};

Und nun Schritt für Schritt:

Mittelwert berechnen (Mean[])

mean=Mean[data]

Bekannte systematische Fehler korrigieren

gSyst=-0.02; (*Wir haben aufgrund von Gründen immer 0.02m/(s^2) zuwenig auf dem Ergebnis*)
meanK=mean-gSyst

Vertauensbreich berechnen

uStat=StandardDeviation[data]/Sqrt[Length[data]] (*statistische Unsicherheit/Vertrauensbereich, t=1*)

Messunsicherheit durch phytagoreisches Addieren aller Unsicherheiten berechnen:

uAblese=0.004*mean; (*es gab einen 0.4% Ablesefehler*)
uGeraet1=0.03+0.004*mean; (*es gab einen Gerätefehler von 0.4%+3dgt.*)
uGeraet2=0.008; (*und einen statitischen Gerätefehler von 0.008 m/(s^2)*)
u=Sqrt[uStat^2+uAblese^2+uGeraet1^2+uGeraet2^2]

Ergebnis runden und mit SI-Präfix und Einheit angeben

Roh: 9.80848 ± 0.081053

Endergebnis :  m/(s^2)

Einführung in Around[]

Seit Version 12 kann Mathematica selbst Fehlerrechnungen machen. Das ist eine erhebliche Erleichterung. Dazu wird der Around[] befehl benutzt.
Dieser braucht den Mittelwert und eine Unsicherheit und verhält sich danach so als wäre es eine Zahl.

Around[wert,unsicherheit]

Mit Zahlenwerten zeigt Around sogar gerundete (nach ISO) Werte an, behält aber im Hintergrund die ungerundeten Werte bei.

value=Around[9.8156173,0.0282579]

Mit “Value” und “Error” kann man die Werte vollständig ungerundet extrahieren:

value["Value"]
value["Error"]

Oder man wendet einfach FullForm[] auf Around an

FullForm[value]

Wie gesagt kann man Around-Objekte wie Zahlen behandeln. Also auch ganz normal mit ihnen Rechnen:

Log[value+5]*2

Dabei arbeitet Mathematica im Hintegrund mit Fehlerfortpflanzung. Dazu in einem späteren Kapitel mehr.

Auch mehrere Around-Objekte können miteinander verrechnet werden:

value1=Around[5.321,0.252];
value2=Around[0.7249,0.0511];
(Pi-value1)/value2+1

Man könnte also zum Beispiel die Fehler einer Datenserie berechnen mit:

data={1.2,3.4,5.2,3.2,2.5,3.3};
value=Around[Mean[data],StandardDeviation[data]/Sqrt[Length[data]]]

Im nächsten Unterkapitel wird sich aber zeigen, dass das noch einfacher geht.

Man kann diese Werte sogar mit Unsicherheitsbalken plotten lassen:

ListPlot[{
{1,value},
{Around[2,0.3],value1},
{3,value2*3}
}]

Messwertberechnung mit Around[]

Bevor ihr hier weiter lest, geht am besten nochmal durch alle Kapitel zurück und lest euch dass durch, dass ihr verpasst habt weil ihr dafür Around[] gebraucht hättet.

Wir starten mit den Daten aus den Manuellen Unterkapitel

data={9.81279,9.83616,9.76557,9.78496,9.84230,9.75096,9.72767,9.84866,9.74724,9.76847};

Wir können direkt mit MeanAround and Mittelwert mit Vertrauensbereich kommen:

value=MeanAround[data]

Dann fügen wir noch die drei zusätzlichen Unsicherheiten durch pythagoreisches addieren hinzu und korrigieren den systematischen Fehler.

gSyst=-0.02; (*systematischer Fehler*)
uAblese=0.004*Mean[data]; (*es gab einen 0.4% Ablesefehler*)
uGeraet1=0.03+0.004*Mean[data]; (*es gab einen Gerätefehler von 0.4%+3dgt.*)
uGeraet2=0.008; (*und einen statitischen Gerätefehler von 0.008 m/(s^2)*)

kValue=value-gSyst (*korrigieren des syst. Fehlers*)
finalValue=kValue+Around[0,uAblese]+Around[0,uGeraet1]+Around[0,uGeraet2] (*phytagoreische Addition der zusätzlichen Unsicherheiten*)

Und fertig [naja denkt euch die Einheiten dazu].
Es sei noch dazu gesagt, das Mathematica nie strikt aufrundet am Fehler. Was dahingehend okay ist dass einige Standards (ISO) das auch so machen. DIN will aber aufrunden.
Wir können uns den Fehler genau angucken um zu sehen ob wir aufrunden müssen:

finalValue["Error"]

Wenn wir also nach DIN angeben wollen, dann müssten wir

schreiben.

Fehlerfortpflanzung

Wir werden ab hier nurnoch Around[] benutzen. Natürlich solltet ihr ein paar mal im Leben manuell eine Fehlerfortpflanzung machen.
Aber wenn man es ein paar mal selbst gemacht hat, ist es wie Fahrrad fahren. Ab da an lassen wir es uns von Mathematica übernehmen.

Fortpflanzung mit Around[]

Wie in der Around[] Einführung bereits erwähnt, pflanzt Around automatisch anhand der Gaußschen Fehlerfortpflanzung fort.

a=Around[5.6789,0.3456];
b=Around[1.2345,0.6789];
result=a/b
result["Error"]

Wir überprüfen einmal durch manuelles rechnen:

Sqrt[(uA*D[A/B,A])^2+(uB*D[A/B,B])^2]/.{
A->5.6789,
uA->0.3456,
B->1.2345,
uB->0.6789
}

Okay cool passt. Ab hier könnt ihr das in beliebig komplizierte Formeln kloppen:

a=Around[5.6789,0.3456];
b=Around[1.2345,0.6789];
c=Around[4.16348,3.1961];
d=Around[0.0156,0.00516];
ArcSin[0.1*a-b/Gamma[c]]/Log[d]

Around kann sogar symbolisch benutzt werden:

Around[x,ux]/Around[y,uy]+Exp[Around[z,uz]]

Das scheint jetzt alles zu einfach, zu cool zu sein. Da muss es doch einen Haken geben?!
Und tatsächlich gibt es einen winzig kleinen.
Mathematica nimmt immer an, dass alle Größen die über Around fortgepflanzt sind, vollständig unkorreliert sind.
Das gilt auch für Variable mit sich selbst. Was heißt das jetzt genau? Naja zwei Sachen:
- Korrelierte Variablen (was auch immer das ist) können nicht mit Around bearbeitet werden.
- Eine Variable darf niemals zweimal in einer Formel stehen wenn sie mit Around benutzt wird.
Wärend uns der erste Punkt noch nicht viel sagt, ist der zweite schon eine Einschränkung.

Man nehme beispielsweise folgendes Beispiel:

a=Around[5,1];
a-a

Nun, eine Variable von sich selbst abzuziehen, muss IMMER  ergeben. Denn egal was man misst, egal wie ungenau
man ist, wenn man eine Variable von sich selbst abzieht, dann ist das Wert-Wert, oder Wahrheitswert-Wahrheitswert.
Das ist immer präzise genau null. Aber Mathematica reported uns ein . Dass passiert da Mathematica annimmt
das beide Seiten am Minuszeichen unkorreliert sind. Also eigentlich berechnet Mathematica das hier:

a=Around[5,1];
b=Around[5,1];
a-b

Und jetzt stimmt das, voraussgesetzt, a und b sind unkorreliert. 
Ihr meint dieser Fall passiert nicht? Das passiert allein im Grundpraktikum 1&2 eine handvoll male.
(Ein Beispiel dazu stelle ich auch einmal im “Reale Beispiel” Kapitel dar)
Wie lösen wir das? Siehe Nächstes Kapitel AroundReplace[]

Fortpflanzung mit AroundReplace[]

Wann immer korrelationen ins Spiel kommen, unabhängig ob zwischen zwei verschiedenen Variablen oder ein und der selben,
muss AroundReplace[] benutzt werden. Dieser Befehl führt die Gaußsche Fehlerfortpflanzung für eine gegebene Formel aus.

?AroundReplace

Ansonsten verhält sie sich genau wie die Rechnungen mit Around.
Formel symbolisch angeben und dann eine Liste an Regeln angeben mit was die Symbole ersetzt werden sollen:

a=Around[5,1];
b=Around[5,1];

AroundReplace[A/B,{A->a,B->b}]

Das führt im unkorrelierten Fall erstmal zum selben Ergebnis wie nur mit Around allein:

a/b

Wie gesagt aber, beachtet sie Selbskorrelationen:

AroundReplace[A-A,{A->a}]

Wie man mit der Formel allgemeine Korrelationen zwischen unterschiedlichen Variablen beachtet, siehe das Hauptkapitel dazu.

Korrelierte Fehlerfortpflanzung

Siehe zwei Hauptkapitel weiter.

Regressionen

Mathematica bietet mehrere Funktionen zum allgemeinen Fitten von Daten an. Davon sind
nur wirklich zwei für uns nützlich:

Information[LinearModelFit]
Information[NonlinearModelFit]

LinearModelFit

Generell

Als erstes räumen wir mal das Gerücht aus dem Weg, das lineare Fits, lineare Funktionen erzeugen. Das ist Blödsinn.
Lineare Fits sind Regressionen in denen die Fitparameter linear kombiniert mit unabhängigen Basisfunktionen sind.
Klingt kompliziert, ist aber eigentlich ganz einfach. Z.B.
 
kann zum Beispiel linear gefittet werden. Warum? Nun, folgende Parameter mit ihren jeweiligen Basisfunktionen Paare hat man:
 mit 
 mit 
 mit 
Und da ,  und  linear unabhängig sind und jeweils mit einen parameter kombiniert sind, können wir das linear fitten.

Um das zu machen, erstellen wir erstmal ein paar zufällige Daten nach diesen Modell (Im Sinne der Transparenz: ):

(*Keine Sorge wenn der code hier zu hart zu verstehen ist. Der ist nur dafür da halbwegs sinnvolle Zufallsdaten zu erzeugen.*)
errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,0.03*Exp[2*x]-3*x-3*2},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];

ListPlot[data,IntervalMarkers->"Fences"]

LinearModelFit nimmt nur drei Funktionargumente entgegen.
1. Die Daten
2. Die Basisfunktionen
3. Die Funktionsvariable(n)

Ein weiteres Wort zuvor: In dieser Mathematica Version (12.3.0), nimmt LinearModelFit, sowie auch NonlinearModelFit KEINE Around Objekte entgegen.
Ich bin mir sicher, dass wird sich in zukünftigen Versionen ändern aber wer weiß ob ich das update hier. Also schaut mal selbst.
Allerdings können wir das Around sehr schnell herauslösen indem wir den Daten eine Regel einsetzen:

dataOhneAround=data/.{Around[v_,_]:>v};
fit=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x]

Nach einem Fit, liefert Mathematica ein solches FittedModel-Objekt. Das tut es, weil man eine MEEEENGEEEE Informationen aus so einem Fit ziehen kann.
Die sind alle in diesem Objekt hinterlegt. Hier mal ein paar wichtige:

fit["BestFit"] (*Das Fitmodell mit den Besten Parametern eingesetzt*)

fit["BestFitParameters"] (*Die besten fit Parameter. Zur Erinnerung, die Wahrheitswerte waren 0.03, -3, 2*)

fit["ParameterErrors"] (*Die Unsicherheiten zu den parameters*)

fit["ParameterTable"] (*Eine nette Zusammenfassung der Parameterergebnisse*)

Wir extrahieren die Parameter mal als Around Objekt

values=fit["BestFitParameters"]
errors=fit["ParameterErrors"]
Around@@@Transpose[{values,errors}]

Weiter, kann das FittedModel-Objekt auch als Funktion benutzt werden. Dabei tut es so, als wäre es die BestFit-Funktion:

fit[1.5]

Somit ist das plotten sehr einfach. (Mehrere Plots kann man mit Show[] verbinden)

Show[
ListPlot[data,IntervalMarkers->"Fences"],
Plot[fit[x],{x,-1,3}]
]

Gewichtungen

Beginnen wir mit den gleichen Daten von oben:

errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,0.03*Exp[2*x]-3*x-3*2},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];
dataOhneAround=data/.{Around[v_,_]:>v};

Gewichtungen kann man über die Weights-Option mit einfügen. Gewichtet wird natürlich wieder instrumentell:

fit=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x,Weights->(1/errors^2)]

fit["ParameterTable"]

Wenn ihr eure daten als ein Around array gegeben habt, und ihr braucht die Unsicherheiten daraus zum gewichten,
dann kann man diese leicht mit Part rausholen:

data[[All,2,2]] (*[[All,2,2]] heißt wir nehmen von Allen Datenpunkten, den zweiten Eintrag (y) und davon den zweiten Eintrag (die Unsicherheit in Around)*)

Alternativ können natürlich wieder Regeln verwendet werden:

data/.{{x_,Around[y_,ye_]}:>ye}

Somit geht dann

fit=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x,Weights->(1/data[[All,2,2]]^2)]

IncludeConstantBasis

Leider ist LinearModelFit schon soooo alt, dass man damals bei der Einführung dieses Befehls eine schlechte Entscheidung getroffen hat,
die wir bis heute noch mitschleppen. Und diese ist IncludeConstantBasis.
Diese Option ist standardmäßig auf True und bedeutet, Mathematica wird immer einen konstanten Term mitfitten, ob man den jetzt angibt oder nicht.

Also wieder mit zufälligen Daten:

errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,2x^2+3*x},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];
dataOhneAround=data/.{Around[v_,_]:>v};

Also probieren wir doch mal die Funktion nur nach  zu fitten:

fit=LinearModelFit[dataOhneAround,{x^2,x},x] (*Wir geben explizit nur den fit nach (x^2) und x an*)
fit["ParameterTable"]

Ja huch?! da ist ja ein konstanter Term aufgetaucht.
Wir müssen also diese Option explizit ausschalten:

fit=LinearModelFit[dataOhneAround,{x^2,x},x,IncludeConstantBasis->False]
fit["ParameterTable"]

Wenn man übrigens selbst einen konstanten Term angibt, merkt Mathematica das und fügt keinen weiteren, eigenen hinzu.

fit=LinearModelFit[dataOhneAround,{x^2,x,2},x,IncludeConstantBasis->False]
fit["ParameterTable"]

VarianceEstimatorFunction

Standardmäßig skaliert Mathematica die Abschätzungen der Unsicherheiten mit der ausgerechneten Varianz der Anpassung.
Das ist in den meisten Fällen vollkommen okay und gerechtfertigt.
ABER
Wenn ihr euch vollkommen, absolut, defintiv, total sicher seid, dass die Unsicherheiten die ihr für die Gewichtungen abgegeben habt,
die einzigen sind die eure Streuungen beeinflusst haben, dann solltet ihr Mathematica das mitteilen. Denn dann sind die Varianzen
statistisch direkt durch eure Gewichte gegeben und müssen nicht von Mathematica abgeschätzt werden.
Das macht ihr indem ihr die Option VarianceEstimatorFunction auf Konstant 1 setzt.

Wir starten wieder mit dem Datenset von oben:

errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,0.03*Exp[2*x]-3*x-3*2},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];
dataOhneAround=data/.{Around[v_,_]:>v};

Nun fitten wir zweimal. Einmal mit Abschätzung und einmal mit der Option auf konstant eins:

fit=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x,Weights->(1/errors^2)]
fitVar1=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x,Weights->(1/errors^2),VarianceEstimatorFunction->(1&)]

Wir sehen ein dass die Fit Ergebnisse die gleichen sind. Die Unsicherheiten sind aber unterschiedlich groß:

fit["ParameterErrors"]
fitVar1["ParameterErrors"]

In diesem Fall sind die Unterschiede nicht groß. Das liegt daran dass ich den Datensatz so konzipiert habe.
In echt ist es aber wirklich selten der Fall das ihr alle Unsicherheiten so richtig eingeschätzt habt.
Eine “Probe” ist übrigens zu gucken wie groß Mathematica die Varianz abschätzt. Ist der Wert nahe an
eins, habt ihr einen guten Job gemacht und könnt überlegen VarianceEstimatorFunction->(1&) zu setzen.
Andernfalls lasst es. Wir schätzen also ab:

fit["EstimatedVariance"]

Passt.

Wir können mal ganz explizit die Abschätzung der Gewichte versauen, indem wir sie mit absicht 10x vergrößern:

fit=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x,Weights->(10*1/errors^2)]
fitVar1=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x,Weights->(10*1/errors^2),VarianceEstimatorFunction->(1&)]

Jetzt sind die Unsicherheiten im Varianz=1 gezwungenen Fall viel zu klein wärend der Fall in dem die Varianz abgeschätzt wird
immernoch das gleiche Ergebnis wie zuvor produziert:

fit["ParameterErrors"]
fitVar1["ParameterErrors"]

Das ist demnach auch in der Varianz zu sehen:

fit["EstimatedVariance"]

Cool. Die ist sogar 10x so groß wie zuvor. Also sehen wir dass die Gewichte (die wir ja mit x10 angegeben haben) auch um diesen Faktor zu groß waren.
Unterm Strich also: wenn ihr euch unsicher seid, fasst VarianceEstimatorFunction einfach nicht an. Easy.
(Die Motivation für die VarianceEstimatorFunction findet ihr unten im Kapitel Distributionstheorie->LogLikelihood)

OffsetFunction

LinearModelFit hat noch ein paar andere Optionen aber LinearOffsetFunction ist so die letzte die wirklich nice to know ist.
Ich habe glaube bis heute noch nie geschafft sie einmal in der realen Welt anzuwenden aber der Tag wird kommen.

Die Motivation ist die folgende: bisher haben wir gesagt ein linearer Fit ist eine Linearkombination von Parametern mit unabhängigen
Basisfunktionen. Dem kann man noch hinzufügen: “Plus eine Funktion unabhängig von den Parametern”. Hintegrund ist einfach, dass
man diese nehmen könnte, um sie direkt von den y Daten abzuziehen um dann normal linear ohne sie zu fitten und das würde keinen Unterschied machen.
Also nehmen wir mal ein Modell der Form
 
Wir sehen dort die Terme  zusätzlich welche unabhängig von den Parametern sind.

Erzeugen wir also Beispieldaten dieses Modells:

errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,0.03*Exp[2*x]-3*x-3*2+Log[x+1.1]+1},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];
dataOhneAround=data/.{Around[v_,_]:>v};

Nun fitten wir das mit  als LinearOffsetFunction:

fit=LinearModelFit[dataOhneAround,{Exp[2x],x,-3},x,LinearOffsetFunction->(Log[x+1.1]+1)]

Done. Nun noch nen hübscher Plot dazu und wir sind fertig mit LinearModelFit

Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

NonlinearModelFit

Generell

Wenn LinearModelFit ein Messer wäre, wäre NonlinearModelFit eine Kettensäge.
Wie der Name schon vermuten lässt, kann man mit NonlinearModelFit auch Funktionen fitten die man nicht als Linearkombination
von Parametern mit Basisfunktionen darstellen kann.
Daher hat diese Funktion vier Argumente:
1. Die Daten
2. Die Fitfunktion
3. Die Fitparameter
4. Die Funktionsvariable(n)

Weil es ironisch ist, erstellen wir mal ein zufälliges Beispieldatenset für eine lineare Funktion:

errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,2x-1},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];
dataOhneAround=data/.{Around[v_,_]:>v};

ListPlot[data]

Und fitten das jetzt mit NonlinearModelFit (und ja klar könnten wir hierfür LinearModelFit benutzen aber dann wäre das nichtmehr ironisch):

fit=NonlinearModelFit[dataOhneAround,A*x+B,{A,B},x]
Show[{Plot[fit[x],{x,-1,3}],ListPlot[data]}]

Und wir bekommen von diesem Objekt natürlich wieder alle Informationen die wir auch bei LinearModelFit bekommen haben.
Hier mal exemplarisch nur ParameterTable und BestFitParameters aufgeführt:

fit["ParameterTable"]
fit["BestFitParameters"]

Wir vergleichen trotzdem mal mit LinearModelFit und hoffen dass das Selbe rauskommt:

LinearModelFit[dataOhneAround,{x,1},x]["ParameterTable"]

Alles gleich. Das ist gut.

Start Parameter Wahl

Jetzt solltet ihr euch fragen wieso es überhaupt ein LinearModelFit gibt, wenn NonlinearModelFit doch auch alles lösen kann und eh das gleiche Ergebnis rauskommt?!
Und die Antwort auf diese Frage ist der Grund warum ich Co-Author auf einem Paper von der Universität Tübingen bin (kein Scherz, spendiert mir ein Bier dann erzähl ich euch die Geschichte).
LinearModelFit kann das Fitsystem deterministisch per Matrixinversionen lösen und kommt so IMMER zum globalen Optimum der Parameter. D.h. wenn ihr mit LinearModelFit arbeitet,
ist das Ergebnis das dort herauskommt immer das Beste das erreichbar ist. NonlinearModelFit on the other hand, kann kein deterministisches System aufstellen. Es arbeitet im Gegenzug
mit einem iterativen Optimierungsalgorithmus den ihr sogar steuern könnt, und manchmal sogar müsst. Dieser Algorithmus kann aber Gefahr laufen, in einem lokalen Minimum stecken
zu bleiben oder garnicht erst ein convexes Minimum zu finden. Das klingt nach einem utopisch seltenen Problem. Wird euch aber bereits im GP2 über den Weg laufen
(Fraunhofer Beugung und Besselfunktionen lassen grüßen - Dafür gibt es nachher im “Reale Beispiele” Kapitel auch genau dieses eine Beispiel).

So ein Beispiel ist auch schnell kreirt:

errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,Sin[5*x]},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];
dataOhneAround=data/.{Around[v_,_]:>v};

Ein einfacher Sinus . Kann ja wohl nicht allzuschwer sein?!

fit=NonlinearModelFit[dataOhneAround,Sin[A*x],{A},x]
Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

wtf? Was ist denn hier passiert? Das passt doch wohl garnicht an die Daten. Man hat eher das Gefühl, dass NonlinearModelFit versucht hat den groben Trend der Funktion nachzuarmen
als tatsächlich die Oszillationen mitzumachen. Wieso das?
Um das zu beantworten schauen wir uns mal an was die Fits tatsächlich optimieren. Das Abstandsquadrat zu allen Datenpunkten:

(*Den Code hier müsst ihr jetzt nicht verstehen*)
Show[Plot[Total[(dataOhneAround[[All,2]]-Sin[A*dataOhneAround[[All,1]]])^2],{A,-1,7},PlotRange->{{-1,7},{0,70}}],ListPlot[{Callout[{A,Total[(dataOhneAround[[All,2]]-Sin[A*dataOhneAround[[All,1]]])^2]}/.fit["BestFitParameters"],Style["NonlinearModelFit's Lösung",11]],Callout[{5,Total[(dataOhneAround[[All,2]]-Sin[5*dataOhneAround[[All,1]]])^2]},Style["Tatsächlich beste's Lösung",11]],Callout[{1,Total[(dataOhneAround[[All,2]]-Sin[dataOhneAround[[All,1]]])^2]},Style["Wo NonlinearModelFit angefangen hat zu optimieren",11]]},PlotStyle->Black],FrameLabel->Thread[Style[{"Fitparameter A","Summe Abstandsquadrat (kleiner=besser)"},16]],ImageSize->Large]

Ahhh jetzt wird das Bild ein bisschen klarer. NonlinearModelFit startet immer am Startpunkt, dass alle Variablen ungefähr 1 sind und optimiert daraufhin los.
Für ihn (um NonlinearModelFit mal zu vermenschlichen) sah es von diesen Standpunkt so aus, als wäre zu kleineren A-Werten zu gehen besser weil das
Abstandsquadrat in diese Richtung abnimmt. Woher sollte es wissen dass auf der anderen Seite des Berges ein Tal wartet das soviel besser ist?
Theoretisch sogar zwei Berge weiter. Es gibt noch ein weiteres Lokales Minimum bei etwa . 
Wir können NonlinearModelFit über mehrere Wege dazu bringen das mitzubekommen. An den Algorithmen spielen wir mal nicht rum weil das wäre eine
viel größere Box als wir jetzt gewillt sind aufzumachen. 
Was wir aber tun können: Wir fragen NonlinearModelFit nicht bei  zu beginnen sondern mal bei  oder so. Das kann man indem man den Parameter als Liste mit seinem Startwert angiebt. Mal sehen:

fit=NonlinearModelFit[dataOhneAround,Sin[A*x],{{A,4}},x]
Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

Ja nice. Und schon haben wir es sofort gefunden. Wir erinnern uns noch daran, dass bei  ein weiteres lokales Minimum war. Bleibt NonlinearModelFit darin stecken?

fit=NonlinearModelFit[dataOhneAround,Sin[A*x],{{A,2.5}},x]
Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

Ja tatsächlich tut es das. Bei .
Kann man irgendwie die Parameter den geschickt wählen um das zu umgehen? Naja das ist leider von Mal zu Mal unterschiedlich.
Ein guter Tipp: plottet einmal selbst das Fitmodel und versucht mit Hand gute Werte für den Fit zu finden. Die müssen nicht präzise sein.
Nur präzise genug um den Fit konvergieren zu lassen.
Ansonsten: Im “Reale Beispiele” Kapitel sind Beispiele gegeben zu diesem Problem.

Gewichtungen & VarianceEstimatorFunction

Es gilt hier das Selbe wie auch bei LinearModelFit. Also für Erklärungen verweise ich auf die vorhergehenden Kapitel bei LinearModelFit zu den beiden Themen.

(χ^2) und (R^2) Tests

Manchmal möchte man als ein Maß wie “gut” eine Regression geworden ist den (χ^2) und/oder (R^2) Test angeben.
Das können wir in Mathematica natürlich schnell bekommen.
Also machen wir mal schnell einen Fit:

errors=RandomVariate[NormalDistribution[0,0.4],41];
data=Table[{x,0.03*Exp[2*x]-3*x-3},{x,-1,3,0.1}];
data[[All,2]]=Around@@@Transpose[{data[[All,2]]+errors,Abs[errors]+0.2}];
dataOhneAround=data/.{Around[v_,_]:>v};
fit=LinearModelFit[dataOhneAround,{Exp[2x],x,1},x]

(R^2) Test

Der (R^2) Test ist sehr schnell bekommen:

fit["RSquared"]

Mehr ist dazu nicht zu sagen.

Adjusted (R^2) Test

Der (R^2) Test hat ein Problem. Und das können wir leicht nachvollziehen:

data=Transpose[{Range[-1,3,0.01],0.001*Range[-1,3,0.01]^2+2*Range[-1,3,0.01]+1+RandomVariate[NormalDistribution[0,0.3],Round[4/0.01+1]]}];
fits=Table[LinearModelFit[data,x^Range[0,endN],x]["RSquared"],{endN,1,30}]

ListLogPlot[fits,Joined->True,FrameLabel->Thread[Style[{"Anzahl an Parametern \!\(\*TemplateBox[<|"boxes" -> FormBox[
StyleBox["p", "TI"], TraditionalForm], "errors" -> {}, "input" -> "p", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)","(R^2) Test"},16]]]

Wir sehen also: jeh mehr Parameter wir in das Modell stecken, desto besser ist unser (R^2) Test. Das sollte natürlich nicht so sein.
Der angepasste (R^2) Test nimmt das unter Berücksichtigung:

afits=Table[LinearModelFit[data,x^Range[0,endN],x]["AdjustedRSquared"],{endN,1,30}]
ListLogPlot[{fits,afits},Joined->True,FrameLabel->Thread[Style[{"Anzahl an Parametern \!\(\*TemplateBox[<|"boxes" -> FormBox[
StyleBox["p", "TI"], TraditionalForm], "errors" -> {}, "input" -> "p", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)","Jeweiliger (R^2) Test"},16]],PlotLegends->{"(R^2) Test","Angepasster (R^2) Test"}]

Da wir in physikalischen Modellen aber eh selten den (R^2) Test benutzen (und der eh nicht allzuviel Unterschied macht) ist das das echt nicht so wichtig ehrlich gesagt.
Also make your own mind about that

fit["RSquared"]
fit["AdjustedRSquared"]

Test

Der (χ^2) Test wiederrum ist relativ nützlich. Der gibt wenigstens konkrete Vorschläge was los ist.
Dafür rufen wir von Mathematica mal etwas auf namens ANOVATable (Analysis of Variance Table):

fit["ANOVATable"]

Diese Tabelle enthält einige Aussagen über Abweichungen/Varianzen eures Modells.
Uns interessiert hierbei nur die “Error”-Zeile.
Da haben wir drei Einträge: DF (degrees of Freedom), SS (Sum of Squares), MS (Mean Sum of Squares).
Interessieren tuen wir uns nur für die MS Spalte. Denn genau das ist unser  Test. Also:

Den können wir sogar ohne Ablesen spezifischer bekommen:

Last[fit["ANOVATableMeanSquares"]]

Korrelation & Korrelierte Fehlerfortpflanzung

Korrelation

Ohne viel Hintergrundwissen gilt allgemein: zwei Größen sind miteinander Korreliert, wenn man die eine nicht ohne die
andere wählen kann. Man kann sich dieses Prinzip sehr leicht an einer Regression einer linearen Funktion  mit den korrelierten Parametern  erklären.

data=Table[{x,2x+1+RandomVariate[NormalDistribution[0,0.3]]},{x,-1,3,0.1}];
fit=LinearModelFit[data,{x,1},x];
Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

Wenn man jetzt zum Beispiel den Anstieg  festhält, und ein kleines bisschen nach unten zwingt, dann muss der
Achsenabschnitt  darauf reagieren und selbst steigen um die Daten noch halbwegs gut zu treffen:

fit=LinearModelFit[data,{1},x,LinearOffsetFunction->0.5x];
Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

Vice versa, wenn den Achsenabschnitt senkt, steigt der Anstieg:

fit=LinearModelFit[data,{x},x,IncludeConstantBasis->False,LinearOffsetFunction->-0.5];
Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

Und das gilt also auch in die andere Richtung. Also wir sehen: für jede Änderung die der Anstieg  macht,
macht der Achsenabschnitt  eine Änderung in die entgegengesetzte Richtung. Die beiden sind also (über die Daten)
miteinander korrliert.

Man führt hierfür einen Korrelationsfaktor  zwischen den Variablen  und  ein. Diese kann Werte annehmen von . 
Dafür gibt jetzt drei Fälle:
1.  - unkorreliert. Die Variablen können sich unabhängig voneinander ändern.
2.  - korreliert. Die Variablen gehen bei Änderungen immer gemeinsam in eine Richtung.
3.  - antikorreliert. Die Variablen gehen bei Änderungen immer in entgegengesetzte Richtungen (siehe Beispiel: lineare Funktion)

Damit können wir aber noch nicht rechnen. Durch Herleitung aus Multinormal-Verteilungen (siehe Distributionstheorie) kann man eine Größe
Gewinnen, die sich Kovarianz  nennt. Definieren kann man diese als

mit  der Korrelationsfaktor, und  die Unsicherheiten von  und .
Damit lässt sich jetzt rechnen. Für jedes Paar an Korrelationen fängt man sich einen weiteren Term  unter der Wurzel der Fehlerfortpflanzung ein:


Also rechnen wir mal ein Beispiel.
Nehmen wir an, wir möchten  berechnen und der Korrelationsfaktor ist gegeben mit .
So gilt also:

Sqrt[
(D[A/B,A]*uA)^2 +(*unkorrelierter Term für A*)
(D[A/B,B]*uB)^2+(*unkorrelierter Term für B*)
2*D[A/B,A]*D[A/B,B]*uAB(*Korrelationsterm*)
]/.{
uAB->uA*uB*(-1) (*Die Kovarianz kann angegeben werden wenn die Korrelation bekannt ist.  in diesem Fall*)
}

Und dass könnte man jetzt Problemlos ausrechnen.

Bleibt nur die Frage, wenn man nicht gerade eine Ungleichung findet wie z.B.  die quasi  impliziert, wie berechnet man dann
die Kovarianz. Nun, da gibt es verschiedenste Möglichkeiten. Aber eine ist eine Regression.
Also machen wir mal eine von einer linearen Funktion:

data=Table[{x,2x+1+RandomVariate[NormalDistribution[0,0.1]]},{x,-1,3,0.1}];
fit=LinearModelFit[data,{x,1},x]

Von dem FittedModel-Objekt können wir tatsächlich explizit und direkt nach den Korrelationsmatrizen / Kovarianzmatrizen fragen:

fit["CorrelationMatrix"]//MatrixForm
fit["CovarianceMatrix"]//MatrixForm

Die Matrizen sind folgendermaßen aufgebaut:
Korrelationmatrix: (({{1, , , ⋯}, {, 1, , ⋯}, {, , 1, ⋱}, {⋮, ⋮, ⋱, ⋱}}))
Kovarianzmatrix:    (({{, , , ⋯}, {, , , ⋯}, {, , , ⋱}, {⋮, ⋮, ⋱, ⋱}}))

Man kann sich aus der Kovarianzmatrix also alles ablesen was man zum korrelierten Fehlerfortpflanzen brauch.
Aber wir sehen gleich dass es einfacher geht.
Ein paar Nebenbemerkungen die der Leser gern selbst beweisen kann (geht alles über die Definition der Kovarianz/Korrelation):
1. Sowohl Korrelationsmatrix als auch Kovarianzmatrix sind immer symmetrisch und somit  und .
2. Die Selbstkovarianz  ist identisch zum Unsicherheitsquadrat 
3. Die Selbstkorrelation  ist immer identisch eins

Herleitung: Korrelierte Fehlerfortpflanzung

Dieses Kapitel hier ist nur für die geeignet, die tiefer in die Materie einsteigen wollen.
Es ist definitv nicht notwendig um die Praktika zu absolvieren.

Es gibt einige Möglichkeiten die erste Ordnung der Gaußschen Fehlerfortpflanzung herzuleiten.
Sobald wir aber Korrelationen betrachten wollen, hören die meisten anschaulichen Herleitungen auf.
Deswegen schreibe ich es hier mal auf für alle die es interessiert.

Wir fangen an mit der Aufgabenstellung. Wir wollen berechnen:

Dabei sind ,   und  .

Wir wiederholen nochmal die Definition für die Standardabweichung (eigentlich Varianz in der Form - fpr Standardabweichung einmal Wurzel ziehen):

 

Dann, wir haben es zwar noch nicht gesagt, führen wir auch die Definition für die Kovarianz ein:

(Der/Dem findingen Leser*in sollte auffallen, dass die Standardabweichung im Quadrat nichts anders ist als die Kovarianz mit sich selbst, also . )

Hier jetzt der Schlachtplan zur Herleitung:
 - Die Unsicherheit  ist uns im Prinzip bekannt. Zumindest als Formel. Was wir nicht kennen ist der  Term.
 - Dieser Term ist im Prinzip eine Differenz des  Wertes: 
 - Im linearen erster Ordnung, können wir erstmal annehmen dass wir  klein machen. Also  .
 - Wir sollten in der Lage sein  schreiben zu können.
 - Mit allen obrigen Definitionen sollten wir das ganze ausschreiben und vereinfachen können.
 Lets try:
 
 Wir beginnen bei der Ausgangssituation  und wenden die Kettenregel an:
 
 Wir erinnern uns daran, dass wir sagen können:
 
 
 In gleicher Richtung wie wir  geschrieben haben, können wir jetzt rückwärts schreiben:  und analog für y.
 
 
Das ganze setzen wir jetzt in die Definition der Standardabweichung für  ein:

Nun ein bisschen 8. Klasse Algebra:

und ein bisschen Verteilen der Summe:


Wir sehen diese  Terme und können diese aus den Definitionen oben als  identifizieren:


Genauso sehen wir aus den Definitionen oben, dass der  Term die Kovarianz ist in der Form   ist.


Schlussendlich ziehen wir die Wurzel und sind fertig:


□

Vermeiden von Korrelierten Rechnungen bei linearen Fits

Bevor wir anfangen überhaupt mit Korrelationen viel zu rechnen, macht es Sinn mal zu gucken wie wir sowas eher vermeiden können.

Explizit schauen wir mal über das fitten einer linearen Funktion .
Es kommt erstaunlich häufig vor, dass wir  oder  ausrechnen wollen. Wir haben ja bereits festgestellt dass die beiden Parameter
einer linearen Funtion stark antikorreliert miteinandern sind. Also hier nicht mit korrelierten Fehlern zu arbeiten ist nicht so toll.
Starten wir also mal ein Beispiel und fitten eine lineare Funktion an ein Beispieldatenset:

data=Table[{x,2x+1+RandomVariate[NormalDistribution[0,0.1]]},{x,-1,3,0.1}];
fit=LinearModelFit[data,{x,1},x]
{a,b}=Around@@@Transpose[{fit["BestFitParameters"],fit["ParameterErrors"]}]

Unkorreliert bekommen wir also für   und  respektiv:

a/b
b/a

Wir wissen dass das “falsch” ist aber machen wir mal weiter.
Wenn wir uns die lineare Funktion  so ansehen, könnten wir mal auf die Idee kommen diese so Umzuformen,
dass die jeweiligen Brüche dastehen würden. Also in etwa so:

und

wenn wir jetzt die Brüche mit Platzhalter-Variablen ersetzen, bekommen wir zwei neue Fitmodelle
 mit 
und
 mit 
Diese Modelle sind zwar jetzt nichtlinear aber irgendeinen Preis müssen wir wohl schon bezahlen.
Wenn die Regressionen jetzt also gut funktionieren bekommen wir als Ergebniss  oder  zurück.
Die Unsicherheit der beiden Fitparameter entspricht genau der korreliert fortgepflanzten Brüche.
Also machen wir:

fitBdurchA=NonlinearModelFit[data,A*(x+U),{A,U},x];
fitBdurchA["ParameterTable"]
fitAdurchB=NonlinearModelFit[data,B*(V*x+1),{B,V},x];
fitAdurchB["ParameterTable"]

Wir bekommen also die voll korreliert ausgerechneten Brüche durch nur Ablesen von  und :

und

Die unkorrelierten Ergebnisse oben hatten also die Unsicherheit der Brüche unterschätzt.

VectorAround[]

Wenn Mathematica mit unkorrelierten Variablen über Around[] rechnen kann, dann wird es wohl auch ein Objekt für
korrelierte Größen geben. Das ist VectorAround[]:

Information[VectorAround]

Benutzen werden wir nur die Definitionen zwei und drei (eigentlich auch wirklich nur zwei).
Aber nehmen wir mal kurz drei. Wir können zwei Größen miteinander korrelieren wenn wir den Korrelationsfaktor kennen.
Also setzen wir mal:

a=5.156; (*Wert für a*)
ua=0.277; (*Unsicherheit für a*)
b=0.248; (*Wert für b*)
ub=0.008; (*Unsicherheit für b*)
Cab=-0.5; (*Korrelationsfaktor zwischen a & b*)

Nach der dritten Definition können wir also schreiben:

value=VectorAround[{a,b},{{ua,ub},Cab}]

Witzigerweise, hat VectorAround bei der Eingabe, das ganze direkt zu einer Kovarianzmatrix umgewandelt. Also das was wir in Definiton zwei eh eingegeben hätten.
Nevermind. Mit dieser Größe könnt ihr zwar auch normal rechnen:

2*value+1

aber das passiert unabhängig miteinander. Und ehrlich gesagt sollten wir das nicht tun.
Stattdessen empfehle ich VectorAround nur mit AroundReplace benutzen.

Schlussendlich überzeugen wir uns davon, dass symbolisch das gleiche rauskommt wie wir bereits in den Kapiteln zuvor geschrieben haben:

Simplify@AroundReplace[F[A,B],{{A,B}->VectorAround[{A,B},{{UA^2,UAB},{UAB,UB^2}}]}]//TraditionalForm

Wir sehen klar unsere beiden unkorrelierten Terme, sowie den einen Korrelationsterm wie wir ihn bereits geschrieben hatten.

Korreliertes AroundReplace[]

Dann benutzen wir mal AroundReplace mit Werten:

a=5.156; (*Wert für a*)
ua=0.277; (*Unsicherheit für a*)
b=0.248; (*Wert für b*)
ub=0.008; (*Unsicherheit für b*)
Cab=-0.5; (*Korrelationsfaktor zwishcen a & b*)
value=VectorAround[{a,b},{{ua,ub},Cab}];

AroundReplace[A/B,{{A,B}->value}]

Bumm. Voll fertig korreliert Fehlerfortgepflanzt. Easy peasy.

Jetzt machen wir das nochmal mit einer Regression:

data=Table[{x,2x+1+RandomVariate[NormalDistribution[0,0.1]]},{x,-1,3,0.1}];
fit=LinearModelFit[data,{x,1},x]

Wir denken daran, die zweite Definition zu nutzen, und einfach nach der Kovarianzmatrix zu fragen:

value=VectorAround[fit["BestFitParameters"],fit["CovarianceMatrix"]]

Cool. Also ab in AroundReplace[] damit:

AroundReplace[A/B,{{A,B}->value}]

Btw. ich will wetten, spätestens in Mathematica version 13 oder so, hat FittedModel-Objekt bestimmt
direkt Properties um direkt VectorAround’s zu bekommen. Wenn nicht wäre ich enttäuscht. Mal sehen.

- Nachtrag: 17.01.2022 - Mathematica 13.0 hat es nicht...bin enttäuscht..

Ein Schlusswort noch. Machen wir nochmal die Regression:

data=Table[{x,2x+1+RandomVariate[NormalDistribution[0,0.1]]},{x,-1,3,0.1}];
fit=LinearModelFit[data,{x,1},x]

Wir denken daran, die zweite Definition zu nutzen, und einfach nach der Kovarianzmatrix zu fragen:

value=VectorAround[fit["BestFitParameters"],fit["CovarianceMatrix"]]

Uhi, was denn hier passiert? Was nerviges, actually! VectorAround erwartet eine perfekt symmetrische Matrix.
LinearModelFit macht, wenn man nach der Kovarianzmatrix fragt, diese Matrix durch eine Matrixinversion. Und die
ist leider nicht numerisch absolut auf alle 16 Nachkommastellen perfekt symmetrisch. Nur fast.
Wir können die aber super leicht symmetrisieren indem wir die kleinen Abweichungen weg mitteln.
Das geht einfach durch

covMat=fit["CovarianceMatrix"];
value=VectorAround[fit["BestFitParameters"],1/2*(covMat+Transpose[covMat])]

Auch hier hoffe ich, dass Mathematica das mal nachbessert.

Plotten

Plot

Information[Plot]

Plot[] ist wohl der wichtigste Plot-Befehl, es nimmt eine Funktion oder eine Liste an Funktionen entgegen, und plottet diese über einen spezifischen bereich.

Plotten fir mal  für :

Plot[x^2,{x,-1,3}]

Oder Jetzt

Plot[{1,x,x^2,x^3,x^4},{x,-1,2}]

Der curious case of HoldAll

Plot hat sogenannte Attributes (wie so ziemlich alle eingebaute Funktionen). Diese verändern das Verhalten dieser Funktion. Plot hat folgende:

Attributes[Plot]

Wärend Protected und ReadProtected nur Attribute sind um das überschreiben von Plot zu verhindern, ist HoldAll ein sehr interessates und am Anfang sehr verwirrendes Attribut.
Aber machen wir erstmal ein Beispiel indem dieses Attribut seine Wirkung zeigt.
Wir plotten wieder einmal, die Funktionen . Nur dieses mal, einmal durch manuelle Eingabe, und einmal als Table:

Plot[{1,x,x^2,x^3,x^4},{x,-1,2}]
Plot[Table[x^n,{n,0,4}],{x,-1,2}]

Mhh. Das ist ja komisch. Die Linien sind zwar die gleichen, aber die Farben nicht. Der Verursacher ist natürlich dieses HoldAll-Attribut.
Es verursacht nämlich, das Argumente der Funktion Plot[] nicht vorher ausgewertet werden, sondern verbatim in die Funktion eingehen.
D.h. unsere Table[] wird nicht ausgeführt sondern geht als EINE Funktion in den Plot ein.
Verwirrt? Machen wir mal ein leichtes Beispiel und definieren unsere eigene Funktion:

FunktionOhneHoldAll[___];

Attributes[FunktionMitHoldAll]={HoldAll};
FunktionMitHoldAll[___];

Wenn wir nun die Funktionen mit der Formel  aufrufen, dann schaut mal was passiert:

FunktionOhneHoldAll[5+5]
FunktionMitHoldAll[5+5]

AHA! Die Funktion die wir mit einem HoldAll-Attribut garniert hatten, weigert sich das  auszuwerten.
Und das passiert mit allen Argumenten:

FunktionOhneHoldAll[5+5,Table[x^n,{n,0,4}]]
FunktionMitHoldAll[5+5,Table[x^n,{n,0,4}]]

Und Plot hat dieses Attribut auch. D.h. als wir Table[x^n,{n,0,4}] als Argument gegeben haben, wurde es NICHT vorher
zu der Liste {1,x,x^2,x^3,x^4} gemacht sondern blieb in der Table Form. Dadurch hat Plot das als EINE Funktion gesehen
und demnach alles in einer Farbe eingefärt.
Oki doki. Aber warum eigentlich? Wieso sollte Plot dieses Attribut brauchen?
Kurze Antwort: es beschützt die Plot Variablen. Man nehme folgendes Beispiel:

a=5
Plot[a^2,{a,-1,3}]

Dadurch das die Definition von  NICHT eingesetzt wird, kann Mathematica das überhaupt plotten. Hätte Mathematica
das Attribut nicht, dann hätten wir der Funktion am Ende Plot[5^2,{5,-1,3}] übergeben. Und damit kann Mathematica nichts anfangen.

Okay letzte Frage, was ist aber wenn ich er vorher ausführen WILL. Z.b. unser Table example um x^0 bis x^4 zu plotten?
Da hilft uns der Evaluate[] Befehl. Dieser zwingt eine Auswertung des Ausdrucks selbst wenn ein HoldAll-Attribut ansteht.
Also finally:

Plot[Evaluate[Table[x^n,{n,0,4}]],{x,-1,2}]

(Letzte Nebenbemerkung: LinearModelFit und NonlinearModelFit haben leider kein HoldAll Argument. Und das ist blöd weil ich jedes Semester
Emails von Studenten bekommen wieso ihr Fit nicht funktioniert und am Ende haben sie einfach nur einen der Parameter oder Fitvariablen
bereits definiert und 5 ist nunmal kein Parameter. Don’t get me wrong, send me emails wenn etwas nicht läuft. Aber der Fehler macht in
irgendeiner Form etwa 50% der emails aus. :D )

ListPlot

Information[ListPlot]

ListPlot ist wohl der zweitwichtigste Befehl. Den er erlaubt uns daten zu plotten.

ListPlot[{1,2,3,Pi,0}]

Wenn wir auch x-Daten angeben wollen, muss das in vielen kleinen Punktelisten der Form {{(Subscript[x, 1]), (Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])}, ...} kommen:

ListPlot[{{1,1},{2,4},{3,Pi},{10,5}}]

Arounds bekommen by default die Fehlerbalken garniert:

ListPlot[{{1,1},{2,Around[4,1]},{Around[5,2],2},{Around[6,1],Around[6,2]}}]

ListLinePlot

ListLinePlot ist quasi ListPlot, bloß mit Linien:

ListLinePlot[{{1,1},{3,Pi},{4,0},{6,10}}]

ListLinePlot ist sogar eigentlich nur Zucker für ListPlot[...,Joined->True]

ListPlot[{{1,1},{2,4},{3,Pi},{10,5}},Joined->True]

ColorData

Teil davon einen hübschen Plot zu machen, ist es die richtige Farbe zu finden.
Mathematica kennt einige Farben sogar bei Namen und PlotStyle ist da um diese einzusetzen:

Plot[Evaluate[Table[x+n,{n,0,5}]],{x,0,1},PlotStyle->{Black,Red,Green,Blue,LightGray,Magenta}]

Die sehen jetzt aber auch nicht so hübsch aus. Stattdessen hat Mathematica eine Sammlung an verschiedenen Farbpaletten.
Diese sind über einen Index abfragbar.

indizes=ColorData["Indexed"]
randomIndizes=RandomChoice[indizes,5]; (*fünf zufällige auswählen zum anzeigen*)
Grid[{#,ColorData[#]}&/@randomIndizes,Frame->All]

Ich persönlich bin ein großer Fan von ColorData[97]. Aber was ihr benutzt bleibt euch überlassen:

Plot[Evaluate[Table[x+n,{n,0,5}]],{x,0,1},PlotStyle->{ColorData[97,1],ColorData[97,2],ColorData[97,3],ColorData[97,4],ColorData[97,5],ColorData[97,6]}]

PlotThemes

Plots hübsch und presentabel zu machen gehört auch zur Show. Hierbei helfen euch PlotThemes. Diese stellen bei Default viel in euren Plots automatisch an.
Hier ein paar Beispiele:

Partition[Plot[{x,x^2,x^3},{x,-1,1.5},PlotLabel->#,PlotTheme->#]&/@{"Web","Minimal","Detailed","Business","Marketing","Scientific","Monochrome","Classic","Default"},3]//Grid[#,Frame->All]&

Ich brauch nicht erwähnen, dass für euch nur das “Scientific” oder “Detailed” Theme in Frage kommt.
Allerdings, könnt ihr einen normalen Plot auch so aussehen lassen wie jeden anderen. Alles was man ändern kann
ist über Optionen änderbar.

Als Beispiel, wandel ich mal das Monochrome-Theme hier so um, dass es so aussieht als wäre es das Detailed Theme.

Plot[{x,x^2,x^3},{x,-1,1.5},PlotTheme->"Monochrome"]
Plot[{x,x^2,x^3},{x,-1,1.5},PlotTheme->"Monochrome",Frame->True,FrameTicks->Automatic,PlotLegends->"Expressions",PlotStyle->(Directive[ColorData[97,#],Dashing[None]]&/@{1,2,3}),GridLines->Automatic,Axes->False]

Plots miteinander verbinden: Show[]

Am Ende muss man mehrere Plots miteinander verbinden um z.B. eine Regession in den Daten darzstellen. Dazu gibt es Show[]

Information[Show]

Also verbinden wir mal einen ListPlot mit einem Plot

Show[ListPlot[Table[{x,x^2+RandomReal[0.2{-1,1}]},{x,-1,3,0.1}]],Plot[x^2,{x,-1,3}]]

Dabei werden immer die Optionen des ersten Plots übernommen. Hier mal anhand des PlotLabels gezeigt:

Show[ListPlot[{{-1,0},{2,5}},PlotLabel->"Plot1"],Plot[x^2,{x,-1,3},PlotLabel->"Plot2"]]

Ein letztes Wort der Warnung noch. Das Show heißt nicht das automatisch die Optionen der einzelnen Plots angepasst werden.
Z.b.:

Show[ListPlot[Table[{x,Exp[x]+RandomReal[{-1,1}]},{x,-1,6,0.1}]],Plot[Exp[x],{x,-1,6}]]

Wir sehen, dass die Linie von Plot nicht ganz bis  durchgeht. Wieso? Nun Plot hat entschieden, dass es mehr von dem Funktionsverlauf zeigen kann wenn es den Plot Bereich ändert.

Plot[Exp[x],{x,-1,6}]

ListPlot hat diese Entscheidung aber nicht getroffen. So sehen wir also ListPlot mit dem vollen Plotbereich wärend Plot begrenzt ist. Wir müssen Plot also zwingen den gesamten Bereich zu zeigen.

Show[ListPlot[Table[{x,Exp[x]+RandomReal[{-1,1}]},{x,-1,6,0.1}]],Plot[Exp[x],{x,-1,6},PlotRange->All]]

DensityPlot

Wenn ihr eine Funktion  habt, dann ist man häufig besser gelegen einen DensityPlot zu benutzen statt eines 3D Plots.
Dieser kodiert den -Wert als Farbe statt wirklicher Höhe.

DensityPlot[Sin[x]*Cos[3y],{x,-4,4},{y,-3,5},PlotLegends->Automatic]

Diese Plots sind für gewöhnlich auf Schnelligkeit ausgelegt statt auf Qualität. Wenn man das ins Protokoll/Paper/BA/MA machen will,
muss man die Qualität erhöhen. Das kann man machen indem man mehr Punkte ausrechnen lässt. Dafür gibt es:
- PlotPoints (Die Anzahl pro Dimension auszurechnen)
- MaxRecursion (Die Anzahl an Verfeinerungen an stark Ändernden Orten)
- PerformanceGoal (Was das Ziel beim plotten ist: Schnelligkeit oder Qualität)

DensityPlot[Sin[x]*Cos[3y],{x,-4,4},{y,-3,5},PerformanceGoal->"Quality",PlotPoints->50,MaxRecursion->2]

Den Unterschied sieht man wenn man sich die Meshes ansieht die kreiert werden:

{DensityPlot[Sin[x]*Cos[3y],{x,-4,4},{y,-3,5},ImageSize->Medium,Mesh->All],DensityPlot[Sin[x]*Cos[3y],{x,-4,4},{y,-3,5},ImageSize->Medium,PerformanceGoal->"Quality",PlotPoints->50,MaxRecursion->2,Mesh->All]}

Das rechte ist viel feiner, demnach ist der Plot auch sehr fein.

Jetzt ist natürlich noch die Farb-Palette zu nennen. Jeh mehr Farben in einer Palette, desto genauer sieht man Details. Das kann aber auch zu Falschdarstellungen kommen und ehrlich gesagt:
es gibt viele Wissenschaftler die hier sehr religiös sind und einige Farbpaletten verabscheuen usw. das ist teilweise etwas ridiculous.
Eine verschmähte Palette sind z.B. Rainbow-Tables:

DensityPlot[Exp[-x^2+y/3],{x,-4,4},{y,-3,5},PlotPoints->100,MaxRecursion->0,PlotRange->All,ColorFunction->"Rainbow"]

Am Ende ist das wieder eure Wahl. Und Mathematica bietet natürlich wieder viele Gradienten-basierte ColorFunctions an:

ColorData["Gradients"]

Ich persönlich mag ja die Temperatur-basierten Schemes und die die Parula ähnlich sehen:

DensityPlot[Exp[-x^2+y/3],{x,-4,4},{y,-3,5},PlotPoints->100,MaxRecursion->0,PlotRange->All,ColorFunction->"LightTemperatureMap"]
DensityPlot[Exp[-x^2+y/3],{x,-4,4},{y,-3,5},PlotPoints->100,MaxRecursion->0,PlotRange->All,ColorFunction->"BlueGreenYellow"]

Es gibt sogar Custom Color Functions vom Wolfram Function Repository:

TurboColor=ResourceFunction["TurboColor"]
ViridisColor=ResourceFunction["ViridisColor"]
DensityPlot[Exp[-x^2+y/3],{x,-4,4},{y,-3,5},PlotPoints->100,MaxRecursion->0,PlotRange->All,ColorFunction->(TurboColor[#]&)]
DensityPlot[Exp[-x^2+y/3],{x,-4,4},{y,-3,5},PlotPoints->100,MaxRecursion->0,PlotRange->All,ColorFunction->(ViridisColor[#]&)]

Plot3D

Sollte man dennoch einen 3d Plot machen wollen, gibt es Plot3D

Plot3D[Sin[x]*Cos[3y],{x,-4,4},{y,-3,5}]

Und auch hier gelten die gleichen Qualitätsunterschiede:

Plot3D[Sin[x]*Cos[3y],{x,-4,4},{y,-3,5},PlotPoints->50,MaxRecursion->2]

Es gibt aber selten einen Grund ein 3D Plot einem DensityPlot vorzuziehen.

LogPlots

Manchmal muss man die Achsen eines Plots auch skalieren.
Für die logarithmischen Fälle gibt es extra Funktionen

LogPlot[x^2,{x,1,10}]
LogLinearPlot[x^2,{x,1,10}]
LogLogPlot[x^2,{x,1,10}]

Bzw. auch die ListPlot Analoga:

ListLogPlot[Table[{x,x^2},{x,1,10,0.5}]]
ListLogLinearPlot[Table[{x,x^2},{x,1,10,0.5}]]
ListLogLogPlot[Table[{x,x^2},{x,1,10,0.5}]]

Allerdings ist per Option jede Skalierung möglich:

Plot[x^2,{x,1,5},ScalingFunctions->{"Reciprocal","Reverse"}]

Legenden

Angehängt

Legenden werden immer häufig dann gebraucht, wenn man multiple Plots in einem hat und bezeichnen muss,
welche Kurve welche ist.
Es gibt die einfache PlotLegends Option welche mit dem Automatic Argument für Plot sogar ziemlich einfach ist:

Plot[{x,x^2,x^3},{x,0,1},PlotLegends->Automatic]

In dieses 1, 2, 3 könnt ihr einfach draufklicken und was eingeben.
Man Kann aber auch direkt eine Liste an Strings übergeben:

Plot[{x,x^2,x^3},{x,0,1},PlotLegends->{"Eine Legende","💩 Emoji","(x^3) Funktion"}]

Oder wenn man super faul ist, kann Mathematica in Plot auch selbst die Ausdrücke hinzufügen:

Plot[{x,x^2,x^3},{x,0,1},PlotLegends->"Expressions"]

Diese Legende wie sie Plot[] macht, nennt sich übrigens LineLegend/PointLegend und diese kann man selbst machen:

llegend=LineLegend[{Black,Red,Blue},{"1.","Zweitens","III"}]

PointLegend[{Black,Red,Blue},{"1.","Zweitens","III"}]

Diese kann man dann auch an einen beliebigen Plot mit PlotLegend anhängen:

Plot[x,{x,0,1},PlotLegends->llegend]

Passt zwar farblich und von der Anzahl her nicht aber es geht.
Wir können die Legende natürlich auch verschieben. Das geht mit Placed:

Grid[{{
Plot[x,{x,0,1},PlotLegends->Placed[llegend,Right]],
Plot[x,{x,0,1},PlotLegends->Placed[llegend,Left]],
Plot[x,{x,0,1},PlotLegends->Placed[llegend,Top]],
Plot[x,{x,0,1},PlotLegends->Placed[llegend,Bottom]]
}},Frame->All]

Finally noch um es nochmal zu erwähnen: das PlotTheme->”Scientific” ist seit v10 bis v13 von Mathematica die Nummer 108:

colorScheme=ColorData[108]
colors=colorScheme/@Range[1,10]
Plot[Evaluate@Table[x^n,{n,1,10}],{x,0,1},PlotLegends->LineLegend[colors,Table["x hoch "<>ToString[n],{n,1,10}]]]

Das ist ja jetzt alles schön und gut. Aber was ist wenn wir (aus irgendeinen Grund) keine
PlotLegends Option für unser Objekt haben?
Das passiert zum Beispiel bei Show. Oder wenn wir mehrere Plots miteinandern Mischen und nicht alle eine LineLegend sind:

Show[Plot[x,{x,0,1}],ListPlot[{#,#+RandomVariate[NormalDistribution[0,0.01]]}&/@Range[0,1,1/100]],PlotLegends->{"Hallo","Das","Geht","Nicht"}]

Wie wir sehen, entsteht keine Legende trotz angegebener Liste. Woher soll Mathematica auch die Farbinformationen holen?
Naja, ihr könntet jetzt sagen “Naja kein Ding, definiere ich halt eine eigene LineLegend”. Und klar, versuchen wir es mal:

Show[Plot[x,{x,0,1}],ListPlot[{#,#+RandomVariate[NormalDistribution[0,0.01]]}&/@Range[0,1,1/100]],PlotLegends->LineLegend[{Black,Red,Blue},{"Voll","Selbst","Definiert"}]]

Pustekuchen. Aber wir erinnern uns mal was wir bei Show gesagt hatten. Show übernimmt die PlotOptionen des ersten Plots in der Liste.
Ahh, also einfach umquartieren?

Show[Plot[x,{x,0,1},PlotLegends->LineLegend[{Black,Red,Blue},{"Voll","Selbst","Definiert"}]],ListPlot[{#,#+RandomVariate[NormalDistribution[0,0.01]]}&/@Range[0,1,1/100]]]

Jup. Aber was ist, wenn das nicht geht (aus irgendwelchen Gründen), oder wir mischen müssen, oder oder oder?
Naja, dann benutzen wir einfach Grid:

Grid[{
{Null,LineLegend[{Magenta,Purple},{"Oben","Top"},LegendLayout->"Row"],Null},
{LineLegend[{Red,Yellow},{"Eins","2"}],Plot[x,{x,0,1},ImageSize->Medium],LineLegend[{Black,Blue},{"1","Zwei"}]},
{Null,LineLegend[{Orange,Cyan},{"Unten","Bottom"},LegendLayout->"Row"],Null}
}]

Damit sind uns natürlich keine Grenzen gesetzt. Man beachte, dass ich die Grid zu einer 3x3 Matrix gemacht habe in der Form:
{
{Null,            Legende,    Null},
{Legende,  Plot,              Null},
{Null,            Legende,    Null}
}
Die Null sind einfach leere Einträge.

Was ist denn aber nun, wenn wir LineLegends und Punkte brauchen?
Also bereiten wir mal den Plot und beide Legenden vor:

plot=Show[Plot[{x,2x},{x,0,1},ImageSize->Medium,PlotStyle->{Black,Red}],ListPlot[Table[{x,a*x+RandomVariate[NormalDistribution[0,0.05]]},{a,1,2},{x,0,1,1/100}],PlotStyle->{Black,Red}]]
legend1=LineLegend[{Red,Black},{"Rote Linie","Schwarze Linie"}]
legend2=PointLegend[{Red,Black},{"Rote Punkte","Schwarze Punkte"}]

Naja, prinzipiell können wir wieder easy eine Grid benutzen. Bloß doppelt:

grid(({{plot, grid(({{legende 1}, {legende 2}}))}}))

Grid[{{plot,Grid[{{legend1},{legend2}}]}}]

Easy. Aber seihen wir mal ehrlich: der Text wäre Aligned schon cooler.
Wir könnten das auch mit Grid machen in der form
grid(({{plot, grid(({{rote linie, text1}, {schwarze linie, text2}, {roter punkt, text3}, {schwarzer punkt, text4}}))}}))
Wie machen wir also diese Punkte und Linien. Mit dem selben Befehl, den Plot, ListPlot, ..., etc. auch verwenden:

Head[plot]

Graphics! (Und Point[], Line[], Rectangle[], Disk[], ..., etc.)

Graphics[{
Red,Point[{0,0}],
Black,Line[{{0,1},{1,1}}],
Blue,Point[{{0.5,0.5},{-0.5,0}}],
Line[{{-0.25,-0.5},{0,0.2},{0.8,0.8}}],
Yellow,Rectangle[{-0.8,0.2},{-0.2,0.5}],
Purple,Disk[{0.7,-0.7},0.2]
}]

Mit ein paar Optionen die wir schon kennen, ist es super leicht die primitiven Punkte/Linien für die Legende zu erzeugen:

roteLinie=Graphics[{Red,Thick,Line[{{0.15,0},{0.85,0}}]},ImageSize->30,AspectRatio->1/2,PlotRange->{{0,1},{-1,1}}]
schwarzeLinie=Graphics[{Black,Thick,Line[{{0.15,0},{0.85,0}}]},ImageSize->30,AspectRatio->1/2,PlotRange->{{0,1},{-1,1}}]
roterPunkt=Graphics[{Red,PointSize[0.2],Point[{{0.5,0.5}}]},ImageSize->30,AspectRatio->1/2,PlotRange->{{0,1},{0,1}}] (*Mit Point[{{x1,y1}}] und PointSize*)
schwarzerPunkt=Graphics[{Black,Disk[{0.5,0.5},{0.1,0.2}]},ImageSize->30,AspectRatio->1/2,PlotRange->{{0,1},{0,1}}] (*Oder Mit Disk[{x,y},radius]*)

Und damit bauen wir dann easy die Legende zusammen:

Grid[{{plot,
Grid[{
{roteLinie,"Rote Linie"},
{schwarzeLinie,"Schwarze Linie"},
{roterPunkt,"Roter Punkt"},
{schwarzerPunkt,"Schwarzer Punkt"}
}]
}}]

Also schonmal aligned. Das ist cool. Aber jetzt noch linksbündiger Text:

Grid[{{plot,
Grid[{
{roteLinie,"Rote Linie"},
{schwarzeLinie,"Schwarze Linie"},
{roterPunkt,"Roter Punkt"},
{schwarzerPunkt,"Schwarzer Punkt"}
},Alignment->Left]
}}]

Nice.

Inline

Legenden außerhalb des Plots zu platzieren ist ja easy mit PlotLegens oder Grid.
Aber was wenn wir die Legende innerhalb des Plots haben wollen?

Dafür benutzen wir die Epilog oder Prolog Option in Graphics[], Plot[], ..., etc.

Plot[x,{x,0,1},Prolog->{Disk[{0.3,0.5},0.2]}]
Plot[x,{x,0,1},Epilog->{Disk[{0.3,0.5},0.2]}]

Cool. Wir sehen also, Prolog und Epilog nehmen einfach Graphics Objekte entgegen und packen sie in den Plot/Graphics.
Dabei unterscheidet sich Prolog zu Epilog, dass der Erstere es hinter den Plot legt, während Epilog es vor den Plot legt.

(Preisfrage, wir haben oben Disk[] aufgerufen, dieser Befehl sollte eigentlich einen ausgefüllten Kreis erzeugen keine Ellipse.
Wieso passiert das? Tipp: GoldenRatio hat damit zu tun.)

Da stellt sich ja jetzt die einfache Frage: ist LineLegend/PointLegend/etc. ein Graphics-Objekt?:

llegend=LineLegend[{Black,Red},{"B","R"}]
Plot[x,{x,0,1},
Epilog->{llegend},
Prolog->{Disk[{0.3,0.5},{0.2,0.2*GoldenRatio}]}
]

Leider nicht. Aber das ist kein Problem. Mathematica hat ein Befehl um nicht-Graphics-Objete als ebend solche einzubinden. Inset[]

Plot[x,{x,0,1},
Epilog->{Inset[llegend]}
]

Inset[] übernimmt, dann im zweiten Argument die Poisitionierung der Mitte:

Plot[x,{x,0,1},
Epilog->{Inset[llegend,{0.1,0.9}]}
]

Und mit dem dritten Argument die relative Positionierung des zweiten Arguments:

Plot[x,{x,0,1},
Epilog->{Inset[llegend,{0,1},{Left,Top}]}
]

Und mit dem vierten Argument die Größe dazu (was mit LineLegend nicht geht also benutze ich hier mal einen anderen Plot):

inlinePlot=ListPlot[RandomReal[{0,1},{50,2}]]
Plot[x,{x,0,1},
Epilog->{Inset[inlinePlot,{0,1},{Left,Top},0.4]}
]

Und mit dem fünften Argument die Rotation:

inlinePlot=ListPlot[RandomReal[{0,1},{50,2}]];
Plot[x,{x,0,1},
Epilog->{Inset[inlinePlot,{0.5,0.5},{Center,Center},0.4,{1,1}]}
]

Distributionstheorie

Dieses Kapitel hier ist nur für die geeignet, die tiefer in die Materie einsteigen wollen.
Es ist definitv nicht notwendig um die Praktika zu absolvieren.

Die Distributionstheorie ist die eigentlich Theorie hinter der ganzen Fehlerrechnung.
Als Physiker bekommen wir das nur meistens nicht beigebracht. Die pythagoreische Addition
ist das Resultat der Faltung zweier Verteilungen, die Fehlerfortpflanzung ist die lineare Ordnung
von Transformationen von Distributionen, Mittelwert, Varianz, Standardabweichungen sind
durch die sogenannten Momente der Verteilungen gegeben und Regressionen sind
Maximierungen von Logarithmischen Wahrscheinlichkeiten einer Normalverteilung. All das
steckt in der Distributionstheorie und noch viel mehr. Ich schneide hier mal ein paar Themen an.
Man muss es nicht kennen oder verstehen. Aber es erlaubt einen besser über die
Fehlerrechnung zu denken/argumentieren.

Pythagoreische Addition

Wir wollen jetzt mal herausfinden, wie wir überhaupt dazu kommen dass unterschiedliche Fehler phytagoreisch addiert werden.
Unsere erste Annahme: alle Unsicherheiten sind Normalverteilt.

Definieren wir uns also mal zwei Unsicherheiten als Normalverteilungen:

dist1=NormalDistribution[μ1,σ1];
dist2=NormalDistribution[μ2,σ2];

Unsere Vorstellung, die beiden Unsicherheiten miteinander zu verbinden, ist jetzt die folgende: an jeder Stelle der Verteilung der ersten Unsicherheit,
breitet sich die zweite mit ihrer Verteilung skaliert aus. Ein bisschen grafisch dargestellt:

(*Ihr müsst den Code hier wieder nicht verstehen. Ebenfalls: die sub-Verteilungen sind neu skaliert um das verständlicher zu machen*)
subPDFs=Table[PDF[dist2/.{μ2->i,σ2->3},x]/PDF[dist2/.{μ2->0,σ2->3},0]*PDF[NormalDistribution[5,Sqrt[2]*3],0]*PDF[dist1/.{μ1->5,σ1->3},i]/PDF[dist1/.{μ1->0,σ1->3},0]*2,{i,-5,15}];
Plot[Evaluate[subPDFs~Join~{PDF[dist1/.{μ1->5,σ1->3},x],PDF[NormalDistribution[5,Sqrt[2]*3],x]}],{x,-10,20},PlotRange->All,PlotStyle->(ConstantArray[LightGray,21]~Join~{Black,Directive[Black,Dashed]}),PlotLegends->LineLegend[{Black,LightGray,Directive[Black,Dashed]},Thread[Style[{"Verteilung 1","Verteilung 2 mehrmals an der ersten verteilt","Entstehende Verteilung"},16]]]]

Im Prinzip kommen wir also an der Ergebniss von diesen verteilten Verteilungen immer genauer ran, wenn man unendlich viele der Zweiten zusammenaddiert.
Ein Integral also. Genauergesagt eine Faltung der beiden Verteilungen. Mathematica kann auch direkt ein Faltungsbefehl Convolve[] (so müssen wir das Integral nicht direkt eingeben):

expr=Convolve[PDF[dist1,x],PDF[dist2,x],x,y]

So. Wir vermuten jetzt mal, dass das wieder eine Normalverteilung ist. Wenn dass so sein sollte, dann sollten wir Problemlos die Parameter dazu finden:

simplifiedExpr=FullSimplify[expr,μ1∈Reals&&μ2∈Reals&&σ1>0&&σ2>0&&y∈Reals];
sol=FullSimplify[PDF[NormalDistribution[μ,σ],y]==simplifiedExpr,μ∈Reals&&μ1∈Reals&&μ2∈Reals&&σ>0&&σ1>0&&σ2>0&&y∈Reals]

Okay cool . Durch Koeffizientvergleich sehen wir sofort:


Das prüfen wir nochmal schnell kurz:

PDF[NormalDistribution[μ1+μ2,Sqrt[σ1^2+σ2^2]],y]==simplifiedExpr//Simplify

Cool. Wir sehen also explizit, die neue Breite der Verteilung ist die phytagoreische Addition der vorangegangenen Breiten.
Wie zu erwarten war. Der Fall für  viele Verteilungen lässt sich btw. sehr leicht durch Induktion erschlagen.

Andere Distributionen

Als Physiker machen wir in 99,9% der Fälle die arrogante Annahme, dass alles Normalverteilt ist.
Ich mein, der zentrale Grenzwertsatz hilft uns zwar bei dieser Annahme aber in der realen Welt ist nunmal
nicht alles so well behaved. Zumindest die Teilchenphysiker unter euch werden mit mindestens mit
der Poissonverteilung zu tun bekommen. Anders als die Normalverteilung, ist die Poisson-Verteilung diskret.
Das hierbei dann die ganze Fehlerrechnung eine andere wird, sollte intuitiv klar sein.
Bspw. ist die Standardabweichung jetzt gegeben über die simple Formel:

StandardDeviation[PoissonDistribution[n]]

Also einfach die Quadratwurzel des Mittelwerts. Also hier einfach die empirische Formel für die
Standardabweichung zu nehmen funktioniert hier nichtmehr. Noch schlimmer eigentlich:
die Poissionverteilung ist asymmetrisch!

DiscretePlot[PDF[PoissonDistribution[5],x],{x,0,20}]

(Wie man mit asymmetrischen Verteilungen umgeht: nächstes Unterkapitel)
Und ihr könnt euch sicher sein, dass die Fehlerfortpflanzung wie ihr sie kennt so auch nichtmehr funktioniert.
(Siehe Unterkapitel Fehlerfortpflanzung in ∞ Ordnung)
Aber das ist noch nicht das schlimmste. Ihr könnt auch Verteilungen haben, die abgeschnitten sind, gemischt
sind zwischen verschiedenen Verteilungen oder am schlimmsten: Verteilung bei denen die Momente nicht
definiert sind. Alle drei Fälle findet man in der Quantenoptik (yey, ich sprech aus leidlicher Erfahrung). Ein typisches
Beispiel ist z.B. die Lorentzverteilung oder wie Mathematiker sie nennen: Cauchy-Verteilung:

dist=CauchyDistribution[μ,δ];
Plot[PDF[dist/.{μ->5,δ->2},x],{x,-5,15}]

Wieso diese Verteilung doof ist? Naja seht selbst:

StandardDeviation[dist]

Ihr denkt das ist schon schlimm? Hehe:

Mean[dist]

Die Momente größer null sind bei dieser Verteilung einfach nicht definiert. Die Integrale dazu konvergieren nicht.
Immerhin, einen Median hat sie:

Median[dist]

Da bei symmetrischen Verteilungen die Mittelwerte und Mediane gleich sind (wenn sie existieren) können wir fast nen
Auge zudrücken.
Das führt dann dazu, dass man sich in sehr sehr tiefes, dunkles Wasser begiebt wenn man anfängt mit solchen Größen
Fehlerrechnung zu machen. Interessiert das die Leute? Nah. Wir machen trotzdem Gaußsche Fehlerrechnung. Aber im
Prinzip ist das falsch.
Behaltet das also im Hinterkopf: es ist nie falsch mal zu Prüfen ob die Größen die man da hat, auch tatsächlich
Gaussverteilt sind.

Asymmetrische Distributionen

Asymmetrisch Gaußverteilt

Asymmetrische Verteilungen kommen häufiger vor als einem lieb sein sollte. Selbst wenn man nur
mit symmetrischen Verteilungen arbeitet, kann es schnell passieren, dass die asymmetrisch werden.
Wir beschäftigen uns hier mal kurz wie man die aus einem Datensatz findet.
Also kreiren wir mal schnell einen:

data=RandomVariate[SplicedDistribution[{1/2,1/2},{-∞,5,∞},{NormalDistribution[5,2],NormalDistribution[5,3]}],100000];
Histogram[data]

Okay Cool. Wie finden wir jetzt die Parameter? Also bevor wir das klären, müssen wir erstmal gucken mit
welcher Verteilung wir überhaupt arbeiten wollen. Dazu gibt es geschnittene Verteilungen. Wir schneiden
also zwei Normalverteilungen in der Hälfte am Mittelwert mit zwei unterschiedlichen Breiten.
Also zum Beispiel hier ein :

dist=SplicedDistribution[{1/2,1/2},{-∞,μ,∞},{NormalDistribution[μ,σ1],NormalDistribution[μ,σ2]}]
Plot[PDF[dist/.{μ->5,σ1->2,σ2->3},x],{x,-5,15}]

Sieht doof aus aber so wird das im allgemeinen modelliert.
Jetzt müssen wir aus den Daten nurnoch die Parameter finden. Mathematica hat natürlich dafür eine Funktion FindDistributionParameters[].
Die ist aber relativ langsam. VORALLEM mit SplicedDistributions. Wir bedienen uns einfach einem anderen Trick:
Am Ende wollen wir ja wieder 1σ Intervalle angeben. Das heißt wir müssen nur die Fläche finden mit der wir 68.3% einschließen.
Die Hälfte davon, also ungefähr 34%, soll zu jeder Seite zentriert rausgehen. Bleibt noch eine Frage:
worum soll diese Fläche zentriert sein. Dem Mittelwert? Nah! Der Mittelwert wird durch die Asymmetrie irgendwo hingeschoben. Ab hier gibt man eher den Median an.
Eine passendes Maß das uns hier hilft ist die quantile, in Mathematica berechenbar mit der sogenannten Quantile[] Funktion.
Hier können wir eine Prozentzahl an Fläche angeben und Mathematica berechnet bei welchen Punkt diese Erreicht ist.
Der Median ist natürlich klarerweise bei 50%, also 0.5.
(Ob wir jetzt nur Quantile[] oder Median[] dafür benutzen ist fast egal. Quantile[] interpoliert aber noch. Das ist vorallem bei wenig Daten vorteilhaft.)

med=Quantile[data,0.5]
Median[data]

Jetzt müssen wir nurnoch an das 1σ kommen. Rufen wir uns in Erinnerung dass %. Das heißt die Hälfte davon geht von 1/2 aus in beide Richtungen:

einSigma=Erf[1/Sqrt[2]];
{s1,s2}=Quantile[data,{1/2-einSigma/2,1/2+einSigma/2}]

Nun brauchen wir nurnoch den Abstand zum Median und schon sind wir fertig:

{sigma1,sigma2}={med-s1,s2-med}

Dass kann man sogar so in Arounds packen!

value=Around[med,{sigma1,sigma2}]

Und das passt auch. Die Daten wurden ja erzeugt mit .
Mit diesen Arounds kann man sogar normal rechnen:

-1/value

Asymmetrisch nach Distribution

Wenn man jedoch eine Variable hat, die nach einer speziellen Verteilung läuft, dann kommt man nicht drumm rum diese zu benutzen.
Nehmen wir als Beispiel mal die -Verteilung:

dist=GammaDistribution[α,β];
Plot[PDF[dist/.{α->5,β->0.75},x],{x,0,12}]

Und erzeugen zufällige Beispieldaten dafür:

data=RandomVariate[dist/.{α->5,β->0.75},100000];
Histogram[data]

Wie bereits erwähnt, benutzen wir FindDistributionParameters[]

foundParameters=FindDistributionParameters[data,dist]

Cool. Aber wir wollen das ganze in der Around Form. Wir könnten jetzt also wieder die ganze Quantile[] Keule auspacken:

{s1,med,s2}=Quantile[dist/.foundParameters,{1/2-Erf[1/Sqrt[2]]/2,1/2,1/2+Erf[1/Sqrt[2]]/2}];
Around[med,{med-s1,s2-med}]

Viel Schreibarbeit. Hätte Around für uns aber auch gemacht:

value=Around[dist/.foundParameters]

Ein Wort der Warnung: Around merkt sich nicht unter welcher Distribution es aufgerufen wurde, es macht jetzt happily
mit dem SplicedDistribution-Normalverteilt Modell weiter:

value["Distribution"]

Fehlerfortpflanzung in ∞ Ordnung

Dieses Kapitel nicht vor alle anderen hier zu ziehen war schwierig. Jedes Beispiel im Distributionskapitel bisher kann man mit der
Technik machen, die wir hier kennenlernen werden.

Fehlerfortpflanzung mit höheren AroundReplace Ordnungen

Bevor wir anfangen, möchte ich erst ein kleines Problem der linearen gaußschen Fehlerfortpflanzung aufzeigen.
Nehmen wir mal an, wir haben einen Winkel gemessen. Der kam raus mit . Also ein rechter Winkel.
Cool. Und jetzt rechnen wir den Sinus von diesem Winkel aus.

Sin[90°]

Easy. Okay und jetzt mit Fehlerrechnung:

value=Around[90°,3°]
Sin[value]

wtf? (10^(-18)) als Unsicherheit? Ja genau, du mich auch.
In fact, was auch immer das Problem ist: ich kann es auch in anderen Fällen erzeugen:

AroundReplace[2x^3+1/2*x^2+2,{x->Around[-1/6,1/24]}]

Ja äh Hallo? Wo ist meine Unsicherheit?!
Irgendwas läuft also bei diesen Beispielen falsch. Bloß was?

Gehen wir der Sache mit dem Sinus mal auf dem Grund. Was machen wir denn symbolisch im Hintegrund?

Sin[Around[μ,σ]]

Ah, jetzt wirds klarer. Der Fehler geht explizit mit . Und diese Funktion kann null werden!
D.h. egal wie unfassbar groß meine Unsicherheit werden kann, der Fehler bliebe null wenn ich 90° messe.
Macht das Sinn? Sicher nicht.
Das Problem hier ist: das die Gaußsche Fehlerfortpflanzung wie wir sie kennen und Around sie benutzt,
ist eine Rechnung erster Ordnung. Und die ist bei  nunmal leider null. Die höheren Ordnungen sind
es nicht aber die Betrachten wir ja gewöhnlich nicht.

Das gleiche beim Polynom-Beispiel. Der Wert  ist gerade eine Nullstelle der Ableitung des Polynoms
. Die höheren Ordnungen müssen hier ins Spiel kommen.
Glücklicherweise kann AroundReplace das als drittes, optionales Argument in der Funktion! Bis zur vierten Ordnung :

?AroundReplace

Also probieren wir das mal mit unseren Problem und gucken uns alle 4 Ordnungen an:

Grid[{
{"Ordnung","Formel d. Fehlerfortpflanzung","Ergebnis"},
{"Linear    ",AroundReplace[Sin[A],{A->Around[α,uα]},1],AroundReplace[Sin[A],{A->Around[90°,3°]},1]},
{"2. Ordnung",AroundReplace[Sin[A],{A->Around[α,uα]},2],AroundReplace[Sin[A],{A->Around[90°,3°]},2]},
{"3. Ordnung",AroundReplace[Sin[A],{A->Around[α,uα]},3],AroundReplace[Sin[A],{A->Around[90°,3°]},3]},
{"4. Ordnung",AroundReplace[Sin[A],{A->Around[α,uα]},4],AroundReplace[Sin[A],{A->Around[90°,3°]},4]}
},Frame->All,Background->{None,{LightGray}}]

Wir sehen hier ein paar Sachen. 
1. die zweite Ordnung reicht. 
2. Die Fehlerfortpflanzungen werden eklig kompliziert
3. Sogar die Schätzer sind betroffen?!
Der dritte Fakt ist natürlich erstmal komisch. Wieso sollte den mein Schätzer etwas anderes als  sein?
Die Antwort gibt es im Unterkapitel “Erzeugen von Asymmetrien durch Fehlerfortpflanzung”

Erzeugen von Asymmetrien durch Fehlerfortpflanzung

Wir beginnen wieder mit dem Sinus Beispeil vom vorhergehenden Unterkapitel. Dort hatten wir uns gefragt
wieso eine Fehlerfortpflanzung höherer Ordnung auch unseren Schätzer/Mittelwert beeinflusst.
Und die Antwort darauf ist leicht zu sehen, wenn man sich mal eine Beispielscharr
(Scharr = Eine Menge an normalverteilten Punkten in diesem Fall) anguckt und durch den Sinus jagt:

SmoothHistogram[{RandomVariate[NormalDistribution[90°,3°],1000000]}]
SmoothHistogram[{Sin[RandomVariate[NormalDistribution[90°,3°],1000000]]},PlotRange->All]

Sinus ist eine nichtlineare Funktion und ändert unsere Distribution von einer symmetrischen, zu einer asymmetrischen
Verteilung. Und in solchen verschiebt sich der Mittelwert. Deswegen schieben die höheren Ordnungen den auch weg.
Wir vergleichen auch hier mal Mittelwert und Median.

Mean[Sin[RandomVariate[NormalDistribution[90°,3°],1000000]]]
Median[Sin[RandomVariate[NormalDistribution[90°,3°],1000000]]]

Und der Median geht immerhin besser gegen 
Beim Sinus ist das ja vielleicht offensichtlich aber bereits ein simples quadrat macht alles kaputt:

SmoothHistogram[{(RandomVariate[NormalDistribution[0.50,0.1],1000000])^2},PlotRange->All]

Wie können wir es also richtig machen? Naja ab zum nächsten Unterkapitel um das rauszufinden.

Fehlerfortpflanzung mit Distributionen

Gehen wir wieder mal mit unserem Sinus Beispiel. Wir haben eine Zufallsvariable α welche Normalverteilt ist.
Nun fragen wir uns wie sieht die Verteilung aus wenn wir darüber einen Sinus legen. Dazu muss man die Erzeugende
der Momente Transformieren. Mathematica kann das mit dem Befehl TransformedDistribution:
(das \[Distributed] zeichen kann man eingeben mit [ESC]dist[ESC]. Oder man benutzt Distributed[α,NormalDistribution[Pi/2,σ]])

tf=TransformedDistribution[Sin[α],{α\[Distributed]NormalDistribution[Pi/2,σ]}]

Das hilft uns jetzt scheinbar nicht viel weiter. Aber tatsächlich können wir Mathematica nach dem Mittelwert und Unsicherheit dieser Verteilung Fragen:

m=Mean[tf]
s=StandardDeviation[tf]

Und so bekommen wir das bekannte Ergebniss von :

Around[m,s]/.{σ->(3°)}

So. Als zweites Beispiel gucken wir uns mal Produkt von zwei unabhängigen Variablen an.

tf=TransformedDistribution[A*B,{A\[Distributed]NormalDistribution[μ1,σ1],B\[Distributed]NormalDistribution[μ2,σ2]}]

Und wir fragen anch Mittelwert und Standardabweichung

Mean[tf]
StandardDeviation[tf]

Unser Schätzer ist genauso wie wir ihn erwartet haben. Aber die Unsicherheit ist um den Term  reicher als erwartet.

Ist das dass gleiche wie das Quadrat einer Variable?

tf=TransformedDistribution[A*A,{A\[Distributed]NormalDistribution[μ,σ]}]
Mean[tf]
StandardDeviation[tf]

Nein! Die Selbstkorrelation zwingt den Mittelwert zu wandern und die Unsicherheit ändert sich ebenfalls.
Wir können das aber wieder nachahmen, indem wir eine Korrelation zwischen den vorher unabhängigen von A und B einführen.
Das machen wir mit einer mehrdimensionalen Verteilung. In dem Fall, eine Multinormal-Verteilung.
Diese nimmt nur den Mittelwert-Vektor und die Kovarianzmatrix. Da wir den Korrelationsfaktor  setzen, ist die leicht ausgerechnet:

tf=TransformedDistribution[A*B,{{A,B}\[Distributed]MultinormalDistribution[{μ1,μ2},{{σ1^2,σ1*σ2},{σ1*σ2,σ2^2}}]}]
Mean[tf]
StandardDeviation[tf]

Wenn man nun  und  setzt, kommt das gleiche Ergebnis von der Quadrat Berechnung an:

Mean[tf]/.{μ1->μ,μ2->μ,σ1->σ,σ2->σ}
StandardDeviation[tf]/.{μ1->μ,μ2->μ,σ1->σ,σ2->σ}

Als letztes, leichteres Beispiel mal das Polynom von vor zwei Subkapiteln:

tf=TransformedDistribution[2x^3+1/2*x^2+2,{x\[Distributed]NormalDistribution[-1/6,1/24]}]
Around[Mean[tf],StandardDeviation[tf]]

Pythagoreisches Addieren

Wir erinnern uns, als wir das erste mal versucht haben das phytagoreische Addieren zu erklären, mussten wir ein Faltungsintegral knacken.
Die Vorstellung war zwar recht anschaulich aber jetzt haben wir ja TransformedDistribution im Gepäck.
An sich, Unsicherheiten zusammenzufügen, ist das selbe als einfach nur Normalverteilungen mit Mittelwert null in einer Transformation zu addieren.
Also ist es wirklich nur das?

tf=TransformedDistribution[A+B+c,{A\[Distributed]NormalDistribution[0,σ1],B\[Distributed]NormalDistribution[0,σ2],c\[Distributed]NormalDistribution[0,σ3]}]
StandardDeviation[tf]

Wow. Okay das war einfach. Mathematica hat sogar die TransformedDistribution in die NormalDistribution vereinfacht die dahinter steht.
Jetzt haben wir also noch einen weiteren Blickwinkel darauf.

Log-Likelihood

Okay, last but not least, gucken wir uns nochmal Regressionen an. Wir hatten es noch nicht direkt erwähnt aber ich sage es hier mal:
Regressionen finden die Parameter durch Minimierungen der Summe der Abstandsquadrate . Produziert Sinn. Irgendwie.
Die Abstände der Fitfunktion zu den Daten minimieren ergiebt Sinn und das Quadrat bestraft weitere Abstände stärker als kleinere.
Macht das  einen Unterschied? Naja für die Optimierung nicht. Aber sehr wohl für die Fehlerberechnung. Woher weiß ich das?
Naja weil die Formel sicher nicht vom Himmel fällt.😅
Fun-Fact Mathematica kann Emoji’s. Ist mir gerade selbst so aufgefallen.
Nachtrag 30.06.2022 - Mathematica kann jetzt sogar Emojis als Variablen definieren und benutzen. The possibilities are endless :D
Anyway, back to topic 💩
Versuchen wir mal den ganzen Spaß herzuleiten. Schauen wir uns mal nur einen Datenpunkt  an.
Als erstes eine Standardannahme . Somit reduzieren wir unser Problem auf effektiv . Jetzt fragen wir uns:
was ist den die Wahrscheinlichkeit, dass unsere Fitfunktion  genau bei diesem  Wert sitzt. Das können wir Direkt
in Form der PDF (probability density function - Wahrscheinlichkeitsdichtefunktion) ausrdücken:

PDF[NormalDistribution[y,σy],F[x,params]]

Oder direkt ohne einen Mittelwert anzugeben kann man auch direkt den Abstand als Argument um null angeben:

PDF[NormalDistribution[0,σy],F[x,params]-y]

Macht keinen Unterschied.
Das gilt jetzt also für einen Punkt. Es ist klar, dass wir die Wahrscheinlichkeit Maximieren wollen.
Also wir suchen unsere Parameter so, dass diese Wahrscheinlichkeit maximal wird.
Super Plan. Wie sieht es denn jetzt aber mit  vielen Punkten aus? Naja wenn man mehrere Wahrscheinlichkeiten
hat dass etwas passiert, dann multiplizieren diese Wahrscheinlichkeiten zur Gesamtwahrscheinlichkeit .

Jetzt könnten wir das so machen, implementieren und uns nen Loch in den Bauch freuen.
Nur leider ist das Problem so relativ schwierig gestellt. Exponentialfunktionen fallen schnell, sind
aufwendig auszurechnen (vorallem wenn man tausende von Datenpunkten oder mehr hat) und generell
wird das Problem sehr teuer und schwer auszurechnen. Es gibt viele Gründe wie ich das jetzt rechtfertigen
kann aber ich sag es einfach: Wir machen das nicht. Stattdessen schauen wir uns das Problem logarithmiert an.

Okay. Erstmal, Summen sind schöner als Produkte und die Exponentiale werden durch den Log geschluckt werden:

expr=FullSimplify[Log[PDF[NormalDistribution[y,σy],F[x,params]]],σy>0&&y∈Reals&&F[__]∈Reals]

Nebenbei, hat Mathematica dafür natürlich direkt einen Befehl parat:

LogLikelihood[NormalDistribution[0,σy],{y-F[x,params]}]

Sieht erstmal komisch aus. Aber die beiden hinteren Terme sind Konstanten und weder für die Optimierung
wichtig, noch für die Fehlerberechnung (diese läuft über Ableitungen der Parameter). Ergo brauchen wir die nicht direkt.

Plus@@Select[List@@expr,Not@*FreeQ[F]]

Okay fast am Ziel. Also was sehen wir: Wir sehen das Abstandsquadrat, check. Wir sehen das 1/2 das ich angesprochen hatte, check.
Wir sehen instrumentelle Gewichtung, ach da kommt die her, check. Das Minus können wir zum Plus machen in dem wir die
Maximierung einfach zu einer Minimierung machen. Und schon haben wir hergeleitet wieso wir fitten wie wir fitten:
Minimierung des Abstandsquadrat mit instrumenteller Gewichtung. (Und nen Faktor 1/2 davor, meh)

Also machen wir das einmal in unserem Leben manuell.
Dafür erstellen wir zuerst ein Beispieldatenset und fitten das auf die herkömmliche Weise:

data=Table[{x,2x-1+RandomVariate[NormalDistribution[0,0.1]]},{x,-1,3,0.1}];
fit=LinearModelFit[data,{x,1},x];
fit["ParameterTable"]
fit["CovarianceMatrix"]//MatrixForm

Diese Parameterabschätzungen und Kovarianzmatrix  wollen wir bekommen.
Also Summieren wir also mal alle LogLikelihood’s auf:

expr=LogLikelihood[NormalDistribution[0,σ],(A*data[[All,1]]+B-data[[All,2]])]

Was ist eigentlich hier unser ? Naja wenn wir keine Gewichte angegeben haben ist das die abgeschätzte Varianz. Siehe VarianceEstimatorFunction vom Regressionskapitel.
Maximimieren wir also mal die LogLikelihood:

sol=Last@NMaximize[{expr,σ>0},{A,B,σ}]

Und schon haben wir unsere Parameter


und die abgeschätzte Varianz:

Die Herleitung für die Kovarianzmatrix mache ich jetzt mal nicht aber erzähle, dass es das negative inverse der Hessematrix skaliert mit dem Erwartungswert der Varianz ist () ist (p = Anzahl an Parametern):

Length[data]/(Length[data]-2)*Inverse[-D[expr,{{A,B},2}]/.sol]//MatrixForm

Cool. Wir haben also das exakt gleiche Ergebnis wie LinearModelFit. Unsere Theorie ging also auf.

Als Schlusswort: das wir das Problem logarithmiert haben, hatte garnicht soviel Einfluss wie man denken würde. 
Da unsere PDF’s einen Wertebereich  einschließen, ist der Logarithmus auf diesem Interval monoton. D.h.
wir finden mit dem Logarithmierten Beispiel das Selbe Ergebniss wie mit den nicht logarithmierten:

expr=Likelihood[NormalDistribution[0,σ],(A*data[[All,1]]+B-data[[All,2]])]
sol=Last@NMaximize[{expr,σ>0},{A,B,σ},Method->"DifferentialEvolution"]

Yey. Das ist das gleiche Ergebnis (ignoriert die Warnung). Allerdings musste ich vorher den Optimierungsalgorithmus ändern, da die Standardwahl nicht konvergiert ist.
Daran sieht man schon, wie schwer wir uns mit dem nicht logarithmierten Problem tun würden.

Reale Beispiele (GP1, GP2, FP, Forschung)

[O10 Linsensysteme] Die Falle des R1=R2

Im Versuch O10-Linsensysteme gibt es eine nette Falle in die man hineintreten kann.
In der Brennweitenbestimmung mit dem Sphärometer muss man eine Fehlerfortpflanzung über folgende Formel durchführen:

expr=1/((n-1)*(1/R1+1/R2))

Die ist natürlich mühsam (wenn man sie per Hand macht :D ) und wird komplex:

AroundReplace[expr,{n->Around[μn,σn],R1->Around[μR1,σR1],R2->Around[μR2,σR2]}]

Aber hier der scheinbare Retter: ihr benutzt eh höchstwahrscheinlich eine symmetrische Linse, in der die beiden Seiten den gleichen gemessenen Krümmungsradius haben.
Also können wir vermeintlich einfach  setzen und vereinfachen:

simplifiedExpr=expr/.{R1->R,R2->R}

Wieso auch nicht. Die Formel ist jetzt super easy und die Fehlerfortpflanzung ist jetzt auch relativ leicht:

AroundReplace[simplifiedExpr,{n->Around[μn,σn],R->Around[μR,σR]}]

Naja jetzt vergleichen wir doch mal. Denn wenn das okay war diese Vereinfachung zu machen, dann muss das gleiche Ergebnis in Abschätzung und Unsicherheit rauskommen:

AroundReplace[expr,{n->Around[1.458,0.005],R1->Around[100.0,2.8],R2->Around[100.0,2.8]}]
AroundReplace[simplifiedExpr,{n->Around[1.458,0.005],R->Around[100.0,2.8]}]

Na das war ja mal nichts. Die Schätzer sind zwar gleich aber die Unsicherheiten nicht.
Was haben wir denn also falsch gemacht? Durch das setzen von  haben wir einfach mal angenommen, dass beide zu 100% korreliert sind. Also .
Oder in anderen Worten geschildert. Wir haben die Wahrheitswerte gleichgesetzt und gesagt, die Abweichungen (Unsicherheiten) sind identisch zwischen  und .
Das ist nun wirklich nicht der Fall. Diese beiden Seiten sind “unabhängig” voneinander hergestellt und existieren auch so ohne aufeinander Einfluss zu haben in der Messung
mit dem Sphärometer.
In diesem Fall der Fehlerpropagation führt eine positive korrelation zwischen den Variablen zu einer Erhöhung der Unsicherheit.
Wir können diesen Fall sogar simulieren indem wir diese Korrelation zwischen  und  einführen:

AroundReplace[expr,{n->Around[1.458,0.005],{R1,R2}->VectorAround[{100,100},{{2.8,2.8},1}]}]

Schön! Die exakt gleiche, falsche Unsicherheit wie oben.
Was könnt ihr also aus diesem Beispiel mitnehmen:
Nur weil zwei Größen im Schätzer gleich sind, sind sie nicht selbst gleich oder haben identische Unsicherheiten.
Oder in kurz: Nur weil zwei Dinge gleich sind, sind sie nicht dasselbe.

[A1: Äußerer Photoeffekt] Bestimmung der maximalen Gegenspannung (Subscript[U, g,max])

ACHTUNG
Durch einen Bug in AroundReplace mit Root Objekten funktioniert dieses Beispiel NICHT unter Mathematica v13.0.0
Das Problem wurde in Mathematica v13.0.1 behoben. Solltet ihr also v13.0.0 haben, dann updated auf  v13.0.1 oder höher.

Im Versuch A1: Äußerer Photoeffekt muss durch Regression der Photokurve innerhalb des Modells
einer Diode im Anlaufbereich gearbeitet werden.
Das Modell ist gleich für beide Kurven und läuft folgendermaßen:

Iph1=Isat1*(Exp[(U-U01)/ηUT1]-1)
Iph2=Isat2*(Exp[(U-U02)/ηUT2]-1)

Dabei fassen wir  zu einem Parameter zusammen.
Die Aufgabe besteht nun darin, für zwei verschiedene Intensitäten diesen Verlauf zu messen,
die Regression nach dem gegebenen Modell durchzuführen um dann den Schnittpunkt der beiden
Verläufe zu bestimmen. Der Punkt an dem sich diese Verläufe schneiden, ist dann (Subscript[U, g,max]).
An sich klingt das nicht zu schwer. Also machen wir mal die Beispieldatensätze und fitten
mal an beide das Modell:

dataI1=Table[{U,76.0*(Exp[(U+1.36)/0.13]-1)+RandomVariate[NormalDistribution[0,1]]},{U,Range[-7,-2,0.5]~Join~Range[-2,-1.25,0.1]}];
dataI2=Table[{U,28.0*(Exp[(U+0.89)/0.57]-1)+RandomVariate[NormalDistribution[0,1]]},{U,Range[-7,-2,0.5]~Join~Range[-2,-1,0.2]}];
fit1=NonlinearModelFit[dataI1,Iph1,{Isat1,U01,ηUT1},U];
fit2=NonlinearModelFit[dataI2,Iph2,{Isat2,U02,ηUT2},U];
Show[ListPlot[{dataI1,dataI2},PlotRange->All],Plot[{fit1[u],fit2[u]},{u,-7,-1}]]

Wir sehen zuerst, dass die Regression nicht konvergiert. Also versuchen wir mal wieder schnell Startwerte zu finden.
Wenn man nun den Limes gegen minus unendlich zieht, sieht man schnell das dieser nur durch -(Subscript[I, sat]) bestimmt ist:

Limit[Iph1,U->-Infinity]

Die Werte lassen sich leicht ablesen und die Angabe als Startwert reicht vollkommen das Modell zum konvergieren zu bringen:

fit1=NonlinearModelFit[dataI1,Iph1,{{Isat1,70},U01,ηUT1},U]
fit2=NonlinearModelFit[dataI2,Iph2,{{Isat2,30},U02,ηUT2},U]
Show[ListPlot[{dataI1,dataI2},PlotRange->All],Plot[{fit1[u],fit2[u]},{u,-7,-1}]]

Soweit so einfach. Jetzt müssen wir also nurnoch gleichsetzen, lösen und mit Fehlerfortpflanzung rechnen. Na auf gehts:

Solve[Iph1==Iph2,U]

Ja super. Das System kann garnicht gelößt werden. Wie sollen wir das also tun? Naja wir könnten numerisch rangehen. Das gäbe uns aber erstmal keine Unsicherheit...
Aber wir machen lieber das, was Physiker am besten können. Taylor Approximationen! (Mathematicas Befehl dafür: Series[])
Wir entwickeln mal gleich um  herum, da dieser Wert nahe an (Subscript[U, g,max]) ist. Und wir starten direkt in die zweite Ordnung.

expr1=Series[Iph1,{U,U01,2}]
expr2=Series[Iph2,{U,U02,2}];
{sol1,sol2}=SolveValues[expr1==expr2//Normal,U]//Simplify

Huiiii. Naja ist nicht hübsch, aber es sind Lösungen. Eine davon ist kein tatsächlicher Schnittpunkt.
Wir werden herausfinden welcher in dem wir mal die Taylor Lösungen gegen die richtigen Verläufe Plotten
und uns anschauen, welcher der beiden Schnittpunkte näher am tatsächlichen Schnittpunkt liegt:

Show[ListPlot[{dataI1,dataI2},PlotRange->{{-2,-1},All}],Plot[{fit1[u],fit2[u]},{u,-2,-1}],Plot[Evaluate@{Normal[expr1]/.fit1["BestFitParameters"],Normal[expr2]/.fit2["BestFitParameters"]},{U,-2,-1},PlotStyle->Dashed]]

Die durchgehenden Linien sind unser fit und die gestrichelten sind die Taylor Approximationen. Wir schauen jetzt welcher Schnittpunkt der Taylor Approximation
ist näher am Schnittpunkt der Fits.
Und wir sehen: die Lösung mit den kleineren Absolutwert stimmt. Also Die Lösung mit dem Minus. Das wäre die zweite Lösung in der Variable mit dem Namen sol2.
Jetzt wo das geklärt ist: ab damit in AroundReplace[] unter voller korrelierter Fehlerfortpflanzung!

AroundReplace[sol2,{
{Isat1,U01,ηUT1}->VectorAround[fit1["BestFitParameters"][[All,2]],fit1["CovarianceMatrix"]],
{Isat2,U02,ηUT2}->VectorAround[fit2["BestFitParameters"][[All,2]],fit2["CovarianceMatrix"]]
}]

Wir überprüfen mal numerisch wo der Schnittpunkt tatsächlich sitzt:

FindRoot[Normal[Iph1/.fit1["BestFitParameters"]]==Normal[Iph2/.fit2["BestFitParameters"]],{U,-1.5}]

Also -1.3844 ist garnicht mal so gut. Aber ein Anfang.
Hilft eine Taylorreihe dritter Ordnung mehr?

expr1=Normal@Series[Iph1,{U,U01,3}];
expr2=Normal@Series[Iph2,{U,U02,3}];
sol=SolveValues[expr1==expr2,U];
Select[AroundReplace[#,{
{Isat1,U01,ηUT1}->VectorAround[fit1["BestFitParameters"][[All,2]],fit1["CovarianceMatrix"]],
{Isat2,U02,ηUT2}->VectorAround[fit2["BestFitParameters"][[All,2]],fit2["CovarianceMatrix"]]
}]&/@sol,FreeQ[_Complex]]

Ja auf jeden Fall (solange wir die unphysikalischen, komplexen Lösungen aussortieren).
Jetzt haben wir schon einen ganz angenehmen Wert dafür.
Können wir das noch weiter treiben in höhere Ordnungen? Jup, können wir!
Und es wird sich herausstellen, dass wir ab der fünften Ordnung innerhalb der signifikanten Stellen
konvergiert sind.

expr1=Normal@Series[Iph1,{U,U01,5}];
expr2=Normal@Series[Iph2,{U,U02,5}];
sol=SolveValues[expr1==expr2,U];
Select[AroundReplace[#,{
{Isat1,U01,ηUT1}->VectorAround[fit1["BestFitParameters"][[All,2]],fit1["CovarianceMatrix"]],
{Isat2,U02,ηUT2}->VectorAround[fit2["BestFitParameters"][[All,2]],fit2["CovarianceMatrix"]]
}]&/@sol,FreeQ[_Complex]]

Also wir finden finally:
(Subscript[U, g,max]) =  V


Nebenbei hier interessant: wir lösen hier ja gerade Polynome fünfter Ordnung. Dafür gibt es keine fertigen Lösungen in form von Radikalen.
Ganz im Gegenteil. Mathematica muss diese Lösung in Form von sogenannten Root-Objekten darstellen. Und jetzt sehen wir auch wieso
das Sinn macht. Mathematica kann im Hintergrund diese Root-Objete ableiten und beliebig genau ausrechnen auch ohne dass eine Form
in Elementarfunktionen bekannt ist. Und dass passiert einfach im Hintergrund ohne das wir es gemerkt haben.
Also nur zum Spaß, wie hätte diese Form den ausgesehen?

solAn=Select[{#,AroundReplace[#,{
{Isat1,U01,ηUT1}->VectorAround[fit1["BestFitParameters"][[All,2]],fit1["CovarianceMatrix"]],
{Isat2,U02,ηUT2}->VectorAround[fit2["BestFitParameters"][[All,2]],fit2["CovarianceMatrix"]]
}]}&/@sol,FreeQ[_Complex]][[1,1]];
solAnExpr=AroundReplace[solAn,{
{Isat1,U01,ηUT1}->VectorAround[{μI1,μU1,μη1},(#+Transpose@#)/2&@{{uI1^2,uIU1,uIη1},{0,uU1^2,uUη1},{0,0,uη1^2}}],
{Isat2,U02,ηUT2}->VectorAround[{μI2,μU2,μη2},(#+Transpose@#)/2&@{{uI2^2,uIU2,uIη2},{0,uU2^2,uUη2},{0,0,uη2^2}}]
}]//Shallow[#,8]&

Und das hier ist eine verkürzte Variante. Hinter jedem <<N>> steht ein weiterer langer Term.
Ich schreib mal nicht den gesamten Term aus, weil der ist mehrere Seiten lang. Aber wir können das mal messen wie lang.
Dafür gibt es den sogenannten LeafCount. Jede Operation und Variable führt hierbei zu einem Leaf.
LeafCount[] zählt die Anzahl der Leaf’s zusammen (das benutzt Simplify nebenbei um zu messen wie kompliziert die Terme sind)
Für unser Modell zum Beispiel haben wir einen Leaf-Count von 15 mit diesem Baum:

LeafCount[Iph1]
TreeForm[Iph1]

So. Dann zählen wir mal die Leafs der Lösung.

LeafCount[solAnExpr]

Über 80 tausend Leafs. Und das ist alles im Hintergrund passiert und das hat Mathematica alles für euch problemlos ausgerechnet!
Probiert das mal mit Python...(mit Julia ginge das sogar mit ein bisschen Mehraufwand)

[O8: Fraunhofer Beugung] Startparameter für einen Fit in einem GP-2 Versuch

Im Versuch O8 des optischen Blockpraktikums (Fraunhofer Beugung), muss man einen Fit einer Besselfunktion machen:

Dabei ist  die Besselfunktion erster Art und erster Ordnung. Die kann Mathematica natürlich auch:

Information[BesselJ]

Also erstellen wir mal das Fitmodell:

model=4*I0*(BesselJ[1,Pi*A*(x-x0)]/(Pi*A*(x-x0)))^2+Ib

Erstellen wir mal schnell Besipieldaten dafür und plotten das in einem LogPlot:

data=Table[{x,model*(1+RandomVariate[NormalDistribution[0,0.05]])/.{I0->9.49*10^-7,A->1.669*10^2,x0->7.732*10^-3,Ib->10^-10}},{x,0,0.04,0.04/100}];
ListLogPlot[data]

Dann probieren wir doch mal direkt die Regression aus:

fit=NonlinearModelFit[data,model,{I0,A,x0,Ib},x]
Show[ListLogPlot[data],LogPlot[fit[x],{x,0,0.04}]]

Joar, dass war ja nichts. Aber NonlinearModelFit meckert ja auch dass da was schief läuft.
Von hier aus, haben wir mehrere Möglichkeiten vorzugehen.
Zum einen, kann man einfach mal die Funktion zu den Daten plotten, und selbst mit den Werten spielen. Kein scheiss.
Man muss ja nur in die richtige Gegend kommen, den Rest mach NonlinearModelFit selbst.
Aber wir sind Physiker, wir gehen da systematisch ran. Die Idee ist leicht, einige Parameter sorgen für spezielle Features im Kurvenverlauf.
Wir versuchen diese Features zu finden aus den Daten zu approximieren und nehmen diese als Startwerte.

Einfache Features sind Maxima, Nullstellen, verhalten in Unendlichen, an speziellen Punkten usw..
Wir gucken mal was die Funktion bei unendlich tut:

Limit[model,x->Infinity,Assumptions->A>0]

Ahhja, also wenn wir ganz weit nach rechts gehen, dann laufen wir in ein Plateau dass nur durch den Parameter Ib bestimmt wird.
Gut, sehen wir in den Daten ein Plateu? Ja! Bei so in etwa :

ListLogPlot[data,Epilog->{InfiniteLine[{{0,Log[10^-10]},{1,Log[10^-10]}}]}]

Denkt darran: wir brauchen nicht den genauen Wert, nur so eine Abschätzung.
Demnach schätzen wir also als erste Variable, dass 

Okay nächste Überlegung.  Wenn immer ich eine Funktion sehe der Form  lohnt es sich zu fragen, was passiert bei . Tun wir das mal:

Limit[model,x->x0]

Uhhh okay. Ib kennen wir schon also steht da im Prinzip: . Cool. Jetzt müssen wir nurnoch herausfinden, wo  ist. Ich hab da
so den Verdacht, dass das, das erste Maximum ist. Dann würde gelten, dass  und . Checken wir mal:

dmodel=D[model,x]
d2model=D[model,{x,2}]
Limit[dmodel,x->x0]
Limit[d2model,x->x0]

YES. Da  und  hatten wir Recht! Das ist nen Maximum und wir nehmen jetzt an dass das Erste stimmt. Also gucken wir wo das passiert. Naja so bei etwa
 &

ListLogPlot[data,Epilog->{InfiniteLine[{{0,Log[10^-6]},{1,Log[10^-6]}}],InfiniteLine[{{0.0075,2},{0.0075,1}}]}]

Daraus können wir also schlussfolgern:


Damit wissen wir also drei Variablen in etwa. Reicht das?

fit=NonlinearModelFit[data,model,{{I0,10^-6},A,{x0,0.0075},{Ib,10^-10}},x]
Show[ListLogPlot[data],LogPlot[fit[x],{x,0,0.04}]]

Nope. Thats sad. ABER. Wir wenden mal einen kleinen Trick an. Wir fitten mal so, als wären unsere Bisherigen Abschätzungen absolut richtig gewesen und fitten nurnoch nach dem einzigen unbekannten Parameter :

reducedModel=model/.{I0->10^-6,x0->0.0075,Ib->10^-10}
fit=NonlinearModelFit[data,reducedModel,{A},x]
Show[ListLogPlot[data],LogPlot[fit[x],{x,0,0.04}]]
fit["ParameterTable"]

Hehe, ja lol, das war einfach. Und somit haben wir die letzte Abschätzung für .
Also probieren wir jetzt:

fit=NonlinearModelFit[data,model,{{I0,10^-6},{A,173},{x0,0.0075},{Ib,10^-10}},x]
Show[ListLogPlot[data],LogPlot[fit[x],{x,0,0.04}]]

Nice. Fertig!
Wäre probieren schneller gewesen? Vielleicht.
Kommen wir somit hier aber auf jeden Fall ans Ziel: vielleicht. In diesem Problem: sicher.

Ein letzten Nebensatz dazu: es macht bei diesem Modell vielleicht Sinn, nicht die direkte Funktion mit den Daten zu optimieren.
Eher, würde man den Log optimieren wollen oder die Gewichte so wählen, dass sie Log-gewichten. Bzw. im Versuch gibt es auch
eine Nette Herleitung, wieso man eher mit  gewichten sollte was einer log-gewichtung sehr nahe kommt. Das macht
auch überhaupt kein Problem für die Unsicherheiten wie es zuerst scheinen sollte, da die VarianceEstimatorFunction diese
zusätzliche Varianz wieder weg kürzt.
Und somit kommen wir finally zu dem sehr sehr schönen fit:

fit=NonlinearModelFit[data,model,{{I0,10^-6},{A,173},{x0,0.0075},{Ib,10^-10}},x,Weights->1/Abs[data[[All,2]]]]
Show[ListLogPlot[data],LogPlot[fit[x],{x,0,0.04}]]

[FP: Einzelphotonen Quantenkryptografie] Startparameter für einen Fit in einem FP Quantenoptik Versuch

Im FPraktikum gibt es den Quantenkryptografie Versuch. In diesem muss ein Fit einer Funktion zu aufgenommenen Daten
einer Korrelationsfunktion gemacht werden.
Diese Funktion ist relativ komplex mit sechs unabhängigen Parametern:

model=C1*(1-(k+1)*Exp[-k1*Abs[x-x0]]+k*Exp[-k2*Abs[x-x0]])+C0

Und natürlich, erstellen wir schnell Beispieldaten und vergewissern uns, dass NonlinearModelFit fehlschlägt:

data=Table[{x,model+RandomVariate[NormalDistribution[0,10]]/.{C1->193,k->0.39,k1->0.048,x0->857.4,k2->0.0033,C0->42}},{x,0,1200}];
fit=NonlinearModelFit[data,model,{C1,k,k1,x0,k2,C0},x]
Show[ListPlot[data,PlotRange->All],Plot[fit[x],{x,0,1200}]]

Oki doki. Dann beginnen wir mal. Wir starten wieder mit unserer  Vermutung:

Limit[model,x->x0]

Ja gut, dass war nicht schwer. Und wir wetten mal ohne nachzuprüfen, dass  dieses Prominente Minimum ist (was sollte es sonst sein?):

Cases[data,{_,Min[data[[All,2]]]}]

Daraus folgt:


Okay, probieren wir doch schonmal:

fit=NonlinearModelFit[data,model,{C1,k,k1,{x0,859},k2,{C0,38.04}},x]
Show[ListPlot[data,PlotRange->All],Plot[fit[x],{x,0,1200},PlotRange->All]]

Uhi, naja das Ging schnell. Tatsächlich reicht sogar schon das  anzugeben. Nevertheless, wer noch Spaß haben will:  ist noch sehr leicht zu finden.
Und  mit ein bisschen Aufwand. Und daraufhin dann schlussendlich auch noch .
Aber nicht notwendig.

Orthogonalregression ohne Gewichtung

Naiv, falsches Beispiel:

Eine jede Regression gehen wir für gewöhnlich mit der Annahme ein, dass die Unsicherheit unser abhängigen Variable  vernachlässigbar klein ist.
Jetzt spielen wir mal Devils Advocate und fragen ganz frech: und wenn das nicht so ist? Was wenn beide signifikant beitragen?!
Dann macht man für gewöhnlich eine sogenannte Orthogonalregression.
Man muss jetzt mal kurz erwähnen, dass es hier kein Richtig und Falsch mehr gibt. Es gibt Annahmen und Ansätze die man Verfolgen kann
und die bringen einen zu einem jeweiligen Ergebnis.
Orthogonalregressionen sind nicht einfach und von Modell zu Modell unterschiedlich. Es gibt auch richtige und professionelle Implementierungen
für Speziallfälle. (Bspw. YorkFit für lineare Funktionen: https://github.com/kmisiunas/mathematica-york-fit )

Wir hier, werden mit unserer LogLikelihood Methode vorgehen.
In 2d stellen wir uns die Unsicherheiten als eine Multinormal verteilung in x und y vor. Der Abstand in  ist leicht gegeben durch , fährend der Abstand in der x-Achse
gegben werden kann als .
Der Einfachheit halber, nehmen wir an dass die Korrelation zwischen  und  vernachlässigbar ist...ja irgendwo gibts immer Approximationen. Probiert es mal ohne (tut es nicht)...

singleTerm=FullSimplify[Log[PDF[MultinormalDistribution[{0,0},{{ux^2,0},{0,uy^2}}],{x-f[y],y-f[x]}]],(x-f[y])^2>0&&(y-f[x])^2>0&&ux>0&&uy>0]

Das ist also der Einzelterm den wir für jeden Datenpunkt  einsetzen werden. Wir sehen auch direkt, dass wir etwas in der Form  optimieren. Das sieht doch schonmal
Orthogonal aus. Das ist gut.
Eine Sache noch dazu, wir brauchen nichtnur  sondern ja auch .
Uff das ist ätzend. Naja willkommen bei Orthogonalregressionen. Da ist nunmal alles ätzend.
Um es nicht ganz so doll ätzend zu machen, entscheiden wir uns auch mal für eine lineare Funktion als Modell.
Also schreiben wir mal den zu optimierenden Ausdruck auf und kreiren zuerst Beispieldaten:

(*Diesmal mit Rauschen in y UND x*)
data=Table[{x+RandomVariate[NormalDistribution[0,0.3]],2x-1+RandomVariate[NormalDistribution[0,0.3]]},{x,-1,3,0.1}];
ListPlot[data]

Super. Nah auf gehts und das LogLikelihood Modell aufstellen:

model=A*x+B;
invModel=First[SolveValues[model==y,x]];
expr=Total[((singleTerm/.{f[x]->model,f[y]->invModel})/.{x->#1,y->#2})&@@@data]

Bevor wir irgendwie weiterrechen, wie performt den eine normale Regression?

fit=LinearModelFit[data,{x,1},x];
fit["ParameterTable"]
Show[ListPlot[data],Plot[fit[x],{x,-1,3}]]

Tatsächlich garnicht mal so gut (natürlich ist dieses Szenario gerade so gewählt).

values=Around@@@Transpose[{fit["BestFitParameters"],fit["ParameterErrors"]}]

Bzw. hält sich die Regression für besser als sie sein sollte! Mit ziemlichen Abstand treffen die angegeben Werte nicht die Wahrheitswerte (2 & -1):

{2,-1}-values

Also machen wir uns mal an den Term den wir hergeleitet haben und optimieren den:

sol=Last[NMaximize[{expr,ux>0&&uy>0},{A,B,ux,uy}]]

Okay, sieht schonmal besser aus auf den ersten Blick. Und die Varianzen scheinen zumindest bei x besser zu passen.
Naja berechnen wir mal die Unsicherheiten:

covMat=Length[data]/(Length[data]-2)*Inverse[-D[expr,{{A,B},2}]/.sol];
MatrixForm[covMat]
solUncert=Sqrt[Diagonal[covMat]]

Coolio. Packen wir das mal in Around’s:

values=Around@@@Transpose[{{A,B}/.sol,solUncert}]

Und wir schauen wie nahe wir am Wahrheitswert vorbei sind

{2,-1}-values

Und hier sehen wir eine tolle Sache. Was unsere relative Genauigkeiten angeht, haben wir uns nicht verbessert. Unsere Unsicherheiten
sind immernoch um in etwa den selben Anteil am Ergebnis falsch.
However, unsere absolute Präzision und Genauigkeit ist besser geworden. Wir sind näher am Wahrheitswert und uns auch sicherer dass wir da sind.
Also, Fazit: ganz nett.

Richtiges, nicht naives Beispiel

Na? Ist euch der Fehler im naiven Beispiel aufgefallen?
Eine der ersten Annahmen die wir gemacht haben schien sinnvoll, unschuldig und trivial richtig.
Namely,
zu sagen der Abstand in  und  ist gegeben durch  und  war zwar richtig aber gibt uns nicht den Orthogonalen Abstand wieder.
Wir sehen das leicht wenn wir uns das mal plotten:

point={2.5,1.5};
func=0.5x+1;
inverseFunc=(y-1)/0.5;
Plot[func,{x,0,3},PlotRange->{{0,3},{0,3}},Epilog->{PointSize[0.02],Point[{point}],
Line[{point,point-{0,y-func}}/.{x->point[[1]],y->point[[2]]}],
Line[{point,point-{x-inverseFunc,0}}/.{x->point[[1]],y->point[[2]]}],Blue,
Line[{point,point-{x-inverseFunc,y-func}/2}/.{x->point[[1]],y->point[[2]]}]
},AspectRatio->1]

Wenn wir uns jetzt mal die blaue Linie ansehen,
dann sollte diese eigentlich in 90° winkel mit der linearen Funktion schneiden damit wir das orthogonal nennen können.
Man braucht jetzt die Augen nicht besonders weit aufmachen um zu sehen, dass das im Leben nicht hinkommt.

Wir müssen also herleiten, wie wir den Orthogonalabstand richtig bestimmen.
Wir beginnen also bei der Standardgeschichte von funktionalen blabla und landen beim normalen phytagoreischen Abstand
eines Punktes  mit der parametrisierten Funktionsgrade

dist=Sqrt[(px-x)^2+(py-f[x])^2]

Orthogonal werden wir im minimalen Abstand also leiten wir nach den freien Parameter  ab:

D[dist,x]==0

Wir sehen, wenn wir annehmen, dass der Punkt nicht genau auf der Linie liegt, dann können wir den Nenner wegmultiplizieren:

Simplify[Numerator[D[dist,x]]==0]

Und hier haben wir es: die Gleichung um eine Orthogonale zu einer beliebigen Funktion zu finden.
(Nebenbei, dass ist eine DGL... die könntet ihr also noch lösen. Aber müssen wir nicht.)
Nur leider hat die bereits bei sehr einfachen Problemen mehrere Lösungen. Aber wir interessieren
uns jetzt mal nur für lineare Funktionen

expr=Simplify[Numerator[D[dist/.f[x]->A*x+B,x]]==0]

Sehr gut, dass können wir vermutlich easy lösen:

sol=Solve[expr,x]

Cool!  haben wir,  ist trivialerweise einsetzen.

deltaPX=Simplify[First[(x/.sol)]]
deltaPY=Simplify[First[(A*x+B/.sol)]]

Checken wir schnell nochmal mit nem Plot:

point={2.5,1.5};
func=0.5x+1;
inverseFunc=(y-1)/0.5;
Plot[func,{x,0,3},PlotRange->{{0,3},{0,3}},Epilog->({PointSize[0.02],Point[{point}],
Line[{point,point+{0,deltaPY-y},point+{0,deltaPY-y}+{deltaPX-x,0}}],
Line[{point,point+{deltaPX-x,0},point+{deltaPX-x,0}+{0,deltaPY-y}}],Blue,
Line[{point,point+{deltaPX-x,deltaPY-y}}]
}//.{x->point[[1]],y->point[[2]],A->0.5,B->1,px->x,py->y}),AspectRatio->1]

Yes. Das sieht orthogonal aus. Wir können das auch schnell durch Rechnung beweisen:

ArcTan[0.5]-Pi/2 (*Anstieg der Gerade um 90° rotiert, also orthogonal zum eigentlichen Anstieg*)
ArcTan[(deltaPY-y)/(deltaPX-x)//.{x->point[[1]],y->point[[2]],A->0.5,B->1,px->x,py->y}]

Gleicher Anstieg. Perfekt.
Jetzt können wir das benutzen, um in unseren Ansatz von oben zu gehen:

singleTerm=FullSimplify[Log[PDF[MultinormalDistribution[{0,0},{{ux^2,0},{0,uy^2}}],{x-deltaPX,y-deltaPY}/.{px->x,py->y}]],x∈Reals&&y∈Reals&&A∈Reals&&B∈Reals&&ux>0&&uy>0]
expr=Total[(singleTerm/.{x->#1,y->#2})&@@@data]//Simplify

Und lösen!

sol=Last[NMaximize[{expr,ux>0&&uy>0},{A,B,ux,uy}]]
covMat=Length[data]/(Length[data]-2)*Inverse[-D[expr,{{A,B},2}]/.sol];
MatrixForm[covMat]
solUncert=Sqrt[Diagonal[covMat]]

Packen wir das wieder in Around’s:

values=Around@@@Transpose[{{A,B}/.sol,solUncert}]

Und wir schauen wie nahe wir am Wahrheitswert vorbei sind

{2,-1}-values

Look! This is perfect! Und so würde man das auch in echt tun. Das ist die richtigste Rangehensweise mit den Randbedingungen bisher bei orthogonal distance regressions, ODR für lineare Funktionen!

Als Sanity-Check können wir mal andere Implementierungen heranziehen. Das Wolfram Function Repository kennt zum Beispiel die OrthogonalLineFit Funktion.

OrthogonalLineFit=ResourceFunction["OrthogonalLineFit"]
Information[OrthogonalLineFit]

Der Output dieser Funktion ist leider ein bisschen doof mit einem Ankerpunkt p0 und einem Steigungsvektor und beides innerhalb einer InfiniteLine.
Aber das können wir schnell benutzen und umrechnen um auf unser Ergebnis zu kommen:

{p0Vec,v0Vec}=List@@OrthogonalLineFit[data];
First[yPt/.Solve[p0Vec+t*v0Vec=={xPt,yPt},{yPt,t}]]/.{bParam_+aParam_*xPt->{aParam,bParam}}

Komischerweise scheint das nicht unseren Ergebnis zu entsprechen. Warum denn?
Wenn wir in die Implementierung von OrthogonalLineFit schauen sehen wir, eine SVD-Zerlegung der Daten.
Ohne großartig ins Detail zu gehen, nimmt eine solche Implementierung an, dass  ist. Inwieweit das sinnvoll
ist, ist sogar gut debatierbar aber wir können beweisen dass wir auch auf das Ergebniss kommen mit dieser Annahme.
Wir müssen eben nur  setzen in unseren Gleichungen:

singleTerm=FullSimplify[Log[PDF[MultinormalDistribution[{0,0},{{uu^2,0},{0,uu^2}}],{x-deltaPX,y-deltaPY}/.{px->x,py->y}]],x∈Reals&&y∈Reals&&A∈Reals&&B∈Reals&&uu>0];
expr=Total[(singleTerm/.{x->#1,y->#2})&@@@data]//Simplify;
sol=Last[NMaximize[{expr,uu>0},{A,B,uu}]]

easy peasy lemon squeezy

Bestimmung von hyperfeinstruktur Dipol und Quadrupolkonstanten durch LineareRegression von

Die Hyperfeinstrukturaufspaltung in einem Atom kann beschrieben werden mit:
⟨(Subscript[H, 2])⟩  ⟨(Subscript[H, 4])⟩
Diese Erwargunswerte können für ein State |n,L,J,I,F⟩ in einem Atom gelößt werden zu:

DipoleH[FF_,II_,JJ_]:=(FF(FF+1)-II(II+1)-JJ(JJ+1))/2 (*⟨(Subscript[H, 2])⟩*)
QuadrupoleH[FF_,II_,JJ_]:=(3/2*(FF(FF+1)-II(II+1)-JJ(JJ+1))((FF(FF+1)-II(II+1)-JJ(JJ+1))+1)-2II(II+1)JJ(JJ+1))/(4II(2II-1)JJ(2JJ-1)) (*⟨(Subscript[H, 4])⟩*)

In einem Experiment kann man nun ausschließlich Frequenzabstände zwischen den verschiedenen States messen.
Ein tatsächlich realer Datensatz für 87Rb in einem 6P3/2 state sei hier mal gegeben:

dataF1={{0,1,1,1,(Around[11.504, 0.029])},{0,1,0,2,(Around[25.612, 0.020])},{0,1,1,2,(Around[37.552, 0.013])},{0,1,2,2,(Around[63.132, 0.016])},{1,1,0,2,(Around[14.108, 0.031])},{1,1,1,2,(Around[26.048, 0.028])},{1,1,2,2,(Around[51.629, 0.029])},{0,2,1,2,(Around[11.940, 0.017])},{0,2,2,2,(Around[37.520, 0.020])},{1,2,2,2,(Around[25.581, 0.013])}};
dataF2={{1,1,1,2,(Around[26.90, 0.12])},{1,1,2,2,(Around[52.82, 0.12])},{1,1,1,3,(Around[70.57, 0.12])},{1,1,2,3,(Around[96.20, 0.12])},{1,1,3,3,(Around[139.82, 0.12])},{1,2,2,2,(Around[25.923, 0.021])},{1,2,1,3,(Around[43.675, 0.018])},{1,2,2,3,(Around[69.296, 0.016])},{1,2,3,3,(Around[112.920, 0.016])},{2,2,1,3,(Around[17.752, 0.016])},{2,2,2,3,(Around[43.373, 0.014])},{2,2,3,3,(Around[86.997, 0.014])},{1,3,2,3,(Around[25.621, 0.009])},{1,3,3,3,(Around[69.245, 0.009])},{2,3,3,3,(Around[43.624, 0.005])}};
data=dataF1~Join~dataF2

Diese Daten kommen jetzt in der Form . Es sind deshalb viermal die F-Quantenzahlen vertreten, da wir jeweils eine Frequenzdifferenz von zwei States  &  nehmen.
Aufgrunddessen wie gemessen wird können zu diesen States aber jeweils der Mittelwert von zwei States  &  beitragen. Warum: ist hier nicht wichtig.
D.h. vereinfacht gesagt, unsere Frequenzdifferenz ist gegeben mit .
Implementieren wir das mal in Mathematica (nebenbei, I=3/2, J=3/2):

dE[F_]:=Ahfs*DipoleH[F,3/2,3/2]+Bhfs*QuadrupoleH[F,3/2,3/2]
dNu=(dE[F1]+dE[F1P])/2-(dE[F2]+dE[F2P])/2;
expr=Expand[FullSimplify[dNu],Ahfs|Bhfs]

So und jetzt sehen wir was schönes. Wir wollen ja mit der Menge an Daten eh eine Regression durchführen um  und  zu bestimmen. Und die gehen linear in unser Problem ein!
D.h. wir können theoretisch eine LinearModelFit benutzen oder? Naja wir haben aber eine Funktion von vier Abhängigen ! Kann LinearModelFit das?
Und die Antwort ist: problemlos!
Wir müssen nurnoch die Basisfunktionen rausholen.

base2=DeleteCases[expr[[1]],Ahfs]
base4=DeleteCases[expr[[2]],Bhfs]

Ab jetzt sollten wir einfach fitten können:

fit=LinearModelFit[data/.{Around[v_,__]:>v},{base2,base4},{F2,F2P,F1,F1P},IncludeConstantBasis->False,Weights->(1/data[[All,5,2]]^2)];
fit["ParameterTable"]
Around@@@Transpose[{fit["BestFitParameters"],fit["ParameterErrors"]}]

Und das sind bis heute competitive state of the art Konstanten für diesen Zustand:

Glaser, C., Karlewski, F., Kluge, J., Grimmel, J., Kaiser, M., Günther, A., Hattermann, H., Krutzik, M. and Fortágh, J., 2020. Absolute frequency measurement of rubidium 5S−6P transitions. Physical Review A, 102(1), p.012804.

Angleichen von  driftenden Spektren

Viele von euch werden in ihrer BA oder MA mit Spektren zu tun bekommen.
D.h. eine Observable die typischerweise von der Frequenz oder Wellenlänge
von eintreffenden Licht/Strahlung/etc. abhängt. Wenn diese Spektren sehr
verrauscht sind, ist es oft schwer gute Werte zu fitten.

Machen wir mal ein Beispiel dreier Cauchyverteilungen ohne Hintergrund.

makeExampleSpectrum[dv_]:=Module[{v},Table[{v-dv,PDF[CauchyDistribution[-0.3,0.2],v]+PDF[CauchyDistribution[1.5,0.3],v]+PDF[CauchyDistribution[0.8,0.4],v]+RandomVariate[NormalDistribution[0,0.1]]},{v,-1+dv,3+dv,0.01}]]
data=makeExampleSpectrum[RandomReal[{-0.5,0.5}]];
ListPlot[data]

Können wir direkt fitten. Als Wort der Warnung hier, wenn ihr sowas gut konvergieren lassen wollt, dann setzt inital die Mittelwerte  ungefähr in die richtige Größenordnung, und die Breiten  sehr dünn.

fit=NonlinearModelFit[data,A*PDF[CauchyDistribution[x1,b1],x]+B*PDF[CauchyDistribution[x2,b2],x]+CC*PDF[CauchyDistribution[x3,b3],x],{A,{x1,-0.5},{b1,0.1},B,{x2,0.5},{b2,0.1},CC,{x3,1.3},{b3,0.1}},x];
g1=Show[ListPlot[data],Plot[fit[x],{x,Min[data[[All,1]]],Max[data[[All,1]]]},PlotStyle->Black]]
fit["ParameterTable"]

Das merken wir uns. Aber was wenn wir besser Präzision für unsere Werte haben wollen.
Naja wir können das Spektrum ja problemlos mehrmals messen (Das macht ja hoffentlich ein Programm für euch).

allData=Table[makeExampleSpectrum[RandomReal[{-0.3,0.3}]],{i,1,15}];
ListLinePlot[allData]

Mhh, ja das ist natürlich die Gefahr in die wir laufen. Einfach so Mittelwert bilden ist nicht drinn wenn die Spektren driften zwischen den Messungen.
Wir müssen die erstmal angleichen.
Nehmen wir mal Spektrum eins und zwei:

data1=allData[[1]];
data2=allData[[2]];
ListLinePlot[{data1,data2}]

Wie können wir die den angleichen? Da gibt es etwas namens Korrelation. Das kennen wir schon. Wir gucken also wie sehr das eine Spektrum sich verändert,
wenn es das andere tut.

Correlation[data1[[All,2]],data2[[All,2]]]

Jetzt könnten wir ja sagen, wenn die Daten perfekt übereinander lägen, dann müsste doch die Korrelation maximal werden. Weil dann macht das eine Spektrum
genau dem anderen nach. Also die Korrelation sollte gegen eins streben. Wir können das herausfinden, indem die Spektren gegeneinander rotieren und in jedem
Schritt die Korrelation berechnen. Dieser Vorgang nennt sich Kreuzkorrelation.
(Es gibt btw. effizientere Varianten die Kreuzkorrelation zu berechnen; die FFT ruft. Aber hier nicht notwendig.)

expr=Table[{dm,Correlation[data1[[All,2]],RotateRight[data2[[All,2]],dm]]},{dm,-50,50}];
ListLinePlot[expr,FrameLabel->Thread[Style[{"Rotation in Datenpunkten","Korrelation"},14]]]
MaximalBy[expr,Last]

Aha, wenn wir also den zweiten Datensatz um 26 Datenpunkte nach Rechts verschieben, wird die Korrelation maximal mit ~0.947.
Schauen wir mal:

ListLinePlot[{data1,Transpose[{data1[[All,1]],RotateRight[data2[[All,2]],26]}]}]

Cool! Ja, genau das wollten wir. Dann richten wir mal alle Datensätze am ersten aus.

referenceData=allData[[1]]; (*Wir referenzieren alles an den ersten Datensatz*)
shiftedData=Table[ (*für jeden anderen Datensatz machen wir folgendes:*)
expr=Table[{dm,Correlation[referenceData[[All,2]],RotateRight[currentDataSet[[All,2]],dm]]},{dm,-50,50}]; (*Die Kreuzkorrelation in einem ±50 Datenpunkte Fenster*)
shifting=MaximalBy[expr,Last][[1,1]]; (*wir finden den Rotationsindex mit maximaler Korrelation*)
{shifting,Transpose[{currentDataSet[[All,1]],RotateRight[currentDataSet[[All,2]],shifting]}]} (*wir speicher ab {roationsindex, neuer-rotierter-datensatz}*)
,{currentDataSet,allData[[2;;-1]]}];

allShifts=MinMax[shiftedData[[All,1]]]; (*was ist die minimale und maximale rotation*)
shiftedData={referenceData}~Join~shiftedData[[All,2]]; (*fügen wir den den Referenzdatenset zu den rotierten hinzu*)
shiftedData=(#[[Last[allShifts];;-Abs[First[allShifts]]]])&/@shiftedData; (*finally, schneiden wir an den Seiten alle Datenpunkte weg die an der Rotation direkt beteiligt waren. Diese sind jetzt nämlich beidseitig gemischt*)

ListLinePlot[shiftedData]

Cool. Finally, verbinden wir alle Datepunkte einzeln mit MeanAround.

data=Transpose[{shiftedData[[1,All,1]],MeanAround/@Transpose[shiftedData[[All,All,2]]]}];
ListPlot[data]

Yes! Dieses Datenset ist viel genauer als die rohen Daten. Also Fitten wir das mal jetzt:

fit2=NonlinearModelFit[data/.{Around[v_,__]:>v},A*PDF[CauchyDistribution[x1,b1],x]+B*PDF[CauchyDistribution[x2,b2],x]+CC*PDF[CauchyDistribution[x3,b3],x],{A,{x1,-0.5},{b1,0.1},B,{x2,0.5},{b2,0.1},CC,{x3,1.3},{b3,0.1}},x,Weights->(1/data[[All,2,2]]^2)];
g2=Show[ListPlot[data],Plot[fit2[x],{x,Min[data[[All,1]]],Max[data[[All,1]]]},PlotStyle->Black]]
Quiet@fit2["ParameterTable"]

Und diese Daten sind natürlich deutlich genauer als die alten!

Grid[{{"Fit mit einem einzigen Datenset","Fit mit allen Datensets zusammen"},{Quiet@fit["ParameterTable"],Quiet@fit2["ParameterTable"]},{g1,g2}},Frame->All]

Denkt allerdings daran, dass  jetzt alle drei um einen beliebigen Wert  gleich verschoben sind. Falls ihr allerdings nur  haben wollt z.B., dann ist das ja egal.

Vergleich von Unabhängigen Messwerten

Man nehme mal folgende Situation:

Ich habe einen Wert einmal gemessen und einmal theoretisch a priori ausgerechnet.
ODER 
Ich habe einen Wert über zwei verschiedene, unabhängige, nicht korrelierte Experimente ausgerechnet.

Beide Werte haben natürliche eine Unsicherheit. Wie prüfen wir ob diese beiden übereinstimmen?
Nehmen wir mal wieder ein Beispiel aus meiner Doktorarbeit. 
Der AC-Stark Energie-shift im Rubidium Atom zwei-photon Übergang  5(Subscript[S, 1/2]) ⟶ 5(Subscript[D, 5/2]) kann sowohl gemessen, als auch durch Theorie ausgerechnet werden.
Dabei bekommen wir die Werte:

dNuIGemessen=Around[96.32114159355, 7.7056913274840015] (*gemessener Wert*)
dNuITheorie=Around[88.80330528173381, 0.34793452017587506] (*Theoretischer Wert aus Störungstheorie*)

Wie wahrscheinlich ist es jetzt, das diese beiden Werte den selben Wahrheitswert teilen?
Dazu müssen wir uns eines statistischen Test bedienen. Den sogenannten Z-Test.
Berechnet wird dabei ein Z-Wert wiefolgt:

Dann kann man im Prinzip die Aussagen treffen: Z<2 die Verteilungen stimmen überein, Z>2 stimmen nicht überein. Also berechnen wir mal:

zWert=(dNuIGemessen["Value"]-dNuITheorie["Value"])/Sqrt[dNuIGemessen["Error"]^2+dNuITheorie["Error"]^2]

Und jetzt könnten wir uns ein Loch in den Bauch freuen und fertig sein. Aber so einfach machen wir uns es noch nicht und wandeln diesen Wert
in einen sogenannten p-Wert. Dieser gibt so etwas ähnliches an, ob die Nullhypothese (die Gegenhypothese zu unseren Testkriterium) abgelehnt werden kann oder nicht.
(Die Nullhypothese die wir uns hier stellen ist: der Mittelwert der Subtraktion unserer Distributionen ist NICHT null und somit sind diese Verteilungen verschieden.)
Bei p<0.05 sagt man, dass die Nullhypothese nicht abgelehnt werden kann.
Wir hoffen also auf p>0.05. Oder eher p>>0.05.
Berechnet werden kann dieser Wert einfach mit der CDF einer Standard Normalverteilung:

2*CDF[NormalDistribution[0,1],-Abs[zWert]]

Year. Nullhypothese abgelehnt. D.h. wir können sagen, diese beiden Werte stimmen überein.
Der genaue Wert bei dem das ganze übrigens umschlägt hängt von dem Signifikanzniveau α ab.
Wie bereits erwähnt wird üblicherweise α=0.05 gewählt.

alpha=0.05;
InverseCDF[NormalDistribution[],1-alpha/2]

Nebenbei, können wir mit unseren vorherigen Wissen jetzt auch die Sicherheit in Standardabweichungen darstellen.

InverseErf[1-CDF[NormalDistribution[0,1],-Abs[zWert]]]*Sqrt[2]

Ich bin mir also zu 1.39σ sicher, dass die beiden übereinstimmen. Naja. Ginge besser.