Post on 22-Feb-2018
7/24/2019 Examen de Estadstica Computacional
1/21
Examen de Estadstica Computacional
Problema 1.
Write R code to compare the performances of de Monte Carlo estimators of
1
2
e
x2
2
2dx
One choice of Monte Carlo estimator should sample from a uniform(1,2)
distribution, while another should sample from a standard normal distribution.
As additional choices, ou can consider importance samplin! ("n#esti!ate$)
usin! an alternati#e to the normal %ernel (remember that it should be hea# &
tailed(e'ample Cauch(1,2))). r two or three additional importance samplin!distributions (the could be within the same famil of distributions).
imulations si*es should be based on 1++ runs and 1++++ runs. ou should
include both !raphical and numerical comparisons.
Parte A. Estimacin usando la distribucin uniforme (1,2) y la N (0,1)
n -(1,2)ormal(+
,1)
1+++.1/0+
2+.0132+4
4
1++++
+.1033+5
1
+.24130/
0
C6di!o en R para simulaci6nsim781++++ 999para n:1++
1++++
u78runif(sim,1,2)
par(mar:c(2,2,2,1),mfrow:c(2,1))
f'78function(');(e'p(8
('rt(2?pi) @
restima78mean(f'(u))
sim781++++ 999para n:1++
1++++
u78rnorm(sim,+,1)
par(mar:c(2,2,2,1),mfrow:c(2,1))
f'78function(');e'p(8
('rt(2?pi)@
restima78mean(f'(u))
Parte B. Estimacin usando importance samplin con tres
distribuciones candidatas o en!ol!entes.
C6di!o en R para simulaci6nsim:1+
7/24/2019 Examen de Estadstica Computacional
2/21
p78e'p((8rt(cumsum((p8estint)
7/24/2019 Examen de Estadstica Computacional
3/21
Como puede obser#arse en cuanto se eli!e una distribuci6n >ue en#uel#a lo
7/24/2019 Examen de Estadstica Computacional
4/21
mDs pr6'imo a la funci6n >ue se tiene interEs en conocer el Drea baFo la
cur#a >ue representa, el muestreo por importancia estima un #alor mDs
cercano al #alor de la inte!ral e'presada >ue deseamos conocer. Ga
distribuci6n lo!normal(+,+.3) se acerca mDs a dicho #alor de +.103+3122,
estD Hltima es de +.1/3 para una tamaIo n:1++++.
7/24/2019 Examen de Estadstica Computacional
5/21
Problema 2. Write a pro!ram to !enerate random de#iates from a normal
distribution usin! the reFection method. -sin! the lo!istic function
f(x )= e
{x}
(1+e {x } )2,plot(sample.',data:,!eom:histo!ram,ll:accept,binwidth:+.+1))
colnames()78c(S#alorS,SdecS)
7/24/2019 Examen de Estadstica Computacional
6/21
r78as.matri'(subset(,dec::SesS))
proporcion78sum(nrow(r))=n
"r#$cos
7/24/2019 Examen de Estadstica Computacional
7/21
n cuanto a un histo!rama apilado de todos los #alores incluidos en la muestra
Funtos, realmente puede #er la cantidad de datos perdidos >ue ha en este
eFemplo. Te hecho, ha 2//4 aceptados 05332 recha*ados. Urobablemente
podrVa haber ele!ido un meFor #alor para la distribuci6n uniforme, pero lo ideal
del muestreo de recha*o es >ue utili*a una distribuci6n conocida >ue es s6lo unpoco diferente de la distribuci6n desconocida >ue estamos tratando de estimar.
Ga proporci6n de datos aceptados es +.2//4.
Problema %.
&ennedy abd "entle(1'0) studied the followin! al!orithm to !enerate
Uoisson de#iates usin! the well8%nown relationship between e'ponential wait8
in! times and Uoisson countsB
1.8 et G : e'p;8@, ':+,p:1
2.8 enerate u from -nif(+,1) and set p:pu.
0.8 "f p G, set ' : '1 and repeat step 2.Otherwise, output '.
'plain wh this al!orithm !enerates Uoisson de#iates. PintB "f Uois(), and
Q1,$ are iid 'p(1=), show that N=ni=1
n
Xi 1 and i=1
n+1
Xi>1 .
Write a second pro!ram to !enerate Uoisson de#iates usin! Xenned andentleYs approach in the abo#e problem. he pro!ram should allow the user to
specif the number of de#iates and the Uoisson mean. he output #ector
should contain the Uoisson de#iates.
dio en *
7/24/2019 Examen de Estadstica Computacional
8/21
poiss78 function(n,lambda)
;
'#78 matri'(+,nrow:n,ncol:1)
l:e'p(8lambda)
for (i in 1Bn);
':+
p:1
u78 runif(1,+,1)
p:p?u
while (pL:l)
;
':'1
u78 runif(1,+,1)
p:p?u
@
'#JiK:'
@par(mfrow:c(1,2))
78rpois(n,lambda)
hist('#)
hist()
return('#)
@
poi78 poiss(1++++,/)
l histo!rama '# es referido a nHmeros aleatorios simulados usando el
al!oritmo del Xenned and entle(14+) , comparando con el de la derecha
>ue es de los nHmeros aleatorios !enerados directamente de la funci6n en R,
rpois, muestran similitud entre ambos procedimientos. Go cual e'plica la
relaci6n e'istente entre la e'ponencial de los tiempos de espera los conteos
poisson. sto demuestra >ue en cuanto mDs cre*ca el #alor de n sur!e el hechode (n:).
Problema +.
uppose we want to simulate from the bi#ariate normal densit.
7/24/2019 Examen de Estadstica Computacional
9/21
(X , Y)'
N((00) ,( 1 0.90.9 1))
-sin! the !ibbs ampler Xy N(y ,(12)) , Yx N(x ,(1
2)) .
a) Zerif that the conditional probabilit of mo#in! from Q:'? to Q:' (a
two8step process) isB
x , x
xy 2
2(12)
yx2
2(12)
e
P
dio en *
!ibbs78function(n,rho); '78matri'(nrow:n,nco:2) ':+ :+ c:s>rt(18rho
7/24/2019 Examen de Estadstica Computacional
10/21
Ga ecuaci6n de arriba nos indica el %ernel de transici6n a usar para la
cadena Qi. Tonde Qies condicionalemente simulada solo en 'i-1 , ;Qi@ es una
cadena de Mar%o#. Ga distribuci6n in#ariante de Qies la distribuci6n mar!inal
fX(x )= f(x , y)dy . Como deseamos simular de una normal bi#ariada con
matri* de #arian*a8co#arian*as e'presada arriba, conocemos la condicionalesusando el ibbs ampler.
-na de las #entaFas importantes del al!oritmo de ibbs ampler es
>ue todas las simulaciones son aceptadas, en cada transici6n se obtiene un
punto diferente de la cadena. sto se debe a >ue la probabilidad de aceptaci6n
basado en el al!oritmo de Metropolis8Pastin!s es 1 en todo momento. AdemDs,
las simulaciones s6lo se reali*an a tra#Es de las densidades condicionales
totales. l hecho de >ue Estas sean densidades unidimensionales representa
una #entaFa computacional.
Como el al!oritmo estD basado en una cadena de Mar%o# , la teorVa menciona>ue dicha Cadena con %ernel de transici6n X satisface la condici6n de balance
si e'iste un funci6n f>ue cumpleB
K(,')f():K(',)f(')
Gue!o si se supone una cadema de Mar%o# con %ernel de transici6n Ksatisface
la condici6n de balance detallado como una funci6n de densidad de
probabilidad f. ntoncesB
Ga densidad fes una densidad in#ariante de la cadena.
Ga cadena es re#ersible.
A continuaci6n se muestra los !rDcos de la con#er!encia de la cadenaB
7/24/2019 Examen de Estadstica Computacional
11/21
7/24/2019 Examen de Estadstica Computacional
12/21
Con#er!encia de la cadena e histo!rama para la condicional de '
7/24/2019 Examen de Estadstica Computacional
13/21
Con#er!encia de la cadena e histo!rama para la condicional de
7/24/2019 Examen de Estadstica Computacional
14/21
b) how that this conditional distribution simplies so that
Xx N(2
x,(14
))
7/24/2019 Examen de Estadstica Computacional
15/21
Problema .
he Michaelis8Menten model is popular for !rowth8cur#e data, where is the
response, and Q, the co#ariate, is some measure of a!e or meB
Yi=
1Xi
2+Xi,i=1, , n
a) Nit the Michaelis8Menten model to the followin! data representin! the
len!th in mm() and a!e in ears(Q) of lar!emouth bassB
Gen!th78c(1+,1,2++,2+,210,213,21,20+,203,204,235,2/,23,252,243,24, 24,23,23,0+1,0+1,0+3,20/,01+,02,02,0/+,0/0,04+,04
1,042,04/, /++,/+1,0/1,053,040,/3,/13,/0,/0,/3,/5+,/5+,3+2,/3
1,/33,30+)A!e78c(rep(1,12),rep(2,1+),rep(0,12),rep(/,/),rep(3,/),,,,5,5,1+)
7/24/2019 Examen de Estadstica Computacional
16/21
NormulaB [ t1 ? Q=(t2 Q)
UarametersB stimate td. rror t #alue Ur(L\t\)t1 32.3++ 23.40 22.442 7 2e81 ???t2 1.1/1 +.21/5 4.1/ 1.04e811 ???888i!nif. codesB + ]???^ +.++1 ]??^ +.+1 ]?^ +.+3 ].^ +.1 ] ^ 1
Residual standard errorB 0/.25 on / de!rees of freedom
umber of iterations to con#er!enceB Achie#ed con#er!ence toleranceB 1.1+4e8+
b) -se bootstrap pro!ram to resample the residuals obtained from ttin!
the Michaelis8Menten model to the data, then construct a +_ pointwise
condence band for the nonlinear re!ression based on `:1+++
bootstrap samples. he !eneral approach that " used for the loess
bootstrap demonstration should wor%. he nlsfunction can be used to
nd a solution to the Michaelis8Menten model. `esides a formula, ou
should onl need to pro#ide start #alues for 1 and 2 the can be
selected b studin! as Q !oes to innit (for 1) , and then studin!
() at Q:1 for (2). ou can retrie#e the residuals from our model
usin! the resid function residuals are not a nls model attribute
otherwise.
Reali*ando el remuestreo de los residuales del modelo se obtiene unabanda de conan*a al +_B J81.245401 1.24/K, esto conFuntamente
con el histo!rama para la medias del remuestreo indican >ue sideseamos probar la hip6tesis de >ue los residuales aFustan a unadistribuci6n normal parece indicar >ue estD no se recha*a debido a >uedentro de la banda de conan*a inclue el cero obser#ando elhisto!rama la distribuci6n marcada dentro del mismo parece indicar >uela media es pr6'ima en cero.
7/24/2019 Examen de Estadstica Computacional
17/21
A continuaci6n se presenta el !rDco para las bandas de conan*a para cadapunto de los residuales del modelo usando bootstrapB
7/24/2019 Examen de Estadstica Computacional
18/21
dio en *78c(1+,1,2++,2+,210,213,21,20+,203,204,235,2/,23,252,243,24,
24,23,23,0+1,0+1,0+3,20/,01+,02,02,0/+,0/0,04+,041,042,04/,
/++,/+1,0/1,053,040,/3,/13,/0,/0,/3,/5+,/5+,3+2,/31,/33,30+)
Q78c(rep(1,12),rep(2,1+),rep(0,12),rep(/,/),rep(3,/),,,,5,5,1+)
n:len!th(Q)
modelo.nl78nls( [ t1?Q=(t2Q),start:list(t1:3,t2:2) )
summar(modelo.nl) 999obtenci6n de los residuales del modelo
residuales78resid(modelo.nl)
9bootstrap samplin! de los residuales
cboot78matri'(+,1,1+++)
boot.mean78function(residuales,nboot)
;
n78len!th(residuales)
cboot78numeric(nboot)
s78matri'(+,nrow:1+++,ncol:1)
for(i in 1Bnboot)
;
id78sample(1Bn,si*e:n,replace:)
cbootJiK78mean(residualesJidK)
sJiK78#ar(residualesJidK)
@
media.ord78sort(cboot)
s.ord78sort(s)
"C"78media.ordJ3+K8>t(+.+3,n81,lower.tail:,lo!.p:N)?s>rt(sJ3+K)=s>rt(n)
"C78media.ordJ3+K>t(+.3,n81,lower.tail:,lo!.p:N)?s>rt(sJ3+K)=s>rt(n)
9int78return(cbind("C","C,n))
dat78return(cbind(media.ord,s.ord))
@
test78boot.mean(residuales,1+++)
m78testJ,1K
med78mean(m)
#ar78testJ,2K#ar78#ar
hist(m,main:Medias usando el bootstrap)
9inter#alo de conan*a al +_
9 Calculo !racacion de bandas de conan*a para la media 9
"CA78 sort(m>t(+.+3, n81,lower.tail : , lo!.p : N)?s>rt(#ar)=s>rt(n)
"C`78 sort(m8>t(+.+3, n81,lower.tail : , lo!.p : N)?s>rt(#ar)=s>rt(n)
999RAN"CATO
plot("CA, tpe:SlS, lt:1, main:S`anda de conan*a al +_ para los residualesS, lab:SResidualS, col:0)
lines("C`, tpe:SlS,lt:1, col:2)
9lines(sort(residuales), tpe:SlS,lt:2, col:/)
lines(m, tpe:SlS,lt:2, col:/)
7/24/2019 Examen de Estadstica Computacional
19/21
Problema -.
uppose ou ha#e n data points Q1,$,Qn that are assumed to come from a
Uoisson() distribution. Carr out a Monte Carlo simulation stud to nd out
about the co#era!e of the nominal asmptotic 3_ condence inter#al (
1.96 s
n, +1.96 s /n ), where is the sample mean and sis the sample
standard de#iation and n is the sample si*e. (A nominal 3_ condence
inter#al is one that is supposedto ha#e 3_ co#era!e. "t ma not necessaril
ha#e 3_ co#era!e in realit which is wh ou are doin! this stud).
"n#esti!ate two escenariesB
a) :0, n:13.
b) :2+ n:1++.
"n each case pic% our Monte Carlo sample si*e Mlar!e enou!h so that our
Monte Carlo standard error for the co#era!e probabilities is small. ou ha#e to
do the followin!B
i) ubmit our R code on an!elii) Report our estimated co#era!e probabilit in each case and whether
ou would use the 3_ condence inter#al for this scenarioiii) "n each case report the M ou used an the Monte Carlo standard error
and condence inter#al for our estimate of the co#era!e probabilit.
Poisson(n,%)
nrrorestandar
GVmite"nferior
Gimitesuperior
1 +.//3442 1.1125/ /.32221%0 +.013513 1.444++0 /.13
0 +.2//0++ 1.4/+/5 0.+52100 +.150511 2.2+5/ 0.5+13
Poisson(n,20)
nrrorestDndar
GVmite"nferior
GVmitesuperior
1 1.13/1/5 13.03+/0 2/.0/25
%0 +.413/+14 1.40540 22.21/3
0 +.0120 15.//+0 22./313410 +.//5/05/ 14.1++3 21.52111
7/24/2019 Examen de Estadstica Computacional
20/21
0
dio en *
poisson78function(n,lambda)
;
9mu.+78+
9s.+78+
mu78matri'(+,nrow:1+++++,ncol:1)
s78matri'(+,nrow:1+++++,ncol:1)
for(i in 1B1+++++)
;
'78rpois(n,lambda)
muJiK78mean(')
sJiK78s>rt(#ar('))
@
mu.+78sort(mu)
s.+78sort(s)
M78s>rt(#ar(mu))
"C"78mu.+J23++K81.?sJ23++K=s>rt(n)"C78mu.+J53++K1.?sJ53++K=s>rt(n)
return(cbind("C","C,n,M))
@
poisson(13,0)
poisson(0+,0)
poisson(3+,0)
poisson(1++,0)
poisson(1++,2+)
Uuede notarse >ue al incrementar el tamaIo de muestra el error estDndar #adisminuendo. lo >ue se espera >ue la cobertura #aa aumentando.
7/24/2019 Examen de Estadstica Computacional
21/21
Gos histo!ramas muestran >ue a medida >ue el tamaIo de muestra crece se
asemeFa mas a la distribuci6n normal a >ue por el teorema central del lVmite,
para #alores !randes de , una #ariable aleatoria Uoisson Q puede apro'imarse
por otra normal dado >ue el cociente
Y=X
Con#er!e a una distribuci6n normal de media nula.
l histo!rama de la i*>uierda es para las medias de una Uoisson(13,0) a la
derecha para las medias de la Uoisson(1++,2+).
ademDs se si!ue cumpliendo >ue el #alor para la media la #arian*a es i!ual
a 2.403 1./1. Go cual indica >ue la simulaci6n por mEtodo MonteCarlo, cumple con una cobertura para funci6n de densidad mu alta.