################################################################# library( moments ) # unbuffer the output ################################################################# set.seed( 1 ) N <- 1000 percent.f <- 0.44 mu.f <- 135 sigma.f <- 23 mu.m <- 176 sigma.m <- 28 ################################################################# f.raw <- rchisq( N * percent.f, 3 ) m.raw <- rchisq( N * ( 1 - percent.f ), 3 ) f.su <- ( f.raw - mean( f.raw ) ) / sd( f.raw ) m.su <- ( m.raw - mean( m.raw ) ) / sd( m.raw ) f.star <- round( mu.f + sigma.f * f.su ) quantile( f.star ) # 0% 25% 50% 75% 100% # 108.0 119.0 127.5 147.0 263.0 c( mean( f.star ), sd( f.star ) ) # [1] 134.99773 22.99356 par( mfrow = c( 2, 1 ) ) plot( density( f.star ), lwd = 2, col = 'red', xlab = 'Weight', main = 'Female Sub-Population PDF', xlim = c( 90, 310 ), ylim = c( 0, 0.025 ) ) ################################################################# m.star <- round( mu.m + sigma.m * m.su ) quantile( m.star ) # 0% 25% 50% 75% 100% # 143.00 154.75 168.00 191.00 310.00 c( mean( m.star ), sd( m.star ) ) # [1] 176.00893 28.01083 plot( density( m.star ), lwd = 2, col = 'blue', xlab = 'Weight', main = 'Male Sub-Population PDF', xlim = c( 90, 310 ), ylim = c( 0, 0.025 ) ) ################################################################# population <- c( f.star, m.star ) c( mean( population ), sd( population ) ) # [1] 157.96000 32.95104 plot( density( population ), main = 'Population PDF', xlab = 'Weight', lwd = 2, xlim = c( 90, 310 ), ylim = c( 0, 0.025 ) ) ################################################################# plot( density( f.star ), lwd = 2, col = 'red', xlab = 'Weight', xlim = c( 90, 310 ), ylim = c( 0, 0.025 ), main = 'Sub-Population PDFs' ) lines( density( m.star ), lwd = 2, col = 'blue' ) text( 100, 0.02, 'F', col = 'red' ) text( 185, 0.015, 'M', col = 'blue' ) ################################################################# seed <- 1 set.seed( seed ) M <- 1000000 S.star.1 <- rep( NA, M ) n <- 192 system.time( for ( m in 1:M ) { S.star.1[ m ] <- sum( sample( population, n, replace = T ) ) if ( ( m %% 100000 ) == 0 ) print( m ) } ) system.time( S.star.2 <- apply( matrix( sample( population, n * M, replace = T ), M, n ), 1, sum ) ) par( mfrow = c( 1, 2 ) ) plot( density( S.star.1 ), main = 'PDF of S', xlab = 'Sum (S)', lwd = 2 ) plot( density( S.star.2 ), main = 'PDF of S', xlab = 'Sum (S)', lwd = 2 ) par( mfrow = c( 1, 1 ) ) qqplot( sample( S.star.1, 1000 ), sample( S.star.2, 1000 ), xlab = 'S.star.1', ylab = 'S.star.2' ) abline( 0, 1, lwd = 2, col = 'red' ) ################################################################# mu.population <- 158 sigma.population <- 33 E.s <- n * mu.population SD.s <- sigma.population * sqrt( n ) plot( density( S.star.2 ), main = 'PDF of S', xlab = 'Sum (S)', lwd = 2 ) s.grid <- seq( E.s - 3 * SD.s, E.s + 3 * SD.s, length = 500 ) lines( s.grid, dnorm( s.grid, E.s, SD.s ), lwd = 2, lty = 2, col = 'red' ) text( E.s, 0.9 * max( dnorm( s.grid, E.s, SD.s ) ), 'CLT', col = 'red' ) ################################################################# tolerance.equivalent <- E.s + ( 31400 - 30336 ) / ( 33 * sqrt( 192 ) ) * SD.s mean( S.star.2 > tolerance.equivalent ) # [1] 0.027859 1 - pnorm( ( tolerance.equivalent - E.s ) / SD.s ) # 0.009985386 ################################################################# skewness( population ) 0.8702212 kurtosis( population ) - 3 0.9603594 c( skewness( S.star.2 ), kurtosis( S.star.2 ) - 3 ) ----- sum S ----- n skewness kurtosis 1 0.870 0.960 2 0.618 0.497 10 0.274 0.089 192 0.061 -0.003 print( n.s <- round( 100 * skewness( population ) ) ) print( n.k <- round( 10 * ( kurtosis( population ) - 3 ) ) ) ################################################################# quantile( S.star.2 ) # 0% 25% 50% 75% 100% # 28250 30018 30324 30635 32647 quantile( S.star.2, probs = 1 - 0.0001 ) # 99.99% # 32088 #################################################################