多元统计分析上机

协方差

'''

data.iris=iris[,-5]

data.iris

cov(data.iris)

cor(data.iris)

'''

offer 数据:

offer.txt

data=read.table("d:/offer.txt",header=TRUE,sep="")

offer=data

M=cor(offer)

M

f=function(x,k1,k2){

 which(x<=k2 & x>=k1)

}

k1=0.8

k2=1

result=apply(M,2,f,k1,k2)

name=names(result) # Extract name from a list

result$SC

names(result$SC)

c("GSP",names(result$GSP))

offer 主成分分析

data=read.table("d:/offer.txt",header=TRUE,sep="")

offer=data

offer=as.matrix(offer)

meanValue=apply(offer,1,mean)

offer=offer-meanValue

M=cov(offer)

N=eigen(M)

value=N$values

vector=N$vectors

value[1:2]

vector[,1:2]

sum(value[1:4]/sum(value))

plot(value,type="l")

Y=offer%*%(vector[,1:4])

Y

工业数据:

industry.txt

data=read.table("d:/industry.txt",header=TRUE,sep="")

industry=data

M=as.data.frame(industry)

industry_object=princomp(M,cor=T,scores=T)

summary(industry_object,loadings=T)

loadings(industry_object)

scores=industry_object$scores

tatal=apply(scores,1,sum)

N=cbind(scores,tatal)

T=round(N,digits=3)

T

对称矩阵谱分解

data=read.table("d:/industry.txt",header=TRUE,sep="")

industry=data

M=as.data.frame(industry)

N=eigen(cor(M))

values=N$values

vectors=N$vectors

accurate=sum(values[1:3])/sum(values)

M=values[1]*vectors[,1]%*%t(vectors[,1])+

  values[2]*vectors[,2]%*%t(vectors[,2])+

  values[3]*vectors[,3]%*%t(vectors[,3])

M

因子分析主成分法

data=read.table("C:/offer.txt",header=TRUE,sep="")

M=as.data.frame(data)

C=cor(M)

N=eigen(cor(M))

values=N$values

vectors=N$vectors

accurate=sum(values[1:4])/sum(values)

A=sqrt(values[1:4])*vectors[,1:4]

D=diag(diag(C-A%*%t(A)))

A%*%t(A)-D

因子分析主因子法

data=read.table("d:/offer.txt",header=TRUE,sep="")

R=cor(data)

invR=solve(R)

D=diag(1/diag(invR))

k=0

repeat{

  R1=R-D

  D1=D

  eig=eigen(R1)

  val=eig$values

  vec=eig$vectors

  A=(vec[,1:3])%*%(diag(sqrt(val[1:3])))

  D=diag(1-diag(A%*%t(A)))

  k=k+1

  if (sqrt(sum(D-D1)^2)<1e-5|k>20){break}

}

A

D

因子旋转

data=read.table("d:/offer.txt",header=TRUE,sep="")

R=cor(data)

invR=solve(R)

D=diag(1/diag(invR))

k=0

repeat{

  R1=R-D

  D1=D

  eig=eigen(R1)

  val=eig$values

  vec=eig$vectors

  A=(vec[,1:3])%*%(diag(sqrt(val[1:3])))

  D=1-diag(A%*%t(A))

  D=diag(D)

  k=k+1

  if(sqrt(sum(D-D1)^2)<1e-5|k>20){break}

}

print(A)

A.v=varimax(A,eps=1e-5)

print(A.v)

马氏距离判别法

#马氏距离

data=as.matrix(iris[,-5])

sample=data[5,]

s_1=data[1:50,]

s_2=data[51:100,]

miu_1=apply(s_1,2,mean)

miu_2=apply(s_2,2,mean)

cor_1=cov(s_1)

cor_2=cov(s_2)

m_1=(sample-miu_1)%*%solve(cor_1)%*%(sample-miu_1)

m_2=(sample-miu_2)%*%solve(cor_2)%*%(sample-miu_2)

m_1

m_2

#欧氏距离

(sample-miu_1)%*%(sample-miu_1)

(sample-miu_2)%*%(sample-miu_2)

方法二:

n1=nrow(data[1:50,])

n2=nrow(data[51:100,])

n3=nrow(data[101:150,])

s1=(n1-1)*cov(data[1:50,])

s2=(n2-1)*cov(data[51:100,])

s3=(n3-1)*cov(data[101:150,])

sigma=(s1+s2+s3)/(n1+n2+n3-3)

miu1=apply(data[1:50,],2,mean)

miu2=apply(data[51:100,],2,mean)

miu3=apply(data[101:150,],2,mean)

d1=mahalanobis(data,miu1,sigma)

d2=mahalanobis(data,miu2,sigma)

d3=mahalanobis(data,miu3,sigma)

M=matrix(0,3,150)

M[1,]=d1

M[2,]=d2

M[3,]=d3

class=apply(M,2,which.min)

x1=class[1:50]

x2=class[51:100]

x3=class[101:150]

length(x1[x1!=1])+length(x2[x2!=2])+length(x3[x3!=3])

Fisher 判别法

data=as.matrix(iris[,-5])

data=data[1:100,]

miu_1=apply(data[1:50,],2,mean)

miu_2=apply(data[51:100,],2,mean)

sigma=cov(data) #假设总体有相同的协方差阵

# 如果未知总体的协方差阵, 则计算总体协方差的无偏估计

s_1=49*cov(data[1:50,])

s_2=49*cov(data[51:100,])

sigma=(s_1+s_2)/(50+50-2)

# if class>=0, then the sample belong to class 1, else class 2.

f=function(x){class=(miu_1-miu_2)%*%solve(sigma)%*%(2*x-miu_1-miu_2)}

apply(data,1,f)

方法二:

library(MASS)

attach(iris)

ld=lda(Species~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width)

pre=predict(ld)

newclass=pre$class

cbind(Species,newclass,pre$x)

table(newclass,Species)

Bayes 判别

library(MASS)

attach(iris)

ld=lda(Species~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width)

prior=c(1/3,1/3,1/3)

pre=predict(ld,newdata)

newclass=pre$class

cbind(Species,newclass,pre$x)

table(newclass,Species)

例 1

a=c(11.2,14.9,14.3,13.5,16.2,14.3,20.0,21.8,19.0,16.0,11.9,8.7,14.3,10.1,9.1,13.8,15.3,11.0,18.0,10.4,8.2,11.4,11.6,8.4,8.2,10.9,15.6)

b=c(57.25,67.19,67.74,55.63,75.51,57.63,83.94,68.03,78.31,57.11,49.97,30.72,37.65,34.63,56.33,65.23,55.62,55.55,62.85,30.01,29.28,62.88,28.57,30.23,15.96,24.75,21.44)

c=c(13.47,7.89,19.41,20.59,11.06,22.51,15.09,39.42,83.13,12.57,30.70,15.41,12.95,7.68,10.30,4.69,6.06,8.02,6.40,4.61,6.11,5.31,5.08,6.03,8.04,8.34,28.62)

d=c(73.41,73.09,72.33,73.33,72.08,77.35,89.50,71.90,89.45,60.51,69.20,60.29,66.42,62.96,66.01,64.24,54.74,67.47,58.83,60.26,50.71,61.49,68.47,55.55,40.26,46.01,46.01)

a1=c(16.5,20.6,8.6)

b1=c(80.05,81.24,42.06)

c1=c(8.81,5.37,8.88)

d1=c(73.04,60.43,56.37)

data=data.frame(a,b,c,d)

sample=data.frame(a1,b1,c1,d1)

sample=as.matrix(sample)

s_1=data[1:10,]

s_2=data[11:27,]

miu_1=apply(s_1,2,mean)

miu_2=apply(s_2,2,mean)

cor_1=cov(s_1)

cor_2=cov(s_2)

sample=sample[3,]

(sample-miu_1)%*%solve(cor_1)%*%(sample-miu_1)

(sample-miu_2)%*%solve(cor_2)%*%(sample-miu_2)

例 2

library(MASS)

a=c(11.2,14.9,14.3,13.5,16.2,14.3,20.0,21.8,19.0,16.0,11.9,8.7,14.3,10.1,9.1,13.8,15.3,11.0,18.0,10.4,8.2,11.4,11.6,8.4,8.2,10.9,15.6,16.5,20.6,8.6)

b=c(57.25,67.19,67.74,55.63,75.51,57.63,83.94,68.03,78.31,57.11,49.97,30.72,37.65,34.63,56.33,65.23,55.62,55.55,62.85,30.01,29.28,62.88,28.57,30.23,15.96,24.75,21.44,80.05,81.24,42.06)

c=c(13.47,7.89,19.41,20.59,11.06,22.51,15.09,39.42,83.13,12.57,30.70,15.41,12.95,7.68,10.30,4.69,6.06,8.02,6.40,4.61,6.11,5.31,5.08,6.03,8.04,8.34,28.62,8.81,5.37,8.88)

d=c(73.41,73.09,72.33,73.33,72.08,77.35,89.50,71.90,89.45,60.51,69.20,60.29,66.42,62.96,66.01,64.24,54.74,67.47,58.83,60.26,50.71,61.49,68.47,55.55,40.26,46.01,46.01,73.04,60.43,56.37)

e=c(rep(c(0),10),rep(c(1),20))

data=data.frame(a,b,c,d,e)

attach(data)

ld=lda(e~a+b+c+d)

pre=predict(ld)

newclass=pre$class

cbind(e,newclass,pre$x)

table(newclass,e)

例 3

library(MASS)

a=c(11.2,14.9,14.3,13.5,16.2,14.3,20.0,21.8,19.0,16.0,11.9,8.7,14.3,10.1,9.1,13.8,15.3,11.0,18.0,10.4,8.2,11.4,11.6,8.4,8.2,10.9,15.6)

b=c(57.25,67.19,67.74,55.63,75.51,57.63,83.94,68.03,78.31,57.11,49.97,30.72,37.65,34.63,56.33,65.23,55.62,55.55,62.85,30.01,29.28,62.88,28.57,30.23,15.96,24.75,21.44)

c=c(13.47,7.89,19.41,20.59,11.06,22.51,15.09,39.42,83.13,12.57,30.70,15.41,12.95,7.68,10.30,4.69,6.06,8.02,6.40,4.61,6.11,5.31,5.08,6.03,8.04,8.34,28.62)

d=c(73.41,73.09,72.33,73.33,72.08,77.35,89.50,71.90,89.45,60.51,69.20,60.29,66.42,62.96,66.01,64.24,54.74,67.47,58.83,60.26,50.71,61.49,68.47,55.55,40.26,46.01,46.01)

e=c(rep(c(0),10),rep(c(1),17))

a1=c(16.5,20.6,8.6)

b1=c(80.05,81.24,42.06)

c1=c(8.81,5.37,8.88)

d1=c(73.04,60.43,56.37)

e1=c(0,0,0)

data=data.frame(a,b,c,d,e)

newdata=data.frame(a1,b1,c1,d1,e1)

attach(data)

ld=lda(e~a+b+c+d)

prior=c(10/27,17/27)

pre=predict(ld,newdata)

newclass=pre$class

cbind(e,newclass,pre$x)

table(newclass,e)

k-means 聚类

1.

trace=function(M){

  miu=apply(M,2,mean)

  r=nrow(M)

  p=ncol(M)

  mur=rep(miu,r)

  Mu=matrix(mur,r,p,byrow=TRUE)

  tr=sum((M-Mu)^2)

  list(mean=miu,trace=tr)

}

y1=c(0,1,2,2,3,4)

y2=c(0,0,0,1,1,1)

df=data.frame(y1,y2)

trace(as.matrix(df[1:5,]))

trace(as.matrix(df[6,]))

2.

data=as.matrix(iris[,-5])

x0=data[3:5,]

dis=NULL

for(i in seq(1,nrow(x0))){

 x=rep(x0[i,],150)

 M=matrix(x,150,4,byrow=T)

 N=apply((data-M)^2,1,sum)

 dis=cbind(dis,N)

}

apply(dis,1,which.max)

3.

nnum=150

d.ir=dist(iris[1:num,-5])

d.m=as.matrix(d.ir)

r=0.55

radi=matrix(rep(r,num*num),num,num)

nei=d.m-radi+diag(rep(r,num))

f=function(x){

 y=x[x<0]

 length(y)

}

n.nei=apply(nei,2,f)

x=nei[,1] # Pick a sample point.

A=x[x<=0]

name.A=names(A)

y=as.numeric(name.A)

y1=as.numeric(names(which(n.nei[y]>=3)))

M=nei[,y1]

g=function(x){which(x<=0)}

lis=apply(M,2,g)

u=NULL

for(i in 1:(length(lis)-1)){

  u=union(u,lis[[i]])

}

#u=sort(u)

z=setdiff(u,y)

all(n.nei[z]>=3)

w=vector(len=0)

x=nei[,136]

A=x[x<=0]

name.A=names(A)

g=function(x){

  which(x<0)

}

y=as.numeric(name.A)

repeat{

  y1=as.numeric(names(which(n.nei[y]>=5)))

  if(length(y1)==0) break else{

    M=nei[,y1]

    lis=apply(M,2,g)

    u=lis[[1]]

    for(i in 1:(length(lis)-1)){

      u=union(u,lis[[i+1]])

    }

    if(length(w)==length(u)) break else{

      w=union(w,u)

      y=w

    }

  }

}

DB-Scan 聚类

num=150

data=iris[1:num,-5]

dis=dist(data)

dis=as.matrix(dis)

cutoff=0.55

cutoffM=matrix(rep(cutoff,num*num),num,num)

dis2=dis-cutoffM+diag(rep(cutoff,num))

f=function(x){

 return(length(x[x<0]))

}

rho=apply(dis2,2,f)

delta=NULL

for(i in as.numeric(names(rho))){


  denseMoreThanI=rho[rho>rho[i]]

  if(length(denseMoreThanI)!=0){

    name=as.numeric(names(denseMoreThanI))

    deltaI=min(dis[i,name])

    delta=c(delta,deltaI)

  }else{

    delta=c(delta,max(dis[i,]))

  }

}

plot(rho,delta)

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 202,529评论 5赞 475
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,015评论 2赞 379
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 149,409评论 0赞 335
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,385评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,387评论 5赞 364
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,466评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,880评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,528评论 0赞 256
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,727评论 1赞 295
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,528评论 2赞 319
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,602评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,302评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,873评论 3赞 306
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,890评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,132评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,777评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,310评论 2赞 342