In the previous STT5100 course, last week, we’ve seen how to use monte carlo simulations. The idea is that we do observe in statistics a sample \{y_1,\cdots,y_n\}, and more generally, in econometrics \{(y_1,\mathbf{x}_1),\cdots,(y_n,\mathbf{x}_n)\}. But let’s get back to statistics (without covariates) to illustrate. We assume that observations y_i are realizations of an underlying random variable Y_i. We assume that Y_i are i.id. random variables, with (unkown) distribution F_{\theta}. Consider here some estimator \widehat{\theta} – which is just a function of our sample \widehat{\theta}=h(y_1,\cdots,y_n). So \widehat{\theta} is a real-valued number like 
. Then, in mathematical statistics, in order to derive properties of the estimator \widehat{\theta}, like a confidence interval, we must define \widehat{\theta}=h(Y_1,\cdots,Y_n), so that now, \widehat{\theta} is a real-valued random variable. What is puzzling for students, is that we use the same notation, and I have to agree, that’s not very clever. So now, \widehat{\theta} is 
.
There are two strategies here. In classical statistics, we use probability theorem, to derive properties of \widehat{\theta} (the random variable) : at least the first two moments, but if possible the distribution. An alternative is to go for computational statistics. We have only one sample, \{y_1,\cdots,y_n\}, and that’s a pity. But maybe we can create another one \{y_1^{(1)},\cdots,y_n^{(1)}\}, as realizations of F_{\theta}, and another one \{y_1^{(2)},\cdots,y_n^{(2)}\}, anoter one \{y_1^{(3)},\cdots,y_n^{(3)}\}, etc. From those counterfactuals, we can now get a collection of estimators, \widehat{\theta}^{(1)},\widehat{\theta}^{(2)}, \widehat{\theta}^{(3)}, etc. Instead of using mathematical tricks to calculate \mathbb{E}(\widehat{\theta}), compute \frac{1}{k}\sum_{s=1}^k\widehat{\theta}^{(s)}That’s what we’ve seen last friday.
I did also mention briefly that looking at densities is lovely, but not very useful to assess goodness of fit, to test for normality, for instance. In this post, I just wanted to illustrate this point. And actually, creating counterfactuals can we a good way to see it. Consider here the height of male students,
Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") Davis[12,c(2,3)]=Davis[12,c(3,2)] X=Davis$height[Davis$sex=="M"] |
We can visualize its distribution (density and cumulative distribution)
u=seq(155,205,by=.5) par(mfrow=c(1,2)) hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) lines(u,dnorm(u,178,6.5),col="black") Xs=sort(X) n=length(X) p=(1:n)/(n+1) plot(Xs,p,type="s",col="blue") lines(u,pnorm(u,178,6.5),col="black") |
Since it looks like a normal distribution, we can add the density a Gaussian distribution on the left, and the cdf on the right. Why not test it properly. To be a little bit more specific, I do not want to test if it’s a Gaussian distribution, but if it’s a \mathcal{N}(178,6.5^2). In order to see if this distribution is relevant, one can use monte carlo simulations to create conterfactuals
hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) Y=rnorm(n,178,6.5) hist(Y,col=rgb(1,0,0,.3)) lines(density(Y),col="red",lwd=2) Ys=sort(Y) plot(Xs,p,type="s",col="white",lwd=2,axes=FALSE,xlab="",ylab="",xlim=c(155,205)) polygon(c(Xs,rev(Ys)),c(p,rev(p)),col="yellow",border=NA) lines(Xs,p,type="s",col="blue",lwd=2) lines(Ys,p,type="s",col="red",lwd=2) |
We can see on the left that it is hard to assess normality from the density (histogram and also kernel based density estimator). One can hardly think of a valid distance, between two densities. But if we look at graph on the right, we can compare the empirical distribution cumulative distribution \widehat{F} obtained from \{y_1,\cdots,y_n\} (the blue curve), and some conterfactual, \widehat{F}^{(s)} obtained from \{y_1^{(s)},\cdots,y_n^{(s)}\} generated from F_{\theta_0} – where \theta_0 is the value we want to test. As suggested above, we can compute the yellow area, as suggest in Cramer-von Mises test, or the Kolmogorov-Smirnov distance.
d=rep(NA,1e5) for(s in 1:1e5){ d[s]=ks.test(rnorm(n,178,6.5),"pnorm",178,6.5)$statistic } ds=density(d) plot(ds,xlab="",ylab="") dks=ks.test(X,"pnorm",178,6.5)$statistic id=which(ds$x>dks) polygon(c(ds$x[id],rev(ds$x[id])),c(ds$y[id],rep(0,length(id))),col=rgb(1,0,0,.4),border=NA) abline(v=dks,col="red") |
If we draw 10,000 counterfactual samples, we can visualize the distribution (here the density) of the distance used a test statistic \widehat{d}^{(1)}, \widehat{d}^{(2)}, etc, and compare it with the one observe on our sample \widehat{d}. The proportion of samples where the test-statistics exceeded the one observed
mean(d>dks) [1] 0.78248 |
is the computational version of the p-value
ks.test(X,"pnorm",178,6.5) One-sample Kolmogorov-Smirnov test data: X D = 0.068182, p-value = 0.8079 alternative hypothesis: two-sided |
I thought about all that a couple of days ago, since I got invited for a panel discussion on “coding”, and why “coding” helped me as professor. And this is precisely why I like coding : in statistics, either manipulate abstract objects, like random variables, or you actually use some lines of code to create counterfactuals, and generate fake samples, to quantify uncertainty. The later is interesting, because it helps to visualize complex quantifies. I do not claim that maths is useless, but coding is really nice, as a starting point, to understand what we talk about (which can be very usefull when there is a lot of confusion on notations).
. Then he flips a coin: tails he has a gain of 1, heads he experiences a loss of 1. At time 
, his wealth is 
where 
is associated to the 
(tails), and -1 with probability 
(heads). It is also possible to write
can be interpreted as the net gain of the player. In order to get a good understanding of results that can be obtained. Assume 
denote the number of heads and 
the number of tails. Then 
, while 
. Let 
denote the number of paths to go from point A (wealth 
at time 
) to point B (wealth 
at time 
denote the number of paths going from 
to 
. And let 
denote the number of paths going from 
denotes the first time 
become null, given 
,
, that we denote 
from now on. The term on the left can be approximated using monte-carlo simulations
, and 
. Then is comes that
for the number of runs, and then letting 
)
.
intégrable, alors pour tout 
et 
, 
. La conséquence de cette propriété (qui est juste) est que 
dès lors que la transformation n’est pas linéaire (ou au moins, il n’y a aucune raison pour qu’il y ait égalité).
(définies sur le même espace), 
, pour tout 
, et que l’on simule les couples 
, alors
sera grande. Donc effectivement, l’approximation pourrait être relativement mauvaise si on est sur des couples présentant une forte dépendance.
Didn’t like the ending, February 10, 2009 By Damien Katz
I found a typo, September 14, 2007 By fanfan
A serious reference work?, October 16, 2006 By BJ
Not Nearly A Million, September 3, 2006 By Liron
on obtient
, et on part de 0.
Et si on regarde l’histogramme et la fonction d’autocorrélation sur les 10 000 premières valeurs, on a gardé l’uniformité (ou tout du moins on ne doit pas en être trop loin), et on n’a plus d’autocorrélation.
Mais cet algorithme présente de très gros défauts numériques. Je n’aborderais pas ici les problèmes de ce générateur (je l’avais fait dans le chapitre sur les simulations dans le tome 2 du livre avec Michel Denuit).
. Par exemple si on prend comme valeur initiale π-3, on obtient, avec une simple fonction R
Plus généralement sur le caractère aléatoire des décimales de π, je renvoie aux publications de David H. Bailey (
Bref, en l’état cette méthode n’est pas vraiment géniale pour générer du “hasard“… à suivre donc.
sur [0,1].
Bref, on a à notre disposition des gens qui ont tiré des boules dans une urne. Alors rigoureusement, il faut faire attention car les tirages sont fait sans remise pour les 7 premières boules, puis on remet les boules pour le tirage suivant… Bref, les tirages ne sont pas rigoureusement indépendant les uns des autres (on ne peut pas avoir 3 fois consécutives le nombre “13”, alors qu’en théorie cela devrait être possible).