\documentclass[12pt,a4paper,twoside]{article}
\usepackage[a4paper,twoside]{geometry}
\usepackage{html}
\usepackage{amsmath}
\usepackage{harvard}
\usepackage[german]{babel}
\usepackage[latin1]{inputenc}
%begin{latexonly}
\usepackage{graphicx}
\usepackage{threeparttable}
\usepackage{footmisc}
\usepackage{endnotes}
\let\footnote=\endnote
\renewcommand{\notesname}{Anmerkungen}
%end{latexonly}

% ensure \bm is defined if not latest LaTeX
%begin{latexonly}
\providecommand{\bm}[1]{\mathbf{#1}}
%end{latexonly}
\begin{imagesonly}
  \providecommand{\bm}[1]{\mathbf{#1}}
\end{imagesonly}


\newcommand*{\mean}[1]{\bar{#1}}

\begin{htmlonly}
  \newcommand*{\footref}[1]{}%
\end{htmlonly}

\hyphenation{Like-li-hood}

\title{Rechnet QtiPlot falsch ?\\ 
\large oder\\ Einige Anmerkungen zur Methode der linearen Regression}

\author{P.~Schäfer\\
Institut für Physik, Humboldt-Universität zu Berlin, \\
Newtonstraße 15, 12489 Berlin}

\date{\today}

\begin{document}

\maketitle

\section{Das Problem}
\label{intro}

Bei der Durchsicht der von den Studenten im Physikalischen
Grundpraktikum angefertigten Protokolle fällt auf, dass die dort mit
der Methode der linearen Regression erhaltenen Ergebnisse zum Teil
erheblich von den mit \citeasnoun{mathematica} berechneten
Kontrollwerten abweichen. 

Das Problem soll an einem einfachen Beispiel aus dem Einführungskurs
verdeutlicht werden. Beim Versuch F4 (siehe hierzu:
\citeasnoun{anleitungEP} Seite 9...11) wird eine Feder schrittweise 
mit Massestücken von ca. 50 g belastet. Die Position $x$ einer an der
Feder befestigten Marke wird mit Hilfe einer Spiegelskala gemessen (erste
Messreihe). Nach Erreichen der Höchstbelastung mit 8 Massestücken wird die
Feder wieder schrittweise entlastet und dabei die Position der Marke erneut
abgelesen (zweite Messreihe). Die Unsicherheit der so bestimmten
Positionen $x$ wird sowohl durch den systematischen Restfehler der
verwendeten Spiegelskala als auch durch Fehler beim Ablesen von der
Spiegelskala bestimmt. Der Mittelwert der Massen $\overline{m}$ aller
im Versuch verwendeten Massestücke beträgt 50.22\,g bei einer
Standardabweichung von 0,36\,g. Auch dieser Wert beeinflusst die
Streuung der gemessenen Positionen. 

Aus dem Anstieg der Geraden
\begin{equation}
  x = x_0 - \frac{\overline{m} g}{k} i
\end{equation}
kann dann die gesuchte Federkonstante bestimmt werden. Hierbei
bedeuten $x_0$ die unbekannte Position der verwendeten Markierung bei
unbelasteter Feder, $g$ die Erdbeschleunigung, $k$ die gesuchte
Federkonstante und $i$ die Anzahl der Massestücke, mit denen die Feder
belastet wird. Das negative Vorzeichen resultiert aus der Tatsache,
dass sich der Nullpunkt der verwendeten Spiegelskala unten befindet,
und daher die Werte von $x$ mit zunehmender Belastung abnehmen. Der
Wert $x_0$ für $i = 0$ kann nicht abgelesen werden, da die Achsen der
unbelasteten Federn  mehr oder weniger stark gebogen
sind. Eine typische Messdatenreihe, aufgenommen am Messplatz 2 am 7.1.2015,
ist in Tabelle \ref{tableMessdaten} aufgelistet. 

\begin{latexonly}
\begin{table}[h]
  \begin{center}
    \begin{tabular}{|c|c|c|}
      \hline
      $i$ & Belastung & Entlastung\\
      \hline
      1&	186.0 & 186.0 \\
      2&	175.5 & 176.0 \\
      3&	165.0 & 165.5 \\
      4&	155.0 & 156.0 \\
      5&	145.5 & 146.0 \\
      6&	135.5 & 136.0 \\
      7&	126.0 & 125.0 \\
      8&	116.0 & 115.5 \\
      \hline
    \end{tabular}
  \end{center}
  \caption{Beispiel einer Messdatenreihe des Versuchs F4}
  \label{tableMessdaten}
\end{table}
\end{latexonly}

\begin{htmlonly}
\begin{minipage}{13cm}
  \begin{center}
  \begin{tabular}{|c|c|c|}
    \hline
    $i$ & Belastung & Entlastung\\
    \hline
    1&	186.0 & 186.0 \\
    2&	175.5 & 176.0 \\
    3&	165.0 & 165.5 \\
    4&	155.0 & 156.0 \\
    5&	145.5 & 146.0 \\
    6&	135.5 & 136.0 \\
    7&	126.0 & 125.0 \\
    8&	116.0 & 115.5 \\
    \hline
  \end{tabular}
  \end{center}
  \vspace{1mm}
    {\bf Tabelle \ref{tableMessdaten}:} Beispiel einer Messdatenreihe des Versuchs F4
\end{minipage}
\vspace{2mm}
\end{htmlonly}

Betrachtet man diese Daten als 16 Wertepaare $i , x$, dann können sie
mit der Methode der linearen Regression, so wie sie zum Beispiel
in \citeasnoun{schenk} (Seiten 9 und 17...18), \citeasnoun{nollau}
Seiten 129...151 , \citeasnoun{fahrmeir} Seiten 59...63 und 90...95
oder im \citeasnoun{einfuehrung} Seite 42 Gleichungen 49, 50 und Seite
42 Gleichung für $s_y$ (nach Gleichung 51) beschrieben ist, ausgewertet
werden. In diesem Fall erhält man unabhängig von der verwendeten
Software \citeasnoun{mathematica} oder \citeasnoun{qtiplot} identische
Resultate (siehe Tabelle \ref{tableResultate}). 

Da mit dem systematischen Restfehler der verwendeten Spiegelskala 
$u_x = 200\mu m +5\;10^{-4}\;x$ (\citeasnoun{anleitungEP} Seite 11) eine
Abschätzung des Fehlers der Position $x$ vorliegt und dieser für die
einzelnen Messpunkte unterschiedlich ist, erscheint die Anwendung der
gewichteten linearen Regression sinnvoll. Diese Methode wird z.B.
von \citeasnoun{nollau} (Seiten 152...156) oder \citeasnoun{fahrmeir}
(Seiten 124...127) ausführlich beschrieben. Bei der instrumentellen
Gewichtung werden die Größen $c \frac{1}{s_i^2}$ als Gewichte
verwendet. $s_i^2$ ist dabei die Varianz, das Fehlerquadrat, der
jeweiligen Position $x_i$  und $c$ eine frei wählbare Konstante, die
oftmals gleich 1 gesetzt wird.

Die Ergebnisse der obigen Messreihe mit der Methode der instrumentell
gewichteten linearen Regression mit den beiden zuvor schon genutzten
Auswerteprogrammen sind ebenfalls in der Tabelle \ref{tableResultate}
aufgelistet.  Zusätzlich erfolgte eine Auswertung der
Daten mit den Formeln aus \citeasnoun{nollau} und
\citeasnoun{fahrmeir} sowie mit den im \citeasnoun{einfuehrung} auf
Seite 41,42 angegebenen Gleichungen 43 und 45. Letztere lässt sich aus
der von \citeasnoun{wolf} auf den Seite 117/118 hergeleiteten
Gleichung 9.20 für $\alpha = 2$ ableiten.

\begin{latexonly}
\begin{table}[h]
  \centering  
  \begin{threeparttable}
    \begin{tabular}{|c|c|c|c|c|}
      \hline
      \multicolumn{5}{|c|}{lineare Regression} \\ \hline
      Software & \multicolumn{2}{|c|}{$x_0$} & \multicolumn{2}{|c|}{Anstieg $A$}\\
      \cline{2-5}
      &  Wert  &  Fehler & Wert & Fehler \\ 
      \hline
      \citename{mathematica} & 195.723 & 0.225292 & -10.0149 & 0.0446145 \\
      \hline
      \citename{qtiplot}\tnote{$\dagger$} & 195.723 &  0.225292 & -10.0149 & 0.0446145\\
      \hline
      \citename{einfuehrung} & 195.723 &  0.225292 & -10.0149 & 0.0446145 \\
      \hline \hline
      \multicolumn{5}{|c|}{gewichtete lineare Regression (instrumentelle
        Gewichtung) }\\ \hline
      Software & \multicolumn{2}{|c|}{$x_0$} & \multicolumn{2}{|c|}{Anstieg $A$}\\
      \cline{2-5}
      &  Wert  &  Fehler & Wert & Fehler \\ 
      \hline
      \citename{mathematica} & 195.707 & 0.233646 & -10.0114 & 0.044768 \\
      \hline
      \citename{nollau}& 195.707 & 0.233646 & -10.0114 & 0.044768 \\
      \hline
      \citename{qtiplot}\tnote{$\dagger$} & 195.707 &  0.156635 & -10.0114 & 0.030012 \\
      \hline
      \citename{einfuehrung} & 195.707 & 0.156635 & -10.0114 & 0.030012 \\
      \hline
      \citename{mathematica}\tnote{$\ddagger$} & 195.707 & 0.156635 & -10.0114 & 0.030012 \\
      \hline
    \end{tabular}
    \vspace{1mm}
    \begin{tablenotes}\footnotesize 
    \item[$\dagger$] LinearFit mit und ohne Option \glqq{}Scale Errors with sqrt(Chi2/doF)\grqq{}
    \item[$\ddagger$] LinearModelFit mit Option \glqq{}VarianceEstimatorFunction $\rightarrow$ (1\&)\grqq{}
    \end{tablenotes}
  \end{threeparttable}
  \caption{Ergebnisse der Auswertung der Messdaten mit verschiedenen
    Methoden und verschiedener Software. Die zugehörigen Ausgaben der
    genutzten Software sind im Anhang enthalten.}
  \label{tableResultate}
\end{table}
\end{latexonly}

\begin{htmlonly}
\begin{minipage}{20cm}
  \begin{center}
  \begin{tabular}{|c|c|c|c|c|}
    \hline
    \multicolumn{5}{|c|}{lineare Regression} \\ \hline
    Software & \multicolumn{2}{|c|}{$x_0$} & \multicolumn{2}{|c|}{Anstieg $A$}\\
    \cline{2-5}
    &  Wert  &  Fehler & Wert & Fehler \\ 
    \hline
    \citename{mathematica} & 195.723 & 0.225292 & -10.0149 & 0.0446145 \\
    \hline
    \citename{qtiplot} $\dagger$ & 195.723 &  0.225292 & -10.0149 & 0.0446145\\
    \hline
    \citename{einfuehrung} & 195.723 &  0.225292 & -10.0149 & 0.0446145 \\
    \hline \hline
    \multicolumn{5}{|c|}{gewichtete lineare Regression (instrumentelle
      Gewichtung) }\\ \hline
    Software & \multicolumn{2}{|c|}{$x_0$} & \multicolumn{2}{|c|}{Anstieg $A$}\\
    \cline{2-5}
    &  Wert  &  Fehler & Wert & Fehler \\ 
    \hline
    \citename{mathematica} & 195.707 & 0.233646 & -10.0114 & 0.044768 \\
    \hline
    \citename{nollau}& 195.707 & 0.233646 & -10.0114 & 0.044768 \\
    \hline
    \citename{qtiplot} $\dagger$ & 195.707 &  0.156635 & -10.0114 & 0.030012 \\
    \hline
    \citename{einfuehrung} & 195.707 & 0.156635 & -10.0114 & 0.030012 \\
    \hline
    \citename{mathematica} $\ddagger$ & 195.707 & 0.156635 & -10.0114 & 0.030012 \\
    \hline
  \end{tabular}
  \vspace{1mm}
  $\dagger$ {\footnotesize LinearFit mit und ohne Option \glqq{}Scale Errors with sqrt(Chi2/doF)}\grqq{}
  \newline
  $\ddagger$ {\footnotesize LinearModelFit mit Option \grqq{}VarianceEstimatorFunction $\rightarrow$ (1\&)}\grqq{}
  \end{center}
{\bf Tabelle \ref{tableResultate}:} Ergebnisse der Auswertung der Messdaten mit verschiedenen
Methoden und verschiedener Software. Die zugehörigen Ausgaben
der genutzten Software sind im Anhang 
(\htmladdnormallink{QtiPlot}{./Ergebnisse-QtiPlot.pdf} und
\htmladdnormallink{Mathematica}{./Ergebnisse-Mathematica.pdf}) enthalten. 
\end{minipage}
\vspace{2mm}
\end{htmlonly}

Sieht man sich die Ergebnisse genauer an, fällt sofort auf, dass die von
\citename{mathematica} und 
\citename{qtiplot} berechneten Standardabweichungen der Parameter
$x_0$ und $A$ nicht übereinstimmen. Demgegenüber liefern sowohl
\citename{mathematica} und die Berechnung mit den Gleichung von
\citeaffixed{fahrmeir,nollau}{z.B.} als auch \citename{qtiplot} und
die Gleichungen aus dem \citename{einfuehrung} jeweils zueinander
passende Ergebnisse. Verwendet man bei der Auswertung mit \citename{mathematica} die
Option \glqq{}VarianceEstimatorFunction $\rightarrow$ (1\&)\grqq{}, so
stimmt das Ergebnis sowohl mit den von \citename{qtiplot} ausgegebenen
Werten als auch mit den mit Hilfe der Gleichungen aus dem
\citename{einfuehrung} berechneten Werten überein. Damit deutet sich
hier ein grundsätzliches Problem an, dass über die eingangs
aufgeworfene Frage, ob die Software QtiPlot richtig rechnet, weit
hinausgeht und auf die Frage hinausläuft, warum stehen die Formeln aus
dem \citename{einfuehrung} im Widerspruch zu gängigen Lehrbüchern
\citeaffixed{fahrmeir,nollau}{z.B.}. 

Ein erster Hinweis darauf findet sich bei \citeasnoun{james} auf Seite 3.   
\begin{itemize}
\item[]\glqq{}Unfortunately, statisticans do not agree on basic
  principles. They can crudely be divided into two schools: Bayesian
  and frequentists (or classical). 

  The Bayesian approach is closer to everyday reasoning, where
  probability is interpreted as a degree of belief that something
  will happen, or that a parameter will have a given value.

  The frequentist approach is closer to scientific reasoning, where 
  probability means the relative frequency of something
  happening. This make it more objective, since it can be determined
  independently of the observer, but restricts its application to
  repeatable phenomina.
\grqq{}
\end{itemize}

Liegt die Ursache für die unterschiedlichen Formeln und Ergebnisse in
den Unterschieden zwischen diesen beiden Richtungen begründet? Dazu
sollen die jeweiligen Aussagen näher betrachtet werden, wobei eine
Beschränkung auf das allgemeine linear Modell erfolgt. Zu den
Grundlagen, zu Herleitungen und Beweisen sowie zu weiterführenden
Fragen sei auf die angegebene Literatur verwiesen. 

\newpage 
\section{Das allgemeine lineare Modell}
Den Ausgangspunkt bildet eine Reihe von $n$ Messwerten, 
\begin{displaymath}
  \bm{y}=\left( \begin{array}{l}
                y_1 \\
                y_2 \\
                \vdots \\
                y_n \\
                \end{array} \right)
\end{displaymath}
die in Abhängigkeit von $r$ nicht stochastischen Einflussgrößen
erhalten wurden. Die $n$ Vektoren $\bm{x}_i = (x_1,x_2,\ldots,x_r)_i$
müssen dabei nicht
alle unterschiedliche Werte besitzen.  Die Abhängigkeit der Messwerte
von den Einflussgrößen kann über 
\begin{equation}
y_i = f(\bm{x}_i,\bm{\theta}) +\epsilon_i \hspace{15mm} i = 1,\dots,n 
\end{equation}
beschrieben werden. Die Spaltenvektoren
\begin{displaymath}
  \bm{\theta}=\left( \begin{array}{l}
                \theta_1 \\
                \theta_2\\
                \vdots \\
                \theta_p \\ 
                \end{array} \right) \hspace{5mm} \mathrm{und} \hspace {5mm}
  \bm{\epsilon}=\left( \begin{array}{l}
                \epsilon_1 \\
                \epsilon_2\\
                \vdots \\
                \epsilon_n \\
                \end{array} \right) \hspace{5mm}
\end{displaymath}
enthalten die $p$ unbekannten Parameter des Modells beziehungsweise 
eine unabhängige Zufallsgröße, die die Abweichungen der Messwerte
$y_i$ von den aus dem Modell berechneten Werten beschreibt. Für den
Erwartungswert $\bm{\mathrm{E}}$ und  die Varianz
$\bm{\mathrm{Var}}$ dieser Größe soll gelten:   
\begin{equation}
  \bm{\mathrm{E}}[\bm{\epsilon}] =  0 \hspace{5mm} \mbox{und} \hspace{5mm}
  \bm{\mathrm{Var}}[\bm{\epsilon}] = \sigma^2 
\end{equation}
wobei $\sigma^2$ im allgemeinen unbekannt ist. Weitere Annahmen über
die Verteilung der $\epsilon_i$, insbesondere die Voraussetzung einer
Normalverteilung, sind an dieser Stelle nicht notwendig.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ Zu den notwendigen Voraussetzungen für ein lineares Modell
  gibt es unterschiedliche Auffassungen:\\ 
  \citeasnoun{martin} Seite 145: Insbesondere ist die Annahme nicht
  notwendig, dass die Fehler der Messwerte normalverteilt sind. \\
  \citeasnoun{james} Seite 186:  Es werden keine Annahmen über die
  Verteilung der Daten $\bm{y}$ gemacht. \\
  \citeasnoun{nollau} setzt schon bei der Herleitung des linearen
  Modells für die $\epsilon_i$ unabhängige normalverteilte
  Zufallsgrößen voraus.  
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Damit gilt für den Erwartungswert
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ Im Unterschied zu \citeasnoun{martin}\,(Seite 145),
  \citeasnoun{fahrmeir}\,(Seite 62) und \citeasnoun{nollau}\,(Seite
  134) verwendet \citeasnoun{james} in diesem Zusammenhang im Kapitel
  8.4 (Seite 183\,ff) anstelle von $\bm{\mathrm{E}}[y_i]$ den
  Ausdruck $\bm{\mathrm{E}}[y_i,\bm{\theta}]$. 
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vom $y_i$:
\begin{equation}
  \bm{\mathrm{E}}[y_i] = f(\bm{x}_i,\bm{\theta}) \hspace{15mm} i = 1,\ldots,n
\end{equation}
Wenn sich die Funktion $f(\bm{x}_i,\bm{\theta})$ als Linearkombination von $p$
beliebigen, nicht notwendigerweise linearen Funktionen $f_j(\bm{x}_i)$
darstellen lässt, kann man dafür schreiben: 
\begin{equation}
 \bm{\mathrm{E}}[y_i] =\sum_{j=1}^p \theta_j f_j(\bm{x}_i) \hspace{15mm} i = 1,\ldots,n 
\end{equation}
In diesem Fall liegt ein in den Parametern $\bm{\theta}$ lineares Modell
vor. Aus diesem lassen sich verschiedene Spezialfälle
ableiten. Wenn zum Beispiel die Messwerte nur von einer Einflussgröße
$\bm{x}_i = (x_1)_i$
linear abhängen, kann  $f_1(\bm{x}_i) = 1$ sowie  $f_2(\bm{x}_i) = x_i$ festgesetzt
werden. In diesem Fall ergibt sich daraus das Modell der linearen Regression
\begin{equation}
  \bm{\mathrm{E}}[y_i] = \theta_1 + \theta_2 x_i \hspace{15mm} i = 1,\ldots,n 
\end{equation}
mit den unbekannten Parametern $\theta_1$ und $\theta_2$. In ähnlicher 
Weise lassen sich auch das Modell für eine Gerade durch den
Koordinatenursprung oder für die Anpassung mit einem Polynom
herleiten. 

Durch die Wahl der $p$ Funktionen $f_j(\bm{x}_i)$ und der $n$
Messstellen $\bm{x}_i$ sind die Werte
\begin{displaymath}
  a_{i,j} =f_j(\bm{x}_i) \hspace{15mm} i = 1,\ldots,n \hspace{5mm} j = 1,\ldots,p
\end{displaymath}
fest vorgeben. Diese $a_{i,j}$ bilden die $n \times p$ Design-Matrix
\begin{displaymath}
  \bm{A}=\left( \begin{array}{llll}
      a_{1,1} & a_{1,2} & \cdots & a_{1,p}\\
      a_{2,1} & a_{2,2} & \cdots & a_{2,p}\\
      \vdots & \vdots & \ddots & \vdots \\
      a_{n,1} & a_{n,2} & \cdots & a_{n,p}\\
                \end{array} \right)
\end{displaymath}
des linearen Modells, das damit als einfache Matrizengleichnung geschrieben
werden kann.
\begin{equation}
  \bm{y} = \bm{A}\;\bm{\theta} + \bm{\epsilon}
  \hspace{5mm} \mathrm{beziehungsweise} \hspace {5mm}
  \bm{\mathrm{E}}[\bm{y}] = \bm{A}\;\bm{\theta}
\label{SimpleLinMod}
\end{equation}Bisher wurde vorausgesetzt, dass alle Messergebnisse die gleiche
Genauigkeit besitzen. Dies ist oftmals nicht der Fall und die
einzelnen Messwerte sollen daher mit unterschiedlicher Gewichtung in das
Modell einfließen. Zu diesem Zweck kann zu jedem Messwert ein entsprechendes
Gewicht $p_i$ festgelegt werden. Damit lassen sich die relativen
Unterscheide in der Genauigkeit der einzelnen Messungen in einfacher
Weise beschreiben. Diese Werte bilden die Gewichtsmatrix
$\bm{P}=diag(p_1,p_2,\ldots,p_n)$. Setzt man für 
\begin{eqnarray}
  \bm{y}^* & = & \bm{P}^{1/2} \bm{y} \nonumber \\
  \bm{A}^* & = & \bm{P}^{1/2} \bm{A} \nonumber \\
  \bm{\epsilon}^* & = & \bm{P}^{1/2} \bm{\epsilon} 
\end{eqnarray}
so ergibt 
\begin{equation}
  \bm{y}^*  =    \bm{A}^* \bm{\theta} + \bm{\epsilon}^*
\end{equation}
wieder ein lineares Modell (Gleichung \ref{SimpleLinMod}).

Liegen Informationen über die Varianz $\sigma_i$ der Fehler der
einzelnen Messwerte $y_i$ vor, so können diese zur Gewichtung
benutzt werden. Dazu geht man von der Vorstellung aus, dass  die
einzelnen Messwerte $y_i$ jeweils Mittelwerte einer Sichprobe des
Umfanges $n_i$ repräsentieren. Die Varianz dieser Mittelwerte folgt
aus der Varianz $\sigma$ der zugrunde liegenden Grundgesamtheit
entsprechend $\sigma_i^2 = \sigma^2/n_i$. Setzt man als Gewichte
den Umfang der jeweiligen Stichproben ein, wie dies auch
von \citeasnoun{fahrmeir} für gruppierte Daten beschrieben wird,
so folgt daraus:     
\begin{displaymath}
  p_i  =  n_i = \frac{\sigma^2}{\sigma_i^2}
\end{displaymath}
Da oftmals die Varianz $\sigma^2$ der Grundgesamtheit vorab nicht
bekannt ist, werden die Gewichte nur proportional zu $n_i$ angesetzt.  
\begin{equation}
    p_i \propto \frac{\sigma^2}{\sigma_i^2}\hspace{10mm} p_i = c\;\frac{1}{\sigma_i^2} 
\label{Gewicht}
\end{equation}
wobei c in diesem Fall, eine  willkürlich festgelegte Konstante
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ \citeasnoun{bevington} wählt im Kapitel 11 für die
  Gewichte (Gleichung (11.2) Seite 194, Gleichung (11.32) Seite
  203 sowie auf Seite 215) den Ausdruck: 
 \begin{displaymath}
   p_i = \frac{1/\sigma_i^2}{(1/n)\sum(1/\sigma_i^2)}
   = \frac{n}{\sum(1/\sigma_i^2)}\; \frac{1}{\sigma_i^2}
 \end{displaymath}
  Daraus folgt für die Konstante c:
  \begin{displaymath}
    c = \frac{n}{\sum(1/\sigma_i^2)}
  \end{displaymath}
  Dies entspricht einer Normalisierung der Gewichtsfaktoren auf deren Mittelwert.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ist, die, wie später zu sehen sein wird, keinen Einfluss auf die
Schätzungen von $\widehat{\bm{\theta}}$
(Gleichung \ref{SolutionAllgLinMod}) und
$\widehat{\bm{\mathrm{Cov}}[\widehat{\bm{\theta}}]}$ (Gleichung \ref{CovAllgLinMod})
hat. 
 
Die Diagonalmatrix $diag(\sigma_1^2,\sigma_2^2\ldots\sigma_n^2)$
entspricht der $n\times n$ Kovarianzmatrix $\bm{\mathrm{Cov}}[\bm{\epsilon}]$ der
bisher als statistisch unabhängig angenommenen Abweichungen
$\epsilon_i$ der Messwerte $y_i$ von den aus dem Modell berechneten
Werten. Diese Annahme ist im weiteren nicht mehr erforderlich. Die
Kovarianzmatrix der Messabweichungen $\bm{\Sigma}$ ist allgemein
gegeben durch:    
\begin{displaymath} 
\bm{\Sigma} =  \bm{\mathrm{Cov}}[\bm{\epsilon}]=\left( \begin{array}{llll}
      \sigma_{1,1} & \sigma_{1,2} & \cdots & \sigma_{1,n}\\
      \sigma_{2,1} & \sigma_{2,2} & \cdots & \sigma_{2,n}\\
      \vdots & \vdots & \ddots & \vdots \\
      \sigma_{n,1} & \sigma_{n,2} & \cdots & \sigma_{n,n}\\
                \end{array} \right)
\end{displaymath}
wobei
\begin{displaymath}
 \sigma_{i,i} = \sigma_i^2= \bm{\mathrm{Var}}[\epsilon_i] \hspace{3mm} \mathrm{und} \hspace{3mm}
 \sigma_{i,j} = \bm{\mathrm{Cov}}[\epsilon_i,\epsilon _j] \hspace{3mm} \mathrm{mit} \hspace{3mm} i \neq j
\end{displaymath}
sind. Die Kovarianzmatrix $\bm{\Sigma}$ ist symmetrisch und positiv
semidefinit. Ist sie positiv definit, kann ihre Inverse, auch als
Präzisionsmatrix bezeichnet, zur Festlegung der Gewichtmatrix genutzt
werden.
\begin{displaymath}
  \bm{P} = \sigma^2\;\bm{\Sigma}^{-1}
\end{displaymath}
$\sigma^2$ ist hierbei der Varianzfaktor, oftmals auch als Varianz der
Gewichtseinheit bezeichnet, der festlegt welchem Messwert das Gewicht
Eins gegeben wird. Bei unbekanntem $\sigma$ wird die Gewichtsmatrix als
proportional zur Präzisionsmatrix angesetzt.
\begin{equation}
  \bm{P} \propto \bm{\Sigma}^{-1} \hspace{10mm} \bm{P} =  c\;\bm{\Sigma}^{-1}
\label{Gewichtsmatrix}
\end{equation}
Für die Konstante $c$ gilt das oben gesagte.

Setzt man für die Gewichtsmatrix die Einheitsmatrix ein,
$\bm{P} =\bm{I}$, so ergibt sich daraus das einfache lineare Modell
(Gleichung \ref{SimpleLinMod}). Eine getrennte Behandlung dieses
Spezialfalles ist im weiteren nicht erforderlich.  

\subsection{Klassische Statistik }

Gesucht sind die Schätzwerte für $\widehat{\bm{\theta}}$
,$\widehat{\sigma^2}$ und  $\widehat{\bm{\mathrm{Cov}}[\widehat{\bm{\theta}}]}$ , die
den wahren Werten der unbekannten Parameter möglichst nahe kommen. Bei
der Methode der kleinsten Quadrate bestimmt das Minimum der Summe der
quadrierten Abweichung 
\begin{eqnarray*}
  Q(\bm{\theta}) & = & \sum_{i=1}^n\sum_{j=1}^n (y_i -\bm{\mathrm{E}}[y_i])p_{i,j} (y_j -\bm{\mathrm{E}}[y_j])\\
                 & = & \bm{\epsilon}^T\bm{P}\bm{\epsilon} \\ 
                 & = & (\bm{y} - \bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta})\\
   & = & \bm{y}^T\bm{P}\bm{y}-2\bm{\theta}^T\bm{A}^T\bm{P}\bm{y}+\bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta}
\end{eqnarray*}
die Schätzwerte der Parameter. Aus
\begin{displaymath}
  \frac{\partial\,Q(\bm{\theta})}{\partial\,\bm{\theta}} = - 2
  \bm{A}^T\bm{P}\bm{y} + 2 \bm{A}^T\bm{P}\bm{A}\;\bm{\theta} = \bm{0}
\end{displaymath}
folgt unmittelbar:
\begin{equation}
  \widehat{\bm{\theta}} = (\bm{A}^T\bm{P}\bm{A})^{-1}\,\bm{A}^T\bm{P}\bm{y} 
\label{SolutionAllgLinMod}
\end{equation}
Der Schätzwert für die Varianz der Messabweichungen kann aus dem Minimum von
$Q(\bm{\theta})$ bestimmt werden. 
\begin{equation}
\widehat{\sigma^2} = \frac{1}{n-p}\,Q(\bm{\theta})_{min}
     = \frac{1}{n-p}\,(\bm{y}-\bm{A}\;\widehat{\bm{\theta}})^T\bm{P}(\bm{y}-\bm{A}\;\widehat{\bm{\theta}})
\label{sigmaKlassisch}
\end{equation}
Für die Kovarianzmatrix von $\widehat{\bm{\theta}}$ ergibt sich damit:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ Nach \citeasnoun{fahrmeir}( Seite 463 ) gilt für die
  Linearkombination eines Zufallsvektors $\bm{X}$ mit 
  einer geeignet dimensionierten Matrix $\bm{A}$ und einem geeignet
  dimensionierten Vektor $\bm{b}$:
  \begin{eqnarray*}
    \bm{\mathrm{E}}[\bm{A}\bm{X} +\bm{b}]&=&\bm{A}\,\bm{\mathrm{E}}[\bm{X}] +\bm{b}\\
    \bm{\mathrm{Cov}}[\bm{A}\bm{X} +\bm{b}]&=&\bm{A}\,\bm{\mathrm{Cov}}[\bm{X}]\,\bm{A}^T
  \end{eqnarray*}
  Damit ergibt sich mit $\bm{\Sigma}=\bm{\mathrm{Cov}}[\bm{y}]$ :
  \begin{eqnarray*}
    \bm{\mathrm{Cov}}[\widehat{\bm{\theta}}]
    &=&\bm{\mathrm{Cov}}[(\bm{A}^T\bm{P}\bm{A})^{-1}\,\bm{A}^T\bm{P}\,\bm{y}]\\
    &=&(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\,
    \bm{\Sigma}\,
    \left((\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right)^T
  \end{eqnarray*}
  Mit $\bm{\Sigma}=\widehat{\sigma^2}\,\bm{P}^{-1} $ folgt daraus:
  \begin{eqnarray*}
    \bm{\mathrm{Cov}}[\widehat{\bm{\theta}}]
    &=&\widehat{\sigma^2}\,(\bm{A}^T\bm{P}\;\bm{A})^{-1}(\bm{A}^T\bm{P}\bm{A})(\bm{A}^T\bm{P}\;\bm{A})^{-1}\\
    &=&\widehat{\sigma^2}\,(\bm{A}^T\bm{P}\;\bm{A})^{-1}
  \end{eqnarray*}
  da die Matrizen $\bm{P}$ und $\bm{A}^T\bm{P}\;\bm{A}$ symmetrisch
  sind.
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{equation}
  \widehat{\bm{\mathrm{Cov}}[\widehat{\bm{\theta}}]} =
  \widehat{\sigma^2}(\bm{A}^T\bm{P} \bm{A})^{-1} 
\label{CovAllgLinMod}
\end{equation}

Bisher wurden keine Annahmen über die Verteilung der Messdaten
zu Grunde gelegt. Die optimalen Eigenschaften der Methode der kleinsten
Quadrate folgen einzig aus der Linearität des Modells. Wird für die
Messabweichungen $\bm{\epsilon}$ eine Normalverteilung angenommen, führt die
Maximum-Likelihood-Schätzung von $\widehat{\bm{\theta}}$
beziehungsweise die Restringierte Maximum-Likelihood-Schätzung 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ Die Maximum-Likelihood-Schätzung
  \begin{displaymath}
    \widehat{\sigma^2}=\frac{1}{n}\,\widehat{\bm{\epsilon}}^T\bm{P}\widehat{\bm{\epsilon}}
  \end{displaymath}  
  ist nicht erwartungstreu. Im allgemeinen wird sie deswegen durch die
  Restringierte Maximum-Likelihood-Schätzung
  \begin{displaymath}
  \widehat{\sigma^2}=\frac{1}{n-p}\,\widehat{\bm{\epsilon}}^T\bm{P}\widehat{\bm{\epsilon}}
  \end{displaymath}
  ersetzt. (\citeasnoun{fahrmeir} Seite 94, \citeasnoun{wakefield}
  Seite 44 ff.)
  \label{RMLS}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
von $\widehat{\sigma^2}$ zu dem gleichen Ergebnis. Unter dieser
Voraussetzung sind sind auch die geschätzten Werte 
$\widehat{\bm{\theta}} \sim \mathrm{N}(\bm{\theta},\sigma^2(\bm{A}^T\bm{P}\bm{A})^{-1})$ 
normalverteilt. %% \citeasnoun{fahrmeir} Seite 103 

Das Konfidenzintervall (Fehlerintervall) eines geschätzten Parameters
$\widehat{\theta_j}$ wird durch die Wahrscheinlichkeit $P$ bestimmt,
mit der dieser im Intervall $d_ {min} < \theta_j < d_{max}$ liegt
\begin{displaymath}
P(d_ {min} < \theta_j < d_{max}) = 1 - \alpha
\end{displaymath}
wobei $1 - \alpha$ das Signifikanzniveau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\footnote{
  Im Praktikum wird meist mit $(1 \sigma)$-Fehlerintervallen
  gearbeitet. Das zugehörige Signifikanzniveau ist $68.2689\%$,
  woraus $\alpha = 31.7311\%$ folgt.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
angibt. Zur Festlegung der Intervallgrenzen $d$ wird der Test der
Hypothese $H_0\mathrm{:}\;\theta_j = d$ verwendet. Dieser basiert auf
der Teststatistik 
\begin{displaymath}
  t = \frac{\widehat{\theta_j} - d}{s_j}
\end{displaymath}
mit
\begin{displaymath}
  s_j = \widehat{\bm{\mathrm{Var}}[\widehat{\theta_j}]}^{1/2}
  = \left(\widehat{\bm{\mathrm{Cov}}[\widehat{\bm{\theta}}]}_{j,j}\right)^{1/2}
  = \left(\widehat{\sigma^2}(\bm{A}^T\bm{P} \bm{A})^{-1}_{j,j}\right)^{1/2}
\end{displaymath}
Die Hypothese $H_0$ wird abgelehnt, wenn sich
\begin{displaymath}
  |t| > t_{n-p}(1-\alpha/2)
\end{displaymath}
ergibt. Hierbei bedeutet $t_{n-p}(1-\alpha/2)$ das
$(1-\alpha/2)$-Quantil der Standard t-Verteilung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ Standard t-Verteilung $x \sim t_n$
  (z.B.\citeasnoun{fahrmeir} Seite 461)
  \begin{displaymath}
    p(x) =
    \frac{\Gamma(\frac{n+1}{2})}{\sqrt{n\pi}\Gamma(\frac{n}{2})}
    \left(1+\frac{x^2}{n}\right)^{-\frac{n+1}{2}}
  \end{displaymath}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mit $n-p$ Freiheitsgraden . Die Grenzen des Konfidenzintervalls zum
Signifikanzniveau  $(1-\alpha)$ werden somit bestimmt durch
\begin{equation}
  d_{min,max}  =  \widehat{\theta_j} \mp t_{n-p}(1-\alpha/2)\; s_j
\label{klassKonfInt}
\end{equation}

Ähnliche Gleichungen kann man in vielen Lehrbüchern zur klassischen
Statistik, z.B in 
\citeasnoun{nollau} (152\,ff),
\citeasnoun{james} (183\,ff),
\citeasnoun{fahrmeir} (124\,ff), 
\citeasnoun{martin} (145\,ff) oder  
\citeasnoun{wakefield} (214\,ff)
finden, wobei teilweise anstelle der Gewichtsmatrix $\bm{P}$ direkt
die inverse Kovarianzmatrix $\bm{\Sigma}^{-1}$ benutzt oder das erste
Element des Parametervektors $\bm{\theta}$ mit $\theta_0$ bezeichnet
wird.

%\newpage
\subsection{Bayes Statistik}

Die Bayes Statistik verfolgt einen grundsätzlich anderen Ansatz als
die Klassische Statistik. Die Klassische Statistik versucht nur aus
den gegebenen Daten, den Messwerten, z.B. mit der Methode der
kleinsten Quadrate die wahrscheinlichsten Werte der unbekannten
Parameter zu schätzen, deren Erwartungswerte als fest angenommenen
werden. Demgegenüber geht die Bayes-Statistik davon aus, das auch die
Parameter selbst Zufallsgrößen sein können und bestimmt  nicht nur die
Parameterwerte, für die die Datenwahrscheinlichkeit am größten ist,
sondern auch deren Wahrscheinlichkeitsverteilung, aus der sich die  
Unsicherheit über den wahren Wert der Parameter ableiten lässt. Dazu
arbeitet die Bayes-Statistik mit bedingten Wahrscheinlichkeiten. Mit
Hilfe des Bayes-Theorem 
\begin{equation}
  p(\bm{\theta}\mid\bm{y}) = \frac{p(\bm{y}\mid\bm{\theta})\;p(\bm{\theta})}{p(\bm{y})}
\label{BayesTheorem}
\end{equation}
wird die gesuchte bedingte Wahrscheinlichkeitsverteilung
$p(\bm{\theta}\mid\bm{y})$ der Parameter, deren
Posteriori-Wahrscheinlichkeit bestimmt.   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ Zu bedingten Wahrscheinlichkeiten und zum Bayes-Theorem
  siehe z.B. \citeasnoun{koch} (Kapitel 2) oder \citeasnoun{wakefield}
  Kapitel 3}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Neben der Likelihood-Funktion $p(\bm{y}\mid\bm{\theta})$ sind
Vorkenntnisse über die Wahrscheinlichkeitsverteilung der
unbekannten Parameter $p(\bm{\theta})$, über die sogenannte
Priori-Verteilung notwendig. Je nach Festlegung beinhaltet diese eine
mehr oder weniger große Subjektivität. Auch wenn die statistischen
Aussagen ausschließlich aus der  Posteriori-Verteilung abgeleitet
werden, bleibt diese subjektive Komponente erhalten. Im Gegensatz dazu
arbeitet die klassische Statistik nur mit der Likelihood-Funktion. 
Da der Ausdruck ${p(\bm{y})}$ nur eine Normierungskonstante
%% (\citeasnoun{Koch} Seite 14, \citeasnoun{wakefield} Seite 86)
ist, reicht es aus, Gleichung \ref{BayesTheorem} in der Form  
\begin{equation}
  p(\bm{\theta}\mid\bm{y}) \propto p(\bm{y}\mid\bm{\theta})\times p(\bm{\theta})
\label{PropBayesTheorem}
\end{equation}
zu verwenden.

In der Schreibweise der Bayes-Statistik lautet das allgemeine lineare Modell:
\begin{eqnarray}
  \bm{y}  & = & \bm{A}\;\bm{\theta} + \bm{\epsilon}
  \hspace{5mm} \mathrm{beziehungsweise} \hspace {5mm}
  \bm{\mathrm{E}}[\bm{y}\mid\bm{\theta}] = \bm{A}\;\bm{\theta}
  \nonumber \\
  & \mathrm{mit} &
  \bm{\mathrm{E}}[\bm{\epsilon}\mid\bm{\theta}] = 0 \nonumber \\
  & \mathrm{und} &
  \bm{\mathrm{Cov}}[\bm{\epsilon}\mid\bm{\theta},\sigma^2] =
  \bm{\mathrm{Cov}}[\bm{y}\mid\sigma^2] = \bm{\Sigma} = \sigma^2\bm{P}^{-1}
\label{BayesAllgLinMod}
\end{eqnarray}
wobei im weiteren für die Verteilung der Messabweichungen eine mehrdimensionale
Normalverteilung vorausgesetzt wird.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{Die in den Gleichungen \ref{Gewicht} und
  \ref{Gewichtsmatrix} genutzte, frei wählbare Konstante $c$
  beeinflusst den  Varianzfaktor $\sigma^2$, der auch als Varianz 
  der Gewichtseinheit bezeichnet wird.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Unter diesen Voraussetzungen ergibt sich die Likelihood-Funktion zu:
\begin{displaymath}
 p(\bm{y}\mid\bm{\theta},\sigma) = \frac{1}{(2\pi)^{n/2}|(\sigma^2\bm{P}^{-1})|^{1/2}}\;
 \exp\left(-\frac{1}{2\sigma^2}(\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta})\right)
 \end{displaymath}

%% Argumentation nach \citeasnoun{wakefield} Seite 221 ff (ohne Gewichtsmatrix)}
Da über die Parameter $\bm{\theta}$ und $\sigma^2$ keine
Vorinformationen vorliegen, wird im weiteren von einem nichtinformativen
Priori-Verteilung ausgegangen. Für $\theta$ wird angenommen, das alle
möglichen Werte gleich wahrscheinlich sind. 
\begin{displaymath}
p(\bm{\theta}) \propto 1 \hspace{5mm} \mbox{für} \hspace{5mm}
-\infty < \theta_i < \infty\;,\; i=1,\ldots,n 
\end{displaymath}
Da $\sigma^2$ nur Werte im Intervall $0 < \sigma^2 < \infty$ annehmen
kann, wird von einer Gleichverteilung für $\ln{\sigma}$
ausgegangen. Dies führt auf die nichtinformative Priori-Verteilung für
die Varianz: \footref{BayesKoch}
\begin{displaymath}
p(\sigma)  \propto \sigma^{-2} \hspace{5mm} \mbox{für} \hspace{5mm} 0 < \sigma < \infty
\end{displaymath}
%% \citeasnoun{koch} Seite 57 Gleichung 2.218
Somit ergibt sich für die Priori-Verteilung $p(\bm{\theta},\sigma^2)$ :
\begin{equation}
p(\bm{\theta},\sigma^2)=p(\bm{\theta})\times p(\sigma^2)\propto\sigma^{-2}
\label{PriorWakefield}
\end{equation}
%% \citeasnoun{wakefield} Seite 221 Gleichung 5.42
Diese Verteilungsfunktion ist nicht normierbar und wird deswegen auch als
uneigentliche Priori-Verteilung bezeichnet. Da die verwendete
Likelihood-Funktion normierbar ist, folgt aus der Verknüpfung nach dem
Bayes-Theorem eine Posteriori-Verteilung, die normierbar ist.
\begin{displaymath}
  p(\bm{\theta},\sigma\mid\bm{y}) \propto (\sigma^2)^{-(n+2)/2}
  \exp\left(-\frac{1}{2\sigma^2}(\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}\,(\bm{y}-\bm{A}\;\bm{\theta})\right)
\end{displaymath}
%%  \citeasnoun{wakefield} Seite 221 Gleichung 5.43
Der Exponenten lässt sich geeignet umformen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% \citeasnoun{wakefield} Seite 44\,ff,  Gleichung 2.23, Kapitel 2
\footnote{
  \label{substitution}
  Die Umformung erfolgt in mehreren Schritten. Dabei wird ausgenutzt,
  dass die Gewichtsmatrix $\bm{P}$ entsprechend ihrer Definition
  symmetrisch ist. Daraus ergibt sich , dass auch die Matrix
  $(\bm{A}^T\bm{P}\;\bm{A})$ symmetrisch ist.
  \begin{latexonly}
  \begin{eqnarray*}
    \lefteqn{\mbox{entwickle:}}\\
    \lefteqn{(\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta})}\\
    & = & \bm{y}^T\bm{P}\;\bm{y}
    - \bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y}
    - \bm{y}^T\bm{P}\;\bm{A}\;\bm{\theta}
    + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta} \\[3mm]
    \lefteqn{\mbox{alle  Terme sind Skalare}\,\rightarrow\,
      \bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y} =
      \bm{y}^T\bm{P}\;\bm{A}\;\bm{\theta} }\\[3mm]
    & = &  \bm{y}^T\bm{P}\;\bm{y} - 2 \bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y}
    + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta}\\[3mm]
    \lefteqn{\mbox{mit:}\hspace{5mm}2\bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y} =
      2\bm{\theta}^T(\bm{A}^T\bm{P}\bm{A})(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}}\\
    \lefteqn{\mbox{und:}\hspace{5mm}\bm{\mu_0}=(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}}\\[3mm]
    & = & \bm{y}^T\bm{P}\;\bm{y} - 2\bm{\theta}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}
    + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta}\\[3mm]
    \lefteqn{\mbox{setze:}}\\
    \lefteqn{\bm{y}^T\bm{P}\,\bm{y}}\\
    & = & \bm{y}^T\bm{P}\,\bm{y}
    -\bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}
    +\bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}\\
    & = &
    \bm{y}^T\left(\bm{P}-\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right)\bm{y} 
    +\bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}
  \end{eqnarray*}
  \begin{eqnarray*}
    \lefteqn{\mbox{aus}}\\
    \lefteqn{(\bm{y}-\bm{A}\;\bm{\mu_0})^T\bm{P}\;(\bm{y}-\bm{A}\;\bm{\mu_0})}\\
    & = &
      (\bm{y}-\bm{A}\;(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y})^T\bm{P}\;
    (\bm{y}-\bm{A}(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y})\\
    & = &
    \bm{y}^T\left(\bm{I}-\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\right)\;\bm{P}\;
    \left(\bm{I}-\bm{A}(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\right)\bm{y}\\
    & = & \bm{y}^T\left(
    \bm{P}-2\;\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right.\\
    & &\left.\hspace{12mm} +\;\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{A}
    (\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}
    \right)\bm{y}\\
    & = & \bm{y}^T\left(
    \bm{P}-\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right)\bm{y}\\
    \lefteqn{\mbox{und}}\\
    \lefteqn{\bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}}\\
    & = &
    \bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}(\bm{A}^T\bm{P}\bm{A})(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}\\
    & = & \bm{\mu_0}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}\\[3mm]
    \lefteqn{\mbox{folgt mit:}\hspace{5mm} (n-p)\;s^2 = (\bm{y}-\bm{A}\;\bm{\mu_0})^T\bm{P}\;(\bm{y}-\bm{A}\;\bm{\mu_0})}\\
    \lefteqn{\bm{y}^T\bm{P}\,\bm{y}}\\
    & = &
    (n-p)\;s^2 + \bm{\mu_0}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}\\
    \lefteqn{\mbox{damit ergibt sich:}}\\
    \lefteqn{(\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta})}\\
    & = & (n-p)\;s^2 + \bm{\mu_0}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}
    - 2\bm{\theta}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0} + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta}\\
    & = & (n-p)\;s^2 +(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})
  \end{eqnarray*}
\end{latexonly}
\begin{htmlonly}
  \begin{eqnarray*}
    \mbox{entwickle:} & & \\
    (\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta})
    & = & \bm{y}^T\bm{P}\;\bm{y}
    - \bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y}
    - \bm{y}^T\bm{P}\;\bm{A}\;\bm{\theta}
    + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta} \\[3mm]
    \mbox{alle  Terme sind Skalare} & \rightarrow &
      \bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y} =
      \bm{y}^T\bm{P}\;\bm{A}\;\bm{\theta} \\[3mm]
    & = &  \bm{y}^T\bm{P}\;\bm{y} - 2 \bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y}
    + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta}\\[3mm]
    \mbox{mit:}& &2\bm{\theta}^T\bm{A}^T\bm{P}\;\bm{y} = 
      2\bm{\theta}^T(\bm{A}^T\bm{P}\bm{A})(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}\\
    \mbox{und:}& &\bm{\mu_0}=(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}\\[3mm]
    & = & \bm{y}^T\bm{P}\;\bm{y} - 2\bm{\theta}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}
    + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta}\\[5mm]
    \mbox{setze:}& & \\
    \bm{y}^T\bm{P}\,\bm{y}
    & = & \bm{y}^T\bm{P}\,\bm{y}
    -\bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}
    +\bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}\\
    & = &
    \bm{y}^T\left(\bm{P}-\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right)\bm{y} 
    +\bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}\\[3mm]
    \mbox{aus}& & \\
    (\bm{y}-\bm{A}\;\bm{\mu_0})^T\bm{P}\;(\bm{y}-\bm{A}\;\bm{\mu_0})
    & = &
      (\bm{y}-\bm{A}\;(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y})^T\bm{P}\;
    (\bm{y}-\bm{A}(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y})\\
    & = &
    \bm{y}^T\left(\bm{I}-\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\right)\;\bm{P}\;
    \left(\bm{I}-\bm{A}(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\right)\bm{y}\\
    & = & \bm{y}^T\left(
    \bm{P}-2\;\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right.\\
    &   &\left.\hspace{12mm} +\;\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{A}
    (\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}
    \right)\bm{y}\\
    & = & \bm{y}^T\left(
    \bm{P}-\bm{P}\;\bm{A}(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right)\bm{y}\\
    \mbox{und}& & \\
    \bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}
    & = &
    \bm{y}^T\bm{P}\bm{A}\;(\bm{A}^T\bm{P}\bm{A})^{-1}(\bm{A}^T\bm{P}\bm{A})(\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\;\bm{y}\\
    & = & \bm{\mu_0}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}\\[3mm]
    \mbox{folgt mit:}& & (n-p)\;s^2 = (\bm{y}-\bm{A}\;\bm{\mu_0})^T\bm{P}\;(\bm{y}-\bm{A}\;\bm{\mu_0})\\
    \bm{y}^T\bm{P}\,\bm{y}
    & = & (n-p)\;s^2 + \bm{\mu_0}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}\\
    \mbox{damit ergibt sich:}& & \\
    (\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta})
    & = & (n-p)\;s^2 + \bm{\mu_0}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0}
    - 2\bm{\theta}^T(\bm{A}^T\bm{P}\bm{A})\bm{\mu_0} + \bm{\theta}^T\bm{A}^T\bm{P}\bm{A}\;\bm{\theta}\\
    & = & (n-p)\;s^2 +(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})
  \end{eqnarray*}
\end{htmlonly}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{displaymath}
  (\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta}) =
  (n-p)\;s^2 +(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})
\end{displaymath}
Darin bedeuten 
\begin{displaymath}
  \bm{\mu_0} = (\bm{A}^T\bm{P}\bm{A})^{-1}\bm{A}^T\bm{P}\bm{y}
\end{displaymath}
und 
\begin{displaymath}
  s^2 =\frac{1}{n-p}\;(\bm{y}-\bm{A}\;\bm{\mu_0})^T\bm{P}\,(\bm{y}-\bm{A}\;\bm{\mu_0})
\end{displaymath}
Damit ergibt sich eine Normal-inverse Gammaverteilung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{
  Normal-inverse Gamma-Verteilung
  $\bm{\theta},\sigma^2 \sim \mathrm{NIG}(\bm{m},\,\bm{\Sigma},\,a,\,b)$\\
  (z.B. \citeasnoun{fahrmeir2013} Seite 652):
\begin{eqnarray*}
  p(\bm{\theta},\sigma^2) 
  & = &\frac{1}{(2\pi)^{p/2}|\bm{\Sigma}|^{1/2}}\frac{b^a}{\Gamma(a)}\\
  &   &\frac{1}{(\sigma^2)^{a+1}}\exp\left( -\frac{b}{\sigma^2} -
  \frac{1}{2} (\bm{\theta} - \bm{m})^T \bm{\Sigma}^{-1}(\bm{\theta} - \bm{m})\right)
\end{eqnarray*}
Für die praktische Anwendung wird oftmals
$\bm{\Sigma}=\sigma^2\,\bm{M}$ gesetzt und\\ 
$ \bm{\theta},\sigma^2 \sim \mathrm{NIG}(\bm{m},\,\bm{M},\,a,\,b)$
verwendet (z.B. \citeasnoun{fahrmeir2013} Seite 227). 
\begin{eqnarray*} 
  p(\bm{\theta},\sigma^2) 
  & = &\frac{1}{(2\pi)^{p/2}|\bm{M}|^{1/2}}\frac{b^a}{\Gamma(a)}\\
  &   &\frac{1}{(\sigma^2)^{p/2+a+1}}\exp\left( -\frac{b}{\sigma^2} -
  \frac{1}{2\sigma^2} (\bm{\theta} - \bm{m})^T \bm{M}^{-1}(\bm{\theta}
  - \bm{m})\right)
\end{eqnarray*}  
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 $\mathrm{NIG}(\bm{m},\,\bm{M},\,a,\,b)$ mit den Parametern
\begin{eqnarray*}
 \bm{m} & \rightarrow & \bm{\mu_0}\,\mathrm{,}\hspace{12mm}
 \bm{M} \rightarrow  (\bm{A}^T\bm{P}\;\bm{A})^{-1}\\
 a & \rightarrow & \frac{n-p}{2}\,\mathrm{,}\hspace{5mm}
 b   \rightarrow  \frac{n-p}{2}\;s^2
\end{eqnarray*}
und der Dichte
\begin{latexonly}
\begin{eqnarray}
  \label{PostWakefield}
  \lefteqn{p(\bm{\theta},\sigma\mid\bm{y}) \propto} \\
  & &(\sigma^2)^{-(n+2)/2} \exp \left\{ -\frac{1}{2\sigma^2}
  \left[s^2(n-p)+(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})\right]\right\}
  \nonumber
\end{eqnarray}
\end{latexonly}
\begin{htmlonly}
  \begin{equation}
    \label{PostWakefield}
    p(\bm{\theta},\sigma\mid\bm{y}) \propto(\sigma^2)^{-(n+2)/2} \exp \left\{ -\frac{1}{2\sigma^2}
  \left[s^2(n-p)+(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})\right]\right\}
  \end{equation}
\end{htmlonly}  
Um daraus die Verteilung $p(\bm{\theta}\mid\bm{y})$ zu erhalten,
wird über $\sigma^2$ integriert.
\begin{latexonly}
\begin{eqnarray*}
  \lefteqn{p(\bm{\theta}\mid\bm{y}) \propto }\\
  & & \int (\sigma^2)^{-(n+2)/2} \exp \left\{ -\frac{1}{2\sigma^2} \left[s^2(n-p)+(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})\right] \right\} d\sigma^2
\end{eqnarray*}
\end{latexonly}
\begin{htmlonly}
  \begin{displaymath}
    p(\bm{\theta}\mid\bm{y}) \propto \int (\sigma^2)^{-(n+2)/2} \exp \left\{ -\frac{1}{2\sigma^2} \left[s^2(n-p)+(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})\right] \right\} d\sigma^2
  \end{displaymath}
\end{htmlonly}  
Der Integrand ist der Kern einer inversen Gamma-Verteilung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{
  \label{InvGa}
  Inverse Gamma-Verteilung $x \sim \mathrm{InvGa}(a,b)$\\
  (z.B. \citeasnoun{fahrmeir} Seite 461):
  \begin{eqnarray*}
    p(x) & = & \frac{b^a}{\Gamma(a)}x^{-(a+1)}\exp(-b/x),\hspace{3mm}x>0\\
    \mbox{mit:}   &    & \bm{\mathrm{E}}[x] = \frac{b}{a-1}\\
    \mbox{es gilt:}   & \rightarrow &\int p(x)\;d x = 1 \\
    \mbox{daraus folgt:} & \rightarrow &\int x^{-(a+1)}\exp(-b/x)\;d x = \Gamma(a)\;b^{-a}
  \end{eqnarray*}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\mathrm{InvGa}(a,b)$ für $\sigma^2$ mit den Parametern
\begin{eqnarray*}
  a & \rightarrow & n/2\\
  b & \rightarrow & \frac{s^2(n-p)+(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})}{2}
\end{eqnarray*}
Aus dem Normierungsparameter der inversen Gamma-Verteilung folgt: 
\begin{eqnarray*}
  p(\bm{\theta}\mid\bm{y}) & \propto &\Gamma(\frac{n}{2})\left[ \frac
    {s^2(n-p)   
  +(\bm{\theta}-\bm{\mu_0})^T\bm{A}^T\bm{P}\bm{A}\,(\bm{\theta}-\bm{\mu_0})}{2}
    \right]^{-\frac{n}{2}}\\
  & \propto & \left[1 + \frac{(\bm{\theta}-\bm{\mu_0})^T\left(s^2(\bm{A}^T\bm{P}\bm{A})^{-1}\right)^{-1}(\bm{\theta}-\bm{\mu_0})}{n-p}\right]^{\frac{n-p+p}{2}}
\end{eqnarray*}
Dieser Ausdruck entspricht dem Kern einer $p$-dimensionalen
t-Verteilung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{$p$-dimensionale  t-Verteilung
  $\bm{x} \sim \mathrm{T}_p(\bm{\mu},\bm{\Sigma},d)$\\  
  (z.B. \citeasnoun{fahrmeir} Seite 467)
%%  \citeasnoun{fahrmeir2013} Seite 652,
%%  \citeasnoun{wakefield} Seite 660, \citeasnoun{koch} Seite 55)  
  \begin{displaymath}
    p(\bm{x}) = 
     \frac{\Gamma[(d+p)/2]}{\Gamma[d/2](d\;\pi)^{p/2}}|\bm{\Sigma}|^{-1/2}\left[1+\frac{(\bm{x}-\bm{\mu})^T\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu})}{d}\right]^{-(d+p)/2}
  \end{displaymath}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newline
$\mathrm{T}_p(\bm{\mu},\bm{D},d)$ mit :
\begin{eqnarray*}
  \bm{\mu} & \rightarrow & \bm{\mu_0} \\
  \bm{D} &\rightarrow & s^2(\bm{A}^T\bm{P}\bm{A})^{-1}\\
  d & \rightarrow & n-p
\end{eqnarray*}
Für den Parametervektor $\bm{\theta}$ ergibt sich damit die $p$-dimensionalen t-Verteilung
\begin{equation}
\bm{\theta} \mid \bm{y} \sim
\mathrm{T}_p\left(\bm{\mu_0},\,s^2(\bm{A}^T\bm{P}\bm{A})^{-1},\,n-p\right)
\label{PosteriorRV}
\end{equation}
als Posteriori-Randverteilung. Zu dem gleichen Ergebnis gelangt man auch mit anderen
nichtinformativen Priori-Verteilungen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{\label{BayesKoch}
Zu der gleichen Verteilung $\bm{\theta} \mid \bm{y}$ gelangt
\citeasnoun{koch}( Seite 111 ff) mit einer anderen
Priori-Verteilung. Er substituiert den Varianzfaktor $\sigma^2$
mit dem Gewichtsparameter  $ \tau = 1 / \sigma^2 $ und benutzt die
Priori-Verteilungen %% Gleichungen 4.113; S.57 Gleichungen 2.216($\theta$), 2.218($\sigma$), 2.220($\tau$)
\begin{eqnarray*}
p(\bm{\theta}) \propto & const & \hspace{5mm} \mbox{für} \hspace{5mm}
-\infty < \theta_i < \infty\;,\; i=1,\ldots,n \\
p(\tau)  \propto & \tau^{-1} & \hspace{5mm} \mbox{für} \hspace{5mm} 0 < \tau < \infty\\
p(\bm{\theta},\tau) & = &  p(\bm{\theta})\times p(\tau)  \propto  \tau^{-1}  
\end{eqnarray*}
Daraus ergibt sich die Posteriori-Verteilung zu: %% Gleichung 4.115
\begin{displaymath}
p(\bm{\theta},\tau\mid\bm{y}) \propto \tau^{n/2-1}
\exp\left(-\frac{\tau}{2}(\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}(\bm{y}-\bm{A}\;\bm{\theta})\right)
\end{displaymath}
Der Exponent wird in der gleichen Weise umgeformt\footref{substitution}.
Damit lässt sich die Posteriori-Verteilung als Normal-Gammaverteilung
darstellen. Aus dieser folgt dann die mehrdimensionale t-Verteilung:
\begin{displaymath}
\bm{\theta} \mid \bm{y} \sim \mathrm{T}_p\left(\bm{\mu_0},\,s^2(\bm{A}^T\bm{P}\bm{A})^{-1}\,,\,n-p\right)
\end{displaymath}
als Posteriori-Randverteilung für $\bm{\theta}$.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textsuperscript{,}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{\label{bayesFahrmeir}
\citeasnoun{fahrmeir2013}(Seite 227\,ff.) verwendete als konjugierte
Priori-Verteilung eine Normal-inverse Gammaverteilung 
$\theta,\sigma^2 \sim \mathrm{NIG}(\bm{m},\,\bm{M},\,a,\,b)$. Aus
dieser folgt als Posteriori-Verteilung wieder eine Normal-inverse
Gammaverteilung. Mit den von \citeasnoun{fahrmeir2013} auf Seite 231
angegebenen Werten $\bm{m} = \bm{0}$, $\bm{M}^{-1}=\bm{0}$, $a=-p$ und
$b=0$ folgt als Priori-Verteilung:  
\begin{displaymath}
  p(\bm{\theta},\sigma^2) \propto \frac{1}{(\sigma^2)^{-p/2+1}} = (\sigma^2)^{p/2-1}
\end{displaymath}
die nicht mit der von \citeasnoun{fahrmeir2013}(Gleichung. 4.16)
angestrebten, nichtinformativen Priori-Verteilung
\begin{displaymath}
  p(\bm{\theta},\sigma^2) \propto \frac{1}{\sigma^2}
\end{displaymath}
übereinstimmt. Diese ist mit der Vorgabe $a=-p/2$ zu erreichen.
Daraus ergibt sich dann als Posteriori-Verteilung 
\begin{displaymath}
  p(\bm{\theta},\sigma \mid \bm{y}) \propto
  \frac{1}{(\sigma^2)^{n/2+1}}
  \exp\left(-\frac{1}{2\sigma^2}(\bm{y}-\bm{A}\;\bm{\theta})^T\bm{P}\,(\bm{y}-\bm{A}\;\bm{\theta})\right)
\end{displaymath}
die sich durch die Normal-inverse Gammaverteilung
$\mathrm{NIG}(\bm{m},\,\bm{M},\,a,\,b)$ mit
\begin{eqnarray*}
  \mbox{mit}\hspace{5mm} \bm{m} & \rightarrow & \bm{\mu_0}\,\mathrm{,}\hspace{12mm}
  \bm{M} \rightarrow  (\bm{A}^T\bm{P}\;\bm{A})^{-1}\\
  \mbox{und}\hspace{7mm} a & \rightarrow & \frac{n-p}{2}\,\mathrm{,}\hspace{5mm}
  b  \rightarrow  \frac{n-p}{2}\;s^2
\end{eqnarray*}
beschreiben lässt.
%% Seite 232 Gleichung 4.17
%% Jeffreys Priori oder Referenz Prioriverteilung $p(\bm{\mu},\sigma^2) \sim  1/\sigma^2$\\
%% Seite 668
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% James Seite 161 
% nutzt $\ln \sigma^2$ gleichverteilt über $0 < \sigma^2 < \infty$
% und kommt damit zu der gleichen Posteriori-Verteilung Gleichung 7.59
%%%%%%%%% Achtung Druckfehler: richtig $ \exp[-\frac{[\ldots]}{2 \theta_2}$ %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
. Aus dem Erwartungswert der mehrdimensionalen t-Verteilung
\begin{displaymath}
 \bm{\mathrm{E}}[\bm{\theta}\mid\bm{y}] = \bm{\mu_0}
\end{displaymath}   
ergibt sich die Bayesschätzung $\widehat{\bm{\theta}}_B$ zu
\begin{equation}
\widehat{\bm{\theta}}_B = (\bm{A}^T\bm{P}\bm{A})^{-1}\,\bm{A}^T\bm{P}\bm{y} 
\label{BayesSolutionAllgLinMod}
\end{equation}   
Diese stimmt mit dem Ergebnis der Klassischen Statistik (\ref{SolutionAllgLinMod}) überein.

Aus der Posteriori-Randverteilung (\ref{PosteriorRV}) lässt sich der
Konfidenzbereich für jeden Teilvektor des Parametervektor $\bm{\theta}$
und damit auch das Konfidenzintervall des einzelnen Parameters
$\theta_j$ ableiten. Es ergibt sich dann eine eindimensionale
t-Verteilung mit $n-p$ Freiheitsgraden, die sich mit der Substitution
\begin{displaymath}
  x = \frac{{\bm{\mu_0}}_j - \theta_j}{s_j}
\end{displaymath}
und dem entsprechenden Hauptdiagonalenelement der Dispersionsmatrix $\bm{D}$
\begin{displaymath}
  s_j = \left(s^2(\bm{A}^T\bm{P} \bm{A})^{-1}_{j,j}\right)^{1/2}
\end{displaymath}
auf eine Standard t-Verteilung mit $n-p$ Freiheitsgraden zurückführen
lässt. Die sich daraus ergebenden Konfidenzintervalle
\begin{equation}  
  {\bm{\mu_0}}_j - t_{n-p}(1-\alpha/2)\; s_j \le \theta_j \le
  {\bm{\mu_0}}_j + t_{n-p}(1-\alpha/2)\; s_j
\label{BayesKonfInt}
\end{equation}
stimmen mit denen der Klassischen Statistik (\ref{klassKonfInt}) überein.

Aus der Posteriori-Verteilung (\ref{PostWakefield}) folgt als
Randverteilung für die Varianz $\sigma^2$ eine inverse
Gammaverteilung\footref{InvGa} $\mathrm{InvGa}(a,b)$ 
%% koch Seite 112 
mit den Parametern
\begin{eqnarray*}
 a & \rightarrow & \frac{n-p}{2}\\
 b & \rightarrow & \frac{n-p}{2}\;s^2
\end{eqnarray*}
Aus dem Erwartungswert der inversen Gammaverteilung ergibt sich damit die
Bayesschätzung des Varianzfaktors $\widehat{\sigma^2}_B$ zu 
\begin{eqnarray*}
  \widehat{\sigma^2}_B & = & \bm{\mathrm{E}}[\sigma^2\mid\bm{y}]\\
  & = &  \frac{n-p}{n-p-2}\;s^2\\
  & = &  \frac{1}{n-p-2}\;(\bm{y}-\bm{A}\;\bm{\mu_0})^T\bm{P}\,(\bm{y}-\bm{A}\;\bm{\mu_0})
\end{eqnarray*}
Das gleiche Ergebnis kann auch auf anderen Weise gefunden werden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ \citeasnoun{wakefield}(Seite 222) kommt auf eine inverse
  Chi Quadrat Verteilung mit $n-p$ Freiheitsgraden
  \begin{displaymath}
    \sigma^2\mid\bm{y} \sim (n-p)\;s^2 \times\chi^{-2}_{n-p} 
  \end{displaymath}
  wobei zu beachten ist, dass \citeasnoun{wakefield} den
  Parametervektor $\bm{\theta}$ von $0 \ldots k = p_W-2 $ indiziert und
  dort in $p_W$, der Anzahl der unbekannten Parameter, $\sigma$
  enthalten ist (Seite 209).\\
  Der Erwartungswert der inversen Chi Quadrat Verteilung mit k
  Freiheitsgraden ist gegeben durch
  \begin{displaymath}
    \bm{\mathrm{E}}[\chi^{-2}_k] =\frac{1}{k-2}
  \end{displaymath}
  Daraus folgt:
  \begin{displaymath}
    \widehat{\sigma^2}_B = \frac{n-p}{n-p-2}\;s^2
  \end{displaymath}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
. Die Bayesschätzung  $\widehat{\sigma^2}_B$ ist um den Faktor
$(n-p)/(n-p-2)$ größer als die Schätzung $\widehat{\sigma^2}$ der
klassischen Statistik(\ref{sigmaKlassisch}). 


%% Aus der Dispersionsmatrix $\bm{\mathrm{N}}^{-1}$ ergibt sich die
%% Kovarianzmatrix $\bm{\mathrm{D}} = \nu/(nu-2)\bm{\mathrm{N}}^{-1}$ 
%% \citeasnoun{koch} (Seite 55 Gleichung: 2.208)
%%
%% \citeasnoun{Bättig} S.193
%% $\mathrm{likelihood} \propto \frac{1}{\sigma^n}\exp(-0.5 \chi^2$
%% mit $\chi^2 = \frac{\sum_{i=1}^n (y_i-\mu_y(x_i))^2}{\sigma^2}$

%% \citeasnoun{Bättig} S.194 $\sigma$ wird als Skalierungsparameter behandelt 
%%        S.172 minimale Vorinformation:
%% Lageparameter $p(\mu) \propto 1
%% Skalierungsparameter $p(1/mu) \propto 1/\mu$

\subsection{Fazit}

Die Klassische Statistik liefert im Rahmen des Allgemeinen Linearen
Modells für den Parametervektor und für die Fehlerintervalle das
gleiche Ergebnis wie die Bayes Statistik bei Verwendung einer
nicht-informativen Priori-Verteilung. 
Für große Werte von $n$ gilt das auch für den Varianzfaktor $\sigma^2$. 
Der Varianzfaktor, der sich aus der Klassischen Statistik ergibt,
stellt die untere Grenze für den Varianzfaktor dar, den man über den
entsprechenden Ansatz der Bayes-Statistik erhält.
%% Dies gilt nicht nur für den Fall, dass für die y-Werte eine
%% Gauß-Verteilung angenommen wird, sondern auch für den Fall einer
%% Poisson-Verteilung (Zählstatistik) \citeasnoun{barnett}. ?????

Ein Vergleich mit Tabelle \ref{tableResultate} zeigt, dass die mit
\citename{qtiplot} und den Formeln im Skript \citename{einfuehrung}
berechneten Werte für die Fehlerintervalle der Parameter signifikant
kleiner sind als die nach der Klassischen Statistik berechneten
Größen. Die grundlegend unterschiedliche Vorgehensweise der beiden 
Richtungen in der Statistik kann nicht die Ursache für die
aufgetretenen Differenzen sein.

\newpage
\section{Fehlerfortpflanzung}
Im \citename{einfuehrung} findet man auf Seite 42 die Aussage:
\glqq{} Die Standardabweichungen für die beiden Parameter ... erhält man
durch Anwendung des Fortpflanzungsgesetzes für die
Unsicherheiten\grqq{}. Als Quelle kann hier das im Vorwort erwähnte
Buch \citeasnoun{bevington} vermutet werden. Dort findet sich auf
Seite 109 die Herleitung der im \citename{einfuehrung} verwendeten
Formeln.  Eine ähnliche Aussage wie im \citename{einfuehrung} ist auch
bei \citeasnoun{wolf} auf Seite 118 zu finden, ohne dass dazu eine
Begründung oder eine Quellenangabe erfolgt. 

Was bedeutet die Anwendung der Fehlerfortpflanzung auf das Allgemeine
Lineare Modell? Dazu wird zuerst eine skalare Größe $\theta$ betrachtet,
die eine beliebige Funktion des $n$ dimensionalen Vektors $\bm{y}$
ist.
\begin{displaymath}
\theta \equiv \theta(\bm{y}) = \theta(y_1,y_2,\ldots,y_n)
\end{displaymath}
%% Fehlerfortpflanzungsgesetz allgemein \citeasnoun{martin} Seite 99 ff
%% %% siehe auch \citeasnoun{bevington} Seite 39 ff
Wenn die Fehler $\epsilon_i=y_i-\mean{y_i}$ mit $\mean{y_i} = \bm{\mathrm{E}}[y_i]$ klein sind, kann
die Taylorentwicklung von  $\theta(\bm{y})$ in der Umgebung des
Punktes $\bm{y}=\mean{\bm{y}}$ nach der ersten Ordnung
abgebrochen werden. 
\begin{equation}
  \theta(\bm{y}) =
  \theta\left(\mean{\bm{y}}\right)+\sum_{i=1}^n\left(y_i-\mean{y_i}\right)
  \left.\frac{\partial \theta(\bm{y})}{\partial y_i} \right|_{\bm{y}=\mean{\bm{y}}}
\label{ThetaTaylor}
\end{equation}
Für die Varianz von $\theta(\bm{y})$ gilt:
\begin{displaymath}
  \bm{\mathrm{Var}}[\theta(\bm{y})]=\bm{\mathrm{E}}\left[\Big(\theta(\bm{y})-\bm{\mathrm{E}}[\theta(\bm{y})]\Big)^2\right]
  \simeq \bm{\mathrm{E}}\left[\Big(\theta(\bm{y})-\theta(\mean{\bm{y}})\Big)^2\right]
\end{displaymath}
Wird die Taylorentwicklung \ref{ThetaTaylor} in diese Gleichung eingesetzt,
ergibt sich:
\begin{eqnarray*}
 \bm{\mathrm{Var}}[\theta(\bm{y})]&=&\bm{\mathrm{E}}\left[\left( \sum_{i=1}^n\left(y_i-\mean{y_i}\right)
  \left.\frac{\partial \theta(\bm{y})}{\partial y_i} \right|_{\bm{y}=\mean{\bm{y}}} \right)^2\right]\\
  &=&\sum_{i=1}^n\sum_{j=1}^n
  \left.\frac{\partial\theta(\bm{y})}{\partial y_i}\right|_{\bm{y}=\mean{\bm{y}}}
  \bm{\mathrm{E}}\left[\Big(y_i-\mean{y_i}\Big)\Big(y_j-\mean{y_j}\Big) \right]
  \left.\frac{\partial\theta(\bm{y})}{\partial y_j}\right|_{\bm{y}=\mean{\bm{y}}}
\end{eqnarray*}
Dabei sind die
%% mit \citeasnoun{martin} Seite 50 Gleichung 3.28
\begin{displaymath}
  \bm{\mathrm{E}}\left[\Big(y_i-\mean{y_i}\Big)\Big(y_j-\mean{y_j}\Big)\right] 
  = \sigma_{i,j} = \bm{\mathrm{Cov}}[y_i,y_j]
\end{displaymath}
die Elemente der Kovarianzmatrix $\bm{\Sigma}$ des Vektors $\bm{y}$.
Wird außerdem
\begin{displaymath}
(\Delta\,\theta)^2 = \bm{\mathrm{Var}}[\theta(\bm{y})]
\end{displaymath}
gesetzt, erhalten wir das allgemeine Fehlerfortpflanzungsgesetz
\begin{equation}
  (\Delta\,\theta)^2 = \sum_{i=1}^n\sum_{j=1}^n
  \left.\frac{\partial\theta(\bm{y})}{\partial y_i}\right|_{\bm{y}=\mean{\bm{y}}} 
  \;\sigma_{i,j}\;
  \left.\frac{\partial\theta(\bm{y})}{\partial y_j}\right|_{\bm{y}=\mean{\bm{y}}}
\label{allgFFG}
\end{equation}
Für den Fall, dass die Werte des Vektors $\bm{y}$ nicht korreliert
sind, wenn für $ i\ne j$ gilt $\bm{\mathrm{Cov}}[y_i,y_j] = 0$,
vereinfacht sich die Gleichung \ref{allgFFG} zu der häufig verwendeten Form:
\begin{displaymath}
  (\Delta\,\theta)^2 = \sum_{i=1}^n\left(
  \left.\frac{\partial\theta(\bm{y})}{\partial y_i}\right|_{\bm{y}=\mean{\bm{y}}}
  \Delta\,y_i \right)^2
\end{displaymath}
mit $(\Delta\,y_i)^2 = \sigma_{i,i} = \bm{\mathrm{Var}}[y_i]$.

%% Erweiterung für Vektor $\bm{\theta}$ \citeasnoun{martin} Seite 101
Das allgemeine Fehlerfortpflanzungsgesetz \ref{allgFFG} für eine Größe
$\theta$ lässt sich auf den $p$-dimensionalen Vektor $\bm{\theta}$ 
erweitern. Dessen Elemente $\theta_k$ sind alle
Funktionen des selben Vektors $\bm{y}$ und daher korreliert. Für die
Elemente der Kovarianzmatrix $\bm{\mathrm{Cov}}[\bm{\theta}]$ gilt dann:
\begin{equation}
  \bm{\mathrm{Cov}}[\theta_k,\theta_l] =  \sum_{i=1}^n\sum_{j=1}^n
  \left.\frac{\partial\theta_k}{\partial y_i}\right|_{\bm{y}=\mean{\bm{y}}} 
  \;\sigma_{i,j}\;
  \left.\frac{\partial\theta_l}{\partial y_j}\right|_{\bm{y}=\mean{\bm{y}}}
\label{FFG}
\end{equation}
%% \citeasnoun{martin} Seite 101 Gleichung 5.47 fehlerhaft. Summation muss von $1\ldots p$ über
%% die $p$ Parameter $\theta_i$ erfolgen. Vergleiche auch mit Gleichung 5.46 für
%% $l = k$
Wird eine $p \times n$ Matrix $\bm{G}$ mit den Elementen 
\begin{displaymath}
  g_{k,i} = \left.\frac{\partial\theta_k}{\partial y_i}\right|_{\bm{y}=\mean{\bm{y}}}
  \hspace{4mm} \mathrm{mit} \hspace{4mm} k = 1\dots p \hspace{4mm}
  \mathrm{und} \hspace{4mm} i = 1\ldots n
\end{displaymath}  
eingeführt, kann Gleichung \ref{FFG} als Matrizengleichung geschrieben
werden.
\begin{equation}
  \bm{\mathrm{Cov}}[\bm{\theta}] = \bm{G} \, \bm{\Sigma} \, \bm{G}^T
\label{MatrixFFG}
\end{equation}

Bezeichnet man die $p \times n $ Matrix
$(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}$ aus Gleichung \ref{SolutionAllgLinMod}
mit $\bm{H}$, 
%% $\bm{H}$ ist eine $p \times n$ Matrix
%% p => Anzahl der Zeilen
%% n => Anzahl der Spalten
dann kann die Gleichung \ref{SolutionAllgLinMod} auch geschrieben
werden als:
\begin{eqnarray*}
  \widehat{\bm{\theta}} & = & \bm{H} \, \bm{y}\\ 
  \widehat{\theta_k} & = &\sum_{i=1}^n h_{k,i}\, y_i 
\end{eqnarray*}
Aus dieser Darstellung wird sofort klar, dass für die einzelne
Elemente der Matrix $\bm{H}$ gilt:
\begin{displaymath}
  h_{k,i} = \frac{\partial \theta_k}{\partial y_i}
\end{displaymath}
Damit kann aus den Gleichungen \ref{MatrixFFG} und
\ref{Gewichtsmatrix} die Beziehung
\begin{eqnarray}
  \label{CovFFG}
  \bm{\mathrm{Cov}}[\bm{\theta}] &=& \bm{H}\,\bm{\Sigma}\,\bm{H}^T\\ \nonumber
  &=& c (\bm{A}^T\bm{P}\;\bm{A})^{-1}
\end{eqnarray} 
abgeleitet 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{ Bei der Herleitung wird ausgenutzt, dass die Matrizen
  $\bm{P}$ und $\bm{A}^T\bm{P}\;\bm{A}$ symmetrisch sind.
  \begin{eqnarray*}
    \bm{\mathrm{Cov}}[\bm{\theta}] &=& \bm{H}\,\bm{\Sigma}\,\bm{H}^T\\
    &=&(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\,\bm{\Sigma}\,
    \left((\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\right)^T\\
    &=&(\bm{A}^T\bm{P}\;\bm{A})^{-1}\bm{A}^T\bm{P}\,\bm{\Sigma}\,
    \bm{P}^T\bm{A}\left((\bm{A}^T\bm{P}\;\bm{A})^{-1}\right)^T\\
    & & \mathrm{mit:} \hspace{4mm}\bm{\Sigma} =  c\;\bm{P}^{-1} \\   
    &=&c\,(\bm{A}^T\bm{P}\;\bm{A})^{-1}(\bm{A}^T\bm{P}\bm{A})(\bm{A}^T\bm{P}\;\bm{A})^{-1}\\
    &=&c\,(\bm{A}^T\bm{P}\;\bm{A})^{-1}
  \end{eqnarray*}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
werden. Die so bestimmte Kovarianzmatrix enthält die, bei der
Aufstellung des Allgemeinen Linearen Modells willkürlich festlegbare
Konstante $c$. Alle aus der Kovarianzmatrix abgeleiteten
Fehlerintervalle sind daher von der Festlegung dieser Konstanten
abhängig. Nur wenn die Kovarianzmatrix $\bm{\Sigma}$ vollständig bekannt ist, oder
aus den vorliegenden Daten bestimmt werden kann, so wie es von
\citeasnoun{grabe} im Kapitel 21 beschrieben wird, führt die Festlegung
$c = 1$ beziehungsweise $\bm{P}=\bm{\Sigma}^{-1}$ zu richtigen
Ergebnissen. Nur in diesem Fall gilt:
\begin{equation}
  \bm{\mathrm{Cov}}[\bm{\theta}] = (\bm{A}^T\bm{\Sigma}^{-1}\bm{A})^{-1}
\label{CovC1}
\end{equation}
Auf diesen Umstand macht auch \citeasnoun{bevington} auf Seite 108
aufmerksam. Er schlägt vor, die bestimmten Varianzen, die Quadrate der
Fehlerintervalle, zu korrigieren, sie mit dem Wert von
$\chi^2/(n-p)$ zu multiplizieren, wenn $\chi^2/(n-p)$ zu stark von 1
abweicht. Genau dann ergibt sich die von von der klassischen Statistik
(Gleichung \ref{CovAllgLinMod}) berechnet Kovarianzmatrix. Auf
die zwingende Notwendigkeit der Berücksichtigung des Varianzfaktors
$\widehat{\sigma^2}$ weist auch \citeasnoun{martin} ausdrücklich im
Abschnitt 8.1.2 hin. Wird dies unterlassen, dann führt die Anwendung
des Fehlerfortpflanzungsgesetzes auf das Allgemeine Lineare Modell zu
falschen Ergebnissen. 

Damit ist die Ursache der Widersprüche in Tabelle \ref{tableResultate}
geklärt. Die Formeln aus dem  \citeasnoun{einfuehrung} liefern das
gleiche falsche Ergebnis wie die \citename{mathematica} Routine
\glqq{}LinearModelFit\grqq{} mit der Option
\glqq{}VarianceEstimatorFunction $\rightarrow$ (1\&)\grqq{}. In beiden
Fällen wird der Varianzfaktor auf den Wert 1 festgelegt. Ob das auch
bei QtiPlot der Fall ist, wird im weiteren noch untersucht.

Da in der im \citeasnoun{einfuehrung} auf Seite 42 angegebenen
Gleichung 45 der Varianzfaktor nicht berücksichtigt wird und dort
auch die obigen Hinweise in keiner Weise erwähnt werden, ist diese
Gleichung und die ihr zu Grunde liegende Vorgehensweise als falsch zu
betrachten.
\newpage

\section{Rechnet QtiPlot falsch ?}
Um diese Frage zu beantworten, ist es notwendig den Quellcode des
Programms genauer zu analysieren, da die zugehörige Dokumentation
keinerlei Informationen zu den genutzten Algorithmen, zu den
verwendeten statistischen Verfahren und den dahinterliegenden
statistischen Theorien enthält. \citeasnoun{qtiplot} nutzt die Routinen der GNU
Scientific Library. Diese sind ausführlich in \citeasnoun{gslref}
beschrieben. Dort bezieht man  sich auf \citeasnoun{barnett}, der
älteren Version von \citeasnoun{behringer}.

Wie die
\begin{latexonly}Quellcodeanalyse\end{latexonly}
\begin{htmlonly}\htmladdnormallink{Quellcodeanalyse}{./Quelltextanalyse-QtiPlot.pdf}
\end{htmlonly}zeigt, werden
für die Fitverfahren, die auf dem Allgemeinen Linearen Modell basieren,
mehrere verschiedene Bibliotheksfunktionen aufgerufen. Für den
Fall ungewichteter Daten verwendet QtiPlot die Bibliotheksfunktionen
{\tt int gsl\_fit\_linear (\ldots)} für den linearen Fit
beziehungsweise die Funktion {\tt int gsl\_fit\_mul (\ldots)} für den
Fit einer Geraden durch den 
Nullpunkt. Beide Routinen geben nach \citeasnoun{gslref}
(Seite 404) die Elemente der \glqq{}Varianz-Kovarianz-Matrix\grqq{}
der Parameter zurück. Diese Matrix entspricht der Gleichung
\ref{CovAllgLinMod}. Für den Fall gewichteter Daten kommen die
Bibliotheksfunktionen {\tt int gsl\_fit\_wlinear (\ldots)}
beziehungsweise die Funktion {\tt int gsl\_fit\_wmul (\ldots)} (\citeasnoun{gslref}
Seite 404 folgende) zur Anwendung. Diese berechnen die Kovarianz-Matrix
der Parameter. Diese entspricht nach \citeasnoun{barnett} (Seite 161,
Gleichung. 28.28) beziehungsweise \citeasnoun{behringer} (Seite 392 Gleichung 36.17)
der Gleichung \ref{CovC1}. Um daraus die
\glqq{}Varianz-Kovarianz-Matrix\grqq{} zu erhalten, ist die
Multiplikation mit $\chi^2/(n-p) = \widehat{\sigma^2}$ notwendig. Der
Wert von  $\chi^2$ wird von allen vier Bibliotheksfunktionen berechnet.
Dass zwischen gewichteten und ungewichteten Daten unterschieden wird,
obwohl dies vom zugrundeliegenden Allgemeinen Linearen Modell nicht
notwendig wäre, hat rein numerische Gründe. Es erspart an mehreren
Stellen die Ausführung der Multiplikation mit der Gewichtsmatrix, die
in diesem Falle die Einheitsmatrix ist. Eine Folge der unterschiedlichen
Bedeutung der von den verwendeten Funktionen zurückgegebenen Werte ist, dass
die Fehler im Fall gewichteter Daten noch mit {\tt sqrt(chiquadrat/doF)}    
multipliziert werden müssen, während dies im Fall von ungewichtete Daten
bereits erfolgt ist.  Die Quelltextanalyse ergab, dass die Option
\glqq{}Gewichtung mit {\tt sqrt(chiquadrat/doF)}\grqq{} bei allen auf
dem linearen Modell basierenden Fitverfahren wirkungslos ist. Die
Berücksichtigung des Varianzfaktors $\widehat{\sigma^2}$ bei der
Berechnung der Fehlerintervalle muss im Nachgang von Hand durch den
Anwender von QtiPlot erfolgen.

Die Frage, ob QtiPlot numerisch falsch rechnet, kann man nun
eindeutig mit nein beantworten. Dies gilt zumindest für alle
Fitroutinen, die auf dem Allgemeinen Linearen Modell basieren und die
obengenannten Funktionen der GNU Scientific Library verwenden. Die
richtige Interpretation der berechneten Werte durch die Anwender des
Programms ist jedoch  auf Grund der mangelhaften Dokumentation der
Bedeutung der ausgegebenen Werte und der verwendeten statistischen
Verfahren nicht möglich. 

\section{Schlussfolgerungen}
\label{con}
Wie sich gezeigt hat, widersprechen die im \citename{einfuehrung}
angegebenen Formeln, den Aussagen der beiden dominierenden Richtungen
der Statistik. Sie stützen sich auf eine nur unter ganz strengen
Rahmenbedingungen anwendbare Vorgehensweise. Die dafür notwendigen
Voraussetzungen bei der Datenerhebung werden jedoch nicht
behandelt. An diesem Punkt ist das \citename{einfuehrung} als
fehlerhaft anzusehen.

\citename{qtiplot} rechnet numerisch richtig. Durch die ungenügende
Dokumentation ist der Anwender des Programms nicht in der Lage die
gelieferten numerischen Werte im Sinne der zu Grunde liegenden
statistischen Methoden, unabhängig davon ob diese auf der Klassischen
Statistik oder Bayes-Statistik beruhen, richtig zu interpretieren. 
Vor diesem Hintergrund ist nicht nur von der Verwendung des noch
weniger dokumentierten Fit-Wizard von \citename{qtiplot} abzuraten, sondern
generell von der Anwendung von \citename{qtiplot} und anderer aus \citename{qtiplot}
hervorgegangenen Programmen wie \citename{scidavis}. Für die Belange der
Ausbildung von Studenten im Rahmen des Praktikums sind gut
dokumentierte, frei zugängliche Werkzeuge wie \citename{gnuplot} wesentlich
sinnvoller. 

\newpage

%begin{latexonly}
\begingroup
     \parindent 0pt
     \parskip 2ex
     \def\enotesize{\normalsize}
     \theendnotes
     \endgroup 
%end{latexonly}
     
\newpage

\begin{thebibliography}{00}

\harvarditem[Ein\-führungs\-skript]{Ein\-führungs\-skript}{2007}{einfuehrung}
Physikalisches Grundpraktikum: 
Einführung in die Messung, Auswertung und Darstellung experimenteller
Ergebnisse in der Physik\\
Mathematisch-Naturwissenschaftliche Fakultät\,I der HUB\\
Institut für Physik {\bf 2007}\\
\htmladdnormallink{\footnotesize
 http://gpr.physik.hu-berlin.de/Skripten/Einfuehrung/PDF-Datei/Einfuehrung.pdf}{ http://gpr.physik.hu-berlin.de/Skripten/Einfuehrung/PDF-Datei/Einfuehrung.pdf}
 {\footnotesize abgerufen am 4.1.2016 10:39 Uhr}
\normalsize 

\harvarditem[Einführungspraktikum]{Einführungspraktikum}{2007}{anleitungEP}
Physikalisches Grundpraktikum: Einführungspraktikum\\
Mathematisch-Naturwissenschaftliche Fakultät\,I der HUB\\
Institut für Physik {\bf 2007}\\
\htmladdnormallink{\footnotesize http://gpr.physik.hu-berlin.de/Skripten/Einfuehrungspraktikum/PDF-Dateien/Einfuehrungspraktikum.pdf}{http://gpr.physik.hu-berlin.de/Skripten/Einfuehrungspraktikum/PDF-Dateien/Einfuehrungspraktikum.pdf}
\normalsize 

\harvarditem[Barnett]{R.M.\,Barnett u.a.}{1996}{barnett}
R.M. Barnett et al.\\
Review of Particle Properties\\
28. STATISTICS page 159 \\
Physical Review D Volume {\bf 54}, 1 {\bf 1996}

\harvarditem[Behringer]{J.\,Beringer u.a.}{2012}{behringer}
J. Beringer et al. (Particle Data Group)\\ 
REVIEW OF PARTICLE PHYSICS\\
36. STATISTICS page 390\\
Physical Review D Volume {\bf 86}, 010001 {\bf 2012}

\harvarditem[Bevington]{P.R.\,Bevington, D.K.\,Robinson}{2003}{bevington}
P.R.\,Bevington and D.K.\,Robinson\\
Data Reduction and Error Analysis for the Physical Sciences\\
McGraw-Hill Book Co., New York, 3.\,edition {\bf 2003}

\harvarditem[James]{F.\,James}{2006}{james}
Frederick James\\
Statistical Methods in Experimental Physics\\
World Scientific Publishing, Singapore, 2nd edition {\bf 2006}

\harvarditem[Fahrmeir]{L.\,Fahrmeier u.a.}{2009}{fahrmeir}
Ludwig Fahrmeir, Thomas Kneib, Stefan Lang\\
Regression \hspace{3mm}
Modelle, Methoden und Anwendungen\\
Springer-Verlag, Berlin Heidelberg 2.\,Auflage {\bf 2009}

\harvarditem[Fahrmeir]{L.\,Fahrmeier u.a.}{2013}{fahrmeir2013}
Ludwig Fahrmeir, Thomas Kneib, Stefan Lang, Brian Marx \\
Regression \hspace{3mm}
Models, Methods and Applications\\
Springer-Verlag, Berlin Heidelberg  First edition {\bf 2013}

\harvarditem[Grabe]{M.\,Grabe}{2011}{grabe}
Michael Grabe\\
Grundriss der Generalisierten Gauß´chen Fehlerrechnung\\
Springer-Verlag, Heidelberg  Dordrecht London {\bf 2011}

\harvarditem[Koch]{K.\,Koch}{2000}{koch}
Karl-Rudolf Koch\\
Einführung in die Bayes-Statistik\\
Springer Verlag Berlin Heidelberg 1.Auflage {\bf 2000}

\harvarditem[Martin]{B.\,R.\,Martin}{2012}{martin}
Brain R. Martin\\
Statistics for physical Sciences\\
Academic Press Elsevier, Amsterdam, First edition {\bf 2012}

\harvarditem[Nollau]{V.\,Nollau}{1975}{nollau}
Volker Nollau\\
Statistische Analysen \\
Mathematische Methoden der Planung und Auswertung von Versuchen\\
Fachbuchverlag Leipzig 1.Auflage {\bf 1975}

\harvarditem[Schenk]{W.\,Schenk u.a.}{2012}{schenk}
Wolfgang Schenk, Friedrich Kremer, Gunter Beddies, Thomas Franke, Petrik Galvosas, Peter Rieger\\
Physikalisches Praktikum  14.,überarbeitete und erweiterte Auflage\\
Springer Fachmedien Wiesbaden {\bf 2014}

\harvarditem[Wakefield]{J.\,Wakefield}{2013}{wakefield}
Jon Wakefield\\
Bayesian and Frequentist Regression Methods\\
Springer Verlag New York Heidelberg 1.Edition {\bf 2013}

\harvarditem[Wolff]{U.\,Wolff}{2014}{wolf}
Ulli Wolf, Peter Galler, Tomasz Korzec\\
Vorlesung Computational Physics I\\ 
Institut fur Physik der HUB {\bf 2014}\\
\htmladdnormallink{\footnotesize
  http://www.physik.hu-berlin.de/com/teachingandseminars/previousCPI}{http://www.physik.hu-berlin.de/com/teachingandseminars/previousCPI}\\
{\footnotesize abgerufen am 4.1.2016 11:31 Uhr}
\normalsize 

\harvarditem[Mathematica]{Mathematica 10.0}{2014}{mathematica} 
Wolfram Research, Inc.\\
Mathematica, Version 10.0, Champaign, IL {\bf 2014}\\
\htmladdnormallink{\footnotesize http://www.wolframalpha.com}{http://www.wolframalpha.com}
\normalsize

\harvarditem[QtiPlot]{QtiPlot 0.9.8.9-9}{2014}{qtiplot}
Ion Vasilief, Scott Howard (Debian Maintainer)\\
QtiPlot Version 0.9.8.9-9+b3\\
Released: 2014/04/29\\
ubuntu package: qtiplot\_0.9.8.9-9+b3\_amd64.deb\\
\htmladdnormallink{\footnotesize http://www.qtiplot.com}{http://www.qtiplot.com}
\normalsize

\harvarditem[SciDAVis]{SciDAVis}{2015}{scidavis}
Russell Standish,  Georges Khaznadar (Debian Maintainer)\\
SciDAVis Version 1.D9-1\\
Released: 2015/12/10\\
Debian package: scidavis\_1.D9-1\_amd64.deb\\
\htmladdnormallink{\footnotesize http://scidavis.sourceforge.net}{http://scidavis.sourceforge.net}
\normalsize 

\harvarditem[GnuPlot]{GnuPlot}{2015}{gnuplot}
Anton Gladky (Debian Maintainer)\\
GnuPlot Version 5.0.1+dfsg1-3\\
Released: 2015/08/11\\
Debian package: gnuplot5\_5.0.1+dfsg1-3\_all.deb\\
\htmladdnormallink{\footnotesize http://gnuplot.sourceforge.net}{http://gnuplot.sourceforge.net}
\normalsize

\harvarditem[GSL Reference Manual]{GNU Scientific Library, Reference Manual}{2011}{gslref}
GNU Scientific Library Reference Manual \\
Edition 1.15 for GSL Version 1.15, 29 April 2011\\
Abschnitt 37 page 403\\
\htmladdnormallink{\footnotesize https://www.gnu.org/s/gsl/manual/gsl-ref.pdf}{https://www.gnu.org/s/gsl/manual/gsl-ref.pdf}
\normalsize

%% freie Software
%% BayesX
%% http://www.statistik.lmu.de/~kneib/regressionsbuch/software_e.html
%% Windows, Linux, Source
%% GNU General Public License
%%
%% The R Project for Statistical Computing
%% https://www.r-project.org/
%% Windows, Mac, Linux, Source
%% GNU General Public License
%%

\end{thebibliography}
%begin{latexonly}
\newpage
\newgeometry{textheight=25cm}
\setlength{\unitlength}{1mm}
\section{Anhang}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Ergebnisse-QtiPlot}
     \end{picture}
     \caption{Ergebnisse QtiPlot}
     \label{ErgebnisseQtiPlot}
\end{figure}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Ergebnisse-Mathematica-1}
     \end{picture}
     \caption{Mathematica Notebook Seite 1}
     \label{ErgebnisseMathematic1}
\end{figure}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Ergebnisse-Mathematica-2}
     \end{picture}
     \caption{Mathematica Notebook Seite 2}
     \label{ErgebnisseMathematic2}
\end{figure}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Ergebnisse-Mathematica-3}
     \end{picture}
     \caption{Mathematica Notebook Seite 3}
     \label{ErgebnisseMathematic3}
\end{figure}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Quelltextanalyse-QtiPlot-1}
     \end{picture}
     \caption{Quellcodeanalyse Seite 1}
     \label{Quellcode1}
\end{figure}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Quelltextanalyse-QtiPlot-2}
     \end{picture}
     \caption{Quellcodeanalyse Seite 2}
     \label{Quellcode2}
\end{figure}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Quelltextanalyse-QtiPlot-3}
     \end{picture}
     \caption{Quellcodeanalyse Seite 3}
     \label{Quellcode3}
\end{figure}
\begin{figure}[htbp]
     \begin{picture}(148,200)
       \includegraphics[scale=0.7]{Quelltextanalyse-QtiPlot-4}
     \end{picture}
     \caption{Quellcodeanalyse Seite 4}
     \label{Quellcode4}
\end{figure}
\clearpage
%end{latexonly}

\end{document}
