# 
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

k = 3 # a constant of H_m
M = 5 # the number of classes

r = 0:20 # a range of possible observations

H1 = dpois( r, 1*k )
H2 = dpois( r, 2*k )
H3 = dpois( r, 3*k )
H4 = dpois( r, 4*k )
H5 = dpois( r, 5*k )

#postscript("../../WriteUp/Graphics/Chapter2/prob_2.3.6.eps", onefile=FALSE, horizontal=FALSE)

plot( r, H1, type="l", col=1 ); text( 0.9, 0.2, "H_1" )
lines( r, H2, col=2 ); text( 7.3, 0.15, "H_2" )
lines( r, H3, col=3 ); text( 10.3, 0.13, "H_3" )
lines( r, H4, col=4 ); text( 13.3, 0.115, "H_4" )
lines( r, H5, col=5 ); text( 18.5, 0.08, "H_5" )

#dev.off()

mprime = 1:M
n_mprime_to_m = k / log( 1 + 1/mprime )

# Lets make sure our formulas for the lower and upper bounds of the hypthesis regions for H_m are correct: 
# 
ms = 1:(M-1)
ubs = floor( k / log( 1 + 1/ms ) )
lbs = c( 1, ubs+1 )
print(lbs)
print(ubs)

# Lets plot the decision region in the l_1 and l_2 plane:
#
P_0 = P_1 = P_2 = 1./3
sigma = 2
m_0 = matrix( c(0,0,0), nrow=3 )
m_1 = matrix( c(1,1,1), nrow=3 )
m_2 = matrix( c(-1,-1,-1), nrow=3 )

nPts = 100
l1_grid = seq( -10, +10, length.out=nPts )
l2_grid = seq( -10, +10, length.out=nPts )

Lambda_1 = exp( -( as.double( t(m_1) %*% m_1 - t(m_0) %*% m_0 )/(2*sigma^2) ) ) * exp( l1_grid )
Lambda_2 = exp( -( as.double( t(m_2) %*% m_2 - t(m_0) %*% m_0 )/(2*sigma^2) ) ) * exp( l2_grid )

#postscript("../../WriteUp/Graphics/Chapter2/prob_2.3.6_decision_region.eps", onefile=FALSE, horizontal=FALSE)

# First decision region :
#
Pt_1 = P_1 * Lambda_1
Pt_2 = - ( P_0 - P_2 * Lambda_2 )
D1 = outer( Pt_1, Pt_2, FUN="+" )
contour( l1_grid, l2_grid, D1, levels=c(0), xlab="l_1", ylab="l_2", col="black" )

# Second decision region :
#
Pt_1 = as.double( - P_0 * ones( nPts, p=1 ) ) # first argument to outer must be a discretization of l_1
Pt_2 = P_2 * Lambda_2 # second argument to outer must be a discretization of l_2
D2 = outer( Pt_1, Pt_2, FUN="+" )
contour( l1_grid, l2_grid, D2, levels=c(0), add=T, xlab="l_1", ylab="l_2", col="red" )

# Third decision region :
#
Pt_1 = - ( P_0 + P_1 * Lambda_1 )
Pt_2 = P_2 * Lambda_2
D3 = outer( Pt_1, Pt_2, FUN="+" )
contour( l1_grid, l2_grid, D3, levels=c(0), add=T, xlab="l_1", ylab="l_2", col="green" )

text(-5,-5,"H_0")
text(5,-5,"H_1")
text(7,5,"H_1")
text(-5,5,"H_2")

#dev.off()