# # 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. # #----- # Ex 1.2: pbinom(14,20,0.8) cs = 10:15 alphas = 1 - pbinom(cs-1,20,0.5) betas = pbinom(cs-1,20,0.8) ab_bind = rbind( alphas, betas ) colnames( ab_bind ) = cs print( ab_bind, digits=4 ) # Ex 1.3: n = 60 p = 0.5 alphas = c(0.01, 0.05, 0.1) cs = qnorm( 1 - alphas ) * sqrt( n * p * (1-p) ) + n * p + 0.5 # these are floating point numbers cs = ceiling( cs ) p_A = 0.7 betas = pnorm( ( cs - 1 + 0.5 - n*p_A ) / sqrt( n * p_A * (1-p_A) ) ) # Ex 1.4: p_A = 0.6 betas = pnorm( ( cs - 1 + 0.5 - n*p_A ) / sqrt( n * p_A * (1-p_A) ) ) # Ex 1.5: n = 60 p_A = 0.4 alphas = c(0.05, 0.1) cs = qnorm( 1 - alphas ) * sqrt( n * p * (1-p) ) + n * p + 0.5 # these are floating point numbers cs = ceiling( cs ) betas = pnorm( ( cs - 1 + 0.5 - n*p_A ) / sqrt( n * p_A * (1-p_A) ) ) # Ex 1.6: dbinom( 10, 10, 0.7 ) 1 - dbinom( 10, 10, 1-c(0.2,0.1,0.05,0.0) ) sum( dbinom( 9:10, 10, 0.7 ) ) pbinom( 8, 10, 1 - c(0.2,0.1,0.05,0.0) ) # Ex 1.7: 1 - pbinom( 47, 60, 0.7 ) # Ex 1.8: n = 80 p_H = 0.3 c = floor( sqrt( n * p_H * (1-p_H) ) * qnorm(0.05) + n * p_H - 0.5 ) beta = 1 - pbinom( c, n, 0.2 ) # find the value of p such that beta = 0.1 beta_fn = function(p){ 1 - pbinom( c, n, p ) - 0.1 } uniroot( beta_fn, c(0.1,0.4) ) # Ex 1.9: pbinom( c, n, 0.35 ) # Ex 1.10: n = 100 c = 11 p = 0.6 sum( dbinom( (-c+0.5*n) : (+c+0.5*n), n, p ) ) # Ex 1.11: n = 500 c = 185 p = 0.4 pnorm( ( c-1+0.5-n*p )/sqrt( n*p*(1-p) ) ) # Ex 1.12: 1 / choose(8,4) # Ex 1.13 (Part i): sum( dhyper( 0:5, 80, 1000-80, 100 ) ) # Ex 1.13 (Part ii): prob_alpha = function(c){ sum( dhyper( c:20, 20, 1000-20, 100 ) ) } prob_beta = function(c){ sum( dhyper( 0:c, 80, 1000-80, 100 ) ) } result = rbind( sapply( 0:8, prob_alpha ), sapply( 0:8, prob_beta ) ) rownames(result) = c("alpha","beta") colnames(result) = seq(0,8) print( result, digits=2 ) # Ex 1.14: root_fn = function(c){ n = 50 alpha = 1 - pnorm( (c-1+0.5 - n*(1/2)) / sqrt(n*(1/4)) ) beta = pnorm( (c-1+0.5 - n*(0.7)) / sqrt(n*0.7*0.3) ) alpha - beta } c = uniroot( root_fn, c(1,100) ) root_fn( floor(c$root) ) root_fn( ceiling(c$root) ) # Ex 2.1: # n = ( ( qnorm(0.05) * sqrt(0.21) - qnorm(0.95) * sqrt(0.25) )/0.2 )^2 c = 0.5 + 0.5 * n + qnorm(0.95) * sqrt( 0.25 * n ) # Ex 2.2: # p_A = 0.6 alpha = 0.05 beta = 0.05 # Ex 2.3: # alpha = 0.05 beta = 0.05 p_A = 0.15 # Ex 2.4: # find_smallest_n_binomial = function(n_max=100,p_H=0.3,p_A=0.8,alpha_max=0.1,beta_max=0.2){ # Start by increasing n looking for the first time we have a c # such that both the alpha value and the beta value are less than # the required threshold: # for( n in 2:n_max ){ cs = 1:n alphas = 1 - pbinom( cs-1, n, p_H ) betas = pbinom( cs-1, n, p_A ) inds = ( alphas < alpha_max ) & ( betas < beta_max ) # both conditions are met if( sum(inds) > 0 ){ print( n ) print( cs[inds] ) print( alphas[inds] ) print( betas[inds] ) return( list( n, cs[inds], alphas[inds], betas[inds] ) ) } } print( sprintf( "No smallest n in search to n_max= %5d. Try increasing n_max and running again.", n_max ) ) } find_smallest_n_binomial(); # Ex 2.5: # s_values = 13:20 P_A_D_eq_0 = choose( 18, s_values ) / choose( 20, s_values ) print( rbind( s_values, P_A_D_eq_0 ), digits=3 ) # Ex 2.6: # find_smallest_n_binomial(n_max=10000,p_H=0.514,p_A=0.55,alpha_max=1/100,beta_max=1/20) # Ex 2.7: # qnorm(0.95) # Ex 2.8: # 0.1 * 0.9 / 499 # Ex 2.9: # alpha_max = 0.1 beta_max = 0.2 coef = qnorm(1-alpha_max) * sqrt( 4 / 600 ) - qnorm(beta_max) * sqrt( 6 / 600 ) # Ex 2.10: # find_smallest_n_binomial(n_max=10000,p_H=1/3,p_A=9/25,alpha_max=0.05,beta_max=0.08) # Ex 3.2: # ps = seq( 0, 1, by=0.05 ) power = 1 - pbinom( 8, 10, ps ) plot( ps, power ) # Ex 3.3: # ps = seq( 0, 1, by=0.05 ) power = 1 - pbinom( 96, 120, ps ) plot( ps, power ) # Ex 3.4: # n = 100 c = 10 ps = seq( 0, 1, by=0.025 ) power = pnorm( ( c + 0.5 * n + 0.5 - n * ps ) / sqrt( n * ps * ( 1 - ps ) ) ) - pnorm( ( -c + 0.5 * n - 0.5 - n * ps ) / sqrt( n * ps * ( 1 - ps ) ) ) #postscript("../../WriteUp/Graphics/Chapter13/chap_13_sect_3_prob_4_plot.eps", onefile=FALSE, horizontal=FALSE) plot( ps, power, type='l' ) #dev.off() # Ex 3.5: # N = 1000 r = 90 s = 100 E_D = r * s / N var_D = ( ( N - s )/( N-1 ) ) * s * ( r / N ) * ( 1 - r / N ) 1 - pnorm( ( 4.5 - E_D ) / sqrt( var_D ) ) # Ex 3.9: # n_poems = 0:4 probs = c( 0.019, 0.08075, 0.2652, 0.633, 1 ) plot( n_poems, probs, type='l' ) # Ex 4.1: # cs = 1:5 pbinom( 5 - cs, 10, 1/2 ) + (1-pbinom( (5 + cs)-1, 10, 1/2 )) c = 4 # Plot the power curve : ps = seq( 0, 1, length.out=100 ) plot( ps, pbinom( 5 - c, 10, ps ) + (1-pbinom( (5 + c)-1, 10, ps )) ) # Ex 4.2: # cs = 1:200 rhs = pbinom( 200 - cs, 400, 1/2 ) + (1-pbinom( (200 + cs)-1, 400, 1/2 )) which.min( abs( rhs - 0.1 ) ) # gives 17 c = 17 pwr = pbinom( 200 - c, 400, 0.55 ) + (1-pbinom( (200 + c)-1, 400, 0.55 )) c = qbinom( 1 - 0.1, 400, 0.5 ) - 199 pwr = 1 - pbinom( 200+c-1, 400, 0.55 ) # Ex 4.3: # pbinom( (175:200)-1, 500, 1/3 ) # find when greater than 0.95 c = 185 1 - pbinom( c-1, 500, 1/3 ) # compute actual alpha 1 - pbinom( c-1, 500, 0.4 ) # compute power # Ex 4.4: # c( 1, 2, 4, 7, 12, 19, 30, 45, 66, 93, 129, 174, 232, 302, 388, 489 ) / 6435 # Ex 4.5: # c = qbinom( 0.1/5, 200, 1/4 ) - 1 5 * pbinom( c, 200, 1/4 ) # check value is less than 0.1 (yes) d_minus_one = qbinom( 1 - 4 * pbinom( c, 200, 1/4 ), 200, 1/4 ) # Ex 5.5: # #postscript("../../WriteUp/Graphics/Chapter13/chap_13_sect_5_prob_5_plot.eps", onefile=FALSE, horizontal=FALSE) Ns = 5:20 pi_1 = c( 1, 1/2, 1/3, 1/4, 1/5, ( Ns-4 ) / Ns ) plot( pi_1, type='b', col='blue', ylim=c(0,1), ylab="power function pi" ) pi_2 = c( 0, 1/2, 1/3, 1/4, 1/5, ( Ns-4 ) / Ns ) lines( pi_2, type='b', col='red' ) pi_3 = c( 0, 0, 1/3, 1/4, 1/5, ( Ns-4 ) / Ns ) lines( pi_3, type='b', col='green' ) pi_4 = c( 0, 0, 0, 1/4, 1/5, ( Ns-4 ) / Ns ) lines( pi_4, type='b', col='cyan' ) pi_5 = c( 0, 0, 0, 0, 1/5, ( Ns-4 ) / Ns ) lines( pi_5, type='b', col='black' ) #dev.off() Ns = 3:10 pi_s = c( 1, 1/choose(Ns,2) ) plot( pi_s ) # Ex 5.6: # p_h = c( 0.05, 0.05, 0.1, 0.8 ) p_a = c( 0.1, 0.2, 0.36, 0.34 ) lambda = p_h / p_a alphas = rev( cumsum( rev(p_h) ) ) pi_s = rev( cumsum( rev(p_a) ) ) # Ex 6.4: # mu_h = c(-2.5,4.01,0.352) sigmas = c(1.3,0.2,0.08) x_bars = c(-1.7,4.26,0.369) ns = c(12,4,20) zs = sqrt( ns ) * ( ( x_bars - mu_h ) / sigmas ) 1 - pnorm(zs) # Ex 6.6: # mu_h = c(0,2.13,0.5) sigmas = c(3,0.03,1) x_bars = c(3,2.12,0.7) ns = c(10,30,20) zs = sqrt( ns ) * ( ( x_bars - mu_h ) / sigmas ) pnorm(zs) # Ex 6.8: # z = 3.1622777 2 * pnorm(-z) # Ex 6.10: # alpha = 0.01 u = qnorm( 1 - alpha ) 1 - pnorm( u - ( sqrt(25) / 2 ) ) # Ex 6.11: # alpha = 0.02 u = qnorm( 1 - alpha ) 1 - pnorm( u - ( sqrt(25) * (1-0.5) / 1 ) ) # Ex 6.12: # alpha = 0.05 u = qnorm( 1 - alpha ) ( 3 * ( u - qnorm( 1-0.95 ) ) / (-0.6+0.8) )^2 # Ex 6.13: # v = -2.328 frac = seq( -5, +5, length.out=100 ) plot( frac, pnorm( v - frac ), main='pi(Delta)' ) # Ex 6.14: # alpha = 0.01 w = -qnorm( alpha / 2 ) n = 25 Delta = 2.5 - 2 sigma = 2 1 - pnorm( w - sqrt(n) * Delta / sigma ) + pnorm( -w - sqrt(n) * Delta / sigma ) # Ex 6.15: # alpha = 0.05 w = -qnorm( alpha / 2 ) n = 2000:3000 Delta = -0.6 - (-0.8) sigma = 3 # We want to find the value of n such that the power is 0.95 : root_fn = function(n){ 1 - pnorm( w - sqrt(n) * Delta / sigma ) + pnorm( -w - sqrt(n) * Delta / sigma ) - 0.95 } uniroot( root_fn, c(1000,5000) ) # Ex 6.17: # data = c( 6.59, 6.62, 6.88, 6.08, 6.26, 6.65, 6.36, 6.29, 6.58, 6.65, 6.55, 6.63, 6.39, 6.44, 6.34, 6.79, 6.20, 6.37, 6.39, 6.43, 6.63, 6.34, 6.46, 6.30, 6.75, 6.68, 6.85, 6.57, 6.49 ) x_bar = mean(data) s = sd(data) n = length(data) mu_h = 6.4 t = sqrt(n) * ( x_bar - mu_h ) / s 1 - pnorm( t - t * ( t^2 + 1 ) / ( 4 * n ) ) # Ex 6.19: # z = ( -1.5 - (-2) ) / sqrt( 4/20 + 1/20 ) 1 - pnorm(z) # Ex 6.22: # measurements = c( 18.83, 19.03, 18.61, 19.46, 18.80, 18.96, 19.37, 19.20, 18.88, 19.34 ) mu_gold = 19.3 x_bar = mean(measurements) S = sd(measurements) n = length(measurements) t = sqrt(n) * ( x_bar - mu_gold ) / S pnorm( t - t * (t^2 + 1)/(4*n) ) m_diff = measurements - mu_gold abs_diff = abs( m_diff ) so = order( abs_diff ) m_diff[so] # Ex 6.23: # pnorm(-3.984)