source('utils.R')
PL = iris$Petal.Length

w = diff( range( PL ) )
bw0 = w/5 # an initial bin size

# Plot the number of modes as a function of h (the bandwidth) for the iris data
#
adjusts = seq( 0.1, 0.23, length.out=75 ) # good for plotting
adjusts = seq( 0.1, 1.4, length.out=200 ) # large enough to find a k=1 mode
n_of_modes = rep( NA, length(adjusts) )
for( ai in 1:length(adjusts) ){
    a = adjusts[ai]
    nm = count_number_of_modes( PL, bw0*a )
    n_of_modes[ai] = nm
}
#postscript("../../WriteUp/Graphics/Chapter10/ex_5_hk_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( bw0*adjusts, n_of_modes, type='b', pch=16, xlab='bandwidth (h)', ylab='estimated number of modes (k)' )
grid()
#dev.off()

# Based on the above plot I'll take:
#
# h_4 = min( adjusts[ n_of_modes == 4 ] ) # how I calculated the numbers below
#
v=c()
for( n in 1:7 ){
  m = min( adjusts[ n_of_modes == n ] )
  v = c( v, m ) 
}
h_k = bw0 * v

s2 = var( PL ) # the sample variance
c_k = 1 / sqrt( 1 + h_k^2 / s2 )

# Now do the Silverman bootstraps:
#
max_n_modes = 5
P_k = rep( NA, max_n_modes ) # the proportion of times where the number of modes exceedes k
nboot = 200
for( k in 1:max_n_modes ){ # how many modes
  nm_from_boot_samples = rep( NA, nboot )
  for( bi in 1:nboot ){
      X_star = sample( PL, replace=TRUE )
      Z = rnorm( length(PL) )
      Y_star = c_k[k] * ( X_star + h_k[k] * Z )
      nm = count_number_of_modes( Y_star, h_k[k] )
      nm_from_boot_samples[bi] = nm
  }
  P_k[k] = sum( nm_from_boot_samples >= k ) / nboot
}

# Make a table of number of modes, the critical widths used, and P_k:
#
#postscript("../../WriteUp/Graphics/Chapter10/ex_5_pk_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( 1:max_n_modes, P_k, type='b', pch=19, cex=1.25, xlab='k', ylab='P_k', main='' )
grid()
#dev.off()