代写 R 


MAT 4376F/5313F

R – Fitting multivariate distributions

library(mvtnorm);
library(MASS);
library(QRM);
library(evir);
# Simulate bivariate normal data
mu <- c(0,0) # Mean Sigma <- matrix(c(1, .9, .9, 1), 2) # Covariance matrix n=1000; Xnorm<- rmvnorm(n, mean=mu, sigma=Sigma, method="chol") par(mfrow=c(1,2)) qqnorm(Xnorm[,1]); qqnorm(Xnorm[,2]); # Marginally data are normal par(mfrow=c(1,2)) hist(Xnorm[,1]); hist(Xnorm[,2]); # scatter plot par(mfrow=c(1,1)) plot(Xnorm[,1],Xnorm[,2]); # Calculate kernel density estimate Bivnorm.kde <- kde2d(Xnorm[,1], Xnorm[,2], n = 50) # from MASS package par(mfrow=c(1,1)) persp(Bivnorm.kde, phi = 45, theta = 30, shade = .1, border = NA); # Contour plot for bivariate normal density x.points <-seq(-3,3,length.out=100); y.points <- x.points; z <-matrix(0,nrow=100,ncol=100); mu <-c(0,0); for (i in 1:100) {for (j in 1:100) {z[i,j] <-dmvnorm(c(x.points[i],y.points[j]),mean=mu,sigma=Sigma)}} contour(x.points,y.points,z) # Contour plot for the bivariate data contour(Bivnorm.kde) # Simulate bivariate t Sigma <- matrix(c(1, .9, .9, 1), 2) # Sigma matrix nu <- 4 # number of degrees of freedom Xt<- rmvt(n, sigma=Sigma, df=nu); par(mfrow=c(1,2)) qqnorm(Xt[,1]); qqnorm(Xt[,2]); # Marginally data are not normal # scatter plot par(mfrow=c(1,1)) plot(Xt[,1],Xt[,2]); # Calculate kernel density estimate Bivt.kde <- kde2d(Xt[,1], Xt[,2], n = 50) # from MASS package par(mfrow=c(1,1)) persp(Bivt.kde, phi = 45, theta = 30, shade = .1, border = NA); # Comparison of bivariate normal and t par(mfrow=c(1,2)) plot(Xnorm[,1],Xnorm[,2]); plot(Xt[,1],Xt[,2]); par(mfrow=c(1,2)) persp(Bivnorm.kde, phi = 45, theta = 30, shade = .1, border = NA); persp(Bivt.kde, phi = 45, theta = 30, shade = .1, border = NA); par(mfrow=c(1,2)) contour(Bivnorm.kde) contour(Bivt.kde) # Fitting bivariate t distribution fit <- fit.mst(Xt) # fit a multivariate t distribution mu.fit <- fit$mu # estimated location vector; Sigma,fit <- as.matrix(fit$Sigma) # estimated scale matrix; nu.fit <- fit$df # estimated degrees of freedom fit <- fit.mst(Xnorm) # fit a multivariate t distribution to normal data mu.fit <- fit$mu # estimated location vector; Sigma.fit <- as.matrix(fit$Sigma) # estimated scale matrix; nu.fit <- fit$df # estimated degrees of freedom # Working with real data data(bmw); data(siemens); X1=-bmw; X2=-siemens; par(mfrow=c(1,2)); plot.ts(X1); plot.ts(X2); par(mfrow=c(1,2)); qqnorm(X1); qqnorm(X2); par(mfrow=c(1,1)); plot(X1,X2); Biv.kde <- kde2d(X1, X2, n = 50) # from MASS package par(mfrow=c(1,1)) persp(Biv.kde, phi = 45, theta = 30, shade = .1, border = NA); par(mfrow=c(1,1)) contour(Biv.kde) X=cbind(X1,X2); fit <- fit.mst(X) # fit a multivariate t distribution mu.fit <- fit$mu # estimated location vector; Sigma.fit <- as.matrix(fit$Sigma) # estimated scale matrix; nu.fit <- fit$df # estimated degrees of freedom # "Goodness of fit test" n=length(X1); Xnorm<- rmvnorm(n, mean=mu.fit, sigma=Sigma.fit, method="chol") Xt<- rmvt(n, sigma=Sigma.fit, df=nu.fit); Bivnorm.kde <- kde2d(Xnorm[,1], Xnorm[,2], n = 50); Bivt.kde <- kde2d(Xt[,1], Xt[,2], n = 50); Biv.kde <- kde2d(X1, X2, n = 50); par(mfrow=c(1,3)); persp(Biv.kde, phi = 45, theta = 30, shade = .1, border = NA); persp(Bivnorm.kde, phi = 45, theta = 30, shade = .1, border = NA); persp(Bivt.kde, phi = 45, theta = 30, shade = .1, border = NA); par(mfrow=c(1,3)); contour(Biv.kde) contour(Bivnorm.kde) contour(Bivt.kde)