macho=c(120,107,110,116, 114, 111, 113,117,114,112)
femea=c(110,111,107, 108,110,105,107,106,111,111)
media.m=mean(macho)
media.m
media.f=mean(femea)
media.f
#vetor.media
chacal=c(macho,femea)
sexo=factor(rep(c("macho","femea"),each=10))

boxplot(chacal~sexo)
plot(1:20,chacal, pch=rep(c(15,16),each=10),col=rep(1:2,each=10))
	for(i in 1:10)
	{
	lines(c(i,i),c(chacal[i],media.m),col=1)
	}
	for(j in 11:20)
	{
	lines(c(j,j),c(chacal[j],media.f),col=2)
	}
lines(c(1,10),c(media.m,media.m),col=1)
lines(c(11,20),c(media.f,media.f),col=2)

mean(chacal)
sd(chacal)
mean(femea)-mean(macho)
mean(macho)-mean(femea)
dif=abs(mean(femea)-mean(macho))
dif

hist(chacal, freq=FALSE,xlim=c(95,125))
curve(exp=dnorm(x, mean=mean(chacal),sd=sd(chacal)),from=95,to=125, col="red", add=T)

rnorm(10,mean=mean(chacal),sd=sd(chacal))
round(rnorm(10,mean=mean(chacal),sd=sd(chacal)))
abs(round(rnorm(10,mean=mean(chacal),sd=sd(chacal))))
abs(round (mean(rnorm(10,mean=mean(chacal),sd=sd(chacal)))-mean(rnorm(10,mean=mean(chacal),sd=sd(chacal)))))

for(i in 1:10)
{
cat("\n\t", i)
}

resulta=rep(NA,10)
for (i in 1:10) resulta[i]=abs(mean(round(rnorm(10,mean=mean(chacal),sd=sd(chacal))))-	mean(round(rnorm(10,mean=mean(chacal),sd=sd(chacal)))))
	
source("simula.r")
dif.chacal=simula(macho,femea)

dif.chacal=simula(macho,femea,nsim=2000)
table(dif.chacal)
n.maior=sum(dif.chacal>=dif)
n.maior
n.maior/length(dif.chacal)

dif.chacal.uni=simula(macho,femea, teste="uni")
dif.chacal.uni=simula(macho,femea, teste= "uni", nsim=2000)
table(dif.chacal.uni)

n.maior=sum(round(dif.chacal.uni,1)>=round(dif,1))
n.menor=sum(round(dif.chacal.uni,1)<=round((dif*-1),1))
n.maior
n.menor
n.sim=2000

p.bi=(n.maior+n.menor)/2000
p.bi

## Qual a probabilidade de erra ao fazer a afirmação que as mandibulas dos machos são maiores?
p.uni=n.maior/n.sim
p.uni

v1=var(macho)
v1
v2=var(femea)
v2
n1=length(macho)
n1
n2=length(femea)
n2

##desvio padrão das diferenças
s12=sqrt((v1/n1)+(v2/n2))
s12
dif
tvalor=dif/s12
tvalor

res.t= simula(macho,femea,nsim=2000, teste="t")

##
## agora vamos calcular as probabilidades
maior.menor.t=sum(res.t>=tvalor | res.t<=-tvalor)
maior.menor.t
prob.t=maior.menor.t/2000
prob.t
## compare com o resultado do teste t da função t.test()
t.test(macho,femea)

are=c(6,10,8,6,14,17,9,11,7,11)
are
arg=c(17,15,3,11,14,12,12,8,10,13)
arg
hum=c(13,16,9,12,15,16,17,13,18,14)
hum
solos=data.frame(are,arg,hum)
solos
str(solos)
boxplot(solos)

media.geral=mean(c(are,arg,hum))
media.geral
dif.geral=solos-media.geral
dif.geral
sum(dif.geral)
round(sum(dif.geral),10)
ss.solos=dif.geral^2
ss.solos
ss.total=sum(ss.solos)
ss.total

vetor.dados<-c(are,arg,hum)
vetor.obs<-1:length(c(are,arg,hum))
vetor.cor<-rep(c(1,2,3),each=10)

plot(vetor.obs,vetor.dados,ylim=c(0,20),pch=(rep(c(15,16,17),each=10)),col=vetor.cor,ylab="Variável Resposta", xlab="Observações")
	for(i in 1:30)
	{
	lines(c(i,i),c(vetor.dados[i],mean(vetor.dados)),col=vetor.cor[i])
	}
	abline(h=media.geral)



vetor.medias<-rep(rep(c(mean(are),mean(arg),mean(hum)),each=10))
media.solos<-vetor.medias[c(1,11,21)]
plot(vetor.obs,vetor.dados,ylim=c(0,20),pch=(rep(c(15,16,17),each=10)),col=vetor.cor,main="Variação Intra Grupos",ylab="Variável Resposta", xlab="Observações")
	for(i in 1:30)
	{
	lines(c(i,i),c(vetor.medias[i],vetor.dados[i]),col=vetor.cor[i])
	}
	lines(c(1,10),c(media.solos[1],media.solos[1]),col=1)
	lines(c(11,20),c(media.solos[2],media.solos[2]),col=2)
	lines(c(21,30),c(media.solos[3],media.solos[3]),col=3)

solos
media.solos
ss.are=sum((are-media.solos["are"])^2)
ss.are
ss.arg=sum((arg-media.solos["arg"])^2)
ss.arg
ss.hum=sum((hum-media.solos["hum"])^2)
ss.hum
ss.intra=ss.are+ss.arg+ss.hum
ss.intra


	plot(vetor.obs,vetor.medias,ylim=c(5,16),pch=(rep(c(15,16,17),each=10)),col=vetor.cor,main="Variação Entre Grupos",ylab="Variável Resposta", xlab="Observações")
		for(i in 1:30)
		{
		lines(c(i,i),c(vetor.medias[i],mean(vetor.medias)),col=vetor.cor[i])
		}
	abline(h=media.geral)
	points(vetor.obs,vetor.dados,ylim=c(0,20),pch=(rep(c(0,1,2),each=10)),col=vetor.cor,cex=0.5)

#### Cálculo dos valores

media.solos=apply(solos,2,mean)
media.solos
media.geral
ss.entre=10*sum((media.solos-media.geral)^2)
ss.entre

ss.intra+ss.entre
ss.total

ms.entre=ss.entre/2
ms.intra=ss.intra/27
ms.entre
ms.intra
F.solos=ms.entre/ms.intra
F.solos

curve(pf(x,2,27, lower.tail=FALSE),from=0,to=10, main="Disribuição F de Fisher", xlab="Valor F",ylab="Probabilidade",xlim=c(0,10))
abline(v=F.solos,col="red")
lines(c(-1,F.solos),c(1-pf(F.solos,2,27),1-pf(F.solos,2,27)),lty=2,col="red")

curve(expr=df(x, 2,27),main="Disribuição F de Fisher", xlab="Valor F",ylab="Densidade",xlim=c(0,10))
abline(v=F.solos,col="red")

p.solos=pf(F.solos,2,27, lower.tail=FALSE)
p.solos

str(solos)
var.resp=c(solos$are,solos$arg,solos$hum)
var.resp
solos.f=factor(rep(c("are", "arg","hum"),each=10))
solos.f
res.anova=aov(var.resp~solos.f)
res.anova
summary(res.anova)