Sveučilište u Zagrebu
Fakultet elektrotehnike i računarstva
http://www.fer.unizg.hr/predmet/su
Ak. god. 2015./2016.
(c) 2015 Jan Šnajder
Verzija: 0.8 (2015-11-01)
import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
from numpy.random import normal
%pylab inline
Populating the interactive namespace from numpy and matplotlib
Vjerojatnost
Očekivanje, varijanca i kovarijanca
Statistička nezavisnost
Matrica kovarijacije
Teorijske razdiobe
Procjena parametara
Procjenitelj MLE
Procjenitelj MAP
$X$ je slučajna varijabla, $\{x_i\}$ su njezine vrijednosti
Pojednostavljenje notacije: $$P(X=x) \equiv P(x)$$
$P(x_i)\geq 0$, $\sum_i P(x_i)=1$
Distribucija (razdioba) vjerojatnosti
Zajednička (engl. joint) distribucija nad $\{X,Y\}$: $$P(X=x,Y=y)\equiv P(x,y)$$
Kontinuirana slučajna varijabla: funkcija gustoće vjerojatnosti (PDF):
(Marginalna vjerojatnost varijable $X$)
$\rho_{X,Y}\in[-1,+1]$
from scipy import stats
X = sp.random.random(100)
Y0 = sp.random.random(100)
noise = stats.norm.rvs(size=100)
Y1 = X + 0.2 * noise
Y2 = 3 * Y1
Y3 = -Y1
Y4 = 1 - (X - 0.5)**2 + 0.05 * noise
for Y in [Y0, Y1, Y2, Y3, Y4]:
plt.scatter(X,Y, label="r = %.3f" % stats.pearsonr(X, Y)[0])
plt.legend()
plt.show()
/usr/local/lib/python2.7/dist-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison if self._edgecolors == str('face'):
ili $$ P(X|Y) = P(X) \qquad \text{i} \qquad P(Y|X) = P(Y) $$
Znanje o ishodu varijable $Y$ ne utječe na vjerojatnost ishoda varijable $X$ (i obrnuto)
Za nezavisne varijable $X$ i $Y$ vrijedi:
Nezavisne varijable su nekorelirane, ali obrat općenito ne vrijedi: nelinarno zavisne varijable mogu imati nisku korelaciju
Varijable $X$ i $Y$ su uvjetno nezavisne uz danu varijablu Z, što označavamo kao $X\bot Y|Z$, akko
ili $$ P(X,Y|Z) = P(X|Z) P(Y|Z) $$
poznat ishod varijable $Z$, znanje o ishodu varijable $Y$ ne utječe na ishod varijable $X$ (i obrnuto)
$\mathbf{X} = (X_1,\dots,X_n)$ je $n$-dimenzijski slučajan vektor
Matrica kovarijacije $\Sigma$:
Simetrična matrica!
Ako su $X_1...X_n$ međusobno nezavisne, onda $\Sigma = \mathrm{diag}(\sigma_i^2)$
Ako $\sigma^2_i = \sigma^2$, onda $\Sigma = \sigma^2 \mathbf{I}$ (izotropna kovarijanca)
mu = 0.3
p = stats.bernoulli(mu)
xs = sp.array([0,1])
for x in xs:
plt.plot(x, p.pmf(x), 'bo', ms=8, label='bernoulli pmf')
plt.vlines(x, 0, p.pmf(x), colors='b', lw=5, alpha=0.5)
plt.xlim(xmin=-1, xmax=2)
plt.ylim(ymax=1)
plt.show()
/usr/local/lib/python2.7/dist-packages/matplotlib/collections.py:650: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison if self._edgecolors_original != str('face'):
X = p.rvs(size=100); X
array([0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0])
sp.mean(X)
0.29999999999999999
sp.var(X)
0.20999999999999994
xs = linspace(0,1)
plt.plot(xs, xs * (1-xs));
Varijabla koja poprima jednu (i samo jednu) od $K$ mogućih vrijednosti
$\mathbf{x}=(x_1,x_2,\dots,x_K)^\mathrm{T}$ je binaran vektor indikatorskih varijabli
Vjerojatnosti pojedinih vrijednosti: $\boldsymbol{\mu}=(\mu_1,\dots,\mu_K)^\mathrm{T}$, $\sum_k \mu_k=1$, $\mu_k\geq 0$
xs = sp.linspace(-5, 5)
for s in range(1, 5):
plt.plot(xs, stats.norm.pdf(xs, 0, s), label='$\sigma=%d$' % s)
plt.legend()
plt.show()
print stats.norm.cdf(1, 0, 1) - stats.norm.cdf(-1, 0, 1)
print stats.norm.cdf(2, 0, 1) - stats.norm.cdf(-2, 0, 1)
print stats.norm.cdf(3, 0, 1) - stats.norm.cdf(-3, 0, 1)
0.682689492137 0.954499736104 0.997300203937
p = stats.norm(loc=5, scale=3)
X = p.rvs(size=30); X
array([ 3.94363136, 9.6884801 , 6.73686697, 2.23918325, -0.22621627, 8.0346425 , 7.50790913, 8.16860545, 6.52408393, 5.5558834 , 3.70190598, 1.00539899, 5.33881857, 4.42657421, 8.26823883, 4.14257769, 10.56845245, 6.31164921, 8.72326689, 6.81670834, 9.6269925 , 9.05832561, 5.89763569, 8.214693 , 5.37577967, 2.88912765, 2.47697932, 2.19650056, 1.90620298, 1.09675323])
sp.mean(X)
5.5405217061730685
sp.var(X)
8.4165690464455274
$\mathbf{\Sigma}$ mora biti pozitivno definitna. Tada (1) matrica je nesingularna i ima inverz te (2) determinanta joj je pozitivna
Kvadratna forma: $\Delta^2 = (\mathbf{x}-\boldsymbol{\mu})^\mathrm{T}\mathbf{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})$ je Mahalanobisova udaljenost između $\mathbf{x}$ i $\boldsymbol{\mu}$.
mu = [0, 1]
covm = sp.array([[1, 1], [1, 3]])
p = stats.multivariate_normal(mu, covm)
print covm
[[1 1] [1 3]]
x = np.linspace(-2, 2)
y = np.linspace(-2, 2)
X, Y = np.meshgrid(x, y)
XY = np.dstack((X,Y))
plt.contour(X, Y, p.pdf(XY));
covm1 = sp.array([[1, 0], [0, 5]])
print covm1
[[1 0] [0 5]]
plt.contour(X, Y, stats.multivariate_normal.pdf(XY, mean=mu, cov=covm1 ));
covm2 = sp.array([[5, 0], [0, 5]])
print covm2
[[5 0] [0 5]]
plt.contour(X, Y, stats.multivariate_normal.pdf(XY, mean=mu, cov=covm2 ));
plt.contour(X, Y, stats.multivariate_normal.pdf(XY, mean=mu, cov=[[1,0],[0,1]] ));
from scipy import linalg
x00 = sp.array([0,0])
x01 = sp.array([0,1])
x10 = sp.array([1,0])
x11 = sp.array([1,1])
linalg.norm(x00 - x01, ord=2)
1.0
linalg.norm(x00 - x10, ord=2)
1.0
linalg.norm(x00 - x11, ord=2)
1.4142135623730951
sqrt(sp.dot((x00 - x11),(x00 - x11)))
1.4142135623730951
def mahalanobis(x1, x2, covm):
return sqrt(sp.dot(sp.dot((x1 - x2), linalg.inv(covm)), (x1 - x2)))
# ili: from scipy.spatial.distance import mahalanobis
covm1 = sp.array([[1, 0], [0, 5]])
mahalanobis(x00, x01, covm1)
0.44721359549995793
mahalanobis(x00, x10, covm1)
1.0
mahalanobis(x00, x11, covm1)
1.0954451150103321
mahalanobis(x00, x11, sp.eye(2))
1.4142135623730951
Ideja: na temelju slučajnog uzorka izračunati procjenu (estimaciju) parametra teorijske razdiobe
Neka je $(X_1,X_2,\dots,X_n)$ uzorak ($n$-torka slučajnih varijabli koje su iid)
Slučajna varijabla $\Theta=g(X_1,X_2,\dots,X_n)$ naziva se statistika
Statistika $\Theta$ je procjenitelj (estimator) parametra populacije $\theta$
Vrijednost procjenitelja $\hat{\theta} = g(x_1,x_2,\dots,x_n)$ naziva se procjena
Procjenitelj je slučajna varijable, dakle ima očekivanje i varijancu
[Slika: pristranost i varijanca procjenitelja]
Procjenitelj $\Theta$ je nepristran procjenitelj (engl. unbiased estimator) parametra $\theta$ akko
$X$ je slučajna varijabla sa $x\in\mathbb{R}$.
Označimo $\mathbb{E}[X] = \mu$ (srednja vrijednost) i
$\mathrm{Var}(X)=\sigma^2$ (varijanca)
$\mu$ i $\sigma^2$ su parametri populacije i oni su nam nepoznati
Parametre $\mu$ i $\sigma^2$ možemo ih procijeniti na temelju uzorka $\{x^{(i)}\}_{i=1}^N$ pomoću procjenitelja
Za procjenitelje možemo upotrijebiti bilo koje statistike. Npr.
Q: Jesu li ovo dobri procjenitelji? (Jesu li nepristrani?)
$\mathbb{E}[\hat{\mu}]=\mu$ ?
$\mathbb{E}[\hat{\sigma}^2] = \sigma^2$ ?
X = stats.norm.rvs(size=10, loc=0, scale=1) # mean=0, stdev=var=1
sp.mean(X)
0.042575400812075635
mean = 0
n = 10
N = 10000
for i in range(N):
X = stats.norm.rvs(size=n)
mean += sp.sum(X) / len(X)
mean / N
0.0014897662030209388
vrijednosti
nepristran** procjenitelj varijance! $$ \mathbb{E}[\hat{\sigma}^2] = \frac{N-1}{N}\sigma^2 $$
Procjenitelj podcjenjuje (engl. underestimates) pravu varijancu!
Nepristran procjenitelj varijance:
def st_dev(X):
n = len(X)
mean = sp.sum(X) / n
s = 0
for i in range(len(X)):
s += (X[i] - mean)**2
return s / n
X = stats.norm.rvs(size=10, loc=0, scale=1) # mean=0, stdev=var=1
st_dev(X)
0.62610770360872903
stdev = 0
n = 10
N = 10000
for i in range(N):
X = stats.norm.rvs(size=n)
stdev += st_dev(X)
stdev / N
0.90251698503102884
stdev = 0
n = 10
N = 10000
for i in range(N):
X = stats.norm.rvs(size=n)
stdev += st_dev(X)
stdev / N
0.90865744695623196
Kako izvesti procjenitelj za neku teorijsku distribuciju (Bernoullijevu, Gaussovu, ...)?
Tri vrste procjenitelja:
izvlačenje uzorka $\mathcal{D}$ čine najvjerojatnijim $$ p(\mathcal{D} | \boldsymbol{\theta}) = p(\mathbf{x}^{(1)},\dots,\mathbf{x}^{(N)} | \mathbf{\theta}) = \prod_{i=1}^N p(\mathbf{x}^{(i)} | \mathbf{\theta})\ \equiv \color{red}{\mathcal{L}(\boldsymbol{\theta} | \mathcal{D})} $$ NB: Druga jednakost vrijedi uz pretpostavku iid
Funkcija izglednosti $\mathcal{L} : \boldsymbol{\theta}\mapsto p(\mathcal{D} | \boldsymbol{\theta})$ parametrima pridjeljuje vjerojatnost
$\mathcal{L}$ nije PDF! Općenito ne vrijedi $\int_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}|\mathcal{D})\,\mathrm{d}\boldsymbol{\theta}=1$.
def likelihood(mu, m, N):
return mu**m * (1 - mu)**(N - m)
xs = linspace(0,1)
plt.plot(xs, likelihood(xs, 8, 10));
xs = linspace(0,1)
plt.plot(xs, likelihood(xs, 5, 10));
xs = linspace(0,1)
plt.plot(xs, likelihood(xs, 10, 10));
MLE za Bernoullijevu razdiobu je ustvari relativna frekvencija
Vrijedi $\mathbb{E}(\mu_\mathrm{ML})=\mathbb{E}[X]=\mu$, pa je ovo je nepristran procjenitelj
Izraz treba maksimizirati prema $\mu_k$ uz ograničenje $\sum_{k=1}^K\mu_k=1$.
Primjenom metode Lagrangeovih multiplikatora dobivamo:
$N_k$ je broj nastupanja k-te vrijednosti
NB: Procjenitelj $\hat{\sigma}^2_\mathrm{ML}$ je pristran!
MLE ne mora nužno biti nepristran!
p = stats.norm(5, 2)
X = sort(p.rvs(30))
plt.scatter(X, sp.zeros(len(X)));
mean_mle = sp.mean(X); mean_mle
4.8057053896461506
var_mle = np.var(X, axis=0, ddof=1); var_mle
3.1111795224611711
p_mle = stats.norm(mean_mle, sqrt(var_mle))
plt.scatter(X, p_mle.pdf(X))
plt.plot(X, p.pdf(X), c='gray');
plt.plot(X, p_mle.pdf(X), c='blue', linewidth=2)
plt.vlines(X, 0, p_mle.pdf(X), colors='b', lw=2, alpha=0.2)
plt.show()
mu = [3, 2]
covm = sp.array([[5, 2], [2, 10]])
p = stats.multivariate_normal(mu, covm)
x = np.linspace(-10, 10)
y = np.linspace(-10, 10)
X, Y = np.meshgrid(x, y)
XY = np.dstack((X,Y))
plt.contour(X, Y, p.pdf(XY))
plt.show()
D = p.rvs(100)
plt.contour(X, Y, p.pdf(XY), cmap='binary', alpha=0.5)
plt.scatter(D[:,0], D[:,1])
plt.show()
mean_mle = sp.mean(D, axis=0); mean_mle
array([ 3.10940042, 2.13853873])
cov_mle = 0
s = 0
for x in D:
s += sp.outer(x - mean_mle, x - mean_mle)
cov_mle = s / len(D)
cov_mle
array([[ 5.36556926, 1.61872328], [ 1.61872328, 9.05445152]])
sp.cov(D, rowvar=0, bias=0)
array([[ 5.41976693, 1.63507402], [ 1.63507402, 9.14591063]])
p_mle = stats.multivariate_normal(mean_mle, cov_mle)
plt.contour(X, Y, p_mle.pdf(XY));
plt.contour(X, Y, p.pdf(XY), cmap='binary', alpha=0.5)
plt.scatter(D[:,0], D[:,1], c='gray', alpha=0.5)
plt.contour(X, Y, p_mle.pdf(XY), cmap='Blues', linewidths=2);
MLE lako dovodi do prenaučenosti modela
Ideja: nisu sve vrijednosti za $\mu$ jednako vjerojatne!
Definiramo apriornu razdiobu parametra $p(\boldsymbol{\theta})$ i zatim maksimiziramo
aposteriornu vjerojatnost: $$ p(\boldsymbol{\theta}|\mathcal{D}) = \frac{p(\mathcal{D}|\boldsymbol{\theta}) P(\boldsymbol{\theta})} {p(\mathcal{D})} $$
TODO
xs = sp.linspace(0,1)
beta = stats.beta(1,1)
plt.plot(xs,stats.beta.pdf(xs,1,1), label='a=1,b=1')
plt.plot(xs,stats.beta.pdf(xs,2,2), label='a=2,b=2')
plt.plot(xs,stats.beta.pdf(xs,4,2), label='a=4,b=2')
plt.plot(xs,stats.beta.pdf(xs,2,4), label='a=2,b=4')
plt.legend()
plt.show()
TODO
Za strojno učenje posebno su važne Bernoullijeva, kategorička i Gaussova razdioba
Procjenitelj je statistika (slučajna varijabla izračunata iz uzorka) kojom se procjenjuju parametri neke teorijske razdiobe
Dobri procjenitelju su nepristrani
Procjenitelj najveće izglednosti (MLE) odabire parametre koji maksimiziraju vjerojatnost realizacije uzorka (tzv. izglednost)
MLE procjenitelj nije uvijek nepristran i sklon je prenaučenosti
MAP-procjenitelj dodatno koristi apriornu razdiobu parametara i maksimizira aposteriornu vjerojatnost parametara
MAP-procjenitelj na taj način ugrađuje apriorno znanje te izbjegava prenaučenost samo na temelju podataka