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

windows(width=8, height=5)

# Ü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))

windows(width=5, height=5)

tabel=read.csv2(url("http://www.ms.ut.ee/mart/R/komb/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
 

windows(width=7.5, height=5)

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)

# 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)