Til.hah.tun.2010
					M.Nuortio
					Pseudo-R-koodia
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.


   #
   # Alustetaan ensin aineisto.
   #


> Snowfall <- c(126.4 , 82.4 , 78.1 , 51.1 , 90.9 , 76.2 , 104.5 ,
    87.4 , 110.5 , 25.0 , 69.3 , 53.5 , 39.8 , 63.6 , 46.7 , 72.9 ,
    79.6 , 83.6 , 80.7 , 60.3 , 79.0 , 74.4 , 49.6 , 54.7 , 71.8 ,
    49.1 , 103.9 , 51.6 , 82.4 , 83.6 , 77.8 , 79.3 , 89.6 , 85.5 ,
    58.0 , 120.7 , 110.5 , 65.4 , 39.9 , 40.1 , 88.7 , 71.4 , 83.0 ,
    55.9 , 89.9 , 84.8 , 105.2 , 113.7 , 124.7 , 114.5 , 115.6 ,
    102.4 , 101.4 , 89.8 , 71.5 , 70.9 , 98.3 , 55.5 , 66.1 , 78.4 ,
    120.5 , 97.0 , 110.0)


   #
   # Ratkaistaan tehtava kahdella eri tavalla. Ensimmaisessa tavassa
   # "huijataan": R:ssa on paketti, jonka nojalla se osaa luoda
   # ydinestimaattorit suoraan. Seuraavanlaiset ohjeet ovat saatavilla
   # nettiosoitteesta
   # http://stat.ethz.ch/R-manual/R-patched/library/stats/html/density.html
   #
   #
   # The (S3) generic function density computes kernel density
   # estimates. Its default method does so with the given kernel and
   # bandwidth for univariate observations.
   #
   # Usage
   #
   # density(x, ...)
   # ## Default S3 method:
   # density(x, bw = "nrd0", adjust = 1,
   #         kernel = c("gaussian", "epanechnikov", "rectangular",
   #                    "triangular", "biweight",
   #                    "cosine", "optcosine"),
   #         weights = NULL, window = kernel, width,
   #         give.Rkern = FALSE,
   #         n = 512, from, to, cut = 3, na.rm = FALSE, ...)
   #
   # Arguments
   # x 	the data from which the estimate is to be computed.
   # bw 	the smoothing bandwidth to be used. The kernels are
   # scaled such that this is the standard deviation of the smoothing
   # kernel. (Note this differs from the reference books cited below,
   # and from S-PLUS.) bw can also be a character string giving a rule
   # to choose the bandwidth. See bw.nrd. The specified (or computed)
   # value of bw is multiplied by adjust.
   # adjust 	the bandwidth used is actually adjust*bw. This makes it
   # easy to specify values like ‘half the default’ bandwidth.
   # kernel, window 	a character string giving the smoothing kernel
   # to be used. This must be one of "gaussian", "rectangular",
   # "triangular", "epanechnikov", "biweight", "cosine" or "optcosine",
   # with default "gaussian", and may be  abbreviated to a unique
   # prefix (single letter).
   # "cosine" is smoother than "optcosine", which is the usual ‘cosine’
   # kernel in the literature and almost MSE-efficient. However, "cosine"
   # is the version used by S.
   # weights 	numeric vector of non-negative observation weights,
   # hence of same length as x. The default NULL is equivalent to
   # weights = rep(1/nx, nx) where nx is the length of (the finite
   # entries of) x[].
   # width 	this exists for compatibility with S; if given, and bw
   # is not, will set bw to width if this is a character string, or to
   # a kernel-dependent multiple of width if this is numeric.
   # give.Rkern 	logical; if true, no density is estimated, and
   # the ‘canonical bandwidth’ of the chosen kernel is returned instead.
   # n 	the number of equally spaced points at which the density is to
   # be estimated. When n > 512, it is rounded up to the next power of
   # 2 for efficiency reasons (fft).
   # from,to 	the left and right-most points of the grid at which the
   # density is to be estimated; the defaults are cut * bw outside of
   # range(x).
   # cut 	by default, the values of from and to are cut
   # bandwidths beyond the extremes of the data. This allows the
   # estimated density to drop to approximately zero at the extremes.
   # na.rm 	logical; if TRUE, missing values are removed from x.
   # If FALSE any missing values cause an error.
   # ... 	further arguments for (non-default) methods.
   #


> plot(density(Snowfall , bw=1 , adjust=1 , kernel="gaussian" ,
    from=0 , to=150),asp=3000)


> plot(density(Snowfall , bw=6 , adjust=1 , kernel="gaussian" ,
    from=0 , to=150),asp=3000)


> plot(density(Snowfall , bw=15 , adjust=1 , kernel="gaussian" ,
    from=0 , to=150),asp=3000)


> plot(density(Snowfall , bw=30 , adjust=1 , kernel="gaussian" ,
    from=0 , to=150),asp=3000)


> plot(density(Snowfall , bw=2 , adjust=1 , kernel="rectangular" ,
    from=0 , to=150),asp=3000)


> plot(density(Snowfall , bw=8 , adjust=1 , kernel="rectangular" ,
    from=0 , to=150),asp=3000)


> plot(density(Snowfall , bw=16 , adjust=1 , kernel="rectangular" ,
    from=0 , to=150),asp=3000)


> plot(density(Snowfall , bw=40 , adjust=1 , kernel="rectangular" ,
    from=0 , to=150),asp=3000)


   #
   # Ylla piirretyt kuvat ovat hyvin samanlaisia kuin ne, joita
   # naytettiin laskuharjoituksissa.
   #


   #
   # Koodataan seuraavaksi ydinestimaattorit rehellisesti kasi-
   # pelilla.
   #


> GaussYdin <- function(x,h){   1 / sqrt(pi) / sqrt(2) / h * exp( -0.5
                                * x^2 / h^2 )   }


> TasYdin <- function(x,h){   if ( abs(x) <= h ) { 1 / 2 / h }
                              else { 0 }  }


> GaussEst <- function(x,h){   PaluuArvo <- 0   ;
    for( i in 1:63 ){
      PaluuArvo <- PaluuArvo + GaussYdin( x - Snowfall[i] , h ) }   ;
    1 / 63 * PaluuArvo   }


> TasEst <- function(x,h){   PaluuArvo <- 0   ;
    for( i in 1:63 ){
      PaluuArvo <- PaluuArvo + TasYdin( x - Snowfall[i] , h ) }   ;
    1/63 * PaluuArvo   }


   #
   # Alustetaan x-arvojen hila kuvaajien piirtamista varten ja
   # piirretaan kuvaajat. Olkoon hila vaikka sellainen, jossa on
   # 1024 pistetta.
   #


> x.hila <- seq( from=0 , to=150 , by=150/1023 )


> plot( x.hila , GaussEst( x.hila , 1 ) , type='l' , asp=3000 )


> plot( x.hila , GaussEst( x.hila , 6 ) , type='l' , asp=3000 )


> plot( x.hila , GaussEst( x.hila , 15 ) , type='l' , asp=3000 )


> plot( x.hila , GaussEst( x.hila , 30 ) , type='l' , asp=3000 )


   #
   # HUOM! Skalaari- ja vektorisuureiden yhteensopimattomuusongelmien
   # takia tasaisen ytimen kuvaajia ei saa piirrettya yhta helposti.
   # Joudutaan menettelemaan vahan hienovaraisemmin.
   #


> y.hila <- c(1:1024)   ;
    for( j in 1:1024 ){ y.hila[j] <- TasEst( x.hila[j] , 2 ) }   ;
    plot( x.hila , y.hila , type='l' , asp=3000)


> y.hila <- c(1:1024)   ;
    for( j in 1:1024 ){ y.hila[j] <- TasEst( x.hila[j] , 8 ) }   ;
    plot( x.hila , y.hila , type='l' , asp=3000)


> y.hila <- c(1:1024)   ;
    for( j in 1:1024 ){ y.hila[j] <- TasEst( x.hila[j] , 16 ) }   ;
    plot( x.hila , y.hila , type='l' , asp=3000)


> y.hila <- c(1:1024)   ;
    for( j in 1:1024 ){ y.hila[j] <- TasEst( x.hila[j] , 40 ) }   ;
    plot( x.hila , y.hila , type='l' , asp=3000)