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

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

# строим модель
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)+'}$'+' НЕ значимо отличается от нуля')

# тест Чоу
# проверим, действительно ли выборки для новостроек и не-новостроек разные
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 'Выборки однородны'