Til.hah.tun. kev2010
						M. Nuortio
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 'license()' or 'licence()' for distribution details.
     R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
     Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

    #
    # Seuraavassa annetaan tarpeelliset esitiedot.
    #

> Sigma1 <- matrix(c(2,1,1,2),nrow=2,byrow=T)

> Sigma2 <- matrix(c(2,-1,-1,2),nrow=2,byrow=T)

> mu1 <- c(0,0)

> mu2 <- c(0,2)

> library(MASS)

    #
    # Seuraavassa ohjelmoidaan itse luokitin.
    # Testien perusteella luokitin toimii (eivat osana
    # tata tulostetta).
    #

> KvadrLause <- function(x){   0.5 * t(x - mu1) %*% ginv(Sigma1) %*% (x - mu1)
                 - 0.5 * t(x - mu2) %*% ginv(Sigma2) %*% (x - mu2)
                 - 0.5 * log( det(Sigma2) / det(Sigma1) )   }

> KvadrLuok <- function(x){
                if ( KvadrLause(x) <= 0 ) { PaluuArvo <- 1 }
                else { PaluuArvo <- 2 }   ;   PaluuArvo   }

    #
    # Luodaan otokset jakaumista 1 ja 2.
    # HUOM! Vertaa aiemmat tietokoneharjoitukset.
    # Kuvien piirtelyn nojalla ovat oikein
    # (eivat osana tata tulostetta).
    #

> U1 <- eigen(Sigma1)$vectors   ;   U2 <- eigen(Sigma2)$vectors

> Lambda1_nj <- diag(sqrt(eigen(Sigma1)$values))   ;
     Lambda2_nj <- diag(sqrt(eigen(Sigma2)$values))

> d <- NROW(mu1)   ;   n <- 5000   ;   N <- (d * n)   ;
     Otos1 <- matrix(c(1:N),nrow=n,byrow=T)   ;
     for(i in 1:n){  x.randv.iid <- rnorm(d,0,1) ;
        y.randv <- (U1 %*% Lambda1_nj %*% x.randv.iid + mu1) ;
        for(j in 1:d){  Otos1[i,j] <- y.randv[j]   }   }   ;
     Otos2 <- matrix(c(1:N),nrow=n,byrow=T)   ;
     for(i in 1:n){  x.randv.iid <- rnorm(d,0,1) ;
        y.randv <- (U2 %*% Lambda2_nj %*% x.randv.iid + mu2) ;
        for(j in 1:d){  Otos2[i,j] <- y.randv[j]   }   }

    #
    # Seuraavassa ohjelmoidaan laskurit virheille.
    #

> Vaarat1 <- 0   ;
     for(i in 1:n){
        if ( KvadrLuok(Otos1[i,]) == 2 ) { Vaarat1 <- Vaarat1 + 1 }
        else { Vaarat1 <- Vaarat1 }   }   ;
     Vaarat2 <- 0   ;
     for(i in 1:n){
        if ( KvadrLuok(Otos2[i,]) == 1 ) { Vaarat2 <- Vaarat2 + 1 }
             else { Vaarat2 <- Vaarat2 }   }

> Vaarat1   ;   Vaarat2
[1] 1113
[1] 1177

> Vaarat1 / n   ;   Vaarat2 / n
[1] 0.2226
[1] 0.2354

> 0.5 * ( Vaarat1 / n  +  Vaarat2 / n )
[1] 0.229