diak=read.table("D:\\OKTATAS\\2016\\inf_bc\\data_16bc.csv",sep=";",header=T,dec=",")

summary(diak)
#alapstat


par(mfrow=c(2,2))

hist(diak[,7],xlab="jegy",ylab="gyakoriság",main="Analízis jegy")
hist(diak[,5],xlab="óra",ylab="gyakoriság",main="Tanulási idő")
hist(diak[,1],xlab="cm",ylab="gyakoriság",main="Testmagasság")
hist(diak[,3],xlab="",ylab="gyakoriság",main="Cipőméret")
#nem jo


par(mfrow=c(2,2))

u1=hist(diak[,7],xlab="jegy",ylab="gyakoriság",main="Analízis jegy",breaks=c(1:6)-1/2)
hist(diak[,5],xlab="óra",ylab="gyakoriság",main="Tan.id\u151",breaks=c(0:8)*5-1/2)
hist(diak[,1],xlab="cm",ylab="gyakoriság",main="Testmagasság",breaks=5*c(1:9)+152.5)
hist(diak[,3],xlab="",ylab="gyakoriság",main="Cip\u151méret",breaks=c(0:13)+69/2)

summary(diak)
#
par(mfrow=c(2,1))

hist(diak[,6],xlab="perc",ylab="gyakoriság",main="Utazási id\u151",breaks=10*c(1:25)-5)
hist(diak[,2],xlab="kg",ylab="gyakoriság",main="Testsúly",breaks=5*c(6:24)+2.5)

par(mfrow=c(2,2))

boxplot(diak[,1]~diak[,4],ylab="cm",xlab="nem",main="Testmagasság")
boxplot(diak[,2]~diak[,4],ylab="kg",xlab="nem",main="Testsúly")
boxplot(diak[,5]~diak[,4],ylab="h",xlab="nem",main="Tan.id\u151")
boxplot(diak[,7]~diak[,4],ylab="jegy",xlab="nem",main="Jegy")





#hist(diak[diak[,8]== "M",5],xlab="cm",ylab="gyakoriság",main="Testmagasság",breaks=5*c(1:10)+142.5)


#ez a jo
#P-R
par(mfrow=c(1,1))
plot(density(diak[,6]),main="ut.id\u151")
par(mfrow=c(1,2))
plot(density(diak[,6],adjust=2),main="ut.id\u151")
plot(density(diak[,6],adjust=0.5),main="ut.id\u151")
##############################
#homerseklet adatok

homb<-read.table("d:\\oktatas\\2016\\inf_a\\nyir-51-88m.hom")
homk<-read.table("d:\\oktatas\\2016\\inf_a\\karc-51-88.hom")

r=c(0:37)

jan1b<-homb[365*r+1,1]
jan1k<-homk[365*r+1,1]
summary(jan1b)

#ez a jan 1 



par(mfrow=c(1,1))
dat<-jan1b
v<-list()
v$" Nyíregyháza "<-jan1b
v$" Karcag "<-jan1k
boxplot(v,at=c(1:2)+0.1,boxwex=0.5,border=1, main="Január 1-i középhomérsékletek")
#PROBALKOZZUNK  mas beallitasokkal


#tapasztalati elo.fv
#eloszor 12 megf
par(mfrow=c(1,1))
x = rnorm(12)
Fn = ecdf(x)
plot(Fn,xlim=c(-3.2,3.2))
x=c(1:200)
x=(x-100)/6
lines(x,pnorm(x),col=2)

#ez 120 elemu
x <- rnorm(120)
Fn <- ecdf(x)
#lines(Fn,col=3,pch=0.5)

lines(Fn, verticals=F, col.points='blue',
     col.hor='red', col.vert='bisque',cex=0.3)

#ez pedig mar 480
x <- rnorm(480)
Fn <- ecdf(x)
#lines(Fn,col=3,pch=0.5)

lines(Fn, verticals=F, col.points='green',
     col.hor='red', col.vert='bisque',cex=0.1)



###
#becslesek vizsgalata
####

n=c(10,40,160,640)
ered=n
est=matrix(0,1000,4)
for(j in 1:1000)
{
for (i in 1:length(n))
{x=rpois(n[i],1)
est[j,i]=mean(x)}
}

plot(n,est[1,],type="l",main="A Poisson elo. paraméter-becslése",ylab="lambda",ylim=c(min(est[1,],0.99),max(est[1,],1.01)))
#ez csak 1 

for (i in 1:4)
ered[i]=mean(est[,i])
lines(n,ered,col=4)
abline(h=1,col=2)

legend("bottomright",c("egy becslés","1000 becslés átlaga","a paraméter értéke"),lty=c(1,1,1),col=c(1,4,2))

#hogyan tudnank jobban merni az elterest?
#szorodasi meroszamokat is nezzunk

n=c(10,40,160,640)
est=matrix(0,1000,4)
for(j in 1:1000)
{
for (i in 1:length(n))
{x=rexp(n[i],1)
est[j,i]=1/mean(x)}
}
ered=n
for (i in 1:4)
ered[i]=mean((est[,i]-1)^2)
#átlagos négyzetes hiba -gyorsan fogy

plot(n,ered,type="l",ylab="Átl.négyzetes hiba")
abline(h=1,col=2)


####
#ML becsles, idén kimaradt
###
norm.mle <- function(dat){
	
	likelihood <- function(x){
		m=x[1]
		sigma=x[2]
		n=length(dat)
		if(sigma<=0)
			10^6
		else
			-(-n*(log(sqrt(2*pi))+log(sigma))-sum((dat-m)^2)/(2*sigma^2))
	}

	x=optim(c(0,1), likelihood)
	
	z=list()
	z$mle=x$par
	z$nllh=x$value
	z
}

d=rnorm(100)

r=norm.mle(d)

library(MASS)
fitdistr(d, "normal")


###
#Cauchy pelda
####

# neg log likelihood sfvbol
mlogl <- function(mu, x) {
sum(-dcauchy(x, location = mu, log = TRUE))
}
# es kepletbol
mlogl2 <- function(mu, x) {
sum(log(1 + (x - mu)^2))
}
#reprodukalhato legyen a szimulacio
n <- 30
set.seed(42)
x <- rcauchy(n)
#µ=0 de mi becsuljuk

mu.start <- median(x)
 mu.start
#[1] -0.1955062
out <- nlm(mlogl, mu.start, x = x)
mu.hat <- out$estimate
mu.hat
#[1] -0.1816501
#masikbol
out2 <- nlm(mlogl2, mu.start, x = x)
mu.hat <- out2$estimate
mu.hat
#[1] -0.1816501
#ugyanaz
############
# Szimulacio: melyik a jobb: a median vagy az ML becsles
 nsim <- 100
mu.hat=rep(0,times=nsim)
mu.twiddle=mu.hat
mu.atl=mu.hat
 mu <- 0
  for (i in 1:nsim) {
 xsim <- rcauchy(n, location = mu)
 mu.start <- median(xsim)
 out <- nlm(mlogl, mu.start, x = xsim)
 mu.hat[i] <- out$estimate
 mu.twiddle[i] <- mu.start
mu.atl[i]=mean(xsim)
 }
 mean((mu.hat - mu)^2)
#
 mean((mu.twiddle - mu)^2)
#
#Az MSE-hanyados:
mean((mu.atl - mu)^2)
mean((mu.hat - mu)^2)/mean((mu.twiddle - mu)^2)

###
#innentől: konf. int
###

a <- 5
sig <- 2
n <- 20
xdat=rnorm(n,a,sig)
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error
hist(xdat)
abline(v=left,col=2,lty=2)
abline(v=right,col=2,lty=2)


a <- 5
sig <- 2
n <- 2000
xdat=rnorm(n,a,sig)
error <- qt(0.995,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error
hist(xdat)
abline(v=left,col=2,lty=2)
abline(v=right,col=2,lty=2)

#ellenorizzuk a lefedesi vszget
#nov. 29

talal=0
m=1000
for (i in 1:m)
{a <- 5
sig <- 2
n <- 20
xdat=rnorm(n,a,sig)
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error
if (left<a && a<right) talal=talal+1
}
talal

talal=0
m=1000
for (i in 1:m)
{a <- 5
sig <- 2
n <- 2000
xdat=rnorm(n,a,sig)
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error
if (left<a && a<right) talal=talal+1
}
talal

#mi van, ha nem normalis az elo

talal=0
m=1000
for (i in 1:m)
{a <- 5
sig <- 2
n <- 20
xdat=rexp(n,1/a)
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error
if (left<a && a<right) talal=talal+1
}
talal

talal=0
m=1000
for (i in 1:m)
{a <- 5
sig <- 2
n <- 2000
xdat=rexp(n,1/a)
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error
if (left<a && a<right) talal=talal+1
}
talal
#n=20-ra van elteres


########
#konf.int
#######
diak=read.table("D:\\OKTATAS\\2016\\inf_bc\\data_16bc.csv",sep=";",header=T,dec=",")
dat=diak
n=dim(dat)[1]
xdat=dat[,1]

error <- qt(0.995,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error

left 
right
error
sd(xdat)


n=dim(dat)[1]
xdat=dat[,2]

error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error

left 
right
error
sd(xdat)

n=dim(dat)[1]
xdat=dat[,3]
#cipo
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error

left 
right
error
sd(xdat)

xdat=dat[!is.na(dat[,5]),5]
n=length(xdat)
#xdat=dat[,5]
#valszam
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error

left 
right
error
sd(xdat)

n=dim(dat)[1]
xdat=dat[,7]
#jegy
error <- qt(0.975,n-1)*sd(xdat)/sqrt(n)
left <- mean(xdat)-error
right <- mean(xdat)+error

left 
right
error
sd(xdat)



###############################
#t-próba
homb<-read.table("d:\\oktatas\\2016\\inf_a\\nyir-51-88jav.hom")
homk<-read.table("d:\\oktatas\\2016\\inf_a\\karc-51-88.hom")

r<-c(0:37)
#############
jan1b<-homb[365*r+1,1]
jan1k<-homk[365*r+1,1]

#ezek a jan. 1-i adatok

t.test(jan1b,alternative="l",mu=0)
#H1:mu<0 less

t.test(jan1b,alternative="t",mu=0)
#H1:mu=0 two sided

t.test(jan1b,alternative="t",mu=-2)
#H1:mu>-2 greater

p.v=rep(0,times=365)
for (i in 1:365)
{jan1b<-homb[365*r+i,1]
jan1k<-homk[365*r+i,1]
p.v[i]=t.test(jan1b,jan1k,paired=T,alternative="t")$p.value #ez a korrekt, hiszen parositott megfigyeleseink vannak

}

plot(p.v,type="l")


############
#nehezebbek-e az informatikusok, mint a survey-esek

dat=read.table("D:\\OKTATAS\\2016\\tatk\\dat16_ss.csv",sep=";",header=T,dec=",")
pdat=dat[,1]

xdat=diak[,2]
t.test(pdat,xdat)

xdat=diak[diak[,4]=="f",2]

pdat=dat[dat[,3]=="f",1]
t.test(pdat,xdat)



###############
#Wilcoxon proba
###############

x=rnorm(50)
y=rnorm(50,0.5)
t.test(x,y)

wilcox.test(x, y)

##de: valtoztassuk meg az adatokat
y[50]=40

t.test(x,y)

wilcox.test(x, y)

#############
#chi-negyzet
#############

diak=read.table("D:\\OKTATAS\\2016\\inf_bc\\data_16bc.csv",sep=";",header=T,dec=",")

#nem
#H0: egyenletes elo.
chisq.test(table(diak[,4]))
#tovabbi komponensek
chisq.test(table(diak[,4]))$observed
chisq.test(table(diak[,4]))$expected
chisq.test(table(diak[,4]))$residuals

#szimulalt krit ertekkel
chisq.test(table(diak[,4]),simulate.p.value=T)

###############
#11.ora
###############

diak=read.table("D:\\OKTATAS\\2016\\inf_bc\\data_16bc.csv",sep=";",header=T,dec=",")

#normalitasvizsgalat
#itt nincs tablazat
par(mfrow=c(2,2))

u1=hist(diak[,7],xlab="jegy",ylab="gyakoriság",main="Analízis jegy",breaks=c(1:6)-1/2)
hist(diak[,5],xlab="óra",ylab="gyakoriság",main="Tan.id\u151",breaks=c(0:8)*5-1/2)
hist(diak[,1],xlab="cm",ylab="gyakoriság",main="Testmagasság",breaks=5*c(1:9)+152.5)
hist(diak[,2],xlab="",ylab="gyakoriság",main="Testsúly",breaks=c(0:9)*10+69/2)


n<-4;j<-2 #testsúly
 #intervallum: 4 db ill. a valtozo sorszama
 q<-c(1:n)/n
 h<-c(-Inf,qnorm(q))
 mu<-mean(diak[,j])
sig=sd(diak[,j])
 h<-h*sig+mu
 #ezek hataroljak 
 nu<-rep(0,times=n)
 for (i in 1:n)
 nu[i]<-
  sum((diak[,j]<h[i+1]) & (diak[,j]>=h[i]))


 #vart gyak:
 np<-sum(nu)*rep(1/n,times=n)
 chi<-sum((nu-np)^2/np)


 #a p-ertek:
 p<-1-pchisq(chi,n-3)
 p

##############

n<-4;j<-1 #magasság
 #intervallum: 4 db ill. a valtozo sorszama
 q<-c(1:n)/n
 h<-c(-Inf,qnorm(q))
 mu<-mean(diak[,j])
sig=sd(diak[,j])
 h<-h*sig+mu
 #ezek hataroljak 
 nu<-rep(0,times=n)
 for (i in 1:n)
 nu[i]<-
  sum((diak[,j]<h[i+1]) & (diak[,j]>=h[i]))


 #vart gyak:
 np<-sum(nu)*rep(1/n,times=n)
 chi<-sum((nu-np)^2/np)


 #a p-ertek:
 p<-1-pchisq(chi,n-3)
 p

###

n<-4;j<-6 #tan.idő
#diak[1,j]=50
 #intervallum: 4 db ill. a valtozo sorszama
 q<-c(1:n)/n
 h<-c(-Inf,qnorm(q))
 mu<-mean(diak[,j])
sig=sd(diak[,j])
 h<-h*sig+mu
 #ezek hataroljak 
 nu<-rep(0,times=n)
 for (i in 1:n)
 nu[i]<-
  sum((diak[,j]<h[i+1]) & (diak[,j]>=h[i]))


 #vart gyak:
 np<-sum(nu)*rep(1/n,times=n)
 chi<-sum((nu-np)^2/np)


 #a p-ertek:
 p<-1-pchisq(chi,n-3)
 p

#####
#ftlensegvizsg: cipő vs nem
#####

table(diak[,4],diak[,3])
 chisq.test(table(diak[,4],diak[,3]))
 chisq.test(table(diak[,4],
            diak[,3]),simulate.p.value=T)

##########
#linearis regresszio
##############

homb<-read.table("d:\\oktatas\\2016\\inf_a\\nyir-51-88jav.hom")
homk<-read.table("d:\\oktatas\\2016\\inf_a\\karc-51-88.hom")
par(mfrow=c(1,1))
r=c(0:37)

#vajon egyenletesen no-e az atlaghom. marciusban?
atl<-rep(0,times=31)
for (i in 1:31) atl[i]=mean(homb[365*r+59+i,1])
t<-c(1:31)
lm1=lm(atl~t)
summary(lm1)

plot(t,atl)
lines(lm1$fitted.values)
#egesz jo

#vajon egyenletesen no-e az atlaghom. aprilisban?

t<-c(1:30)
atl=t
for (i in 1:30) atl[i]=mean(homb[365*r+90+i,1])

lm1=lm(atl~t)
summary(lm1)

plot(t,atl)
lines(lm1$fitted.values)
#sokkal valtozekonyabb

#majus
t<-c(1:31)
atl=t
for (i in 1:31) atl[i]=mean(homb[365*r+120+i,1])

lm1=lm(atl~t)
summary(lm1)

plot(t,atl)
lines(lm1$fitted.values)

#no es az egesz egyben (julius is):
t<-c(1:153)
atl=t
for (i in 1:153) atl[i]=mean(homb[365*r+59+i,1])

lm1=lm(atl~t)
summary(lm1)

plot(t,atl)
lines(lm1$fitted.values)

#jobb lenne kvadratikus

#lehet nemparameteresen is simitani (Nadarajah-módszer)
lines(loess.smooth(t, atl, span = 1/10, degree = 1,
    family = "gaussian", evaluation = 150),col=2)
lines(loess.smooth(t, atl, span = 1/3, degree = 1,
    family = "gaussian", evaluation = 150),col=4)


t<-c(1:153)
t2=t^2
atl=t
for (i in 1:153) atl[i]=mean(homb[365*r+59+i,1])

lm1=lm(atl~t+t2)
summary(lm1)

plot(t,atl)
lines(lm1$fitted.values)
#ez az igazi

lm1=lm(atl~t2)
summary(lm1)
plot(t,atl)
lines(lm1$fitted.values)


#sajat adataink
#elobb csinaljunk szamot a nembol

#diak=read.table("d:\\oktatas\\2016\\tatk\\dat16_ss.csv",header=T,sep=";")


#probaljuk a testmagassagot elore jelezni:
nem=rep(0,times=dim(diak)[1])
nem[diak[,4]=="f"]=1

lm1=lm(diak[,1]~diak[,2]+nem+diak[,3]+diak[,5]+diak[,6]+diak[,7])
summary(lm1)
#hagyjuk el a legrosszabbat

lm1=lm(diak[,1]~diak[,2]+nem+diak[,3]+diak[,7]+diak[,6])
summary(lm1)

lm1=lm(diak[,1]~nem+diak[,3]+diak[,7])
summary(lm1)

lm1=lm(diak[,1]~diak[,3]+diak[,7])
summary(lm1)

lm1=lm(diak[,2]~diak[,4])
summary(lm1)

plot(diak[,4],diak[,2])
lines(diak[,4],lm1$fitted.values)

#ha nem nezzuk a cipomeretet

lm1=lm(diak[,2]~diak[,1]+nem+diak[,5]+diak[,6]+diak[,7])
summary(lm1)

lm1=lm(diak[,2]~diak[,1]+nem+diak[,5]+diak[,6])
summary(lm1)

lm1=lm(diak[,2]~diak[,1]+nem+diak[,6])
summary(lm1)


#nem igazan normalis, tehat a p-ertekek sem pontosak

cor(cbind(diak[,1:2],nem,diak[,4:7]))


pairs(diak)

###
#idosorok
###

homb<-read.table("d:\\oktatas\\2016\\inf_a\\nyir-51-88jav.hom")

ido<-c(1:(38*365))
plot(ido/365,homb[,1],type="l")
 
ev<-c(1:365)

hombs365<-rep(0,times=365*37+1)
idos365<-c(1:(365*37+1))
idos365<-idos365+366/2-1

#MOZGÓÁTLAG 365 NAPRA
for (i in (0:(37*365)))
{
hombs365[i+1]<-sum(homb[ev+i,1])/365
#idos365[i+1]<-sum(ido[ev+i])/365
}

plot(idos365/365,hombs365,type="l")


acf(homb[,1],lag=730)
#nem stacionarius!

#valasszuk le az eves ciklust
#minden napra a napi atlaggal
r=c(0:37)

meanb<-rep(0, times=365)
for (i in 1:365)
meanb[i]<-mean(homb[365*r+i,1])

plot(meanb,type="l")
#simitsuk egy kicsit
lines(loess.smooth(c(1:365), meanb, span = 1/10, degree = 1,   family = "gaussian", evaluation = 150),col=4)
meanb2=loess.smooth(c(1:365), meanb, span = 1/10, degree = 1,   family = "gaussian", evaluation = 150)


#vonjuk le minden napbol a megfelelo atlagot
meanb<-rep(meanb,times=38)
hombudk<-homb[,1]-meanb
plot(hombudk,type="l")

acf(hombudk,lag=730)
#itt mar nincs ciklikussag

acf(hombudk)
arbud<-ar(hombudk)
#autoregresszios modell illesztes
arbud
acf(arbud$resid[8:length(arbud$resid)])
#ez mar iid-nek tekintheto