options(error=recover)

# See Chapter 14 Devore:
# 
expected_count_under_null = function( T ){
    n = sum( T )
    I = dim(T)[1]
    J = dim(T)[2]
    n_i_dot = rowSums( T )
    n_dot_j = colSums( T )
    E = matrix( n_i_dot, nrow=I, ncol=1 ) %*% matrix( n_dot_j, nrow=1, ncol=J ) / n
    return(E)
}


# Section 5; Question 1:
#
I = 2
J = 2
T = matrix( c( 628, 146, 172, 54 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 2:
#
I = 2
J = 2
T = matrix( c( 168, 73, 42, 26 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 3:
#
I = 4
J = 2
T = matrix( c( 24, 450, 29, 312, 35, 211, 23, 70 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 4:
#
I = 2
J = 2
T = matrix( c( 59, 27, 143, 349 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 5:
#
I = 2
J = 2
T = matrix( c( 301, 1442-301, 345, 1420-345 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 6:
#
I = 3
J = 3
T = matrix( c( 14, 11, 8, 11, 11, 9, 6, 10, 12 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 7:
#
I = 2
J = 2
T = matrix( c( 27, 47, 140, 93 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 8:
#
I = 4
J = 2
T = matrix( c( 2592, 1481, 159, 78, 800, 375, 667, 399 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 9:
#
I = 2
J = 2
T = matrix( c( 11, 10, 55, 24 ), nrow=I, ncol=J, byrow=TRUE )
E = expected_count_under_null( T )
x = sum( ( T - E )^2 / E )
df = (I-1)*(J-1)
p_value = 1-pchisq( x, df)
print( sprintf('d2= %f; df= %d; p_value= %f', x, df, p_value) )


# Section 5; Question 10:
#
n = 141
x = 14
m = 44
y = 5
p = (x+y)/(n+m)
z = ( (x/n) - (y/m) ) /sqrt( (p*(1-p))/n + (p*(1-p))/m) # following Theorem 9.4.1
print(z)