source('../../Data/data_loaders.R')

DF = load_particle_board_data()
DF$temperature = as.factor( DF$temperature )

# Implement Bartlett's test on this data:
#
# Following notes here: https://en.wikipedia.org/wiki/Bartlett%275 test
N = dim(DF)[1]
k = length( unique( DF$temperature ) )
n_i = 3 # 3 samples in each class

# The variance in each class:
#
si2 = tapply( DF$strength, DF$temperature, FUN=var )

# The pooled variance estimate:
#
sp2 = sum( (n_i - 1) * si2 ) / (N-k)

stat_numerator = (N-k) * log(sp2) - sum( (n_i-1) * log(si2) )
stat_denominator = 1 + ( k/(n_i-1) - 1/(N-k) ) / (3*(k-1))
B = stat_numerator / stat_denominator

P_value = 1 - pchisq( B, k-1 )
print( sprintf("B= %10.6f; p_value= %0.6f", B, P_value) )

# Using the canned R procedure:
# 
x = DF$strength
g = as.factor( DF$temperature )
print(bartlett.test(x, g))