LMM_variable_structure.sagews

Построение ЛММ с переменной структурой

def MNK (X, Y) :
    '''
    Построение линейной многофакторной модели методом наименьших квадратов.
    
    Аргументы:
    X, Y - исходные данные
    
    Результат:
    beta_hat - вектор коэффициентов ЛММ
    u_hat - вектор остатков ЛММ
    '''
    beta_hat = (X.transpose()*X)^(-1)*(X.transpose()*Y)
    Y_hat = X * beta_hat
    u_hat = Y - Y_hat
    return beta_hat, u_hat

# чтение из файла (внимание - текстовые переменные! кириллица!!)
X = []
Y = []
f = open('apartment.tsv', 'r')
head = f.readline()
for line in f :
    line = line.decode('utf-8').replace(',', '.').split()
    row = []
    for cell in line :
        if cell.isalpha() : # если ячейка - текст
            row.append(cell)
        else :
            row.append(float(cell))
    Y.append(row[1])
    X.append([1]+row[2:])
Y = vector(Y)
T = len(Y)
# вводим dummy-переменные
# район города - Х2
z21 = vector([int(X[i][2]==u'М') for i in xrange(T)])
z22 = vector([int(X[i][2]==u'П') for i in xrange(T)])
z23 = vector([int(X[i][2]==u'С') for i in xrange(T)])
# новострой/не новострой - X8
z8 = vector([int(X[i][8]==u'Н') for i in xrange(T)])
# удаляем старые текстовые переменные
for row in X :
    del row[2]
    del row[-1]
# делаем X матрицей
X = matrix(X)
# добавляем dummy-переменные
X = X.augment(matrix([z21, z22, z23, z8]).transpose())
show(X[:10])

$\displaystyle \left(\begin{array}{rrrrrrrrrrr} 1.0 & 1.0 & 39.0 & 20.0 & 8.2 & 0.0 & 1.0 & 0.0 & 0.0 & 1.0 & 1.0 \\ 1.0 & 3.0 & 68.4 & 40.5 & 10.7 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0 & 1.0 \\ 1.0 & 2.0 & 54.7 & 28.0 & 10.7 & 0.0 & 1.0 & 0.0 & 0.0 & 1.0 & 1.0 \\ 1.0 & 4.0 & 90.0 & 64.0 & 15.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 0.0 \\ 1.0 & 1.0 & 34.8 & 16.0 & 10.7 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 \\ 1.0 & 2.0 & 48.0 & 29.0 & 8.0 & 1.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 1.0 & 3.0 & 68.1 & 44.4 & 7.2 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 0.0 \\ 1.0 & 4.0 & 132.0 & 89.6 & 11.0 & 1.0 & 1.0 & 0.0 & 1.0 & 0.0 & 1.0 \\ 1.0 & 1.0 & 39.0 & 20.0 & 8.5 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0 & 1.0 \\ 1.0 & 2.0 & 55.5 & 35.0 & 8.0 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \end{array}\right)$

# строим модель
beta, u = MNK(X, Y)
show(r'$\beta='+latex(beta.n(digits=5))+'$')
Sigma_beta = variance(list(u)) * (X.transpose()*X)^(-1)
# найдём коэффициент детерминации
RSS = sum([r^2 for r in u])
TSS = sum([(Y[i]-mean(Y))^2 for i in xrange(T)])
R2 = 1-RSS/TSS
show('$R^2='+latex(R2)+'$')
# проверяем значимость коэффициентов
N = len(beta)
import scipy.stats
t_crit = scipy.stats.t.isf(0.1, T-N)
for k in xrange(N) :
    t_stat = beta[k]/sqrt(Sigma_beta[k,k])
    if t_stat > t_crit :
        show(r'$\beta_{'+str(k)+'}$'+' значимо отличается от нуля')
    else :
        show(r'$\beta_{'+str(k)+'}$'+' НЕ значимо отличается от нуля')

$\beta= \left(-3.8717,\,0.36317,\,0.31476,\,0.081335,\,-0.18805,\,3.5400,\,1.7450,\,5.1333,\,7.0231,\,-0.10794,\,5.4253\right)$

$R^2= 0.902992219201$

$\beta_{0}$ НЕ значимо отличается от нуля

$\beta_{1}$ НЕ значимо отличается от нуля

$\beta_{2}$ значимо отличается от нуля

$\beta_{3}$ НЕ значимо отличается от нуля

$\beta_{4}$ НЕ значимо отличается от нуля

$\beta_{5}$ значимо отличается от нуля

$\beta_{6}$ значимо отличается от нуля

$\beta_{7}$ значимо отличается от нуля

$\beta_{8}$ значимо отличается от нуля

$\beta_{9}$ НЕ значимо отличается от нуля

$\beta_{10}$ значимо отличается от нуля

# тест Чоу
# проверим, действительно ли выборки для новостроек и не-новостроек разные
XX = X[:, :5].augment(X.column(-1))
Y1 = vector([Y[i] for i in xrange(len(Y)) if X[i,-1]==1])
Y2 = vector([Y[i] for i in xrange(len(Y)) if X[i,-1]==0])
X1 = matrix([row[:-1] for row in XX if row[-1]==1])
X2 = matrix([row[:-1] for row in XX if row[-1]==0])
beta, eps = MNK(XX, Y)
beta1, eps1 = MNK(X1, Y1)
beta2, eps2 = MNK(X2, Y2)
p = XX.ncols()
n1 = len(Y1)
n2 = len(Y2)
# используем замечание 10.4
F_stat = (eps*eps-eps1*eps1-eps2*eps2)*(n1+n2-2*p)/p/(eps1*eps1+eps2*eps2)
F_crit = scipy.stats.f.isf(0.1, p, n1+n2-2*p)
if F_stat > F_crit :
    print 'Гипотеза об однородности отвергнута'
else :
    print 'Выборки однородны'

Гипотеза об однородности отвергнута