# Time-stamp: <2012-11-02 13:40 Petri Koistinen>

#
# Paras lineaarinen ennuste sellaisessa tapauksessa, jossa regressiofunktio
# antaa paremman ennusteen.  Yhteisjakauma on kuvan 7.1 jatkuva jakauma.
#

# Ytf
f <- function(x, y) { 
  dgamma(x, 2, 1) * 
  dnorm(y, mean = exp(x/3), sd = sqrt(abs(x)))
}

# Simuloidaan jakaumasta pisteitä visualisointia ja analysointia varten

N <- 1000000
N.small <- 1000
x <- rgamma(N, 2, 1)
y <- rnorm(N, mean = exp(x/3), sd = sqrt(abs(x)))
xvis <- x[1:N.small]
yvis <- y[1:N.small]

xg <- seq(0.0001, 8, len = 61)
yg <- seq(-4, 15, len = 61)
z <- outer(xg, yg, f) 

maxf <- max(z)

# Ytf:n tasa-arvokäyriä ja siitä simuloituja xy-pisteitä

levs <- c(0.9, 0.5, 0.1, 0.01) * maxf
contour(xg, yg, z, levels = levs, drawlabels = F)
points(xvis, yvis, pch = '.')

# Lasketaan Monte Carlo estimaattien avulla paras lineaarinen ennuste,
# ja piirretään se

a.opt <- mean(y)
b.opt <- cov(x, y) / var(x)

mx <- mean(x)

lines(xg, a.opt + b.opt * (xg - mx), col = 'blue')

# Lasketaan regressiofunktio (jakauman konstruktion perusteella),
# ja pirretään se

lines(xg, exp(xg / 3), col = 'red')

# Selite:

legend('topleft', lty = 1,
       c('Best linear predictor', 'Regression function'),
       col = c('blue', 'red'))