diff --git a/statistics/lecture/Makefile b/statistics/lecture/Makefile index 59963db..6f5131d 100644 --- a/statistics/lecture/Makefile +++ b/statistics/lecture/Makefile @@ -1,4 +1,4 @@ -TEXFILES=$(wildcard *.tex) +TEXFILES=descriptivestatistics.tex linear_regression.tex #$(wildcard *.tex) PDFFILES=$(TEXFILES:.tex=.pdf) PYFILES=$(wildcard *.py) PYPDFFILES=$(PYFILES:.py=.pdf) diff --git a/statistics/lecture/descriptivestatistics.tex b/statistics/lecture/descriptivestatistics.tex index 836a1d5..4888ada 100644 --- a/statistics/lecture/descriptivestatistics.tex +++ b/statistics/lecture/descriptivestatistics.tex @@ -293,7 +293,7 @@ \begin{figure}[t] \includegraphics[width=1\textwidth]{quartile} - \caption{\label{quartilefig} Median und Quartile.} + \caption{\label{quartilefig} Median und Quartile einer Normalverteilung.} \end{figure} \begin{definition}[\tr{quartile}{Quartile}] @@ -640,6 +640,204 @@ Stichprobe vorkommen. Was hat das mit dem Standardfehler zu tun? \end{exercise} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\chapter{\tr{Maximum likelihood estimation}{Maximum-Likelihood Methode}} + +In vielen Situationen wollen wir einen oder mehrere Parameter $\theta$ +einer Wahrscheinlichkeitsverteilung sch\"atzen, so dass die Verteilung +die Daten $x_1, x_2, \ldots x_n$ am besten beschreibt. Bei der +Maximum-Likelihood-Methode w\"ahlen wir die Parameter so, dass die +Wahrscheinlichkeit, dass die Daten aus der Verteilung stammen, am +gr\"o{\ss}ten ist. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Maximum Likelihood} +Sei $p(x|\theta)$ (lies ``Wahrscheinlichkeit(sdichte) von $x$ gegeben +$\theta$'') die Wahrscheinlichkeits(dichte)verteilung von $x$ mit dem +Parameter(n) $\theta$. Das k\"onnte die Normalverteilung +\begin{equation} + \label{normpdfmean} + p(x|\theta) = \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x-\theta)^2}{2\sigma^2}} +\end{equation} +sein mit +fester Standardverteilung $\sigma$ und dem Mittelwert $\mu$ als +Parameter $\theta$. + +Wenn nun den $n$ unabh\"angigen Beobachtungen $x_1, x_2, \ldots x_n$ +die Wahrscheinlichkeitsverteilung $p(x|\theta)$ zugrundeliegt, dann +ist die Verbundwahrscheinlichkeit $p(x_1,x_2, \ldots x_n|\theta)$ des +Auftretens der Werte $x_1, x_2, \ldots x_n$ gegeben ein bestimmtes $\theta$ +\[ p(x_1,x_2, \ldots x_n|\theta) = p(x_1|\theta) \cdot p(x_2|\theta) +\ldots p(x_n|\theta) = \prod_{i=1}^n p(x_i|\theta) \; .\] +Andersherum gesehen ist das die Likelihood (deutsch immer noch ``Wahrscheinlichleit'') +den Parameter $\theta$ zu haben, gegeben die Me{\ss}werte $x_1, x_2, \ldots x_n$, +\[ {\cal L}(\theta|x_1,x_2, \ldots x_n) = p(x_1,x_2, \ldots x_n|\theta) \] + +Wir sind nun an dem Wert des Parameters $\theta_{mle}$ interessiert, der die +Likelihood maximiert (``mle'': Maximum-Likelihood Estimate): +\[ \theta_{mle} = \text{argmax}_{\theta} {\cal L}(\theta|x_1,x_2, +\ldots x_n) \] +$\text{argmax}_xf(x)$ bezeichnet den Wert des Arguments $x$ der Funktion $f(x)$, bei +dem $f(x)$ ihr globales Maximum annimmt. Wir suchen also den Wert von $\theta$ +bei dem die Likelihood ${\cal L}(\theta)$ ihr Maximum hat. + +An der Stelle eines Maximums einer Funktion \"andert sich nichts, wenn +man die Funktionswerte mit einer streng monoton steigenden Funktion +transformiert. Aus gleich ersichtlichen mathematischen Gr\"unden wird meistens +das Maximum der logarithmierten Likelihood (``Log-Likelihood'') gesucht: +\begin{eqnarray} + \theta_{mle} & = & \text{argmax}_{\theta}\; {\cal L}(\theta|x_1,x_2, \ldots x_n) \nonumber \\ + & = & \text{argmax}_{\theta}\; \log {\cal L}(\theta|x_1,x_2, \ldots x_n) \nonumber \\ + & = & \text{argmax}_{\theta}\; \log \prod_{i=1}^n p(x_i|\theta) \nonumber \\ + & = & \text{argmax}_{\theta}\; \sum_{i=1}^n \log p(x_i|\theta) \label{loglikelihood} +\end{eqnarray} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{Beispiel: Das arithmetische Mittel} + +Wenn die Me{\ss}daten $x_1, x_2, \ldots x_n$ der Normalverteilung \eqnref{normpdfmean} +entstammen, und wir den Mittelwert $\mu$ als einzigen Parameter der Verteilung betrachten, +welcher Wert von $\theta$ maximiert dessen Likelhood? + +\begin{figure}[t] + \includegraphics[width=1\textwidth]{mlemean} + \caption{\label{mlemeanfig} Maximum Likelihood Estimation des + Mittelwerts. Oben: Die Daten zusammen mit drei m\"oglichen + Normalverteilungen mit unterschiedlichen Mittelwerten (Pfeile) aus + denen die Daten stammen k\"onnten. Unteln links: Die Likelihood + in Abh\"angigkeit des Mittelwerts als Parameter der + Normalverteilungen. Unten rechts: die entsprechende + Log-Likelihood. An der Position des Maximums bei $\theta=2$ + \"andert sich nichts (Pfeil).} +\end{figure} + +Die Log-Likelihood \eqnref{loglikelihood} ist +\begin{eqnarray*} + \log {\cal L}(\theta|x_1,x_2, \ldots x_n) + & = & \sum_{i=1}^n \log \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x_i-\theta)^2}{2\sigma^2}} \\ + & = & \sum_{i=1}^n - \log \sqrt{2\pi \sigma^2} -\frac{(x_i-\theta)^2}{2\sigma^2} +\end{eqnarray*} +Zur Bestimmung des Maximums der Log-Likelihood berechnen wir deren Ableitung +nach dem Parameter $\theta$ und setzen diese gleich Null: +\begin{eqnarray*} + \frac{\text{d}}{\text{d}\theta} \log {\cal L}(\theta|x_1,x_2, \ldots x_n) & = & \sum_{i=1}^n \frac{2(x_i-\theta)}{2\sigma^2} \;\; = \;\; 0 \\ + \Leftrightarrow \quad \sum_{i=1}^n x_i - \sum_{i=1}^n x_i \theta & = & 0 \\ + \Leftrightarrow \quad n \theta & = & \sum_{i=1}^n x_i \\ + \Leftrightarrow \quad \theta & = & \frac{1}{n} \sum_{i=1}^n x_i +\end{eqnarray*} +Der Maximum-Likelihood-Estimator ist das arithmetische Mittel der Daten. D.h. +das arithmetische Mittel maximiert die Wahrscheinlichkeit, dass die Daten aus einer +Normalverteilung mit diesem Mittelwert gezogen worden sind. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Kurvenfit als Maximum Likelihood Estimation} +Beim Kurvenfit soll eine Funktion $f(x;\theta)$ mit den Parametern +$\theta$ an die Datenpaare $(x_i|y_i)$ durch Anpassung der Parameter +$\theta$ gefittet werden. Wenn wir annehmen, dass die $y_i$ um die +entsprechenden Funktionswerte $f(x_i;\theta)$ mit einer +Standardabweichung $\sigma_i$ normalverteilt streuen, dann lautet die +Log-Likelihood +\begin{eqnarray*} + \log {\cal L}(\theta|x_1,x_2, \ldots x_n) + & = & \sum_{i=1}^n \log \frac{1}{\sqrt{2\pi \sigma_i^2}}e^{-\frac{(y_i-f(x_i;\theta))^2}{2\sigma_i^2}} \\ + & = & \sum_{i=1}^n - \log \sqrt{2\pi \sigma_i^2} -\frac{(x_i-f(y_i;\theta))^2}{2\sigma_i^2} \\ +\end{eqnarray*} +Der einzige Unterschied zum vorherigen Beispiel ist, dass die +Mittelwerte der Normalverteilungen nun durch die Funktionswerte +gegeben sind. + +Der Parameter $\theta$ soll so gew\"ahlt werden, dass die +Log-Likelihood maximal wird. Der erste Term der Summe ist +unabh\"angig von $\theta$ und kann deshalb bei der Suche nach dem +Maximum weggelassen werden. +\begin{eqnarray*} + & = & - \frac{1}{2} \sum_{i=1}^n \left( \frac{y_i-f(x_i;\theta)}{\sigma_i} \right)^2 +\end{eqnarray*} +Anstatt nach dem Maximum zu suchen, k\"onnen wir auch das Vorzeichen der Log-Likelihood +umdrehen und nach dem Minimum suchen. Dabei k\"onnen wir auch den Faktor $1/2$ vor der Summe vernachl\"assigen --- auch das \"andert nichts an der Position des Minimums. +\begin{eqnarray*} + \theta_{mle} & = & \text{argmin}_{\theta} \; \sum_{i=1}^n \left( \frac{y_i-f(x_i;\theta)}{\sigma_i} \right)^2 \;\; = \;\; \text{argmin}_{\theta} \; \chi^2 +\end{eqnarray*} +Die Summer der quadratischen Abst\"ande normiert auf die jeweiligen +Standardabweichungen wird auch mit $\chi^2$ bezeichnet. Der Wert des +Parameters $\theta$ welcher den quadratischen Abstand minimiert ist +also identisch mit der Maximierung der Wahrscheinlichkeit, dass die +Daten tats\"achlich aus der Funktion stammen k\"onnen. Minimierung des +$\chi^2$ ist also ein Maximum-Likelihood Estimate. + +\begin{figure}[t] + \includegraphics[width=1\textwidth]{mlepropline} + \caption{\label{mleproplinefig} Maximum Likelihood Estimation der + Steigung einer Ursprungsgeraden.} +\end{figure} + + +\subsection{Beispiel: einfache Proportionalit\"at} +Als Funktion nehmen wir die Ursprungsgerade +\[ f(x) = \theta x \] +mit Steigung $\theta$. Die $\chi^2$-Summe lautet damit +\[ \chi^2 = \sum_{i=1}^n \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \; . \] +Zur Bestimmung des Minimums berechnen wir wieder die erste Ableitung nach $\theta$ +und setzen diese gleich Null: +\begin{eqnarray*} + \frac{\text{d}}{\text{d}\theta}\chi^2 & = & \frac{\text{d}}{\text{d}\theta} \sum_{i=1}^n \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \\ + & = & \sum_{i=1}^n \frac{\text{d}}{\text{d}\theta} \left( \frac{y_i-\theta x_i}{\sigma_i} \right)^2 \\ + & = & -2 \sum_{i=1}^n \frac{x_i}{\sigma_i} \left( \frac{y_i-\theta x_i}{\sigma_i} \right) \\ + & = & -2 \sum_{i=1}^n \left( \frac{x_iy_i}{\sigma_i^2} - \theta \frac{x_i^2}{\sigma_i^2} \right) \;\; = \;\; 0 \\ +\Leftrightarrow \quad \theta \sum_{i=1}^n \frac{x_i^2}{\sigma_i^2} & = & \sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2} \\ +\Leftrightarrow \quad \theta & = & \frac{\sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2}}{ \sum_{i=1}^n \frac{x_i^2}{\sigma_i^2}} +\end{eqnarray*} +Damit haben wir nun einen anlytischen Ausdruck f\"ur die Bestimmung +der Steigung $\theta$ des Regressionsgeraden gewonnen. Ein +Gradientenabstieg ist f\"ur das Fitten der Geradensteigung also gar nicht +n\"otig. Das gilt allgemein f\"ur das fitten von Koeffizienten von +linear kombinierten Basisfunktionen. Parameter die nichtlinear in +einer Funktion enthalten sind k\"onnen aber nicht analytisch aus den +Daten berechnet werden. Da bleibt dann nur auf numerische Verfahren +zur Optimierung der Kostenfunktion, wie z.B. der Gradientenabstieg, +zur\"uckzugreifen. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Fits von Wahrscheinlichkeitsverteilungen} +Zum Abschluss betrachten wir noch den Fall, bei dem wir die Parameter +einer Wahrscheinlichkeitsdichtefunktion (z.B. Mittelwert und +Standardabweichung der Normalverteilung) an ein Datenset fitten wolle. + +Ein erster Gedanke k\"onnte sein, die +Wahrscheinlichkeitsdichtefunktion durch Minimierung des quadratischen +Abstands an ein Histogram der Daten zu fitten. Das ist aber aus +folgenden Gr\"unden nicht die Methode der Wahl: (i) +Wahrscheinlichkeitsdichten k\"onnen nur positiv sein. Darum k\"onnen +insbesondere bei kleinen Werten die Daten nicht symmetrisch streuen, +wie es normalverteilte Daten machen sollten. (ii) Die Datenwerte sind +nicht unabh\"angig, da das normierte Histogram sich zu Eins +aufintegriert. Die beiden Annahmen normalverteilte und unabh\"angige Daten +die die Minimierung des quadratischen Abstands zu einem Maximum +Likelihood Estimator machen sind also verletzt. + +Den direkten Weg, eine Wahrscheinlichkeitsdichtefunktion an ein +Datenset zu fitten, haben wir oben schon bei dem Beispiel zur +Absch\"atzung des Mittelwertes einer Normalverteilung gesehen --- +Maximum Likelihood! Wir suchen einfach die Parameter $\theta$ der +gesuchten Wahrscheinlichkeitsdichtefunktion bei der die Log-Likelihood +\eqnref{loglikelihood} maximal wird. Das ist im allgemeinen ein +nichtlinieares Optimierungsproblem, das mit numerischen Verfahren, wie +z.B. dem Gradientenabstieg, gel\"ost wird. + +\begin{figure}[t] + \includegraphics[width=1\textwidth]{mlepdf} + \caption{\label{mlepdffig} Maximum Likelihood Estimation einer + Wahrscheinlichkeitsdichtefunktion. Links: die 100 Datenpunkte, die aus der Gammaverteilung + 2. Ordnung (rot) gezogen worden sind. Der Maximum-Likelihood-Fit ist orange dargestellt. + Rechts: das normierte Histogramm der Daten zusammen mit der \"uber Minimierung + des quadratischen Abstands zum Histogramm berechneten Fits.} +\end{figure} + + \end{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -653,3 +851,4 @@ What is "a statistic"? % dt. Sch\"atzfunktion \source{http://en.wikipedia.org/wiki/Statistic} \end{definition} + diff --git a/statistics/lecture/mlemean.py b/statistics/lecture/mlemean.py new file mode 100644 index 0000000..892e663 --- /dev/null +++ b/statistics/lecture/mlemean.py @@ -0,0 +1,99 @@ +import numpy as np +import matplotlib.pyplot as plt + +plt.xkcd() +fig = plt.figure( figsize=(6,5) ) + +# the data: +n = 40 +rng = np.random.RandomState(54637281) +sigma = 0.5 +rmu = 2.0 +xd = rng.randn(n)*sigma+rmu +# and possible pdfs: +x = np.arange( 0.0, 4.0, 0.01 ) +mus = [1.5, 2.0, 2.5] +g=np.zeros((len(x), len(mus))) +for k, mu in enumerate(mus) : + g[:,k] = np.exp(-0.5*((x-mu)/sigma)**2.0)/np.sqrt(2.0*np.pi)/sigma +# plot it: +ax = fig.add_subplot( 2, 1, 1 ) +ax.spines['right'].set_visible(False) +ax.spines['top'].set_visible(False) +ax.yaxis.set_ticks_position('left') +ax.xaxis.set_ticks_position('bottom') +ax.set_xlim(0.5, 3.5) +ax.set_ylim(-0.02, 0.85) +ax.set_xticks( np.arange(0, 5)) +ax.set_yticks( np.arange(0, 0.9, 0.2)) +ax.set_xlabel('x') +ax.set_ylabel('Probability density') +s = 1 +for mu in mus : + r = 5.0*rng.rand()+2.0 + cs = 'angle3,angleA={:.0f},angleB={:.0f}'.format(90+s*r, 90-s*r) + s *= -1 + ax.annotate('', xy=(mu, 0.02), xycoords='data', + xytext=(mu, 0.75), textcoords='data', + arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5), + connectionstyle=cs), zorder=1 ) + if mu > rmu : + ax.text(mu-0.1, 0.04, '?', zorder=1, ha='right') + else : + ax.text(mu+0.1, 0.04, '?', zorder=1) +for k in xrange(len(mus)) : + ax.plot(x, g[:,k], zorder=5) +ax.scatter(xd, 0.05*rng.rand(len(xd))+0.2, s=30, zorder=10) + +# likelihood: +thetas=np.arange(1.5, 2.6, 0.01) +ps=np.zeros((len(xd),len(thetas))); +for i, theta in enumerate(thetas) : + ps[:,i]=np.exp(-0.5*((xd-theta)/sigma)**2.0)/np.sqrt(2.0*np.pi)/sigma +p=np.prod(ps,axis=0) +# plot it: +ax = fig.add_subplot( 2, 2, 3 ) +ax.spines['right'].set_visible(False) +ax.spines['top'].set_visible(False) +ax.yaxis.set_ticks_position('left') +ax.xaxis.set_ticks_position('bottom') +ax.set_xlabel(r'Parameter $\theta$') +ax.set_ylabel('Likelihood') +ax.set_xticks( np.arange(1.6, 2.5, 0.4)) +ax.annotate('Maximum', + xy=(2.0, 5.5e-11), xycoords='data', + xytext=(1.0, 1.1), textcoords='axes fraction', ha='right', + arrowprops=dict(arrowstyle="->", relpos=(0.0,0.5), + connectionstyle="angle3,angleA=10,angleB=70") ) +ax.annotate('', + xy=(2.0, 0), xycoords='data', + xytext=(2.0, 5e-11), textcoords='data', + arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5), + connectionstyle="angle3,angleA=90,angleB=80") ) +ax.plot(thetas,p) + +ax = fig.add_subplot( 2, 2, 4 ) +ax.spines['right'].set_visible(False) +ax.spines['top'].set_visible(False) +ax.yaxis.set_ticks_position('left') +ax.xaxis.set_ticks_position('bottom') +ax.set_xlabel(r'Parameter $\theta$') +ax.set_ylabel('Log-Likelihood') +ax.set_ylim(-50,-20) +ax.set_xticks( np.arange(1.6, 2.5, 0.4)) +ax.set_yticks( np.arange(-50, -19, 10.0)) +ax.annotate('Maximum', + xy=(2.0, -23), xycoords='data', + xytext=(1.0, 1.1), textcoords='axes fraction', ha='right', + arrowprops=dict(arrowstyle="->", relpos=(0.0,0.5), + connectionstyle="angle3,angleA=10,angleB=70") ) +ax.annotate('', + xy=(2.0, -50), xycoords='data', + xytext=(2.0, -26), textcoords='data', + arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5), + connectionstyle="angle3,angleA=80,angleB=100") ) +ax.plot(thetas,np.log(p)) + +plt.tight_layout(); +plt.savefig('mlemean.pdf') +#plt.show(); diff --git a/statistics/lecture/mlepdf.py b/statistics/lecture/mlepdf.py new file mode 100644 index 0000000..22f89bf --- /dev/null +++ b/statistics/lecture/mlepdf.py @@ -0,0 +1,70 @@ +import numpy as np +import scipy.stats as st +import scipy.optimize as opt +import matplotlib.pyplot as plt + +plt.xkcd() +fig = plt.figure( figsize=(6,3) ) + +# the data: +n = 100 +shape = 2.0 +scale = 1.0 +rng = np.random.RandomState(4637281) +xd = rng.gamma(shape, scale, n) + +# true pdf: +xx = np.arange(0.0, 10.1, 0.01) +rv = st.gamma(shape) +yy = rv.pdf(xx) + +# mle fit: +a = st.gamma.fit(xd, 5.0) +yf = st.gamma.pdf(xx, *a) + +# plot it: +ax = fig.add_subplot( 1, 2, 1 ) +ax.spines['right'].set_visible(False) +ax.spines['top'].set_visible(False) +ax.yaxis.set_ticks_position('left') +ax.xaxis.set_ticks_position('bottom') +ax.set_xlim(0, 10.0) +ax.set_ylim(0.0, 0.42) +ax.set_xticks( np.arange(0, 11, 2)) +ax.set_yticks( np.arange(0, 0.42, 0.1)) +ax.set_xlabel('x') +ax.set_ylabel('Probability density') +ax.plot(xx, yy, '-', lw=5, color='#ff0000', label='pdf') +ax.plot(xx, yf, '-', lw=2, color='#ffcc00', label='mle') +ax.scatter(xd, 0.025*rng.rand(len(xd))+0.05, s=30, zorder=10) +ax.legend(loc='upper right', frameon=False) + +# histogram: +h,b = np.histogram(xd, np.arange(0, 8.5, 1), density=True) + +# fit histogram: +def gammapdf(x, n, l, s) : + return st.gamma.pdf(x, n, l, s) +popt, pcov = opt.curve_fit(gammapdf, b[:-1]+0.5*(b[1]-b[0]), h) +yc = st.gamma.pdf(xx, *popt) + +# plot it: +ax = fig.add_subplot( 1, 2, 2 ) +ax.spines['right'].set_visible(False) +ax.spines['top'].set_visible(False) +ax.yaxis.set_ticks_position('left') +ax.xaxis.set_ticks_position('bottom') +ax.set_xlim(0, 10.0) +ax.set_xticks( np.arange(0, 11, 2)) +ax.set_xlabel('x') +ax.set_ylim(0.0, 0.42) +ax.set_yticks( np.arange(0, 0.42, 0.1)) +ax.set_ylabel('Probability density') +ax.plot(xx, yy, '-', lw=5, color='#ff0000', label='pdf') +ax.plot(xx, yc, '-', lw=2, color='#ffcc00', label='fit') +ax.bar(b[:-1], h, np.diff(b)) +ax.legend(loc='upper right', frameon=False) + +plt.tight_layout(); +plt.savefig('mlepdf.pdf') +#plt.show(); diff --git a/statistics/lecture/mlepropline.py b/statistics/lecture/mlepropline.py new file mode 100644 index 0000000..fc28d48 --- /dev/null +++ b/statistics/lecture/mlepropline.py @@ -0,0 +1,49 @@ +import numpy as np +import matplotlib.pyplot as plt + +plt.xkcd() +fig = plt.figure( figsize=(6,4) ) + +# the line: +slope = 2.0 +xx = np.arange(0.0, 4.1, 0.1) +yy = slope*xx +# the data: +n = 80 +rng = np.random.RandomState(218) +sigma = 1.5 +x = 4.0*rng.rand(n) +y = slope*x+rng.randn(n)*sigma +# fit: +slopef = np.sum(x*y)/np.sum(x*x) +yf = slopef*xx + +# plot it: +ax = fig.add_subplot( 1, 1, 1 ) +ax.spines['left'].set_position('zero') +ax.spines['bottom'].set_position('zero') +ax.spines['right'].set_visible(False) +ax.spines['top'].set_visible(False) +ax.get_xaxis().set_tick_params(direction='inout', length=10, width=2) +ax.get_yaxis().set_tick_params(direction='inout', length=10, width=2) +ax.yaxis.set_ticks_position('left') +ax.xaxis.set_ticks_position('bottom') +ax.set_xticks(np.arange(0.0, 4.1)) +ax.set_xlim(0.0, 4.2) +#ax.set_ylim(-1, 5) +#ax.set_xticks( np.arange(0, 5)) +#ax.set_yticks( np.arange(0, 0.9, 0.2)) +ax.set_xlabel('x') +ax.set_ylabel('y') +#ax.annotate('', xy=(mu, 0.02), xycoords='data', +# xytext=(mu, 0.75), textcoords='data', +# arrowprops=dict(arrowstyle="->", relpos=(0.5,0.5), +# connectionstyle=cs), zorder=1 ) +ax.scatter(x, y, label='data', s=50, zorder=10) +ax.plot(xx, yy, 'r', lw=6.0, color='#ff0000', label='original', zorder=5) +ax.plot(xx, yf, '--', lw=2.0, color='#ffcc00', label='fit', zorder=7) +ax.legend(loc='upper left', frameon=False) + +plt.tight_layout(); +plt.savefig('mlepropline.pdf') +#plt.show();