note_1_4
Заметка 4. Классы в питоне и собственное распределение.
курса Введение в машинное обучение.
Шокуров Антон В.
shokurov.anton.v@yandex.ru
http://машинноезрение.рф
Версия 0.11

Аннотация

В данной заметке рассматривются такие понятия как классы и иные библиотек относящихся к анализу данных.

Это предварительная версия! Любые замечания приветсвуются.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

import scipy.stats as models

Классы в Python

In [2]:
a = np.random.normal(size=100)
In [3]:
a_mean = np.mean( a )
In [4]:
a_std = np.std( a )
In [5]:
a_mean, a_std
Out[5]:
(-0.0957556268072998, 1.0141719112419765)
In [6]:
b = np.random.normal(size=100)
In [7]:
b_mean = np.mean( b )
b_std = np.std( b )
b_mean, b_std
Out[7]:
(-0.013992204252215933, 1.0328967540508232)
In [8]:
class stats:
    def calc( self, x ):
        self.mean = np.mean( x )
        self.std = np.std( x )
In [9]:
first = stats()
stats.calc( first, a )
In [10]:
print( first.mean, first.std ) # Обе статистики теперь вместе.
-0.0957556268072998 1.0141719112419765
In [11]:
second = stats()
stats.calc( second, b )
In [12]:
print( second.mean, second.std ) # Обе статистики теперь вместе.
-0.013992204252215933 1.0328967540508232

Классы на то и классы, что вызов можно сократить:

In [13]:
first.calc( a )
second.calc( b )

Печать тоже можно сделать методом

In [14]:
class stats:
    def calc( self, x ):
        self.mean = np.mean( x )
        self.std = np.std( x )
        
    def pri( self ):
        print( self.mean, self.std)
In [15]:
first.pri()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-15-a9d512d727bb> in <module>
----> 1 first.pri()

AttributeError: 'stats' object has no attribute 'pri'
In [16]:
first = stats()
stats.calc( first, a )
first.pri()
-0.0957556268072998 1.0141719112419765

Добавление элемента в объект классв.

In [17]:
first.qq = 3 # Добавили переменную в класс.
In [18]:
first.qq
Out[18]:
3
In [19]:
second.qq
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-19-589ee73ebfb8> in <module>
----> 1 second.qq

AttributeError: 'stats' object has no attribute 'qq'
In [20]:
third = stats()
In [21]:
third.pri() # Отсутвуют члены класса, так как они не были ещё созданы.
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-21-062232672ed7> in <module>
----> 1 third.pri() # Отсутвуют члены класса, так как они не были ещё созданы.

<ipython-input-14-e071d05fdfcb> in pri(self)
      5 
      6     def pri( self ):
----> 7         print( self.mean, self.std)

AttributeError: 'stats' object has no attribute 'mean'
In [22]:
stats.__dict__
Out[22]:
mappingproxy({'__module__': '__main__',
              'calc': <function __main__.stats.calc(self, x)>,
              'pri': <function __main__.stats.pri(self)>,
              '__dict__': <attribute '__dict__' of 'stats' objects>,
              '__weakref__': <attribute '__weakref__' of 'stats' objects>,
              '__doc__': None})
In [23]:
first.__dict__
Out[23]:
{'mean': -0.0957556268072998, 'std': 1.0141719112419765, 'qq': 3}
In [24]:
second.__dict__
Out[24]:
{'mean': -0.013992204252215933, 'std': 1.0328967540508232}
In [25]:
from types import MethodType
In [26]:
def get_mean( self ):
    return self.mean
In [27]:
first.getmm=get_mean
In [28]:
first.getmm()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-28-cf0089ee5115> in <module>
----> 1 first.getmm()

TypeError: get_mean() missing 1 required positional argument: 'self'
In [29]:
first.__dict__
Out[29]:
{'mean': -0.0957556268072998,
 'std': 1.0141719112419765,
 'qq': 3,
 'getmm': <function __main__.get_mean(self)>}
In [30]:
first.getmm(first)
Out[30]:
-0.0957556268072998
In [31]:
first.getmm = MethodType(get_mean, first)
In [32]:
first.getmm()
Out[32]:
-0.0957556268072998
In [33]:
first.__dict__
Out[33]:
{'mean': -0.0957556268072998,
 'std': 1.0141719112419765,
 'qq': 3,
 'getmm': <bound method get_mean of <__main__.stats object at 0x7fe75db8c9b0>>}
In [34]:
second.getmm()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-34-fe7cb4023c2d> in <module>
----> 1 second.getmm()

AttributeError: 'stats' object has no attribute 'getmm'
In [35]:
stats.get_mm = MethodType(get_mean, stats)
In [36]:
first.__dict__
Out[36]:
{'mean': -0.0957556268072998,
 'std': 1.0141719112419765,
 'qq': 3,
 'getmm': <bound method get_mean of <__main__.stats object at 0x7fe75db8c9b0>>}
In [37]:
first.get_mm()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-37-73bbf147423b> in <module>
----> 1 first.get_mm()

<ipython-input-26-732b4fea4940> in get_mean(self)
      1 def get_mean( self ):
----> 2     return self.mean

AttributeError: type object 'stats' has no attribute 'mean'
In [38]:
first.getmm()
Out[38]:
-0.0957556268072998

конструктор

In [39]:
class stats:
    def __init__( self ):# Конструктор
        print('first call!')
        self.mean = 0
        self.std = 0

    def calc( self, x ):
        self.mean = np.mean( x )
        self.std = np.std( x )
        
    def pri( self ):
        print( self.mean, self.std)
In [40]:
forth = stats()
forth.pri()
first call!
0 0
In [41]:
class stats:    
    def __init__( self, x = None ):# Конструкотор, расчитан на два случая
        if x is None:
            print('simple call!')
            self.mean = 0
            self.std = 0
            return
        print('data call!')
        self.calc( x )

    def calc( self, x ):
        self.mean = np.mean( x )
        self.std = np.std( x )
        
    def pri( self ):
        print( self.mean, self.std)
In [42]:
fifth = stats()
fifth.pri()
simple call!
0 0
In [43]:
sixth = stats( a )
sixth.pri()
data call!
-0.0957556268072998 1.0141719112419765

Переменная класса

Важно, что ранее объявленные переменные были имено частью self, а не класса. Иначе будет следующий эффект

In [44]:
class stats:
    all_means = []
    def __init__( self, x = None ):# Конструкотор, расчитан на два случая
        if x is None:
            print('simple call!')
            self.mean = 0
            self.std = 0
            return
        print('data call!')
        self.calc( x )

    def calc( self, x ):
        self.mean = np.mean( x )
        self.std = np.std( x )
        self.all_means.append( self.mean )

    def pri( self ):
        print( self.mean, self.std)
In [45]:
first = stats( a )
second = stats( b )
print( first.all_means, second.all_means ) # Переменная all_means одна и таже для всех объектов.
data call!
data call!
[-0.0957556268072998, -0.013992204252215933] [-0.0957556268072998, -0.013992204252215933]

Наследование

In [46]:
class collect( stats ):
    def __init__(self):
        self.data = []
        
    def add( self, x ):
        self.data.append( x )
    
    def calc_stats( self ):
        self.calc( self.data )
        #super( collect, self).calc( self.data )#Можно и так.
In [47]:
acum = collect()
In [48]:
acum.add( 1. )
acum.add( 2. )
In [49]:
acum.data
Out[49]:
[1.0, 2.0]
In [50]:
acum.calc_stats()
acum.pri()
1.5 0.5
In [51]:
acum.add( -2. )
acum.calc_stats()
acum.pri()
0.3333333333333333 1.699673171197595

Подмена вызова

In [52]:
class stats:
    all_means = []
    def __init__( self, x = None ):# Конструкотор, расчитан на два случая
        if x is None:
            print('simple call!')
            self.mean = 0
            self.std = 0
            return
        print('data call!')
        self.calc( x )
    
    def filt( self, x):
        return x

    def calc( self, x ):
        x = self.filt( x ) # Использует метод filt
        self.mean = np.mean( x )
        self.std = np.std( x )
        self.all_means.append( self.mean )

    def pri( self ):
        print( self.mean, self.std)
In [53]:
class collect( stats ):
    def __init__(self):
        self.data = []
    
    def filt( self, x): # Новый фильтр. Он подменяет из родительского.
        xx = [ qq for qq in x if qq > 0]
        return xx
    
    def add( self, x ):
        self.data.append( x )
    
    def calc_stats( self ):
        self.calc( self.data )
In [54]:
acum_pos = collect()
In [55]:
acum_pos.add( 1. )
acum_pos.add( 2. )
acum_pos.add( -2. )
In [56]:
acum_pos.calc_stats() # При подсчете статистик будет вызван новый фильтр.
acum_pos.pri() # Учитывались только положительные!
1.5 0.5
In [ ]:
 

Собственное распределение

In [57]:
x = models.norm.rvs(2, 3, size = 100)
In [58]:
models.norm.fit( x )
Out[58]:
(1.6837935086418192, 3.062994167934736)
In [59]:
np.mean( x )
Out[59]:
1.6837935086418192
In [60]:
np.std( x ) # Вот и все, т.е. по формулам.
Out[60]:
3.062994167934736
In [61]:
pi = 0.3
lmbd = 3
In [62]:
def zip_pmf( x, pi = pi, lmbd = lmbd ):
    if pi < 0 or pi > 1 or lmbd <= 0:
        return np.zeros_like( x )
    return (x == 0) * pi + (1 - pi) * models.poisson.pmf( x, lmbd)
In [63]:
zip_pmf( 0 )
Out[63]:
0.33485094785750474
In [64]:
zip_pmf( 1 )
Out[64]:
0.10455284357251429
In [65]:
models.poisson.pmf( 0, lmbd)
Out[65]:
0.049787068367863944
In [66]:
fig, ax = plt.subplots(figsize=(8, 6))
xs = np.arange(0, 10)

palette = ['g','b']
ax.bar( 2.5*xs, models.poisson.pmf( xs, lmbd), width = 0.9, color = palette[0], label='Poisson')
ax.bar( 2.5*xs + 1, zip_pmf( xs ), width = 0.9, color = palette[1], label='Zero - Poisson')

ax.set_xticks( 2.5 * xs + 1 -0.6*0.9)
ax.set_xticklabels( xs )
ax.set_xlabel('$x$')

ax.set_ylabel('$P(X=x)$')
ax.legend()
Out[66]:
<matplotlib.legend.Legend at 0x7fe75b79d438>
In [67]:
N = 1000
iz = models.bernoulli.rvs( pi, size = N)
x = (1 - iz) * models.poisson.rvs( lmbd, size = N)
print(x[::50])
x.shape
[3 3 5 3 5 1 0 0 0 3 0 2 3 0 1 5 0 3 3 3]
Out[67]:
(1000,)
In [68]:
fig, ax = plt.subplots(figsize=(8, 6))
ax.hist( x, width=0.8,  bins = np.arange( x.max() + 1 ), normed = True)

ax.set_xticks( np.arange( x.max() + 1 ) + 0.4 )
ax.set_xticklabels( np.arange( x.max() + 1 ) )
ax.set_xlabel('$x$')

ax.set_ylabel('Частота')
/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "
Out[68]:
Text(0,0.5,'Частота')
In [69]:
from statsmodels.base.model import GenericLikelihoodModel

class ZeroInflated(GenericLikelihoodModel):
    def __init__(self, endog, exog = None, **kwds):
        if exog == None:
            exog = np.zeros_like( endog )
        print("init !!")
        super( ZeroInflated, self).__init__( endog, exog, **kwds)
    
    def nloglikeobs( self, params):
        pi = params[0]
        lmbd = params[1]
        value = -np.log( zip_pmf( self.endog, pi = pi, lmbd = lmbd) )
        print("log, ", np.sum(value))
        return value
    
    def fit(self, start_params = None, maxiter = 10000, maxfun = 5000, **kwds):
        if start_params is None:
            lmbd_start = self.endog.mean()
            extra_zeros = ( self.endog == 0 ).mean() - models.poisson.pmf(0, lmbd_start )
            start_params = np.array([extra_zeros, lmbd_start])
            
        return super( ZeroInflated, self ).fit(start_params = start_params, maxiter = maxiter, maxfun = maxfun, **kwds)
     
In [70]:
model = ZeroInflated( x )
init !!
In [71]:
#model.fit()
ZeroInflated.fit(model)
log,  1936.0033536207616
log,  1935.1948959409383
log,  1914.6210026045528
log,  1913.3734180367712
log,  1903.4805482814577
log,  1889.1880757580072
log,  1875.3544419883983
log,  1867.626353528197
log,  1855.9819723635474
log,  1857.168579827365
log,  1858.8874815663771
log,  1853.9895642333402
log,  1851.7150738992132
log,  1854.1989648481651
log,  1852.2373701507156
log,  1848.8001724846563
log,  1851.1606879069332
log,  1864.8925049985437
log,  1849.236144820666
log,  1847.987329911325
log,  1849.776874965937
log,  1847.2203342522714
log,  1847.655006445384
log,  1848.0792038482177
log,  1847.3382613047813
log,  1847.1882105013992
log,  1848.6169775963824
log,  1849.7467545079262
log,  1846.9085396092837
log,  1846.8957346869704
log,  1847.1303010103306
log,  1847.181329050878
log,  1846.9323610367828
log,  1846.9329806779233
log,  1846.8696770621464
log,  1847.022201803337
log,  1846.856362402773
log,  1847.0112580523155
log,  1846.8550000115424
log,  1846.8553375302968
log,  1846.8947803583485
log,  1846.848634465715
log,  1846.860466900547
log,  1846.8506095971652
log,  1846.8526559800982
log,  1846.8494182410702
log,  1846.8527188612193
log,  1846.848913273402
log,  1846.8491700343639
log,  1846.8486155992282
log,  1846.8497343236818
log,  1846.8485171631864
log,  1846.8492115961817
log,  1846.8484826260474
log,  1846.8485151597913
log,  1846.8486019562815
log,  1846.8484642990716
log,  1846.8485734324443
log,  1846.8484638268933
log,  1846.8485611258432
log,  1846.8484560334552
log,  1846.8484735236696
log,  1846.8484565186363
log,  1846.8484639439318
log,  1846.848456432212
log,  1846.8484594195666
log,  1846.8484545155802
log,  1846.8484581051703
log,  1846.8484544886471
log,  1846.8484608924514
log,  1846.8484542230276
log,  1846.8484550006056
log,  1846.8484541805365
Optimization terminated successfully.
         Current function value: 1.846848
         Iterations: 38
         Function evaluations: 73
log,  1846.8484629010156
log,  1846.8484541805365
log,  1846.8484541805365
log,  1846.8484653004666
log,  1846.8484652681118
log,  1846.8484747825114
log,  1846.848470397777
log,  1846.848472050082
log,  1846.8485044462695
log,  1846.8484541805365
log,  1846.8484541805365
log,  1846.8485156270704
Out[71]:
<statsmodels.base.model.GenericLikelihoodModelResults at 0x7fe75b6e2e10>
In [72]:
res = model.fit()
log,  1936.0033536207616
log,  1935.1948959409383
log,  1914.6210026045528
log,  1913.3734180367712
log,  1903.4805482814577
log,  1889.1880757580072
log,  1875.3544419883983
log,  1867.626353528197
log,  1855.9819723635474
log,  1857.168579827365
log,  1858.8874815663771
log,  1853.9895642333402
log,  1851.7150738992132
log,  1854.1989648481651
log,  1852.2373701507156
log,  1848.8001724846563
log,  1851.1606879069332
log,  1864.8925049985437
log,  1849.236144820666
log,  1847.987329911325
log,  1849.776874965937
log,  1847.2203342522714
log,  1847.655006445384
log,  1848.0792038482177
log,  1847.3382613047813
log,  1847.1882105013992
log,  1848.6169775963824
log,  1849.7467545079262
log,  1846.9085396092837
log,  1846.8957346869704
log,  1847.1303010103306
log,  1847.181329050878
log,  1846.9323610367828
log,  1846.9329806779233
log,  1846.8696770621464
log,  1847.022201803337
log,  1846.856362402773
log,  1847.0112580523155
log,  1846.8550000115424
log,  1846.8553375302968
log,  1846.8947803583485
log,  1846.848634465715
log,  1846.860466900547
log,  1846.8506095971652
log,  1846.8526559800982
log,  1846.8494182410702
log,  1846.8527188612193
log,  1846.848913273402
log,  1846.8491700343639
log,  1846.8486155992282
log,  1846.8497343236818
log,  1846.8485171631864
log,  1846.8492115961817
log,  1846.8484826260474
log,  1846.8485151597913
log,  1846.8486019562815
log,  1846.8484642990716
log,  1846.8485734324443
log,  1846.8484638268933
log,  1846.8485611258432
log,  1846.8484560334552
log,  1846.8484735236696
log,  1846.8484565186363
log,  1846.8484639439318
log,  1846.848456432212
log,  1846.8484594195666
log,  1846.8484545155802
log,  1846.8484581051703
log,  1846.8484544886471
log,  1846.8484608924514
log,  1846.8484542230276
log,  1846.8484550006056
log,  1846.8484541805365
Optimization terminated successfully.
         Current function value: 1.846848
         Iterations: 38
         Function evaluations: 73
log,  1846.8484629010156
log,  1846.8484541805365
log,  1846.8484541805365
log,  1846.8484653004666
log,  1846.8484652681118
log,  1846.8484747825114
log,  1846.848470397777
log,  1846.848472050082
log,  1846.8485044462695
log,  1846.8484541805365
log,  1846.8484541805365
log,  1846.8485156270704
In [73]:
res.params # Profit!
Out[73]:
array([0.28261887, 3.00954076])

Измерение соотношения количеств

Используем распределение Бернули

Задается вероятностью p происхождения события. Например, вероятность кликнуть на банер, зайти на сайт, прийти в магазин или купить товар. Это все разные вероятности. Важное, что некое событие может произойти и неоднократно.

In [74]:
models.bernoulli.rvs( 0.7, 10) # Так нельзя, так как 10 будет использоватся как сдвиг.
Out[74]:
11

В питоне функции могут использовать два типа аргументов или иначе параметров, которые ей передаются. Первый тип позиционый, который является стандартным для большинства языков программиования. Суть в том, что параметры присваиваются в порядке и задания при вызове. А бывают аргументы именные, которые можно присвоить вне очереди.

Рассмотреть документацию по данной функции https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.bernoulli.html

Конкретно изучим описание функции rvs: rvs(p, loc=0, size=1, random_state=None)

Последнее к сожалению позволит понять какие значения были присвоены аргументам по умолчанию. часть аргументов просто перечислена (а данном случае только p), а часть через знак '='.

In [75]:
# Вообще у любого объекта, функции по соглашению может быть поле __doc__ которое дает краткую документацию по нему.
# Так как это поле является строчкой, то для аккуратного вывода лучше применить функцию print
# Подробнее о ней будет позже. Сейчас она используется для вывода самоописания функции.
print( models.bernoulli.rvs.__doc__ )
        Random variates of given type.

        Parameters
        ----------
        arg1, arg2, arg3,... : array_like
            The shape parameter(s) for the distribution (see docstring of the
            instance object for more information).
        loc : array_like, optional
            Location parameter (default=0).
        size : int or tuple of ints, optional
            Defining number of random variates (Default is 1).  Note that `size`
            has to be given as keyword, not as positional argument.
        random_state : None or int or ``np.random.RandomState`` instance, optional
            If int or RandomState, use it for drawing the random variates.
            If None, rely on ``self.random_state``.
            Default is None.

        Returns
        -------
        rvs : ndarray or scalar
            Random variates of given `size`.

        

Из данного описания следует, что сначала идут позиционные параметры (arg1, arg2, ...), а потом именные (loc, size, random_state). Тем не менее, если они идут подряд, то присвоение идет подряд.

In [76]:
# В данном случае:
models.bernoulli.rvs( 0.7, 10) # 0.7 это параметр распределения, а 10 присвоилось аргументу loc.
Out[76]:
11
In [77]:
models.bernoulli.rvs( 0.7, 10, 5) # Вот теперь очередь дошла до размера, которому было присвоено значение 5.
Out[77]:
array([10, 10, 11, 11, 11])
In [78]:
# Но можно было бы и явно присвоить переменной отвечащей за размер размер:
models.bernoulli.rvs( 0.7, size = 10 ) # Аргументу size присвоено значение вне очереди.
Out[78]:
array([1, 1, 1, 1, 0, 1, 0, 1, 1, 1])
In [79]:
# Нельзя позиционный аргумент подавать после именных. Последнее выдаст соответствующую ошибку:
models.bernoulli.rvs( size = 10, 0.7 )
  File "<ipython-input-79-5776a8475623>", line 2
    models.bernoulli.rvs( size = 10, 0.7 )
                                    ^
SyntaxError: positional argument follows keyword argument
In [80]:
# Но если аргументы заданы по имени, то их можно подавать функции в произвольном порядке.
models.bernoulli.rvs( size = 10, p = 0.7 )
Out[80]:
array([1, 1, 1, 1, 0, 1, 1, 1, 1, 1])
In [81]:
models.bernoulli.rvs( size = 10, p = 0.7, loc = 5 ) # В данном случае было задано и значение параметра loc, положения.
Out[81]:
array([5, 5, 6, 5, 6, 6, 5, 6, 6, 5])

Закон больших чисел

Нам нужно взять среднее количество произошедших событий. Точнее понять как это среднее рапределено. Теория говрит, что оно должно стремится к нормальному распределению.

In [82]:
# Имитация закона больших чисел:
data = [] # Содали пустой список
for i in range( 1000 ): # Мы явно указали размер нешей выборки
    lap = models.bernoulli.rvs( 0.7, size=100 ) # а также количество элементов в устреднение.
    data.append( np.sum( lap )/100 )
In [83]:
plt.hist( data );
In [84]:
# Хотим цикл ещё раз запустить с другим количеством эллементов в сумме.

Продолжаем с Бернули

In [85]:
# Вводим понятие функции, т.е. сами пишем.
def buildBernSum( p, cycles = 1000, elems = 100 ): # Два аргумента. Оба необязательны. Имеют значение по умолчанию.
    data = [] # Содали пустой список
    for i in range( cycles ):
        lap = models.bernoulli.rvs( p, size=elems ) # 0.7
        data.append( np.sum( lap )/elems )
    return data # Вернули результат.
In [86]:
bern = buildBernSum( 0.7 )
plt.hist( bern );
In [87]:
bern = buildBernSum( 0.7, elems = 10000 )
plt.hist( bern );
In [88]:
bern = buildBernSum( elems = 10000, p = 0.7)
plt.hist( bern );
In [89]:
bern = buildBernSum( 0.7, 100000, 10000 ) # Позиционное присвоение параметров
plt.hist( bern );
In [90]:
bern = buildBernSum( 0.7, 4000, 3000 )
plt.hist( bern );

Проверка того, что распределение Нормально

In [91]:
#С увеличением параметров распределение будет стремится к нормальному

Anderson-Darling

Предыдущая статистика работает только для нормального распределения. Для данного же теста можно задать используемое распределение: norm, expon, logistic, gumbel, gumbel_l, gumbel_r, extreme1. Пока будем оставатся в рамках нормального.

В Шапиро-Уилка граница была единственной. Той что обычно принято пользоваться. В методе же Андерсона таких границ считается несколько для различных уровней значимости.

In [92]:
dat = models.uniform.rvs( size = 100 ) # Пвторно проверим на равномером распределении.
plt.hist( dat  );
models.anderson( dat, dist='norm' ) # Тест отвергает данные как принадлежащие равномерному распределению.
# Статистику нужно возвести в квадрат и сравнить с критерием. Если больше, то гипотеза отвергается.
Out[92]:
AndersonResult(statistic=1.4262844372071442, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))

Проверка того что распределение то что нужно.

In [93]:
dat = models.norm.rvs( 2, 3, 100 ) # Для начала прверим для нормального распределения.
res = models.anderson( dat )
In [94]:
res
Out[94]:
AndersonResult(statistic=0.18360768648014414, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
In [95]:
res[1] # Значение значимости критических уровней.
Out[95]:
array([0.555, 0.632, 0.759, 0.885, 1.053])
In [96]:
res[2] # Сами уровни.
Out[96]:
array([15. , 10. ,  5. ,  2.5,  1. ])
In [97]:
# Можно представить их матрицей.
table = np.array([res[1], res[2]])
table
Out[97]:
array([[ 0.555,  0.632,  0.759,  0.885,  1.053],
       [15.   , 10.   ,  5.   ,  2.5  ,  1.   ]])
In [98]:
# Хотим красивое представление результата/критерия в виде таблицы
import pandas as pd # Очередная библиотека полезная для представления данных в таблицах, и их считывания.
In [99]:
pd.DataFrame( table, columns=['Уровни', 'Значимость'] ) # Система ругнется так как матрица не правильного размера.
# Данная функция ожидала матрицу состоящую из двух столбцов. 
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/pandas/core/internals/managers.py in create_block_manager_from_blocks(blocks, axes)
   1650                 blocks = [make_block(values=blocks[0],
-> 1651                                      placement=slice(0, len(axes[0])))]
   1652 

/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/pandas/core/internals/blocks.py in make_block(values, placement, klass, ndim, dtype, fastpath)
   3088 
-> 3089     return klass(values, ndim=ndim, placement=placement)
   3090 

/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/pandas/core/internals/blocks.py in __init__(self, values, placement, ndim)
     86                 'Wrong number of items passed {val}, placement implies '
---> 87                 '{mgr}'.format(val=len(self.values), mgr=len(self.mgr_locs)))
     88 

ValueError: Wrong number of items passed 5, placement implies 2

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-99-bb9297dfa31f> in <module>
----> 1 pd.DataFrame( table, columns=['Уровни', 'Значимость'] ) # Система ругнется так как матрица не правильного размера.
      2 # Данная функция ожидала матрицу состоящую из двух столбцов.

/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/pandas/core/frame.py in __init__(self, data, index, columns, dtype, copy)
    420             else:
    421                 mgr = init_ndarray(data, index, columns, dtype=dtype,
--> 422                                    copy=copy)
    423 
    424         # For data is list-like, or Iterable (will consume into list)

/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/pandas/core/internals/construction.py in init_ndarray(values, index, columns, dtype, copy)
    165         values = maybe_infer_to_datetimelike(values)
    166 
--> 167     return create_block_manager_from_blocks([values], [columns, index])
    168 
    169 

/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/pandas/core/internals/managers.py in create_block_manager_from_blocks(blocks, axes)
   1658         blocks = [getattr(b, 'values', b) for b in blocks]
   1659         tot_items = sum(b.shape[0] for b in blocks)
-> 1660         construction_error(tot_items, blocks[0].shape[1:], axes, e)
   1661 
   1662 

/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/pandas/core/internals/managers.py in construction_error(tot_items, block_shape, axes, e)
   1681         raise ValueError("Empty data passed with indices specified.")
   1682     raise ValueError("Shape of passed values is {0}, indices imply {1}".format(
-> 1683         passed, implied))
   1684 
   1685 

ValueError: Shape of passed values is (5, 2), indices imply (2, 2)
In [100]:
# В Numpy естественно есть функция для транспонирования матрицы.
np.transpose( table )
Out[100]:
array([[ 0.555, 15.   ],
       [ 0.632, 10.   ],
       [ 0.759,  5.   ],
       [ 0.885,  2.5  ],
       [ 1.053,  1.   ]])
In [101]:
# Транспонироваие дает результат и в общем случае: она переставляет оси.
x = np.ones((1, 2, 3)) # ones вспомогательная функция составляющая многомерную матрицу из единиц. В скобках указан размер.
x
Out[101]:
array([[[1., 1., 1.],
        [1., 1., 1.]]])
In [102]:
xT = np.transpose(x, (1, 2, 0)) # В скобках указана перестановка исходных осей, т.е. оси будут идти в порядке 1, 2, 0.
xT
Out[102]:
array([[[1.],
        [1.],
        [1.]],

       [[1.],
        [1.],
        [1.]]])
In [103]:
x.shape, xT.shape # Если оси переставлены, то мы это увидем на соответствующем образом изменившимся порядке их размере.
# Вначале должна идти первая компонента, потом вторая, а уже потом нулевая векто рараеров (1, 2, 3).
# Тогда получим 2, 3, 1
Out[103]:
((1, 2, 3), (2, 3, 1))
In [104]:
# Тогда:
pd.DataFrame( np.transpose( table ), columns=['Значение статистики', 'Уровень значимости'] )
# Получаем красивую таблицу.
Out[104]:
Значение статистики Уровень значимости
0 0.555 15.0
1 0.632 10.0
2 0.759 5.0
3 0.885 2.5
4 1.053 1.0
In [105]:
# Но можно иначе.
table = np.dstack( (res[1], res[2]) )
table, table.shape # Но появится лишняя размерность.
Out[105]:
(array([[[ 0.555, 15.   ],
         [ 0.632, 10.   ],
         [ 0.759,  5.   ],
         [ 0.885,  2.5  ],
         [ 1.053,  1.   ]]]), (1, 5, 2))
In [106]:
# Её можно убрать взяв первый элемент:
table[0], table[0].shape
Out[106]:
(array([[ 0.555, 15.   ],
        [ 0.632, 10.   ],
        [ 0.759,  5.   ],
        [ 0.885,  2.5  ],
        [ 1.053,  1.   ]]), (5, 2))
In [107]:
# тогда опять можем получить таблицу.
pd.DataFrame( table[0], columns=['Значение статистики', 'Уровень значимости'] )
Out[107]:
Значение статистики Уровень значимости
0 0.555 15.0
1 0.632 10.0
2 0.759 5.0
3 0.885 2.5
4 1.053 1.0
In [108]:
# Скорее всего ошибка в коде. Значимости нужны в обратном порядке.
table = np.dstack( (res[1], res[2][::-1] ) )[0]
pd.DataFrame( table, columns=['Значение статистики', 'Уровень значимости'] )
Out[108]:
Значение статистики Уровень значимости
0 0.555 1.0
1 0.632 2.5
2 0.759 5.0
3 0.885 10.0
4 1.053 15.0
In [109]:
# В res[0] у нас значение статистики.
res[0] # Для проверки гипотезы нужно сравнить это число с критическими значениями.
Out[109]:
0.18360768648014414

Бернули

In [110]:
bern = buildBernSum( 0.7, 4000, 3000 )
res = models.anderson( bern ) # Второй параметр указывает на распределение. Поумолчанию проверка на нормальность.
res # Дествительно стремимся к нормальному распределению.
Out[110]:
AndersonResult(statistic=1.0435210798109438, critical_values=array([0.575, 0.655, 0.786, 0.917, 1.091]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
In [111]:
params = models.norm.fit( bern )
params
Out[111]:
(0.6999524166666667, 0.008396367418496462)
In [112]:
models.norm.interval( 0.95, params[0], params[1] )
Out[112]:
(0.6834958389254482, 0.7164089944078853)
In [113]:
models.bernoulli.interval( 0.7, 0.7 ) # Системе плохо.
Out[113]:
(0.0, 1.0)
In [114]:
models.norm.interval( 0.95 )
Out[114]:
(-1.959963984540054, 1.959963984540054)
In [115]:
pp = np.mean( bern )
pp
Out[115]:
0.6999524166666667
In [116]:
dd = np.sqrt(pp*(1-pp)/3000)
dd
Out[116]:
0.00836697936472416
In [117]:
lim = models.norm.interval( 0.95, 0, 1 )
lim
Out[117]:
(-1.959963984540054, 1.959963984540054)
In [118]:
pp + dd*lim[0], pp + dd*lim[1]
Out[118]:
(0.6835534384524176, 0.7163513948809159)
In [119]:
models.laplace.rvs()
Out[119]:
-0.4713776441012143

statmodels

In [120]:
import statsmodels.api as sm
import statsmodels.sandbox.regression.predstd as predstd
In [121]:
nsample = 100
x = np.linspace(-5, 5, 100)
X = np.column_stack((x, x, x**2)) # Отмечу, что переменная x повторена. т.е. есть избыток в данных.
beta = np.array([1, 0.5, 0.5, 0.5])
e = np.random.normal(size=nsample)
In [122]:
#X
sm.
  File "<ipython-input-122-b9fe52f512b9>", line 2
    sm.
       ^
SyntaxError: invalid syntax
In [123]:
X = sm.add_constant(X)
y = np.dot(X, beta) + e
In [124]:
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.954
Model:                            OLS   Adj. R-squared:                  0.954
Method:                 Least Squares   F-statistic:                     1016.
Date:                Sun, 27 Jan 2019   Prob (F-statistic):           8.71e-66
Time:                        10:01:58   Log-Likelihood:                -145.73
No. Observations:                 100   AIC:                             297.5
Df Residuals:                      97   BIC:                             305.3
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.0990      0.158      6.943      0.000       0.785       1.413
x1             0.4876      0.018     26.950      0.000       0.452       0.524
x2             0.4876      0.018     26.950      0.000       0.452       0.524
x3             0.5015      0.014     36.137      0.000       0.474       0.529
==============================================================================
Omnibus:                        3.932   Durbin-Watson:                   1.639
Prob(Omnibus):                  0.140   Jarque-Bera (JB):                2.381
Skew:                          -0.152   Prob(JB):                        0.304
Kurtosis:                       2.308   Cond. No.                          inf
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is      0. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/statsmodels/regression/linear_model.py:1633: RuntimeWarning: divide by zero encountered in double_scalars
  return np.sqrt(eigvals[0]/eigvals[-1])
In [125]:
print('Parameters: ', results.params)
print('R2: ', results.rsquared)
Parameters:  [1.09900322 0.48760878 0.48760878 0.50149395]
R2:  0.9544432141870732
In [126]:
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.954
Model:                            OLS   Adj. R-squared:                  0.954
Method:                 Least Squares   F-statistic:                     1016.
Date:                Sun, 27 Jan 2019   Prob (F-statistic):           8.71e-66
Time:                        10:02:00   Log-Likelihood:                -145.73
No. Observations:                 100   AIC:                             297.5
Df Residuals:                      97   BIC:                             305.3
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.0990      0.158      6.943      0.000       0.785       1.413
x1             0.4876      0.018     26.950      0.000       0.452       0.524
x2             0.4876      0.018     26.950      0.000       0.452       0.524
x3             0.5015      0.014     36.137      0.000       0.474       0.529
==============================================================================
Omnibus:                        3.932   Durbin-Watson:                   1.639
Prob(Omnibus):                  0.140   Jarque-Bera (JB):                2.381
Skew:                          -0.152   Prob(JB):                        0.304
Kurtosis:                       2.308   Cond. No.                          inf
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is      0. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
In [127]:
yy = np.dot( X, results.params)
In [128]:
plt.plot( x, yy)
plt.plot( x, y, ".")
Out[128]:
[<matplotlib.lines.Line2D at 0x7fe74f311160>]
In [129]:
results = sm.OLS(y, sm.add_constant(x)).fit() # Все внимание на F-statistic, Porb(...), R-squared.
print(results.summary()) #  Они показывают, что мы плохо приблизили. "Слуйчайная прямая" дала бы схожий результат.
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.341
Model:                            OLS   Adj. R-squared:                  0.334
Method:                 Least Squares   F-statistic:                     50.74
Date:                Sun, 27 Jan 2019   Prob (F-statistic):           1.79e-10
Time:                        10:02:03   Log-Likelihood:                -279.31
No. Observations:                 100   AIC:                             562.6
Df Residuals:                      98   BIC:                             567.8
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.3625      0.399     13.433      0.000       4.570       6.155
x1             0.9752      0.137      7.123      0.000       0.704       1.247
==============================================================================
Omnibus:                        8.724   Durbin-Watson:                   0.121
Prob(Omnibus):                  0.013   Jarque-Bera (JB):                9.067
Skew:                           0.702   Prob(JB):                       0.0107
Kurtosis:                       2.546   Cond. No.                         2.92
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [130]:
yy = np.dot( sm.add_constant(x), results.params)
plt.plot( x, yy)
plt.plot( x, y, ".")
Out[130]:
[<matplotlib.lines.Line2D at 0x7fe74f2d65c0>]

Доверительный интервал

In [131]:
model = sm.OLS(y, X) # Возвращаеся к прошлой модели.
results = model.fit()
yy = np.dot( X, results.params)

prstd, iv_l, iv_u = predstd.wls_prediction_std( results )
In [132]:
plt.plot(x, yy, 'g-', label="data")
plt.plot(x, y, 'b.', label="True")
plt.plot(x, results.fittedvalues, 'r--.', label="OLS")
plt.plot(x, iv_u, 'r--')
plt.plot(x, iv_l, 'r--')
plt.legend(loc='best');
In [133]:
R = [[0, 1, 0, 0], [0, 0, 1, 0]]
print(np.array(R))
print(results.f_test(R)) # т.е. нестоит отбрасывать гипотезу о том, что 2 и 3 параметр идентичны.
[[0 1 0 0]
 [0 0 1 0]]
<F test: F=array([[726.30483789]]), p=7.723400427367783e-47, df_denom=97, df_num=1>
/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/statsmodels/base/model.py:1532: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 2, but rank is 1
  'rank is %d' % (J, J_), ValueWarning)
In [134]:
R = [[0, 0, 0, 1], [0, 1, 0, 0]]
print(np.array(R))
print(results.f_test(R)) # первая переменная и предпоследняя явно разные.
[[0 0 0 1]
 [0 1 0 0]]
<F test: F=array([[1016.10539598]]), p=8.706110259817051e-66, df_denom=97, df_num=2>
In [ ]: