knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(MATH4793LABStayl0048)
data <- read.delim("C://Users//micha//OneDrive//Desktop//Stats 2//MVN Consistency//T4-3.DAT", header = FALSE, sep="", skip=0, as.is=TRUE)
data=data[,1:4]
normal_mvn=function(x, sigma)
{
obs_count=nrow(x)
columns=ncol(x)
datamat=matrix(nrow=1, ncol=columns)
column_count=0
resultmat=matrix(nrow=1, ncol=columns)

empirical_rule_generator=function(sd_count, iter=20) {

sigma=sd_count 
count=0
data=matrix(nrow=iter, ncol=1)

while(count <= iter)
{
x=rnorm(1000000)
low=-sigma*sd(x)
high=sigma*sd(x)
in_interval=length(which(x > low & x < high))
ratio=in_interval/1000000
data[count,]=ratio
count=count+1
}
mean(data)
}
p=empirical_rule_generator(sd_count=sigma, iter=30)
q=1-p

while(column_count<=columns)
{
dat=x[,column_count]
dat=as.numeric(dat)
mean_adjusted=dat-mean(dat)
sdev=sd(mean_adjusted)

low=-sigma*sd(mean_adjusted)
high=sigma*sd(mean_adjusted)
interval=c(low,high)
in_interval=length(which(mean_adjusted > low & mean_adjusted < high))

ratio=(in_interval)/(obs_count)-p
norm_ind=3*sqrt((p*q)/obs_count)
indicator=ratio-norm_ind
resultmat[1,column_count]=ratio
column_count=column_count+1
}

trigger=any(resultmat < 0)

if(trigger == "TRUE")
{
barplot(resultmat, main="The data is not consistent with multivariate normality")
}

if(trigger == "FALSE")
{
barplot(resultmat, main="The data is consistent with multivariate normality")
}

}
normal_mvn(data, 1)


miketay123/MATH4793tayl0048 documentation built on April 19, 2022, 1:27 a.m.