set.seed( 1357908642 )

mu <- 70.0

sigma.H <- 1.0

sigma.P <- 25.0

I <- 10

J <- 150

y <- matrix( NA, I, J )

a.H <- rep( NA, I )

a.P <- matrix( NA, I, J )

for ( i in 1:I ) {

  a.H[ i ] <- rnorm( 1, 0, sigma.H )

  for ( j in 1:J ) {

    a.P[ i, j ] <- rnorm( 1, 0, sigma.P )

    y[ i, j ] <- round( mu + a.H[ i ] + a.P[ i, j ], 1 )

  }

}

print( y )

print( y.dot.dot <- mean( y ) )

# 68.81193

print( y.i.dot <- apply( y, 1, mean ) )

# [1] 68.18733 70.54733 67.33733 70.44800 67.28933 66.21200 
# [7] 70.24000 71.21533 66.43467 70.20800

print( SS.H <- J * sum( ( y.i.dot - y.dot.dot )^2 ) )

# 4912.045

print( df.H <- I - 1 )

# 9

print( MS.H <- SS.H / df.H )

# 545.7828

SS.P <- 0

for ( i in 1:I ) {

  for ( j in 1:J ) {

    SS.P <- SS.P + ( y[ i , j ] - y.i.dot[ i ] )^2

  }

}

print( SS.P )

# 909749.6

print( df.P <- I * ( J - 1 ) )

# 1490

print( MS.P <- SS.P / df.P )

# 610.5702

print( sigma.hat.P.2 <- MS.P )

# 610.5702

print( sigma.hat.H.2 <- ( MS.H - MS.P ) / J )

# -0.431916

dput( y, "variance-components-example-data-1.txt" )