stat_data_analysis_rcode

Statistical Data Analysis - Fall 2018
Log | Files | Refs | README

R-IC04.R (3764B)

      1 #*******************************************************#
      2 #*******************************************************#
      3 # 		         STAT320 - IC04 Solutions	      	  	    #
      4 #*******************************************************#
      5 #*******************************************************#
      6 #***************************************************#
      7 #          CI for Proportions                    #  
      8 #***************************************************#
      9 
     10 #***************************************************#
     11 #***************************************************#
     12 #***************************************************#
     13 #***************************************************#
     14 #***************************************************#
     15 
     16 #***************************************************#
     17 #   Fair Coin
     18 #***************************************************#
     19 set.seed(257)
     20 #***************************************************#
     21 
     22 #***************************************************#
     23 #***************************************************#
     24 
     25 prop.CI.f = function(a=.05,p0=0.5,NSize=100,NSim=100000)
     26 {
     27   POP = c("T","H")
     28   Succ = "T"
     29   ONSuccs = numeric(NSim)
     30   
     31   ONSuccs = replicate(NSim,
     32             {
     33               Sample=sample(x=POP, prob =c(p0,1-p0),size = NSize,replace = TRUE) 
     34               sum(Sample==Succ)
     35             })
     36    x = quantile(ONSuccs,probs = c(a/2,1-a/2))/NSize
     37   return(x) 
     38 }
     39 
     40 
     41 #***************************************************#
     42 #***************************************************#
     43 NSim = 100000 # Number of simulated samples
     44 a = 0.05 # Alpha level for CI
     45 
     46 
     47 # Sample Estimate
     48 N = 564    # Sample size
     49 nS = 51    # Number of Successes in sample 
     50 phat = nS/N # Estimated sample proportion
     51 
     52 # Set True value to be the estimate from sample
     53 pT = phat
     54 
     55 # simulate 100,000 samples, compute their sample proportions, 
     56 # and store them in the vector SimProps
     57 
     58 POP = c(1,0)
     59 SimProps = numeric(NSim)
     60   
     61 for(i in 1:NSim)
     62 {
     63   # Generate a sample using the given pT and the sample size N
     64   Sample=sample(x=POP, prob =c(pT,1-pT),size = N,replace = TRUE) 
     65 
     66   # Count the number of successes (1's) in the simulated sample
     67   xS = sum(Sample==1)
     68                     
     69   # compute sample proportion of current sample 
     70   phat = xS/N
     71   
     72   # save current sample proportion
     73   SimProps[i] = phat
     74   
     75   # go back and repeat
     76 }
     77 
     78 
     79 # Compute the (1-a)*100% confidence based on the simulated phat's
     80 
     81 CI = quantile(SimProps,probs = c(a/2,1-a/2))
     82 
     83 CI 
     84 
     85 
     86 prop.CI.f(a,p0=phat,NSize=N,NSim=NSim)
     87 
     88 #***************************************************#
     89 
     90 
     91 #***************************************************#
     92 #***************************************************#
     93 # Program to demonstrate what confidence means 
     94 # in Confidence intervals
     95 #***************************************************#
     96 #***************************************************#
     97 
     98 N.CIs = 20   # Number of simulated samples
     99 N = 564      # Sample size
    100 
    101 pT = .15      # True percantage of successes in population
    102 a = 0.01
    103 
    104 NCIIn = 0     # Counter of the samples the produced a CI that
    105               # Included the true value
    106 
    107 POP = c(1,0)
    108 for(i in 1:N.CIs)
    109 {
    110   # Randomly select a sample from the Real Population to get 
    111   # an estmate phat of the true pT 
    112     CurSample = sample(x=POP, prob =c(pT,1-pT),size = N,replace = TRUE)
    113     phat = sum(CurSample)/N
    114     
    115   # we use the phat to compute a confidence interval 
    116     ci = prop.CI.f(a,p0=phat,NSize = N,NSim=NSim)
    117     
    118   # Checking if the particular CI contains the true value pT 
    119     if(ci[1] <= pT && pT <= ci[2])
    120       NCIIn = NCIIn + 1
    121 }
    122                   
    123 # Number and Proportion of CIs who actually contained the true value of pT
    124 NCIIn 
    125 NCIIn/N.CIs
    126 
    127