Examen Parcial (2016)

Examen 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 9
Fecha de subida 27/04/2016
Descargas 5
Subido por

Vista previa del texto

Parcial Anna Olmo, 1363013 April 27, 2016 1 Ejercicio 1 Uso el m´etodo inverso para simular una Pareto sin usar una Pareto.
Pareto: F (x) = 1 − ( x5 )1.05 C´ alculo de la inversa: x = 1 − ( y5 )1.05 ( y5 )1.05 = 1 − x 1 5 y = (1 − x) 1.05 5 y= 1 (1−x) 1.05 Para simular: 5 F (x) = 1 (1−u) 1.05 > > > > + + + + + + + > n<-1000 X<- 5/( (1 - runif(n))^(1/1.05) ) Y<-X for(i in 1:n){ if(X[i]<50){ pd<-max(1-(X[i]/50),0) if(runif(1)<=pd){ Y[i]<-0 } } } Y<-Y[Y>0] 1.1 (a) Ser´ an todas aquellas X que sean menores a 10, entre el total de X.
> length(X[X<10])/n [1] 0.504 1.2 (b) Los cuantiles los calcula la funci´ on quantile. El m´aximo es el cuantil 100.
> quantile(X,probs=c(0.50,0.90,0.95,0.99,1)) 50% 9.892432 1.3 90% 50.532266 95% 95.305920 99% 100% 465.853072 4994.523209 (c) Ser´ an todas aquellas Y que sean mayores a 65, entre el total de Y.
1 > length(Y[Y>65])/length(Y) [1] 0.2341772 1.4 (d) La esperanza se aproxima por la media.
> mean(Y) [1] 100.7185 1.5 (e) Los cuantiles los calcula la funci´ on quantile. El m´aximo es el cuantil 100.
> quantile(Y,probs=c(0.50,0.90,0.95,0.99,1)) 50% 29.62722 1.6 90% 145.07875 95% 99% 100% 270.49000 1489.65256 4994.52321 (f ) > hist(X,breaks=50,col="blue", + main="Histograma de los meteoritos", + xlab="Di´ ametro de los meteoritos", + xlim=c(0,3500),ylim=c(0,1000)) > hist(Y,breaks=50,col="light blue",add=T) > legend("topright",legend=c("Entran","Caen"), + fill=c("blue","light blue")) 2 1000 Histograma de los meteoritos 600 400 0 200 Frequency 800 Entran Caen 0 500 1000 1500 2000 2500 3000 3500 Diámetro de los meteoritos 2 > > > > > > Ejercicio 2 datos<-read.csv2("http://mat.uab.cat/~acabana/data/profundidad.csv") X1<-c(100,75,65,40,73,65,50,30,45,50,45,55) X2<-c(85,45,80,28,75,70,65,55,50,40) X3<-c(96,58,95,90,65,80,85,95,82) X4<-c(93,120,65,105,115,82,99,87,100,90,78,95,93,88,110) mean(X1);mean(X2);mean(X3);mean(X4) [1] 57.75 [1] 59.3 [1] 82.88889 [1] 94.66667 > Z<-c(X1,X2,X3,X4) # datos conjuntos > names(Z)<-c(rep("X1",12),rep("X2",10),rep("X3",9),rep("X4",15)) # a que X pertenece cada dato > Z<-sort(Z) # ordenar los datos 3 > > + + + + + + + + + + + + + > > > R<-rep(0,4) # calculo de las medias de las posiciones for(i in 1:length(Z)){ if(names(Z)[i]=="X1"){ # posicion de los que vienen de R[1]<-R[1]+i/12 } if(names(Z)[i]=="X2"){ # posicion de los que vienen de R[2]<-R[2]+i/10 } if(names(Z)[i]=="X3"){ # posicion de los que vienen de R[3]<-R[3]+i/9 } if(names(Z)[i]=="X4"){ # posicion de los que vienen de R[4]<-R[4]+i/15 } } names(R)<-c("X1","X2","X3","X4") estadT<-3*R[4]+R[3]-3*R[1]-R[2] # estad´ ıstico del test names(estadT)<-NULL 2.1 X1 X2 X3 X4 (a) Porque las X no tienen el mismo n´ umero de datos.
2.2 (b) Porque la hip´ otesi nula es que no hay diferencia entre las excavaciones y, por tanto, no importa de donde proviene cada valor, es decir, se pueden permutar.
2.3 (c) > estadT [1] 80.97222 2.4 (d) > permZ<-function(Z){ + nuevaZ<-sample(Z) + + R<-rep(0,4) + for(i in 1:length(Z)){ + if(names(nuevaZ)[i]=="X1"){ + R[1]<-R[1]+i/12 + } + if(names(nuevaZ)[i]=="X2"){ + R[2]<-R[2]+i/10 + } + if(names(nuevaZ)[i]=="X3"){ + R[3]<-R[3]+i/9 + } + if(names(nuevaZ)[i]=="X4"){ 4 + + + + + + + + + > > + > > > + R[4]<-R[4]+i/15 } } names(R)<-c("X1","X2","X3","X4") estadT<-3*R[4]+R[3]-3*R[1]-R[2] names(estadT)<-NULL estadT } estadTperm<-replicate(1000,permZ(Z)) hist(estadTperm,freq=F,xlim=c(-60,90),col="grey", main="Histograma de la distribuci´ on del estad´ ıstico T") abline(v=quantile(estadTperm,probs=0.95),col="red") points(estadT,0,col="blue") legend("topleft",legend=c("punto cr´ ıtico","test original"), fill=c("red","blue")) 2.5 (e) El pvalor es la probabilidad de que salga ”peor” de lo que ha salido, es decir, que los valores de las permutaciones sean mayores que el valor original. Por tanto, ser´a la cantidad de veces que eso pasa entre el total de permutaciones.
> pvalor<-sum(estadTperm>=estadT)/1000;pvalor [1] 0 2.6 (f ) Tanto en el gr` afico como en el pvalor se ve que no podemos aceptar la hip´otesis nula ya que hay pruebas estad´ısticamente significativas de que las medias de los yacimientos no son iguales.
2.7 (g) Para la muestra original: > > > > > > > estad1 <-R[1]-R[2] estad1[2]<-R[1]-R[3] estad1[3]<-R[1]-R[4] estad1[4]<-R[2]-R[3] estad1[5]<-R[2]-R[4] estad1[6]<-R[3]-R[4] names(estad1)<-NULL Con permutaciones: > permZsep<-function(Z){ + nuevaZ<-sample(Z) + + R<-rep(0,4) 5 + + + + + + + + + + + + + + + + + + + + + + + + + + > > > + + + + + + > for(i in 1:length(Z)){ if(names(nuevaZ)[i]=="X1"){ R[1]<-R[1]+i/12 } if(names(nuevaZ)[i]=="X2"){ R[2]<-R[2]+i/10 } if(names(nuevaZ)[i]=="X3"){ R[3]<-R[3]+i/9 } if(names(nuevaZ)[i]=="X4"){ R[4]<-R[4]+i/15 } } names(R)<-c("X1","X2","X3","X4") estad1<-R[1]-R[2] estad1[2]<-R[1]-R[3] estad1[3]<-R[1]-R[4] estad1[4]<-R[2]-R[3] estad1[5]<-R[2]-R[4] estad1[6]<-R[3]-R[4] names(estad1)<-NULL estad1 } estadperm<-t(replicate(1000,permZsep(Z))) par(mfrow=c(2,3)) for(i in 1:6){ hist(estadperm[,i],xlim=c(-25,25),col="grey", main="Histograma igualdad") points(estad1[i],0,col="blue") abline(v=quantile(estadperm[,i],probs=c(0.025,0.975)), col="red") } layout(1) 6 Frequency 10 20 50 0 0 0 150 250 300 200 50 100 Frequency 200 100 0 50 −20 0 10 20 −20 0 10 20 estadperm[, i] estadperm[, i] Histograma igualdad Histograma igualdad Histograma igualdad −20 0 10 20 estadperm[, i] 300 200 0 50 100 Frequency 200 0 0 50 50 100 150 Frequency 200 250 300 estadperm[, i] 100 Frequency −20 Frequency Histograma igualdad 350 Histograma igualdad 300 Histograma igualdad −20 0 10 20 −20 estadperm[, i] 0 10 20 estadperm[, i] Seg´ un se ve en los gr´ aficos, los u ´nicos que tienen medias iguales son X1 = X2 y X3 = X4.
3 3.1 Ejercicio 3 (a) > inventari<-function(N,lambda,mu,a){ + act<-0 + for(i in 2:1000){ + if(act==0){ # si la producci´ on no est´ a activada + if(N[i-1]<a){ # si el stock es inferior a a + act<-1 # se activa la producci´ on + prod<-rpois(1,mu) + venta<-rpois(1,lambda) + N[i]<-N[i-1]-venta+prod + } + if(N[i-1]>=a){ # si el stock es superiora a a + venta<-rpois(1,lambda) + N[i]<-N[i-1]-venta 7 + } + } + if(act==1){ # si la producci´ on est´ a activada + if(N[i-1]>=10){ # si el stock es superior a 10 + act<-0 # se detiene la producci´ on + venta<-rpois(1,lambda) + N[i]<-N[i-1]-venta + } + if(N[i-1]<10){ # si el stock es inferior a 10 + prod<-rpois(1,mu) + venta<-rpois(1,lambda) + N[i]<-N[i-1]-venta+prod + } + } + if(N[i]<=0){ + t<-i + break + } + } + t + } > tiempo<-inventari(10,0.9,1,5); tiempo [1] 9 3.2 (b) > tempsrep<-replicate(1000,inventari(10,0.9,1,5)) > mean(tempsrep) [1] 44.321 8 ...