Entrega 1 (2016)

Trabajo Español
Universidad Universidad Autónoma de Barcelona (UAB)
Grado Estadística Aplicada - 3º curso
Asignatura Simulació, Remostreig i Aplicacions
Año del apunte 2016
Páginas 16
Fecha de subida 27/04/2016
Descargas 8
Subido por

Vista previa del texto

Primer trabajo evaluable Anna Olmo, 1363013 March 28, 2016 1 Ejercicio 1 Considera un empresario incompetente. Su compa˜ n´ıa se inicia con un capital de 10000 euros, pero tiene perdidas, en media, cada d´ıa. Concretamente, la p´erdida el d´ıa i-´esimo se puede describir seg´ un una variable aleatoria Yi ∼ N (−30, 10000). Si la compa˜ nia dispone de 11000 euros de capital en el banco, este empresario podr´ıa vender su compa˜ n´ıa a un competidor. Sin embargo, si la cuenta corriente del banco se sit´ ua en n´ umeros negativos, ´el se declarar´ a en bancarota.
1.1 a.
Crea una rutina en R que sea capaz de generar, usando alg´ un m´etodo cl´asico de simulaci´on, la distribuci´ on de Yi . Genera una muestra de Yi de tama˜ no N = 104 .
Para el m´etodo cl´ asico es necesario encontrar tribuci´ on Cauchy.
f (x) g(x) e = 1 √ 10000 2π π(1+x2 ) = π(1+x2 )e 40000 √ 10000 2π √ df (x) dg(x) = −(x−30)2 pi/2 40000 10000 [2xe + (1 + x2 )e donde f (x) es la distribuci´on normal y g(x) es la dis- √ −(x−30)2 −1 (x−30)2 2 2∗10000 f (x) g(x) −(x−30)2 40000 π/2 10000 (1 −(x−30)2 + x2 )e 40000 √ pi/2 −(x−30)2 −(x−30) x−30 ] = 40000[2x − (1 + x2 ) 20000 ] 20000 10000 e = df (x) dg(x) = 2x − (1 + x2 ) (x−30) 20000 = 0 2 40000x = (1 + x )(x − 30) 40000x = x3 − 30x2 + x − 30 x3 − 30x2 − 39999x − 30 = 0 El discriminante de la ecuaci´ on es negativo, por tanto no se puede resolver. As´ı que usar´e R para la programaci´ on.
> Y<-rnorm(1e4,mean=-30,sd=sqrt(10000)) > Y[1:12] [1] -104.63830 [7] -38.07838 1.2 427.89395 18.31115 24.56211 -112.13161 -103.57867 -122.31147 167.41446 -59.83138 135.27197 9.81169 b.
Escribe un programa en R para estimar la probabilidad de que este empresario venda su compa˜ n´ıa antes de declararse en bancarota. Calcula tambi´en el error est´andar del estimador. Considera N = 104 y usa el simulador construido en a para contestar el ejercicio (no uses rnorm).
> empresa<-function(C0,N){ + Y<-rnorm(N,mean=-30,sd=sqrt(10000)) + C<-C0+Y[1] + estado<-0 + for(i in 2:N){ + C[i]<-C[i-1]+Y[i] + estado[i]<-0 1 + if(C[i]>=11000){ estado[i] <- 1; break } + if(C[i]<0){ estado[i] <- -1; break } + } + estado[length(estado)] + } > final<-replicate(1e4,empresa(10000,1e4)) La probabilidad de vender antes de banca rota es: [1] 0.0012 La probabilidad de banca rota antes de vender es: [1] 0.9988 2 Ejercicio 2 Considera dos jugadores (A y B) jugando a un juego en el que A gana con probabilidad q en (0,1) y B gana con probabilidad 1 − q. El jugador A empieza con eA euros y el B con eB euros. El jugador que gana le da un euro a su oponente. Ellos repiten el juego de manera infinita. Sup´on, ademas, que un jugador puede tener capital negativo. Sea Xnn>0 un proceso estoc´astico donde Xn denota el capital en euros del jugador A despu´es de n juegos.
2.1 a.
Escribe un programa en R para simular los primeros N valores de Xnn>0 para q = 0.6, eA = 6 y eB = 4.
Dibuja una realizaci´ on para cada N = 10, 100, 10000.
Si la binomial da 1, la partida la gana el jugador A, por tanto ´este recibe una moneda de manos del jugador B. Si la binomial da 0, ocurrir´ a lo contrario. Sabiendo esto: • N=10 > > > + + + + + + + + + eA<-6; eB<-4; q<-0.6; n<-10 partidas<-rbinom(n,size=1,prob=0.6) for(i in 2:(n+1)){ if(partidas[i-1]==1){ eA[i]<-eA[i-1]+1 eB[i]<-eB[i-1]-1 } if(partidas[i-1]==0){ eA[i]<-eA[i-1]-1 eB[i]<-eB[i-1]+1 } } 2 15 5 −5 0 Monedas 10 Jugador A Jugador B 0 2 4 6 8 Partidas • N=100 > > > + + + + + + + + + eA<-6; eB<-4; q<-0.6; n<-100 partidas<-rbinom(n,size=1,prob=q) for(i in 2:(n+1)){ if(partidas[i-1]==1){ eA[i]<-eA[i-1]+1 eB[i]<-eB[i-1]-1 } if(partidas[i-1]==0){ eA[i]<-eA[i-1]-1 eB[i]<-eB[i-1]+1 } } 3 10 −10 0 Monedas 10 20 Jugador A Jugador B 0 20 40 60 Partidas • N=10000 > > > + + + + + + + + + eA<-6; eB<-4; q<-0.6; n<-10000 partidas<-rbinom(n,size=1,prob=q) for(i in 2:(n+1)){ if(partidas[i-1]==1){ eA[i]<-eA[i-1]+1 eB[i]<-eB[i-1]-1 } if(partidas[i-1]==0){ eA[i]<-eA[i-1]-1 eB[i]<-eB[i-1]+1 } } 4 80 100 2000 −2000 −1000 0 Monedas 1000 Jugador A Jugador B 0 2000 4000 6000 8000 10000 Partidas 2.2 b.
Sea τ = inf n > 0 : Xn = 0, eA + eB el n´ umero aleatorio de juegos despu´es de que el juego termine, si un capital negativo no es posible. Escribe un programa en R para estimar P (τ = 8) basado en 1000 realizaciones de Xnn>0 (usa los par´ ametros de a y N = 1000).
Es la misma funci´ on que antes, a˜ nadiendo que si uno de los jugadores se queda con 0 monedas terminan de jugar (break).
> SoloPositivo<-function(n){ + eA<-6; eB<-4; q<-0.6 + partidas<-rbinom(n,size=1,prob=q) + for(i in 2:(n+1)){ + if(partidas[i-1]==1){ + eA[i]<-eA[i-1]+1 + eB[i]<-eB[i-1]-1 + } + if(partidas[i-1]==0){ 5 + + + + + + + > > > eA[i]<-eA[i-1]-1 eB[i]<-eB[i-1]+1 } if(eA[i]==0||eB[i]==0){break} } length(eA) } partidas.jugadas<-replicate(10000,SoloPositivo(1000)) partidas.no.jugadas<-1000-partidas.jugadas summary(partidas.no.jugadas) Min. 1st Qu.
861.0 977.0 Median 987.0 Mean 3rd Qu.
982.5 993.0 Max.
995.0 Si no es posible capital negativo, segun estas muestras, solo pueden jugar aproximadamente de 5 a 180 partidas. Como deber´ıan llegar a las 1000 partidas, no jugaran de 820 a 995 partidas. Por tanto, la probabilidad de que en 1000 partidas u ´nicamente se dejen por jugar 8, es nula.
2.3 c.
Considera ahora que eA = eB = 5 y q = 0.5. Adem´as, si uno de los jugadores se arruina, el otro le cede un euro y contin´ uan con el juego. Sea entonces Xnn>0 el proceso estoc´astico donde Xn es el capital en euros del jugador A despu´es de n juegos.
Escribe un programa en R para simular 1000 muestras de Xnn>0 y dibuja un histograma de Xn0 para n0 = 5, 50, 100, 1000.
Es la funci´ on inicial con el a˜ nadido de que si unos de los jugadores se queda sin monedas, el otro le da una de las suyas.
• F´ ormula general: > Prestando<-function(n){ + eA<-5; eB<-5; q<-0.5 + partidas<-rbinom(n,size=1,prob=q) + for(i in 2:(n+1)){ + if(partidas[i-1]==1){ + eA[i]<-eA[i-1]+1 + eB[i]<-eB[i-1]-1 + } + if(partidas[i-1]==0){ + eA[i]<-eA[i-1]-1 + eB[i]<-eB[i-1]+1 + } + if(eA[i]==0){ + eA[i]<-eA[i-1]+1 + eB[i]<-eB[i-1]-1 + } + if(eB[i]==0){ + eA[i]<-eA[i-1]-1 6 + eB[i]<-eB[i-1]+1 + } + } + eA[n] + } • N=5 > eA.5muestras<-replicate(1000,Prestando(5)) > hist(eA.5muestras,freq=F,main="Dinero final de A para n=5") 0.2 0.0 0.1 Density 0.3 Dinero final de A para n=5 2 4 6 8 eA.5muestras • N=50 > eA.50muestras<-replicate(1000,Prestando(50)) > hist(eA.50muestras,freq=F,main="Dinero final de A para n=50") 7 0.3 0.0 0.1 0.2 Density 0.4 0.5 Dinero final de A para n=50 2 3 4 5 6 7 8 eA.50muestras • N=100 > eA.100muestras<-replicate(1000,Prestando(100)) > hist(eA.100muestras,freq=F,main="Dinero final de A para n=100") 8 0.3 0.0 0.1 0.2 Density 0.4 0.5 Dinero final de A para n=100 2 3 4 5 6 7 8 eA.100muestras • N=1000 > eA.1000muestras<-replicate(1000,Prestando(1000)) > hist(eA.1000muestras,freq=F,main="Dinero final de A para n=1000") 9 0.3 0.0 0.1 0.2 Density 0.4 0.5 Dinero final de A para n=1000 2 3 4 5 6 7 8 eA.1000muestras 2.4 d.
Escribe un programa en R para estimar pˆt = P (X1 000 = i) para i = 0, ..., 10 considerando N = 5000.
Tras hacer la simulaci´ on, simplemente hay que sumar todas las veces que el jugador A ha tenido cierta cantidad de monedas, y luego dividirlo entre el total de juegos simulados.
> > > > + + + + + + eA.5000muestras<-replicate(5000,Prestando(1000)) probs<-rep(0,11) j<-1 while(j<=length(eA.5000muestras)){ for(i in 0:10){ if(eA.5000muestras[j]==i){ probs[i+1]<-probs[i+1]+1 } } j<-j+1 10 + } > probs<-probs/5000 i P(X1000=i) [1,] 0 0.000e+00 [2,] 1 0.000e+00 [3,] 2 4.904e-05 [4,] 3 0.000e+00 [5,] 4 4.932e-05 [6,] 5 0.000e+00 [7,] 6 5.004e-05 [8,] 7 0.000e+00 [9,] 8 5.160e-05 [10,] 9 0.000e+00 [11,] 10 0.000e+00 3 Ejercicio 3 Considera una casa con cuatro habitaciones P1 , P2 , P3 y P4 . Hay puertas entre las habitaciones P1 y P2 , P2 y P3 , P3 y P4 y entre P1 y P4 . Asume que hay un gato en la habitaci´on P1 y un rat´on en la habitaci´ on P3 . Los dos van cambiando independientemente de habitaci´on, cada uno con probabilidad 0.9. En este caso el gato y el rat´ on se mueven hacia las habitaciones vecinas, con la misma probabilidad. Este proceso puede modelarse con una cadena de Markov.
Se mueven con probabilidad 0.9, por tanto, con una binomial deciden si moverse o quedarse donde est´ an.
Luego, desde cada habitaci´ on se puede ir a dos m´as (por ejemplo, de P2 se puede ir a P1 o P3 ), as´ı que con otra binomial de probabilidad 0.5 se decide a cual ir (si da 1 ir´a a la siguiente, si da -1 ir´a a la anterior).
Como de la siguiente de la 4 es la 1 i viceversa, se tiene que poner a mano dentro del if.
Esto se hace tanto para el gato como para el rat´on (por separado ya que son independientes). Por u ´ltimo, se mira en qu´e habitaci´ on se encuentra cada uno, y si est´an en la misma el gato se come al rat´on y termina el juego.
El ejercicio pide encontrar la esperanza de vida y la probabilidad de terminar en una cierta habitaci´ on, asi que la funci´ on dar´ a esos dos valores (cu´ anto han tardado en encontrarse y d´onde ha sido).
> vidaRaton<-function(n){ + pasosGato<-1 + pasosRaton<-3 + for(i in 2:n){ + if(rbinom(1,size=1,prob=0.9)==1){ + siguiente<-2*rbinom(1,size=1,prob=0.5)-1 + pasosGato[i]<-pasosGato[i-1]+siguiente + if(pasosGato[i]==0){pasosGato[i]<-4} + if(pasosGato[i]==5){pasosGato[i]<-1} + } + else{ + pasosGato[i]<-pasosGato[i-1] + } + if(rbinom(1,size=1,prob=0.9)==1){ + siguiente<-2*rbinom(1,size=1,prob=0.5)-1 + pasosRaton[i]<-pasosRaton[i-1]+siguiente 11 + if(pasosRaton[i]==0){pasosRaton[i]<-4} + if(pasosRaton[i]==5){pasosRaton[i]<-1} + } + else{ + pasosRaton[i]<-pasosRaton[i-1] + } + if(pasosGato[length(pasosGato)]==pasosRaton[length(pasosRaton)]){ + vida<-length(pasosRaton) + habitacion<-pasosRaton[vida] + break + } + } + c(vida,habitacion) + } La funci´ on da una matriz donde la primera fila contiene la vida del rat´on, y la segunda fila contiene la habitaci´ on en la que se encontraron.
3.1 a.
Si el gato y el rat´ on est´ an en la misma habitaci´on, el gato se come al rat´on. Sin embargo, el gato no se comer´ a al rat´ on si ellos se encuentran cuando est´an cambiando de habitaci´on. Escribe un programa de R para estimar el tiempo esperado de vida del rat´on. Simula la cadena de Markov 10000 veces.
La esperanza de vida del rat´ on ser´ a la media de las vidas simuladas.
> vida<-replicate(10000,vidaRaton(10000)) > mean(vida[1,]) [1] 5.0163 3.2 b.
Escribe un programa de R para estimar la probabilidad de que el gato se coma al rat´on en la habitaci´ on P3 bas´ andote en 10000 realizaciones de la cadena de Markov.
La probabilidad de que se encuentren en la 3a habitaci´on, es la suma de todas las veces que ha ocurrido en la simulaci´ on entre el total de simulaciones realizadas.
> habitacion<-vida[2,] > sum(habitacion[habitacion==3])/10000 [1] 0.5298 4 Ejercicio 4 Considera una poblaci´ on que inicialmente tiene 5 individuos. Los nuevos individuos crecen con tasa 0.5 y mueren con tasa 0.25. Adem´ as, con una tasa de 0.01 una pandemia mata a muchos de los individuos de la 12 poblaci´ on. En este caso, cada individuo muere independientemente del otro con probabilidad 0.5. Sea Xtt>0 el proceso que describe el n´ umero de individuos en la poblaci´on a tiempo t.
4.1 a.
Programa una rutina en R que sea capaz de simular la distribuci´on de los tiempos entre nacimientos, muertes y ocurrencia de pandemias.
Para resolver el ejercicio sup´ on que el n´ umero de nacimientos sigue un proceso de Poisson(λ1 ), el n´ umero de muertes sigue otro proceso de Poisson(λ2 ) independiente del primero y que el n´ umero de pand´emias sigue tambi´en otro proceso de Poisson(λ3 ) independiente a los dos primeros.
Dada una poblaci´ on, pueden pasar 3 cosas: nacen, mueren o pandemia. Lo que afecta primero es la pandemia, as´ı que con una poisson se mira si ocurre o no. En caso de que suceda, cada individuo tiene una probabilidad de 0.5 de morir. La poblaci´ on superviviente ser´a la inicial menos los muertos.
Para los supervivientes hay 2 opciones: nacer o morir. Cada individuo tiene 0.25 de probabilidad de morir, y 0.5 de tener descendencia. Por lo tanto, al final de t, la poblaci´on ser´a los que hayan sobrevivido a la pandemia y no hayan muerto de forma natural, m´as la descendencia.
> poblacionT<-function(n){ + poblacion<-5 + for(i in 2:(n+1)){ + pandemia<-rpois(1,0.01) + if(pandemia==1){ + muerte<-rpois(1,0.5) + poblacion[i]<-poblacion[i-1]-muerte + muertos<-rpois(1,0.25) + nacidos<-rpois(1,0.5) + poblacion[i]<-poblacion[i]-muertos+nacidos + } + else{ + muertos<-rpois(1,0.25) + nacidos<-rpois(1,0.5) + poblacion[i]<-poblacion[i-1]-muertos+nacidos + } + } + poblacion + } 4.2 b.
Escribe un programa en R para simular Xtt>0 para [0, 1000]. Representa una de las muestras simuladas y utiliza el simulador programado en a.
> simu<-poblacionT(1000) > plot(1:1001,simu,type="l",xlab="tiempo",ylab="individuos") 13 250 200 150 100 0 50 individuos 0 200 400 600 800 1000 tiempo 4.3 c.
Escribe un programa en R para estimar el tiempo esperado hasta que la poblaci´on tenga m´as de 50 individuos.
Usa N = 10000 y el simulador programado en a, si se precisa.
La funci´ on que se necesita es la misma que en a, pero esta vez el output ser´a el momento de t en que la poblaci´ on supera los 50 individuos.
> poblacion50<-function(n){ + poblacion<-5 + tiempo<-0 + for(i in 2:(n+1)){ + pandemia<-rpois(1,0.01) + if(pandemia==1){ + muerte<-rpois(1,0.5) + poblacion[i]<-poblacion[i-1]-muerte + muertos<-rpois(1,0.25) + nacidos<-rpois(1,0.5) 14 + poblacion[i]<-poblacion[i]-muertos+nacidos + } + else{ + muertos<-rpois(1,0.25) + nacidos<-rpois(1,0.5) + poblacion[i]<-poblacion[i-1]-muertos+nacidos + } + if(poblacion[i]>50){ + tiempo<-i + break + } + } + tiempo + } > T50<-replicate(10000,poblacion50(1000)) > mean(T50) [1] 189.9382 El tiempo medio para que la poblaci´ on tenga m´as de 50 individuos es aproximadamente de 190.
15 ...

Tags: