In agracy2246/MATH4753grac0009: Adam Gracy MATH4753

knitr::opts_chunk$set(echo = TRUE)

Task 1

getwd()

mychisim<-function(n1=10,sigma1=3,mean1=5,iter=1000,ymax=0.1,x=20, y=0.1){    # adjust ymax to make graph fit
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1

data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1

ssq1=apply(data1.mat,2,var) # ssq1 is s squared

w=(n1-1)*ssq1/sigma1^2      #chi-sq stat

hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]," = ",n1," statistic = ",chi^2)),
xlab=expression(paste(chi^2, "Statistic",sep=" ")), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(dchisq(x,n1-1),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(chi^2==frac((n[1]-1)*s^2,sigma^2)) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title) # Legend #
invisible(list(w=w,summary=summary(w),sd=sd(w),fun="Chi-sq")) # some output to use if needed
}

Task 2

n_1=10,iter=1000,μ_1=10,σ_1=4

mychisim(n1 = 10, iter = 1000, mean1 = 10, sigma1 = 4,ymax = .12)

n_1=20,iter=1000,μ_1=10,σ_1=4

mychisim(n1 = 20, iter = 1000, mean1 = 10, sigma1 = 4, ymax = .11)

n_1=100,iter=1000,μ_1=10,σ_1=4

mychisim(n1 = 100, iter = 1000, mean1 = 10, sigma1 = 4, ymax = .04)

n_1=200,iter=1000,μ_1=10,σ_1=4

mychisim(n1 = 200, iter = 1000, mean1 = 10, sigma1 = 4, ymax = .03)

chisq hist

chisq <- mychisim(n1 = 10, iter = 1500, mean1 = 20, sigma1 = 10, ymax = .11)
hist(chisq$w)

Task 3

myTsim<-function(n1=10,sigma1=3,mean1=5,iter=1000,ymax=0.1,x=2,y=0.3,...){    # adjust ymax to make graph fit
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1

data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1

sd1=apply(data1.mat,2,sd) # sd
ybar=apply(data1.mat,2,mean)  # mean

w=(ybar-mean1)/(sd1/sqrt(n1))      #T stat

hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]," = ",n1," statistic = ",T," iterations= ",iter)),
xlab=expression(paste(T, "Statistic",sep=" ")), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(dt(x,n1-1),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(T==frac((bar(y)-mu),s/sqrt(n1))) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title) # Legend #
invisible(list(w=w,summary=summary(w),sd=sd(w),fun="T")) # some output to use if needed
}

n_1=10,iter=1000,μ_1=10,σ_1=4

myTsim(n1 = 10, iter = 1000, mean1 = 10, sigma1 = 4, ymax = .5)

n_1=20,iter=1000,μ_1=10,σ_1=4

myTsim(n1 = 20, iter = 1000, mean1 = 10, sigma1 = 4, ymax = .5)

n_1=100,iter=1000,μ_1=10,σ_1=4

myTsim(n1 = 100, iter = 1000, mean1 = 10, sigma1 = 4, ymax = .5)

n_1=200,iter=1000,μ_1=10,σ_1=4

myTsim(n1 = 200, iter = 1000, mean1 = 10, sigma1 = 4, ymax = .5)

T hist

T <- myTsim(n1 = 10, iter = 1500, mean1 = 20, sigma1 = 10)
hist(T$w)

Task 4

mychisim2<-function(n1=10,n2=14,sigma1=3,sigma2=3,mean1=5,mean2=10,iter=1000,ymax=0.07,x=40,y=0.04,...){    # adjust ymax to make graph fit
  y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
  y2=rnorm(n2*iter,mean=mean2,sd=sigma2)
  data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
  data2.mat=matrix(y2,nrow=n2,ncol=iter,byrow=TRUE)
  ssq1=apply(data1.mat,2,var) # ssq1 is s squared
  ssq2=apply(data2.mat,2,var)
  spsq=((n1-1)*ssq1 + (n2-1)*ssq2)/(n1+n2-2) # pooled s squared
  w=(n1+n2-2)*spsq/(sigma1^2)#sigma1=sigma2,  Chi square stat
  hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
  main=substitute(paste("Sample size = ",n[1]+n[2]," = ",n1+n2," statistic = ",chi^2)),
  xlab=expression(paste(chi^2, "Statistic",sep=" ")), las=1)
  lines(density(w),col="Blue",lwd=3) # add a density plot
  curve(dchisq(x,n1+n2-2),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
  title=expression(chi^2==frac((n[1]+n[2]-2)*S[p]^2,sigma^2)) #mathematical annotation -see ?plotmath
  legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title) # Legend #
  invisible(list(w=w,summary=summary(w),sd=sd(w),fun="Chi-sq")) # some output to use if needed
}

n_1=10,n_2=10,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=1000

mychisim2(n1 = 10, n2 = 10, mean1 = 5, mean2 = 10, sigma1=4,sigma2 = 4, iter=1000,x =30)

n_1=20,n_2=10,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=1000

mychisim2(n1 = 20, n2 = 10, mean1 = 3, mean2 = 5,sigma1 = 10,sigma2 = 10, iter=1000)

n_1=50,n_2=50,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=10000

mychisim2(n1 = 50, n2 = 50, mean1 = 5, mean2 = 10, sigma1=4,sigma2 = 4, iter=10000)

n_1=80,n_2=50,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=10000

mychisim2(n1 = 80, n2 = 50, mean1 = 3, mean2 = 5, sigma1=10,sigma2 = 10, iter=10000)

hist

chisq2 <-mychisim2(iter=10000)
hist(chisq2$w)

Task 5

Alter the function to place the legend on the mouse click

myTsim2<-function(n1=10,n2=14,sigma1=3,sigma2=3,mean1=5,mean2=10,iter=1000,ymax=0.5,...){
  y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
  y2=rnorm(n2*iter,mean=mean2,sd=sigma2)
  data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
  data2.mat=matrix(y2,nrow=n2,ncol=iter,byrow=TRUE)
  ssq1=apply(data1.mat,2,var) # ssq1 is s squared
  ybar1= apply(data1.mat,2,mean)
  ssq2=apply(data2.mat,2,var)
  ybar2=apply(data2.mat,2,mean)
  spsq=((n1-1)*ssq1 + (n2-1)*ssq2)/(n1+n2-2) # pooled s squared
  w=((ybar1-ybar2)-(mean1-mean2))/sqrt(spsq*(1/n1+1/n2))#sigma1=sigma2,  Chi square stat
  hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
       main=substitute(paste("Sample size = ",n[1]+n[2]," = ",n1+n2," statistic = ",T)),
       xlab=paste(" T Statistic",sep=""), las=1)
  lines(density(w),col="Blue",lwd=3) # add a density plot
  curve(dt(x,n1+n2-2),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
  title=expression(T==frac((bar(Y)[1]-bar(Y)[2])-(mu[1]-mu[2]),S[p]*sqrt(frac(1,n[1])+frac(1,n[2])))) #mathematical annotation -see ?plotmath
  legend(locator(1),c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title)# Legend #
  return(list(w=w,summary=summary(w),sdw=sd(w),fun="T")) # some output to use if needed
}

T stat

$\frac{(\bar Y_1 - \bar Y_2) - (\mu_1 - \mu_2)}{S_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}$

Where $\bar Y_1 ,\bar Y_2$ are the means of $Y_1,Y_2$.

$\mu_1,\mu_2$ being the sample means

$n_1, n_2$ being the size of samples 1 and 2;

$S_p$ being the deviation of the population

copy and paste

w=((ybar1-ybar2)-(mean1-mean2))/sqrt(spsq*(1/n1+1/n2))

myTsim2<-function(n1=10,n2=14,sigma1=3,sigma2=3,mean1=5,mean2=10,iter=1000,ymax=0.5,...){
  y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
  y2=rnorm(n2*iter,mean=mean2,sd=sigma2)
  data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
  data2.mat=matrix(y2,nrow=n2,ncol=iter,byrow=TRUE)
  ssq1=apply(data1.mat,2,var) # ssq1 is s squared
  ybar1= apply(data1.mat,2,mean)
  ssq2=apply(data2.mat,2,var)
  ybar2=apply(data2.mat,2,mean)
  spsq=((n1-1)*ssq1 + (n2-1)*ssq2)/(n1+n2-2) # pooled s squared
  w=((ybar1-ybar2)-(mean1-mean2))/sqrt(spsq*(1/n1+1/n2))#sigma1=sigma2,  Chi square stat
  hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
       main=substitute(paste("Sample size = ",n[1]+n[2]," = ",n1+n2," statistic = ",T)),
       xlab=paste(" T Statistic",sep=""), las=1)
  lines(density(w),col="Blue",lwd=3) # add a density plot
  curve(dt(x,n1+n2-2),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
  title=expression(T==frac((bar(Y)[1]-bar(Y)[2])-(mu[1]-mu[2]),S[p]*sqrt(frac(1,n[1])+frac(1,n[2])))) #mathematical annotation -see ?plotmath
  legend(2,.3,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title)# Legend #
  return(list(w=w,summary=summary(w),sdw=sd(w),fun="T")) # some output to use if needed
}

n_1=10,n_2=10,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=1000

myTsim2(n1=10,n2=10,sigma1=4,sigma2=4,mean1=5,mean2=10,iter=1000,ymax=0.5)

n_1=20,n_2=10,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=1000

myTsim2(n1=20,n2=10,sigma1=10,sigma2=10,mean1=3,mean2=5,iter=1000,ymax=0.5)

n_1=50,n_2=50,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=10000

myTsim2(n1=50,n2=50,sigma1=4,sigma2=4,mean1=5,mean2=10,iter=10000,ymax=0.5)

n_1=80,n_2=50,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=10000

myTsim2(n1=80,n2=50,sigma1=10,sigma2=10,mean1=3,mean2=5,iter=10000,ymax=0.5)

Tsim2 hist

tsimtwo<-myTsim2()
hist(tsimtwo$w)

Task 6

F Statistic

$F = (\frac{S^2_1}{S^2_2})(\frac{\sigma^2_2}{\sigma^2_1})$

Assumptions

No assumptions are made.

myFsim2<-function(n1=10,n2=14,sigma1=3,sigma2=2,mean1=5,mean2=10,iter=1000,ymax=0.9,x=6,y=0.5,...){
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
y2=rnorm(n2*iter,mean=mean2,sd=sigma2)
data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
data2.mat=matrix(y2,nrow=n2,ncol=iter,byrow=TRUE)
ssq1=apply(data1.mat,2,var) # ssq1 is s squared
ssq2=apply(data2.mat,2,var)
#spsq=((n1-1)*ssq1 + (n2-1)*ssq2)/(n1+n2-2) # pooled s squared
w=ssq1*sigma2^2/(ssq2*sigma1^2) #
hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]+n[2]," = ",n1+n2," statistic = ",F)),
xlab=paste("F Statistic",sep=""), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(df(x,n1-1,n2-1),xlim=c(0,6),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(F==frac(s[1]^2,s[2]^2)*frac(sigma[2]^2,sigma[1]^2)) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title)# Legend #
invisible(list(w=w,summary=summary(w),sd=sd(w),fun="F")) # some output to use if needed
}

Plot 1

myFsim2(n1 = 10, n2 = 10,sigma1 = 4, sigma2 =4, mean1 = 5, mean2 = 10 )

Plot 2

myFsim2(n1 = 20, n2 = 10,sigma1 = 4, sigma2 =4, mean1 = 10, mean2 = 10 )

Plot 3

myFsim2(n1 = 30, n2 = 10,sigma1 = 4, sigma2 =7, mean1 = 10, mean2 = 5 )

Plot 4

myFsim2(n1 = 20, n2 = 50,sigma1 = 3, sigma2 =4, mean1 = 14, mean2 = 10 )

F hist

F <- myFsim2()
hist(F$w)

Task 7

library(grac0009MATH4753)
data("fire")
knitr::kable(head(fire))

agracy2246/MATH4753grac0009 documentation built on April 26, 2020, 9:39 a.m.