Une version un peu « bricolée » du développement de la page 28.

In [50]: 
import sympy as sy
In [51]: 
x_sy, theta = sy.symbols('x_sy, theta', real=True, positive=True)
In [52]: 
F = 1 - sy.exp(-x_sy/theta)
In [53]: 
sy.init_printing()
In [54]: 
F
Out[54]:
$$1 - e^{- \frac{x_{sy}}{\theta}}$$
In [55]: 
 sy.factor(F.subs(x_sy,x_sy+sy.Rational(1,2))-F.subs(x_sy,x_sy-sy.Rational(1,2)))
Out[55]:
$$- \frac{1}{e^{\frac{x_{sy}}{\theta}}} \left(- e^{\frac{1}{2 \theta}} + e^{- \frac{1}{2 \theta}}\right)$$
In [56]: 
L=sy.prod([(sy.exp(sy.Rational(1,2)/theta)-sy.exp(-sy.Rational(1,2)/theta))*sy.exp(-x/theta) for x in sy.symbols('x:10',positive=True,real=True)]);L
Out[56]:
$$\frac{\left(e^{\frac{1}{2 \theta}} - e^{- \frac{1}{2 \theta}}\right)^{10}}{e^{\frac{x_{0}}{\theta}} e^{\frac{x_{1}}{\theta}} e^{\frac{x_{2}}{\theta}} e^{\frac{x_{3}}{\theta}} e^{\frac{x_{4}}{\theta}} e^{\frac{x_{5}}{\theta}} e^{\frac{x_{6}}{\theta}} e^{\frac{x_{7}}{\theta}} e^{\frac{x_{8}}{\theta}} e^{\frac{x_{9}}{\theta}}}$$
In [57]: 
sy.expand_log(sy.log(L),force=True)
Out[57]:
$$10 \log{\left (e^{\frac{1}{2 \theta}} - e^{- \frac{1}{2 \theta}} \right )} - \frac{x_{0}}{\theta} - \frac{x_{1}}{\theta} - \frac{x_{2}}{\theta} - \frac{x_{3}}{\theta} - \frac{x_{4}}{\theta} - \frac{x_{5}}{\theta} - \frac{x_{6}}{\theta} - \frac{x_{7}}{\theta} - \frac{x_{8}}{\theta} - \frac{x_{9}}{\theta}$$
In [58]: 
l = sy.collect(sy.expand_log(sy.log(L),force=True),1/theta); l
Out[58]:
$$10 \log{\left (e^{\frac{1}{2 \theta}} - e^{- \frac{1}{2 \theta}} \right )} + \frac{1}{\theta} \left(- x_{0} - x_{1} - x_{2} - x_{3} - x_{4} - x_{5} - x_{6} - x_{7} - x_{8} - x_{9}\right)$$
In [59]: 
sy.solve(sy.diff(l,theta),theta)
Out[59]:
$$\begin{bmatrix}\frac{1}{\log{\left (\frac{x_{0} + x_{1} + x_{2} + x_{3} + x_{4} + x_{5} + x_{6} + x_{7} + x_{8} + x_{9} + 5}{x_{0} + x_{1} + x_{2} + x_{3} + x_{4} + x_{5} + x_{6} + x_{7} + x_{8} + x_{9} - 5} \right )}}\end{bmatrix}$$

Retour à l'example 9.4.1 page 26

Nous chargeons les modules « indispensables » : numpy, scipy et matplotlib avec la commandes magique %pylab

In [60]: 
%pylab

Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

Nous assignons à une variable donnees une liste avec les données :

In [61]: 
donnees = [70,11,66,5,20,4,35,40,29,8]

La valeur moyenne est obtenue avec la fonction mean de scipy :

In [62]: 
mean(donnees)
Out[62]:
$$28.8$$

Nous définissons la fonction de vraisemblance (symbolique) :

In [63]: 
L = sy.prod([sy.exp(-x/theta)/theta for x in sy.symbols('x:10',real=True,positive=True)]); L
Out[63]:
$$\frac{1}{\theta^{10} e^{\frac{x_{0}}{\theta}} e^{\frac{x_{1}}{\theta}} e^{\frac{x_{2}}{\theta}} e^{\frac{x_{3}}{\theta}} e^{\frac{x_{4}}{\theta}} e^{\frac{x_{5}}{\theta}} e^{\frac{x_{6}}{\theta}} e^{\frac{x_{7}}{\theta}} e^{\frac{x_{8}}{\theta}} e^{\frac{x_{9}}{\theta}}}$$

Nous obtenons ensuite la log-vraisemblance :

In [64]: 
l = sy.collect(sy.expand_log(sy.log(L),force=True),1/theta); l
Out[64]:
$$- 10 \log{\left (\theta \right )} + \frac{1}{\theta} \left(- x_{0} - x_{1} - x_{2} - x_{3} - x_{4} - x_{5} - x_{6} - x_{7} - x_{8} - x_{9}\right)$$
In [65]: 
x_dot = sy.symbols('x_dot',real=True,positive=True)
l = l.subs(sum([x for x in sy.symbols('x:10',real=True,positive=True)]), x_dot)
l
Out[65]:
$$- 10 \log{\left (\theta \right )} - \frac{x_{dot}}{\theta}$$

Nous définissons la fonction score :

In [66]: 
S = sy.diff(l,theta)
S.subs(sum([x for x in sy.symbols('x:10',real=True,positive=True)]), x_dot)
Out[66]:
$$- \frac{10}{\theta} + \frac{x_{dot}}{\theta^{2}}$$

Nous définissons la fonction d'information :

In [67]: 
I = - sy.diff(S,theta)
I.subs(sum([x for x in sy.symbols('x:10',real=True,positive=True)]), x_dot)
Out[67]:
$$- \frac{10}{\theta^{2}} + \frac{2 x_{dot}}{\theta^{3}}$$

Nous obtenons l'EMV comme l'argument qui annulle la fonction score :

In [71]: 
theta_chapeau = sy.solve(S,theta)[0]
theta_chapeau
Out[71]:
$$\frac{x_{dot}}{10}$$

Nous vérifions que l'information est positive en theta_chapeau :

In [72]: 
I.subs(theta,theta_chapeau)
Out[72]:
$$\frac{1000}{x_{dot}^{2}}$$

Nous définissons une version numérique de la fonction de vraisemblance :

In [89]: 
def l_fct(x, x_dot = sum(donnees), n=len(donnees)):
    return -n*log(x) - float(x_dot)/x
In [90]: 
l_fct(sum(donnees)/float(len(donnees)))
Out[90]:
$$-43.6037538714$$

Nous définissons la fonction de log-vraisemblance relative :

In [91]: 
def r(x, l_max = l_fct(sum(donnees)/float(len(donnees)))):
    return l_fct(x) - l_max

Nous la traçons après avoir utiliser une autre commande magique qui nous intégre les figure dans notre carnet de notes :

In [95]: 
%matplotlib inline
In [96]: 
xx = arange(10,100)
plot(xx,r(xx))
grid()
ylim([-5,0])
axhline(y=log(0.5),color='red')
axhline(y=log(0.1),color='red')
axhline(y=log(0.01),color='red')
Out[96]:
<matplotlib.lines.Line2D at 0x7f50538a6a10>
In []: