Невозможно получить коэффициенты линейной регрессии в R после успешного поиска вопросов и ответов через Householder

Я пытаюсь вручную рассчитать коэффициенты регрессии, а не использовать данные по умолчанию http://people.sc.fsu.edu/~jburkardt/datasets/regression/x31.txt

Вот мой код, который правильно создает Q&R, удовлетворяющий A=QR. Но я не могу найти коэффициенты как измерения вопросов и ответов, создающих проблемы. Могут ли мне здесь помочь специалисты? Когда у меня есть правильные вопросы и ответы, как поиск коэффициентов может пойти не так?

library(xlsx)
data.df<-read.xlsx("regression.xlsx",2,header = F)
#Remove unneccesary Index position
data.df<-data.df[2:5]

#Decomposing
#coefficients [b]=inv(t(X)(Matrix Multiplication)(X))(Matrix Multiplication)t(X)(Matrix Multiplication)y
Y<-as.matrix(data.df[4])
#But note that if you need to express Y=Xb, the coefficient of b_0 must be X_0 
#which is 1
X_0<-as.data.frame(c(1,1,1,1,1,1,1,1,1,1))
X<-(cbind(X_0,data.df[1:3]))
names(X)<-c("X1","X2","X3","X4")
X<-as.matrix(X)
#Create copy for final evaluvations
A<-X  


x1<-as.matrix(X[,1])
deter_x<-sqrt(sum(x1^2))
n=dim(x1)[1]
deter_e1<-as.matrix(c(deter_x,rep(0,n-1)))
v1=x1+deter_e1
#c_i=2/(transpose(v_i)%*%(v_i))
c1<-as.numeric(2/(t(v1)%*%v1))
#H_i = I - c_i*v_i%*%transpose(v_i)
I<-diag(n)
H1<-I-c1*(v1%*%t(v1))
R1<-H1%*%X
#Check R1 and see if it is Upper Triangle Matrix
R1
#We will take rest of the interesting portion of matrix R1.
n=dim(R1)[1]
X<-as.matrix(as.data.frame(cbind(R1[2:n,2],R1[2:n,3],R1[2:n,4])))
x1<-as.matrix(X[,1])
deter_x<-sqrt(sum(x1^2))
n=dim(x1)[1]
deter_e1<-as.matrix(c(deter_x,rep(0,n-1)))
v1=x1-deter_e1
#c_i=2/(transpose(v_i)%*%(v_i))
c1<-as.numeric(2/(t(v1)%*%v1))
#H_i = I - c_i*v_i%*%transpose(v_i)
I<-diag(n)
H2<-I-c1*(v1%*%t(v1))
R2<-H2%*%X

#Check R2 and see if it is Upper Triangle Matrix, if no go for R3
n=dim(R2)[1]
X<-as.matrix(as.data.frame(cbind(R2[2:n,2],R2[2:n,3])))
x1<-as.matrix(X[,1])
deter_x<-sqrt(sum(x1^2))
n=dim(x1)[1]
deter_e1<-as.matrix(c(deter_x,rep(0,n-1)))
v1=x1+deter_e1
#c_i=2/(transpose(v_i)%*%(v_i))
c1<-as.numeric(2/(t(v1)%*%v1))
#H_i = I - c_i*v_i%*%transpose(v_i)
I<-diag(n)
H3<-I-c1*(v1%*%t(v1))
R3<-H3%*%X
R3

#Check R3 and see if it is Upper Triangle Matrix, if no go for R4
n=dim(R3)[1]
X<-as.matrix(as.data.frame(cbind(R3[2:n,2])))
x1<-as.matrix(X[,1])
deter_x<-sqrt(sum(x1^2))
n=dim(x1)[1]
deter_e1<-as.matrix(c(deter_x,rep(0,n-1)))
v1=x1-deter_e1
#c_i=2/(transpose(v_i)%*%(v_i))
c1<-as.numeric(2/(t(v1)%*%v1))
#H_i = I - c_i*v_i%*%transpose(v_i)
I<-diag(n)
H4<-I-c1*(v1%*%t(v1))
R4<-H4%*%X
R4
#As we can see R4 has all values except first element as zero
#Let us replace Matrices iteratively in R1 from R2 to R4 and round it of
R1[2:10,2:4]<-R2
R1[3:10,3:4]<-R3
R1[4:10,4]<-R4
R<-round(R1,5)
R
#Find Complete H1
#Q=H1%*%H2%*%H3%*%H4
H1_COM<-H1
# 
H_temp<-diag(10)
n=dim(H_temp)[1]
dim(H_temp[2:n,2:n])
dim(H2)
H_temp[2:n,2:n]<-H2
H2_COM<-H_temp
H2_COM
H_temp<-diag(10)
n=dim(H_temp)[1]
dim(H_temp[3:n,3:n])
dim(H3)
H_temp[3:n,3:n]<-H3
H3_COM<-H_temp
H3_COM

H_temp<-diag(10)
n=dim(H_temp)[1]
dim(H_temp[4:n,4:n])
dim(H4)
H_temp[4:n,4:n]<-H4
H4_COM<-H_temp
Q=H1_COM%*%H2_COM%*%H3_COM%*%H4_COM
# The following code properly reconstructs A Matrix proving proper Q & R
A=round(Q%*%R)
# When you try to find coefficients using Q&R you will end up in error.
solve(R)%*%t(Q)%*%Y
#Error in solve.default(R) : 'a' (10 x 4) must be square
#So trying to get matrix R without all 0 rows R[1:4,1:4]
solve(R[1:4,1:4])%*%t(Q)%*%Y
#Error in solve(R[1:4, 1:4]) %*% t(Q) : non-conformable arguments
dim(solve(R[1:4,1:4]))
dim(solve(R[1:4,1:4]))
#4 4
dim(t(Q))
#[1] 10 10
dim(Y)
#10  1

person Hari Prasad    schedule 16.03.2017    source источник


Ответы (1)


Я хотел бы указать вам на эту тему, на которую я ответил (довольно исчерпывающе): Множественный регрессионный анализ в R с использованием декомпозиции QR . Сообщество решит, будет ли ваш вопрос рассматриваться как дубликат. Обратите внимание, что среди вариантов последний использует факторизацию QR, написанную в простом коде R.

Учитывая, что ваш игрушечный код для QR-факторизации верен (как вы прокомментировали в своем коде), основная проблема связана с вашими последними несколькими строками.

Решение простое:

solve(R) %*% (t(Q) %*% Y)[1:4,]

1:4 выбирает первые 4 элемента в одностолбцовой матрице t(Q) %*% Y.

Если вы сверитесь с моим связанным ответом, вы увидите, что вместо solve я использую backsolve, поскольку это треугольная система уравнений.

Вы также заметите, что я использую crossprod, когда могу, вместо t и %*%. Мой недавний ответ Когда "кросспрод" предпочтительнее "%*%", а когда нет? содержит всестороннее обсуждение на этих двух модах.

person Zheyuan Li    schedule 16.03.2017
comment
Вот Это Да! Это было быстро. Спасибо. Просто чтобы изменить, ваш ответ правильный, если я изменяю свою R на квадратную матрицу 4X4 для данного примера: решить (R [1: 4, 1: 4]) % *% (t (Q) % *% Y) [1 :4,]. Единственное, что я хочу понять, это насколько правильно с математической точки зрения выбрать первые 4 элемента в одностолбцовой матрице из t(Q) %*% Y, дав правильные результаты? Я определенно собираюсь прочитать ваши другие ответы, но если вы можете помочь мне дать ответ, когда вы свободны, это большая помощь! Спасибо еще раз! Ваше здоровье! - person Hari Prasad; 16.03.2017
comment
Знаешь что! Ты великолепен! Спасибо за помощь. Я отметил это как ответ и проголосовал за него. Я думаю, что ключевой вывод из вашего ответа: он показывает, как используется ортогональное преобразование к Y. . Спасибо еще раз. :-) - person Hari Prasad; 16.03.2017
comment
Привет Чжэюань Ли, я извиняюсь за недостаток знаний. Тем не менее я не могу понять, как правильно брать только 4 верхние матрицы единственного столбца? Я прочитал всю литературу, которую вы указали. Но как оставление остальных значений даст мне правильные результаты, я не могу понять :-( - person Hari Prasad; 20.03.2017