Untitled

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

#Exercice 0

#1.Simuler une dizaine de réalisations indépendantes de variables aléatoires de loi uniforme sur [0,1].

runif(10)

##  [1] 0.18311902 0.53963300 0.22735962 0.99112732 0.79571575 0.45036542
##  [7] 0.72153543 0.62566786 0.04908972 0.08636519

#2. Changer la graine (seed) et recommencer.

set.seed(1)

runif(10)

##  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968
##  [7] 0.94467527 0.66079779 0.62911404 0.06178627

set.seed(2)

runif(10)

##  [1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393 0.9434750 0.1291590
##  [8] 0.8334488 0.4680185 0.5499837

#On peut voir que mettre une certaine graine(seed) nous renvoie la même distribution uniforme chaque fois qu'on relance. Ça se passe parce que le générateur n'est pas vraiment aléatoire, mais pseudoaléatoire.

#3. Simuler alors un échantillon de taille entre 1000 et 10 000, choisie individuellement et donnez une représentation de la distribution empirique de votre échantillon.

x<-runif(10000)

hist(x)

#Nous pouvons voir que nous avons presque une distribution uniforme, avec chaque valeur ayant une fréquence similaire.

#Exercice 1

#1. Simulez Nmc réalisations indépendantes d'une variable aléatoire Z. Vous choisirez individuellement les valeurs de p1,p2,...,p7 ainsi que la valeur de Nmc.

nfamilles=7

  


jeu_de_familles<-function(ntrials){
  trials = rep(0,ntrials)
  for (j in 1:ntrials)
{
  trials[j] = sample(1:nfamilles,1)
}
  return(trials)
}


trials<-jeu_de_familles(1000)

#Choisissez individuellement un k E {1,...,7} donnez un estimateur pk du parametre pk construit à partir de vos simulations. Justifiez la convergence asymptotique de cet estimateur.

k=1

pk<-sum(trials==k)/length(trials)

#La vraie valeur pk doit être 1/7, alors que l'estimateur pk va se rapprocher de 1/7 en augmentant le n. Ça se vérifie avec la loi forte des grands nombres.

#3. Donnez une représentation de la distribution empirique de pk.

n<-c(1:1000)

pk_vec<-rep(0,1000)

for(i in 1:1000){
  trials<-jeu_de_familles(i)
  pk_vec[i]<-sum(trials==k)/length(trials)
}

plot(n,pk_vec)

mean(pk_vec)

## [1] 0.1437907

#Nous pouvons voir que l'estimateur de pk va se rapprocher de l'esperance de pk, égale à 1/7, avec un plus grand nombre de simulations choisies.

#Exercice 2 

#1. Simulez une variable aleatoire reelle de loi N(0,1) suivant la methode de Box-Muller puis la methode du rejet polaire (Marsaglia).

#La methode de Box Muller

box<-function(mean,variance,n)
{
    u1<-runif(n)
    r<--2*log(u1)
    u2<-runif(n)
    th<-2*pi*u2
    z1<-sqrt(r)*cos(th)
    z2<-sqrt(r)*sin(th)
    z<-vector(length=n*2)
    i=1
    j=1
    while (i<=n*2)
    {
        z[i]=z1[j]
        i=i+1
        z[i]=z2[j]
        j=j+1
        i=i+1
    }
    z<-z*variance+mean
    return(z)
}
z<-box(0,1,1000)
hist(z,main="box-mueller",breaks=50,freq=FALSE)

mean(z)

## [1] -0.02580264

sd(z)

## [1] 1.008534

#La methode Marsaglia
marsaglia<-function(mean,variance,n){
  a<-2*runif(n)-1
  b<-2*runif(n)-1

  

w<-a^2+b^2
r<-(-2*log(w)/w)

#z_1<-a*sqrt(-2*log(w)/w)
#z_2<-b*sqrt(-2*log(w)/w)

z<-vector(length=n)
    i=1

    while (i<=n)
    {
      if(r[i]>=0){
        z[i]=a[i]*sqrt(r[i])
        
      }
        
        
        i=i+1
    }
  z<-   z[z!=0]
    return(mean+variance*z)
}


z<-marsaglia(0,1,5000)
hist(z,breaks=50,freq=FALSE)

mean(z)

## [1] -0.0183878

sd(z)

## [1] 0.9916237

#2. Comparez les temps de simulation des deux methodes en fonction du nombre de gausiennes a simuler. Expliquer en quelques mots les raisons d'eventuelles differences.



system.time(z<-marsaglia(0,1,10000))

##    user  system elapsed 
##   0.003   0.000   0.002

system.time(z<-box(0,1,10000))

##    user  system elapsed 
##   0.004   0.000   0.003

#La methode de Marsaglia va avoir un temps d'execution moindre, car elle n'a pas besoin de calculer les cos et les sin.

#3


install.packages("fitdistrplus")

## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6'
## (as 'lib' is unspecified)

library(fitdistrplus)

## Loading required package: MASS

## Loading required package: survival

## Loading required package: npsurv

## Loading required package: lsei

plotdist(z,"norm",para=list(mean=5, sd=10))

#4. E[X] Monte Carlo. 

n<-rep(0,1000)

marsaglia_mean<-rep(0,1000)

for(i in 1:1000){
  marsaglia_mean[i]=mean(marsaglia(5,10,i))
}

plot(1:1000,marsaglia_mean,type="l")

#Exercice 3

#1.

inverse_exp<-function(u,mu){
    return(-log(1-u)/mu)
}

e<-inverse_exp(runif(25000),1)
z<-box(0,1,12500)

hist(e,freq=FALSE,breaks=500)

install.packages("Rlab")

## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6'
## (as 'lib' is unspecified)

library(Rlab)

## Rlab 2.15.1 attached.

## 
## Attaching package: 'Rlab'

## The following object is masked from 'package:survival':
## 
##     cancer

## The following object is masked from 'package:MASS':
## 
##     michelson

## The following objects are masked from 'package:stats':
## 
##     dexp, dgamma, dweibull, pexp, pgamma, pweibull, qexp, qgamma,
##     qweibull, rexp, rgamma, rweibull

## The following object is masked from 'package:datasets':
## 
##     precip

bernoulli<-rbern(25000, 0.2)


A=(z*((bernoulli-1)^2))+e*bernoulli

hist(A,breaks=500,freq=FALSE)

mean(A)

## [1] 0.2033087

#Exercice 4

N = 500
k_mean = 0
k_var = 0.5
h = 0.3
X = sample(x = c(0:9),size = N,replac=TRUE, prob = c(0.1,0.05,0.15,0.06,0.14,0.25,0.02,0.02,0.02,0.19))
#On choisit comme fonction la densité de la loi normale 
Kh<- function(x,h){
  return (dnorm(x/h,mean = k_mean,sd = k_var))/h
}
f<-function(x,h){
  sum(Kh(x-X,h))/N
}
empiric<-numeric(N)
for (i in 1:N){
  empiric[i]=f(10*(i-1)/N,h)
} 
hist(x = X,freq=FALSE)
lines(seq(from=0,to=9,length.out = N),empiric,col='red')

Untitled

R Markdown

Including Plots