M. Nuortio    Til.hah.tun.2010
					Harjoitus 4.5
R version 2.8.1 (2008-12-22)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
      R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type licenseO or licenceO for distribution details.
      R is a collaborative project with many contributors.
Type contributorsO for more information and
   citationO on how to cite R or R packages in publications.
      Type demo(?) for some demos, helpO for on-line help, or
   help.startO for an HTML browser interface to help.
Type qO to quit R.


> Cyksi <- matrix(c(1.75,-1.299,-1.299,3.25),nrow=2,byrow=T)

> Ckaksi <- matrix(c(5.625,3.375,3.375,5.625),nrow=2,byrow=T)

> byksi <- c(4,0)

> bkaksi <- c(-5,0)

> Uyksi <- eigen(Cyksi)$vectors

> Ukaksi <- eigen(Ckaksi)$vectors

> Lambda_yksi <- diag(eigen(Cyksi)$values)

> Lambda_kaksi <- diag(eigen(Ckaksi)$values)

> Lambda_yksi_nj <- diag(sqrt(eigen(Cyksi)$values))

> Lambda_kaksi_nj <- diag(sqrt(eigen(Ckaksi)$values))

> Cyksi_nj <- Uyksi %*% Lambda_yksi_nj %*% t(Uyksi)

> Ckaksi_nj <- Ukaksi %*% Lambda_kaksi_nj %*% t(Ukaksi)

> Lambda_yksi_nj_k <- diag(1/sqrt(eigen(Cyksi)$values))

> Lambda_kaksi_nj_k <- diag(1/sqrt(eigen(Ckaksi)$values))

    # HUOM! Kaanteismatriisit voitaisiin laskea jollakin
    # kaanteismatriisifunktiollakin suoraan, mutta
    # varmuuden vuoksi ohjelmoin sen noin "matalalla tasolla".
    # Kaanteismatriisifunktio voisi olla esimerkiksi
    # numeerisesti epastabiilimpi - en tosin tieda,
    # onko oikeasti.

> Cyksi_nj_k <- Uyksi %*% Lambda_yksi_nj_k %*% t(Uyksi)

> Ckaksi_nj_k <- Ukaksi %*% Lambda_kaksi_nj_k %*% t(Ukaksi)

> V <- eigen(Cyksi_nj_k %*% Ckaksi %*% Cyksi_nj_k)$vectors

> A <- t(V) %*% Cyksi_nj_k

> Ayksi <- Uyksi %*% Lambda_yksi_nj

> Akaksi <- Ukaksi %*% Lambda_kaksi_nj

> d <- NROW(byksi)   ;   n <- 1200   ;   N <- (d * n)   ;

> Pistematriisi <- matrix(c(1:N),nrow=n,byrow=T)   ;
   for(i in 1:(n/2)){   x.randv.iid <- rnorm(d,0,1) ;
     y.randv <- (Ayksi %*% x.randv.iid + byksi) ;
     for(j in 1:d){   Pistematriisi[i,j] <- y.randv[j]   }   }
   for(i in (n/2 + 1):n){   x.randv.iid <- rnorm(d,0,1) ;
     y.randv <- (Akaksi %*% x.randv.iid + bkaksi) ;
     for(j in 1:d){   Pistematriisi[i,j] <- y.randv[j]   }   }

> plot(Pistematriisi[,1],Pistematriisi[,2])

    # HUOM! Tuossa yllakin voisi olla plot(...,...,asp=1)
    # kuten alempana. Talloin kuvia voisi vertailla
    # samassa mittakaavassa. Voisi myos antaa talle
    # pistematriisille nimen Pistematriisi_1 ja
    # alla olevalle Pistematriisi_2, niin ne sailyisivat
    # yhta aikaa muistissa.

> for(i in 1:(n/2)){   x.randv.iid <- rnorm(d,0,1) ;
     y.randv <- A %*% (Ayksi %*% x.randv.iid + byksi) ;
     for(j in 1:d){    Pistematriisi[i,j] <- y.randv[j]   }   }
   for(i in (n/2 + 1):n){   x.randv.iid <- rnorm(d,0,1) ;
     y.randv <- A %*% (Akaksi %*% x.randv.iid + bkaksi) ;
     for(j in 1:d){   Pistematriisi[i,j] <- y.randv[j]   }   }

> plot(Pistematriisi[,1],Pistematriisi[,2],asp=1)