setwd("C:/Users/yetian/Desktop")
X<-read.table("d9.1.txt",header=T)
X<-as.matrix(scale(X))
X
## x1 x2 x3 x4 x5
## 冀东水泥 1.07204470 1.08614241 -0.496387668 0.7825346 0.15724597
## 大同水泥 0.02961712 -0.01168649 2.034761517 -1.2541471 -1.09517509
## 四川双马 -0.74970573 -0.75938814 -0.586250952 -0.6268359 -0.18669737
## 牡丹江 -1.38415446 -1.33529838 -0.002139602 -0.1542725 -0.79830492
## 西水股份 -1.08608012 -1.13569313 0.117678111 -0.1592643 1.86150069
## 狮头股份 0.12953031 0.11920220 2.259419729 -2.1501872 -0.59592467
## 太行股份 0.51419608 0.49387108 0.027814826 -0.3514512 -0.14924500
## 海螺水泥 1.48668442 1.44772242 -1.095476232 0.6610658 2.05941414
## 尖峰集团 -1.73551583 -1.77541161 -0.586250952 0.1394157 -0.06471947
## 四川金顶 0.22944350 0.21082429 -0.601228167 1.2559299 -1.54219836
## 祁连山 0.99044893 1.01251752 0.252473038 -0.3797384 0.04351445
## 华新水泥 -0.38002694 -0.32581935 -0.541319310 1.5795693 0.39398619
## 福建水泥 1.19027530 1.19903390 -0.930726877 0.9356185 0.44655695
## 天鹅股份 -0.30675727 -0.22601672 0.147632539 -0.2782372 -0.52995352
## x6
## 冀东水泥 0.6043139
## 大同水泥 -0.3009265
## 四川双马 -0.3408551
## 牡丹江 -0.5825889
## 西水股份 1.1587015
## 狮头股份 0.1388646
## 太行股份 0.1294951
## 海螺水泥 2.0636536
## 尖峰集团 -0.9870642
## 四川金顶 -0.4244601
## 祁连山 0.9374366
## 华新水泥 0.2216047
## 福建水泥 -1.9952252
## 天鹅股份 -0.6229499
## attr(,"scaled:center")
## x1 x2 x3 x4 x5 x6
## 27.362143 28.111429 1.001429 50.364286 10.913571 -25.573571
## attr(,"scaled:scale")
## x1 x2 x3 x4 x5 x6
## 6.0052133 6.1120636 0.6676809 12.0195512 29.1036307 69.3738381
Fr1<-factanal(X,3,scores="regression")
loadr<-Fr1$loadings
ldr<-cbind(loadr[,1],loadr[,2],loadr[,3])
R<-solve(t(X)%*%X)
Fr2<-X%*%R%*%ldr
Fr2
## [,1] [,2] [,3]
## 冀东水泥 0.08131519 0.038352531 -0.001485989
## 大同水泥 0.01929590 -0.151678115 -0.042355297
## 四川双马 -0.05860636 0.047643195 -0.027417705
## 牡丹江 -0.09709158 0.008331800 -0.063453607
## 西水股份 -0.10864694 -0.028092126 0.161415013
## 狮头股份 0.02302039 -0.175697402 0.005030787
## 太行股份 0.04128859 -0.001326998 -0.012729006
## 海螺水泥 0.08756259 0.066222223 0.142729640
## 尖峰集团 -0.13838833 0.047801980 -0.015565896
## 四川金顶 0.03382008 0.064542688 -0.144247035
## 祁连山 0.07861897 -0.021350910 0.007874588
## 华新水泥 -0.03370105 0.041013082 0.020009985
## 福建水泥 0.08571936 0.070760223 0.010431470
## 天鹅股份 -0.01420680 -0.006522171 -0.040236949
Fr1$scores
## Factor1 Factor2 Factor3
## 冀东水泥 1.0570974 0.49858291 -0.01931786
## 大同水泥 0.2508467 -1.97181549 -0.55061886
## 四川双马 -0.7618827 0.61936154 -0.35643016
## 牡丹江 -1.2621905 0.10831340 -0.82489689
## 西水股份 -1.4124102 -0.36519764 2.09839517
## 狮头股份 0.2992651 -2.28406622 0.06540023
## 太行股份 0.5367516 -0.01725097 -0.16547708
## 海螺水泥 1.1383136 0.86088890 1.85548532
## 尖峰集团 -1.7990483 0.62142573 -0.20235665
## 四川金顶 0.4396610 0.83905494 -1.87521145
## 祁连山 1.0220466 -0.27756182 0.10236965
## 华新水泥 -0.4381137 0.53317006 0.26012980
## 福建水泥 1.1143517 0.91988289 0.13560911
## 天鹅股份 -0.1846884 -0.08478822 -0.52308034
用教材公式带入进行矩阵计算,结果与函数计算结果不同。
Fr2c<-X%*%solve(cor(X))%*%ldr
Fr2c
## [,1] [,2] [,3]
## 冀东水泥 1.0570974 0.49858291 -0.01931786
## 大同水泥 0.2508467 -1.97181549 -0.55061886
## 四川双马 -0.7618827 0.61936154 -0.35643016
## 牡丹江 -1.2621905 0.10831340 -0.82489689
## 西水股份 -1.4124102 -0.36519764 2.09839517
## 狮头股份 0.2992651 -2.28406622 0.06540023
## 太行股份 0.5367516 -0.01725097 -0.16547708
## 海螺水泥 1.1383136 0.86088890 1.85548532
## 尖峰集团 -1.7990483 0.62142573 -0.20235665
## 四川金顶 0.4396610 0.83905494 -1.87521145
## 祁连山 1.0220466 -0.27756182 0.10236965
## 华新水泥 -0.4381137 0.53317006 0.26012980
## 福建水泥 1.1143517 0.91988289 0.13560911
## 天鹅股份 -0.1846884 -0.08478822 -0.52308034
将X的相关系数矩阵代替公式中的R,计算的结果与函数计算结果相同。
Fb1<-factanal(X,3,scores="Bartlett")
loadb<-Fb1$loadings
ldb<-cbind(loadb[,1],loadb[,2],loadb[,3])
omega<-diag(Fr1$uniquenesses)
Fb2<-t(solve(t(ldb)%*%solve(omega)%*%ldb)%*%t(ldb)%*%solve(omega)%*%t(X))
Fb2
## [,1] [,2] [,3]
## 冀东水泥 1.0599581 0.50128979 -0.02155828
## 大同水泥 0.2521911 -1.98142175 -0.54969418
## 四川双马 -0.7637094 0.62370640 -0.35955744
## 牡丹江 -1.2649294 0.11099850 -0.82940358
## 西水股份 -1.4180692 -0.37210624 2.11411057
## 狮头股份 0.3002321 -2.29687590 0.07102437
## 太行股份 0.5383645 -0.01699915 -0.16697768
## 海螺水泥 1.1397265 0.86110853 1.86446655
## 尖峰集团 -1.8038441 0.62550728 -0.20358403
## 四川金顶 0.4424241 0.84813722 -1.88976307
## 祁连山 1.0247754 -0.27943325 0.10280481
## 华新水泥 -0.4395863 0.53553500 0.26092456
## 福建水泥 1.1171955 0.92454310 0.13331611
## 天鹅股份 -0.1847288 -0.08398954 -0.52610873
Fb1$scores
## Factor1 Factor2 Factor3
## 冀东水泥 1.0599581 0.50128979 -0.02155828
## 大同水泥 0.2521911 -1.98142175 -0.54969418
## 四川双马 -0.7637094 0.62370640 -0.35955744
## 牡丹江 -1.2649294 0.11099850 -0.82940358
## 西水股份 -1.4180692 -0.37210624 2.11411057
## 狮头股份 0.3002321 -2.29687590 0.07102437
## 太行股份 0.5383645 -0.01699915 -0.16697768
## 海螺水泥 1.1397265 0.86110853 1.86446655
## 尖峰集团 -1.8038441 0.62550728 -0.20358403
## 四川金顶 0.4424241 0.84813722 -1.88976307
## 祁连山 1.0247754 -0.27943325 0.10280481
## 华新水泥 -0.4395863 0.53553500 0.26092456
## 福建水泥 1.1171955 0.92454310 0.13331611
## 天鹅股份 -0.1847288 -0.08398954 -0.52610873
计算的结果与函数计算结果相同。
library(mvstats)
factanal.rank(Fr1,plot=T)
## $Fs
## Factor1 Factor2 Factor3
## 冀东水泥 1.0570974 0.49858291 -0.01931786
## 大同水泥 0.2508467 -1.97181549 -0.55061886
## 四川双马 -0.7618827 0.61936154 -0.35643016
## 牡丹江 -1.2621905 0.10831340 -0.82489689
## 西水股份 -1.4124102 -0.36519764 2.09839517
## 狮头股份 0.2992651 -2.28406622 0.06540023
## 太行股份 0.5367516 -0.01725097 -0.16547708
## 海螺水泥 1.1383136 0.86088890 1.85548532
## 尖峰集团 -1.7990483 0.62142573 -0.20235665
## 四川金顶 0.4396610 0.83905494 -1.87521145
## 祁连山 1.0220466 -0.27756182 0.10236965
## 华新水泥 -0.4381137 0.53317006 0.26012980
## 福建水泥 1.1143517 0.91988289 0.13560911
## 天鹅股份 -0.1846884 -0.08478822 -0.52308034
##
## $Ri
## F rank
## 冀东水泥 0.57762286 3
## 大同水泥 -0.73582797 14
## 四川双马 -0.17323177 9
## 牡丹江 -0.66885707 13
## 西水股份 -0.11845569 8
## 狮头股份 -0.66290875 12
## 太行股份 0.15786371 5
## 海螺水泥 1.23140339 1
## 尖峰集团 -0.53300988 11
## 四川金顶 -0.03369328 7
## 祁连山 0.32577583 4
## 华新水泥 0.08514207 6
## 福建水泥 0.78759236 2
## 天鹅股份 -0.23941579 10
Fr1$loadings
##
## Loadings:
## Factor1 Factor2 Factor3
## x1 0.983 0.155
## x2 0.985 0.142
## x3 -0.990 -0.124
## x4 0.127 0.844
## x5 0.293 0.953
## x6 0.210 0.631
##
## Factor1 Factor2 Factor3
## SS loadings 1.998 1.800 1.367
## Proportion Var 0.333 0.300 0.228
## Cumulative Var 0.333 0.633 0.861
自编函数是将得分由大到小排列的。而从因子载荷不难看出,方差贡献率高的变量与原变量正相关的多些,而且第一个贡献率最高的是正相关的。所以对综合得分由大到小排列较好。