teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added maximum likelihood chapter

Browse Source
This commit is contained in:
Jan Benda 2015-10-24 23:05:44 +02:00
parent ad10b61ab7
commit 93089b4be2
5 changed files with 419 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
statistics/lecture/Makefile
View File

@ -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)

201
statistics/lecture/descriptivestatistics.tex
View File

@ -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}

99
statistics/lecture/mlemean.py Normal file
View File

@ -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();

70
statistics/lecture/mlepdf.py Normal file
View File

@ -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();

49
statistics/lecture/mlepropline.py Normal file
View File

@ -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();