Generate 100 Multinomial Random variable having n= 30, respective probabilities of (0.3, 0.2, 0.4, 0.1) using method from Textbook Chapter 4.7. Derive mean and variance of X1, and draw histogram. Diagnose the result if it has np1 = 9, np1(1-p1)= 0.63 which corresponds to theoretic conjectures.

prob<-c(0.3,0.2,0.4,0.1)
q<-2013122032
myfunc<-function(n, p, k, q){
  set.seed(q)
  pr<-(1-prob[k])^n
  F<- pr
  c<-prob[k]/(1-prob[k])
  i<-0
    repeat{
    U<-runif(1)
    if(U < F) break
    else{
    pr<-(c*(n-i)/(i+1))*pr
    F<-F+pr
    i<-i+1
    }
    X<-i
    }
return(X)
}
myfunc(30,0.3,1, 2013122032)

## [1] 8

myfunc(22,2/7,2, 2013122032)

## [1] 3

myfunc(19,0.5,3, 2013122032)

## [1] 7

this.is.the.function<-function(n){
e<-matrix(0,4,n)
q<-2013122032
for(i in 1:n){
q<-q+1
a<-myfunc(30,0.3,1, q)
b<-myfunc(30-a,2/7,2, q)
c<-myfunc(30-a-b,0.8,3, q)
d<-30-a-b-c
e[,i]<-c(a,b,c,d)
}
return(e)
}
this.is.the.function(100)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,]    9   10    8    8    7    5    9    7    3     7     8     9     9
## [2,]    2    3    3    4    3    5    4    4    3     3     5     3     3
## [3,]    9    3    6    4    7    5    4    7   10     7     8     8     4
## [4,]   10   14   13   14   13   15   13   12   14    13     9    10    14
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,]     9     4     6     9     9     6     8     7     7     8     7
## [2,]     3     3     5     4     3     3     5     5     5     2     4
## [3,]     7     4     6     8     7     6     7     7     5     8     7
## [4,]    11    19    13     9    11    15    10    11    13    12    12
##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,]     6     7     8     7     6     6     7     5    10     6     8
## [2,]     3     3     3     4     6     5     3     3     4     5     7
## [3,]     6     7     8     7     6     6     7     5     4     6     7
## [4,]    15    13    11    12    12    13    13    17    12    13     8
##      [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
## [1,]     6    10     9     7     6     9     6     7     8     8     7
## [2,]     3     3     4     3     3     6     5     5     2     4     5
## [3,]     6     5     7     7     6     6     6     6     8     6     7
## [4,]    15    12    10    13    15     9    13    12    12    12    11
##      [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57]
## [1,]     6     6     8     7     8     3    10     6     6     8     9
## [2,]     5     4     4     3     1     3     2     6     5     4     3
## [3,]     5     6     7     7     8    10     6     6     6     4     9
## [4,]    14    14    11    13    13    14    12    12    13    14     9
##      [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [1,]     6     7     9     8     8     6     8     7     5     6     7
## [2,]     3     1     4     5     5     5     4     4     1     3     5
## [3,]     6     7     4     5     5     5     5     7     5     6     7
## [4,]    15    15    13    12    12    14    13    12    19    15    11
##      [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79]
## [1,]     8     8     6     9     9     7     7    10     6    10     6
## [2,]     4     4     6     4     1     4     5     3     4     3     4
## [3,]     5     7     6     7     9     4     7     7     6     8     6
## [4,]    13    11    12    10    11    15    11    10    14     9    14
##      [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
## [1,]     7     6     7     9     7     5     7     7     7     9    10
## [2,]     5     3     4     5     2     5     3     3     2     3     1
## [3,]     6     6     6     5     7     5     6     7     7     7     7
## [4,]    12    15    13    11    14    15    14    13    14    11    12
##      [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## [1,]    10     4     9     8     8    11     9     6     6      4
## [2,]     4     4     2     4     5     3     4     2     4      4
## [3,]     5     4     9     7     8     5     7     6     6      4
## [4,]    11    18    10    11     9    11    10    16    14     18

First row indicates X1

a<-this.is.the.function(100)
a[1,1:10]

##  [1]  9 10  8  8  7  5  9  7  3  7

mean(a[1,])

## [1] 7.34

var(a[1,])

## [1] 2.630707

Create a function to generate sequence of random variables.

myfunc<-function(n, p){
  i<-0
    repeat{
    U<-runif(1)
    if(U < F) break
    else{
    pr<-(c*(n-i)/(i+1))*pr
    F<-F+pr
    i<-i+1
    }
    X<-i
    }
return(X)
}

2. Generate 200 random variable following pdf of Exp(0.2). Use methods from Exampe 5b. Draw histograms to compare

1) First we generate 200 exponential random varible using basic R function : rexp()

method1<-rexp(200, 0.2)

Error with : hist(method1, main=“Generating exponential random variable using basic R function”)

set.seed(2013122032)
unif.v<-runif(200,0,1)
X.rv<-log(1-unif.v)/(-0.2)

Error with : hist(X.rv, main=“Generating exponential random variable using inverse transformation method”)

par(mfrow=c(2,1))
hist(method1, main="Generating exponential random variable using basic R function")
hist(X.rv, main="Generating exponential random variable using inverse transformation method")

3. Generate 200 standard normal random variables using multiple methods; rejection method using exponential distribution, polar method and basic r function of rnorm(). Draw histograms for each of the results to diagnose their distributions.

A. Rejection method using exponential distribution

myfuncA<-function(n){
set.seed(2013122032)
pocket<-vector(length=30)
c<-sqrt(2*exp(1)/pi)
for (i in 1:n){
  repeat{
Y<-rexp(1,1)
U<-runif(1)
if(U<=exp(-(Y-1)^2/2)) break
}
  U<-runif(1)
  if(U>1/2) pocket[i]<-(-1)*Y else pocket[i]<-Y
}
return(pocket)
}
hist(myfuncA(200))

fn<-fivenum(myfuncA(200))
((fn[4]-fn[3])-(fn[3]-fn[2]))/((fn[4]-fn[3])+(fn[3]-fn[2]))

## [1] -0.008768848

sort.data<-sort(myfuncA(200))
q.125<-sort.data[25]
q.625<-sort.data[75]
kurtosis<-(q.625-q.125)/(fn[4]-fn[2])-1.705
kurtosis

## [1] -1.056089

B. Polar method

myfuncB<-function(n){
  set.seed(2013122032)
  X.vector<-c()
  Y.vector<-c()
  for(i in 1:n){
    Rsqr<-(-2*log(runif(1)))
    theta<-(2*pi*runif(1))
    X.vector[i]<-sqrt(Rsqr)*cos(theta)
    Y.vector[i]<-sqrt(Rsqr)*sin(theta)
  }
normal.rv<-c(X.vector, Y.vector)
return(normal.rv)
  }
hist(myfuncB(200))

fn<-fivenum(myfuncB(200))
((fn[4]-fn[3])-(fn[3]-fn[2]))/((fn[4]-fn[3])+(fn[3]-fn[2]))

## [1] 0.07524688

sort.data<-sort(myfuncB(200))
q.125<-sort.data[25]
q.625<-sort.data[75]
kurtosis<-(q.625-q.125)/(fn[4]-fn[2])-1.705
kurtosis

## [1] -1.175469

C. Using rnorm(200,0,1)

myfuncC<-function(n){
set.seed(2013122032)
rnorm(n,0,1)
}
myfuncC(200)

##   [1]  0.07795212 -0.52452241 -1.32165596  1.16614873 -1.10271469
##   [6] -0.97968308 -0.09317199 -0.08189487 -0.98270656 -0.46371171
##  [11] -0.55730870 -1.63390528  2.00617285 -0.32841597  0.00525315
##  [16]  1.91808160 -0.69127812 -1.33037185  1.50923898 -0.41746466
##  [21]  1.06306149 -0.20239130 -0.61009713  0.09601668 -0.29651834
##  [26]  0.92698219  0.59782695 -0.91463477 -1.62027297 -1.27918525
##  [31] -2.08031614  0.52870506  0.86789293 -1.15706861  1.29978153
##  [36]  0.70722466  1.75992443 -0.78665185 -1.01974308 -1.27569937
##  [41]  0.70965527 -0.69440695 -0.56346838  2.36968573  0.03235621
##  [46]  0.50300192 -0.07959173  1.34337061  0.39264119 -0.27559279
##  [51] -1.01218398 -0.96351026  0.74737771  0.78186053 -0.31989315
##  [56]  0.01359434  0.75577754 -0.11053702 -0.63467925  0.55735073
##  [61]  0.74787950  0.99722903 -0.94828725  1.96999052 -1.97973918
##  [66]  1.28714485 -0.18988412 -1.49492216  0.03067364 -1.87689611
##  [71]  0.67414661  1.40878445 -1.01104172 -1.14357890  0.03440163
##  [76] -0.57275227  1.04026604  1.92491882 -0.23704671  0.16918144
##  [81] -0.25118576  0.82453811  0.11530490 -0.33634047 -0.81818068
##  [86]  0.61249933  0.86483678 -0.97449908  0.10152273 -0.18136230
##  [91]  1.43219115  0.10376441 -0.93860223 -1.22544012  1.52238262
##  [96]  1.58296771  2.40752189 -0.43835338 -1.67467288  0.81790992
## [101]  0.37600210 -1.11229260 -1.30964311 -0.29312600  0.88798011
## [106] -0.20487181 -0.76220488 -1.09756313  1.53984457 -0.21946630
## [111] -0.12283314 -0.41950418 -1.65493903  1.55839433  0.91136729
## [116] -1.46363377  1.52897618 -1.09146609 -0.44306551 -2.15315023
## [121]  0.51619279  1.49893399 -0.71584246  1.77208159 -1.80191432
## [126] -0.14394629 -1.66717941  0.40994538 -0.08673844  1.92917139
## [131]  1.17149080  0.73387495  0.06137114  0.10410468  0.84157441
## [136] -0.39938728  1.89658537 -0.66495074  0.06280490  1.04824072
## [141]  0.25882294 -1.39301275 -0.68733310  1.27547688 -0.62712467
## [146] -0.16047208 -0.49254450 -0.61210499  0.35616835 -1.40951305
## [151]  1.35672672 -1.28364006 -0.24006383 -0.35218234 -0.37421434
## [156] -0.64126630  0.30215577 -0.18062520  1.19437775 -1.53658579
## [161]  1.17415410  0.74130413  0.91792560 -0.21092488  0.13070972
## [166] -1.10252369 -0.19892235 -0.04274412 -2.31008392  1.32559926
## [171]  0.84998396  0.37327917 -1.71329563  0.48346149 -0.58757907
## [176] -0.30309707  1.00342287  1.75259680  0.26864600  0.71381811
## [181]  2.75507843 -1.39568806  0.42230280  0.11179428 -0.16084027
## [186]  0.97648962  0.87727586  0.08122942 -0.49708558 -0.18457833
## [191] -0.87598820 -0.18748795 -1.23581530  2.05591354  1.28014885
## [196]  0.58569160 -0.96867825  0.96777037 -1.69228109 -0.69660209

hist(myfuncC(200))

fn<-fivenum(rnorm(200,0,1))
((fn[4]-fn[3])-(fn[3]-fn[2]))/((fn[4]-fn[3])+(fn[3]-fn[2]))

## [1] 0.1010391

set.seed(2013122032)
sort.data<-sort(rnorm(200,0,1))
q.125<-sort.data[25]
q.625<-sort.data[75]
kurtosis<-(q.625-q.125)/(fn[4]-fn[2])-1.705
kurtosis

## [1] -1.028629

Calculate duration of each of random variable methods used in question number 3.

system.time(myfuncA(20000))

##    user  system elapsed 
##    0.89    0.00    0.89

system.time(myfuncB(20000))

##    user  system elapsed 
##    1.48    0.03    1.51

system.time(myfuncC(200000))

##    user  system elapsed 
##    0.02    0.00    0.02

Statistical Computing_Chapter 4,5 Assignment

SeongJin Kim

April 1, 2018

2. Generate 200 random variable following pdf of Exp(0.2). Use methods from Exampe 5b. Draw histograms to compare

1) First we generate 200 exponential random varible using basic R function : rexp()

3. Generate 200 standard normal random variables using multiple methods; rejection method using exponential distribution, polar method and basic r function of rnorm(). Draw histograms for each of the results to diagnose their distributions.