# 
# layout-käsu võimaluste demonstreerimine
#
# Viimati muudetud - 18. juulil 2005
# Koostanud: Märt Möls
#


# Millega png-pilte vaadata
vaataja="C:\\Program Files\\Mozilla Firefox\\firefox.exe"

# Kuhu pildfailid salvestada
tee="F:\\kursus\\R\\temp\\komb\\"

algus=function(pilt, xs, ys){
assign("teenimi", paste(tee,pilt,sep=""), env=.GlobalEnv)
assign("teenimi2", paste('"',vaataja,'" ',tee,pilt,sep=""), env=.GlobalEnv)
png(filename = teenimi, width = xs, height = ys, pointsize = 12, bg = "white")
}



algus("pilt20.png", 300, 300)
mat=diag(1:3)
layout(mat=mat, widths=rep(1,3), heights=rep(1,3))
par(mar=c(2,2,1,1))
x=rnorm(100); y1=rnorm(100); y2=rnorm(100)+1*x; y3=rnorm(100)+3*x;
plot(x,y1, pch=21, col="darkblue", bg=rgb(0,0,runif(100)), cex=1.5)
plot(x,y2, pch=21, col="darkblue", bg=rgb(0,0,runif(100)), cex=1.5)
plot(x,y3, pch=21, col="darkblue", bg=rgb(0,0,runif(100)), cex=1.5)
dev.off()
system(teenimi2)


algus("pilt21.png", 300, 300)
mat=rbind(
c(1,2,3),
c(4,5,6))
layout(mat=mat, widths=c(1,1,1), heights=c(1,2))
par(mar=c(3,3,1,1))
x1=rnorm(100)
x2=rexp(100)
x3=runif(100)
qqnorm(x1, main="", ylab=""); qqline(x1);
qqnorm(x2, main="", ylab=""); qqline(x2);
qqnorm(x3, main="", ylab=""); qqline(x3);
hist(x1, main="", col=terrain.colors(18), ylab="")
hist(x2, main="", col=terrain.colors(18), ylab="")
hist(x3, main="", col=terrain.colors(18), ylab="")
dev.off()
system(teenimi2)




algus("pilt22.png", 300, 300)
mat=rbind(
c(1,1,1),
c(1,2,1),
c(1,1,1))
layout(mat=mat, widths=c(4,4,1), heights=c(1,3,3))
x=rexp(100)
hist(x, col=topo.colors(12), xlab="", ylab="", main="")
par(mar=c(3,3,1,1))
plot(ecdf(x), xlab="", ylab="", main="")
dev.off()
system(teenimi2)




algus("pilt23.png", 800, 500)

# Ühisjaotuse hindamiseks vajame moodulit KernSmooth
library(KernSmooth)

# Algandmete genereerimine:
n=10000; a=matrix(c(rchisq(n, df=2),rchisq(n, df=2)),n ,2)

# Ekraani jaotuse defineerimine:
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(5,1), c(1,3), TRUE)

# Ühisjaotuse hindamine ja jooniste tegemine
est<-bkde2D(a,bandwidth=c(0.2,0.2),range.x=list(c(-1,6),c(-1,6)))

# 1. joonis (ühisjaotus)
par(mar=c(5,4,1,1))
contour(est$x1,est$x2,est$fhat,xlim=c(-1,6),ylim=c(-1,6), xlab=expression(x[1]), ylab=expression(x[2]))

# 2. joonis (x1 marginaaljaotus)
par(mar=c(0,4,1,1)) 
plot(density(a[,1]), xlim=c(-1,6), main=expression(x[1]), xaxt="n", ylim=c(0,0.5))

# 3. joonis
par(mar=c(5,0,1,1)); abi=density(a[,2])

plot(abi$y, abi$x,
main="", 
xlab="Density",
ylab=expression(x[2]),
type="l", 
ylim=c(-1,6), 
xlim=c(0,0.5))

dev.off()
system(teenimi2)






algus("pilt24.png", 500, 500)

tabel=read.csv2("E:/konsult/down/emavanus.csv")
tabel2=tabel[2:39, 2:16]

DEN=tabel2/t(matrix(rep(colSums(tabel2),38),15,38))

layout(matrix(rbind(
c(1,1,1),
c(0,2,0),
c(0,3,0)),3,3), c(1,3,1),c(3,1,1))

image(1989:2003,13:50,t( as.matrix(DEN)), col=terrain.colors(50), xlab="Aasta", ylab="Sünnitaja vanus", breaks=seq(0,0.085,length=51))
contour(1989:2003,13:50,t( as.matrix(DEN)), add=TRUE, col="brown")
par(xpd=NA)
lines(c(1992, 1992), c(13,50), lwd=1)
text(1992.25, 14, "A", cex=1.5)
text(1992.25, 49, "B", cex=1.5)

text(2002.3, 14, "A'", cex=1.5)
text(2002.3, 49, "B'", cex=1.5)

lines(c(2002, 2002), c(13,50), lwd=1)
# arrows(1992, 12, 1993.5, 4, length=0.1)
# lines(c(2002, 2002), c(13,13-30), lwd=1)
# arrows(2002, 13-30, 1999.5, 13-36, length=0.1)
par(xpd=FALSE)

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

bar2=function(korgus, varv, piirid, xlim, ylim, a, ...){
for(i in 1:length(korgus)){
for(j in 1:(length(piirid[piirid<korgus[i]]))){
  rect(i-a, piirid[j], i-a+1, min(piirid[j+1],korgus[i]) , col=varv[j], border=varv[j]) }
  rect(i-a, 0, i-a+1, korgus[i], border="gray40")
}}

# varv=rep(NA,length(DEN[,4]))
# for(i in 1:length(DEN[,4])){
# abi=(seq(0,0.085,length=51)+0.5*0.085/51-DEN[i,4])^2
# varv[i]=which(abi==min(abi))
# }

par(mar=c(3,1,3,0))
abi=barplot(DEN[,4], main="Sünnitajate vanusejaotus 1992", xlab="Vanus", ylim=c(0,0.1), col=terrain.colors(50), space=0, axes=FALSE)
axis(side=1, at=abi[seq(3, 38, 5)], label=seq(15,50,5), las=3)

bar2(DEN[,4], terrain.colors(50), seq(-0.0001,0.085,length=51), a=1)
text(0, 0.02, "A", cex=1.5)
text(38, 0.02, "B", cex=1.5)

abi=barplot(DEN[,14], main="Sünnitajate vanusejaotus 2002", xlab="Vanus", ylim=c(0,0.1), col=terrain.colors(50), space=0, axes=FALSE)
axis(side=1, at=abi[seq(3, 38, 5)], label=seq(15,50,5), las=3)

bar2(DEN[,14], terrain.colors(50), seq(-0.0001,0.085,length=51), a=1)

text(0, 0.02, "A'", cex=1.5)
text(38, 0.02, "B'", cex=1.5)

title("\n Sünnitajate vanusejaotuse muutumine", outer=TRUE, cex.main=1.4)

dev.off()
system(teenimi2)





algus("pilt25.png", 750, 500)

x=rchisq(100,4)*2+35

a1=dnorm(x1, mean=pi, sd=sqrt(2*pi))

toep=function(x,keskv,disp){
exp(colSums(log(sapply(keskv,dnorm,x=x,sd=sqrt(disp)))))
}

d=seq(30,50,0.1)
y=toep(x,d,30)
plot(d, y,type="l", main="Tõepärafunktsioon",xlab=expression(theta),
ylab=expression(L(theta)))
dmax=d[which(y==max(y))]
arrows(dmax-10,max(y)*0.5,dmax,max(y)*0.55)
text(dmax-10,max(y)*0.48,expression(paste("Suurima tõepära hinnang ", hat(theta))))
lines(c(dmax,dmax),c(0,max(y)),lty=2)

# Eeljaotus: mu=30, tau=3

# Järeljaotuse parameetrid:
jarel=function(mu,tau,x){
n=length(x)
sig2=var(x)
tau2=tau*tau
taup2=1/(n/sig2+1/tau2)
mup=taup2*(n*mean(x)/sig2+mu/tau2)
vastus=list(mup=mup,tau2p=taup2)
class(vastus)=c("Posterior","numeric")
return(vastus)
}
print.Posterior=function(x){
cat("Järeljaotuse dispersioon:",x$tau2p,"\n")
cat("             standardhälve:",sqrt(x$tau2p),"\n")
cat("             täpsus:",1/x$tau2p,"\n")
cat("             keskmine:",x$mup,"\n")
cat("             95%-tol. intervall: ",qnorm(0.025,mean=x$mup,sd=sqrt(x$tau2p)),"...",
qnorm(0.975,mean=x$mup,sd=sqrt(x$tau2p)),"\n")
}

jarel(30,4,x)

# Teeme A4-trükiks sobiva graafikaakna

# Joonistame tõepärafunktsiooni

mat=
cbind(
c(1,1,1,1,1,1,1),
c(1,2,1,3,1,4,1),
c(1,1,1,1,1,1,1))
nf=layout(mat,widths=c(1,4,7),heights=c(2,4,1,4,1,4,3))

d=seq(25,55,0.1)
y=toep(x,d,var(x))
plot(d, y,type="l", main="Tõepärafunktsioon ja järeljaotused erinevate eeljaotuste korral",xlab=expression(theta),
ylab=expression(L(theta)),lty=2)

legend(48,max(y)*0.9,c("tõepärafunktsioon","eeljaotus","järeljaotus"),lty=c(2,1,1),lwd=c(1,1,2),cex=1.5)


a=jarel(30,4,x)
y1=dnorm(d,mean=30,sd=4)
y2=dnorm(d,mean=a$mup,sd=sqrt(a$tau2p))

par(mar=c(2,1,2,1))
plot(d,y1, main="Minu eeljaotus",xlab=expression(paste("Reieümbermõõdu keskväärtus - ", theta)),ylab=expression(p(theta)), ylim=c(0,max(c(y1,y2))),type="l")
lines(d,dnorm(d,mean=a$mup,sd=sqrt(a$tau2p)), lwd=2)


a=jarel(30,100,x)
y1=dnorm(d,mean=30,sd=100)
y2=dnorm(d,mean=a$mup,sd=sqrt(a$tau2p))

plot(d,y1, main="Mitteinformatiivne eeljaotus",xlab=expression(paste("Reieümbermõõdu keskväärtus - ", theta)),ylab=expression(p(theta)), ylim=c(0,max(c(y1,y2))),type="l")
lines(d,dnorm(d,mean=a$mup,sd=sqrt(a$tau2p)), lwd=2)


a=jarel(30,0.5,x)
y1=dnorm(d,mean=30,sd=0.5)
y2=dnorm(d,mean=a$mup,sd=sqrt(a$tau2p))

plot(d,y1, main="Informatiivne ja vale eeljaotus",xlab=expression(paste("Reieümbermõõdu keskväärtus - ", theta)),ylab=expression(p(theta)), ylim=c(0,max(c(y1,y2))),type="l")
lines(d,dnorm(d,mean=a$mup,sd=sqrt(a$tau2p)), lwd=2)

dev.off()
system(teenimi2)