R.matlab尾部相关性分析

library("R.matlab")
path<-('D:/M/mat/(你的mat文件所在的路径)')
pathname<-file.path(path,'s.mat')
testdata<-readMat(pathname)
print(testdata)

**

library("R.matlab")
path<-('D:/Alice/研/任务2021.9/')
rainfall<-file.path(path,'rainfall.mat')
runoff<-file.path(path,'runoff.mat')
X_1<-readMat(rainfall)
X_2<-readMat(runoff)
unlist(X_1) ->X_1
unlist(X_2) ->X_2

**

library(copula)
library(psych)
library(VineCopula)
plot(X_1,X_2)
abline(lm(X_2~X_1),col='red ',lwd=1)
cor(X_1,X_2,method='spearman')##皮尔森相关性系数
#以下4行代码是Copula函数的非参数估计，在不知道哪一个Copula函数下使用。
#u<- pobs(as.matrix(cbind(X_1,X_2)))[,1]
#v<- pobs( as.matrix(cbind(x_1,x_2)))[,2]
#selectedCopula <- BiCopSelect(u, v,familyset=NA)
#selectedCopula

#将数据映射到累计概率分布值上，因为copula函数的自变量是对应边际分布的累计概率分布值。这个量是服从[0,1]上的均匀分布的。

##################################################
# 二维高斯cCopula函数
gaussian.cop <- normalCopula(dim=2)
# 生成二维高斯cCopula函数
set.seed (500)


#使用pobs函数将数据映射到累计概率分布值上，因为copula函数的自变量是对应边际分布的累计概率分布值。这个量服从[0,1]上的均匀分布的。
m <- pobs(as.matrix(cbind(X_1,X_2)))
fit <- fitCopula(gaussian.cop,m,method='ml')#用数据集拟合高斯Copula函数，对其参数进行估计
#获得参数估计
coef(fit)
rho <- coef(fit)[1]
rhogaussain <- rho
#画出拟合后的联合概率密度函数
persp(normalCopula(dim=2,rho) ,dCopula)
#画出拟合后的联合概率密度函数在二维的投影。
u<- rCopula(3965,normalCopula(dim=2,rho))
plot(u[,1],u[,2],pch= '.',col='blue ')
#使用以上高斯copula函数，对数据的分布进行联合（这里假设原始数据的分布是付出正太分布的)。
cor(u ,method= 'spearman' )
cor(u ,method= 'kendall' )

X_1_mu <- mean(X_1)
X_1_sd <- sd(X_1)
X_2_mu <- mean(X_2)
X_2_sd <- sd(X_2)
copula_dist <- mvdc(copula=normalCopula(rho,dim=2),margins=c("norm" , "norm"),paramMargins=list( list(mean=X_1_mu, sd=X_1_sd),list(mean=X_2_mu,sd=X_2_sd)))
#对联合后的分布，进行模拟，生成模拟数据，并与远数据比较。
sim <- rMvdc(3965, copula_dist)
plot(X_1,X_2,main= 'relation')
points(sim[,1],sim[,2],col='red' ,pch='.')
legend( 'bottomright',c( 'Observed' , 'Simulated'),col=c( 'black' , 'red'),pch=21)

#####################################
#t-coupla
t.cop<- tCopula(dim=2)
set.seed(500)
# 生成t-Copula函数

#使用pobs函数将数据映射到累计概率分布值上，因为copula函数的自变量是对应边际分布的累计概率分布值。
m <- pobs(as.matrix(cbind(X_1,X_2)))
fit <- fitCopula(t.cop,m,method='ml') #拟合进行参数估计
#获得参数估计
coef(fit)
rho <- coef(fit)[1]
df <- coef(fit)[2]
#画出拟合后的联合概率密度函数
persp(normalCopula(dim=2,rho) ,dCopula)
#画出拟合后的联合概率密度函数在二维的投影。
u <- rCopula(3965, tCopula(dim=2,rho,df=df))
plot(u[,1],u[,2],pch='.',col='blue')
#对数据的分布进行联合（这里假设原始数据的分布是付出正太分布的)
cor(u ,method= 'spearman' )
cor(u ,method= 'kendall' )

X_1_mu <- mean(X_1)
X_1_sd<- sd(X_1)
X_2_mu <- mean(X_2)
X_2_sd<- sd(X_2)

copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df),margins=c("norm","norm"),paramMargins=list(list(mean=X_1_mu,sd=X_1_sd),list(mean=X_2_mu,sd=X_2_sd)))
#对联合后的分布，进行模拟，生成模拟数据，并与远数据比较。
sim <- rMvdc(3965 , copula_dist)
plot(X_1,X_2, main= 'relation')
points(sim[,1],sim[,2],col='red',pch='.')
legend('bottomright',c('Observed','Simulated'),col=c('black','red' ),pch=21)

clayton.cop <- claytonCopula(dim=2)
set.seed( 500)
m <- pobs(as.matrix(cbind(X_1,X_2)))
fit <- fitCopula(clayton.cop , m ,method='ml')
coef(fit)

alpha <- coef(fit)[1]
persp( claytonCopula(dim=2, alpha),dCopula)

u <- rCopula(3965 , claytonCopula(dim=2, alpha))
plot(u[,1],u[,2],pch='.' ,col='blue')
cor(u, method= 'spearman' )
cor(u ,method= 'kendall' )

X_1_mu <- mean(X_1)
X_1_sd<- sd(X_1)
X_2_mu <- mean(X_2)
X_2_sd<- sd(X_2)
copula_dist<-mvdc(copula=claytonCopula(dim=2,alpha),margins=c("norm","norm"),paramMargins=list(list(mean=X_1_mu,sd=X_1_sd),list(mean=X_2_mu,sd=X_2_sd)))

sim <- rMvdc (3965, copula_dist)
plot(X_1,X_2,main= 'relation')
points(sim[,1],sim[,2],col='red' ,pch='.')
legend('bottomright' ,c('Observed' ,'Simulated'),col=c('black','red' ),pch=21)

gumbel.cop <- gumbelCopula(dim=2)
set.seed(500)
m<- pobs(as.matrix(cbind(X_1,X_2)))
fit<- fitCopula(gumbel.cop,m,method='ml')
coef(fit)
alpha <- coef(fit)[1]

persp(gumbelCopula(dim=2,alpha) ,dCopula)
u <-rCopula(3965,gumbelCopula(dim=2,alpha))
plot(u[,1],u[ ,2],pch='.' ,col='blue')
cor(u,method='spearman')
cor(u ,method= 'kendall' )

X_1_mu <- mean(X_1)
X_1_sd <- sd(X_1)
X_2_mu <- mean(X_2)
X_2_sd <- sd(X_2)

copula_dist<-mvdc(copula=gumbelCopula(dim=2,alpha),margins=c("norm","norm"),paramMargins=list(list(mean=X_1_mu, sd=X_1_sd),list(mean=X_2_mu, sd=X_2_sd)))

sim <- rMvdc(3965,copula_dist)
plot(X_1,X_2,main='relation')
points(sim[,1],sim[,2],col='red' ,pch='.')
legend('bottomright',c('Observed','Simulated'),col=c( 'black' , 'red'),pch=21)

frank.cop <- frankCopula(dim=2)
set.seed(500)
m<- pobs(as.matrix(cbind(X_1,X_2)))
fit<- fitCopula(frank.cop,m,method='ml')
coef(fit)
alpha <- coef(fit)[1]

persp(frankCopula(dim=2,alpha) ,dCopula)
u <-rCopula(3965,frankCopula(dim=2,alpha))
plot(u[,1],u[ ,2],pch='.' ,col='blue')
cor(u,method='spearman')
cor(u ,method= 'kendall' )