About Me

My photo
Just wasting time sharing knowledge about, Big Data and Analytics

Jul 16, 2013

What are my chances to talk to this girl? Fisher or Bayes

Robert Mathews said that : "Ronald Fisher gave scientists a mathematical machine for turning baloney into breakthroughs, and ukes into funding. It is time to pull the plug.". He's right.
In one previous life, I wrote a thesis in Philosophy. But, a specific area, Epistemology also called
theory of knowledge, because, It questions what knowledge is and how it can be acquired,
and the extent to which any given subject or entity can be known.

My thesis deal about : The tradition, since Cournot in applying mathematical modelling to social sphere and more specific, how the climate modelling interact with interdiscplinary source of knowledge (mathematics, physics,geography, philosophy).

After reading this : http://www.academia.edu/1075253/Climate_Change_Epistemic_Trust_and_Expert_Trustworthiness 
It seems that the use of bayesian statistics is misunderstood.

Assume two thesis :
  • According to the Bayesian's statistic is the science who deal about the degree of of proofs in the observations. That means that, Bayesian statistics self-contained paradigm providing tools and techniques for all statistical problems.
  • In the classical frequentist view point of statistical theory, a statistical  procedure is judged by averaging its performance over all possible data
However, the bayesian approach gives prime importance to how a given procedure performs for
the actual data observed in a given situation.

The core of this theory have been formalised by Popper (Karl) with two main principles :
  1. Knowledge cannot start from nothing — from a tabula rasa – nor yet from observation. The advance of knowledge consists, mainly, in the modification of earlier knowledge. Although we may sometimes, for example in archaeology, advance through a chance observation,the significance of the discovery will usually depend upon its power to modify our earlier theories.
  2. Any probability is a degree of belief about something; It's not a property.That means that any scientific model produce data absolutely or conditionnaly on probability. It's possible to measure how the data can modify the degrees of belief (baye's principle).
One of the result of frequentist theory, may be the most criticized, is the fisher s p-value.


A problem can be to evaluate how the diploma is important to get in the first job

Hypothese :

- The diploma has no effect on the first job. In frequentist, we compute the p-value that can be
interpreted as the probability to  observe a difference at least as  important observed in the data if our hypothesis is true.

Read more: http://www.answers.com/topic/bayesian-statistics#ixzz2UKZDkq3v

So, let us talk about the p-value  and this problem :
If you cross seven times a beautiful girl each day for  10 days at the same place, can you conclude that she is always there? And then, tomorrow, youalways have a chance to speak with her.
We're going to examine both approach : Fisher (based on p_value) or Jeffreys (Bayesian)

  • With Fisher, we have :
H0 : p=0.5 and H1 :p=0.5
where p is "the probability to cross the girl". If we reject H0 (fisher) we can conclude that the girl is always there.
the p_value is :


> 2*pbinom(3,10,0.5)
0.34375
which is greater than 5% (the arbitrary threshold that everybody likes). indeed, we cannot conclude... Thank's Fisher; I can't trust you. 

With Bayes, we take our two hypothesis again and we suppose p(H0)=p(H1)=0.5. We have the same chance to cross her every day. With bayesian approach, wa need to compute the likelihood,  wich require, a priori distribution, that can be sum up (in our example) by "your a priori belief of cross again and to talk to her"


Let s compute this :
> (p <- c(0.5,0.6,0.7))
[1] 0.5 0.6 0.7
> (apriori <- c(0.5,0.3,0.2))
[1] 0.5 0.3 0.2
> (vraisemblance <- dbinom(7,10,p))
[1] 0.1171875 0.2149908 0.2668279
> (loi.jointe <- apriori*vraisemblance)
[1] 0.05859375 0.06449725 0.05336559
> (p.y <- sum(loi.jointe))
[1] 0.1764566
> (aposteriori <- loi.jointe/p.y)
[1] 0.3320576 0.3655134 0.3024290
 
So, the conditionnal probability p(H0/y)=is 33% and P(H1/y) = 66% .
There are two out of three chance that the girl is always there. In other words;
I have 66% of chances to cross her tomorrow.
Bayesian approach is hopeful. I like it; I can take my time with the girl next door.
 
 
 
 

Jul 10, 2013

Analyse discriminante linéaire ou Regression logistique

Supposons que l'on dispose d'iris de Paris (en population >100khabts) et qu'on veuille pouvoir les classer selon leurs caractéristiques sociodémos :
  • Population
  • taux de chômage
  • Etudiants
  • CSP
  • etc...
Une fois, les iris classés, on se demande si l'on peut transporter cette typologie à une autre grande ville (Lyon) par exemple : Il faudrait alors pouvoir utiliser un modèle d'affectation des iris selon leurs caractéristiques respectives à des classes prédéfinies. Le choix de la méthode dépend en général dans ce cas de ... pas grand chose comme indiqué ici. L'analyse discriminante ou une régression logistique, ou un SVM, pour ne citer que ces trois méthodes, ferait très bien l'affaire; Mais les résultats sont-ils identiques ou mieux, quelle est la méthode qui commet le moins d'erreur et pourquoi?

On va tester deux méthodes ici :
  • L'analyse discriminante 
  • La régression logistique multinomiale
Le programme :
# Etape 1 : Réduction de la dimension par une ACP pour diminuer le nombre de variables
# Etape 2 : CAH sur les facteurs de  l'ACP --> Définition des classes
# Etape 3: Analyse discriminante pour trouver les règles d'affectation aux classes
# Etape 4: Reg. log. pour trouver les mêmes règles d'affectation
# Etape 5: Comparaison des résultats des méthodes sur un échantillon test et mesure du taux de mal classés




1. Réduction de la dimension

# Etape 1 : Réduction de la dimension
require(ade4)
Vm.acp <-dudi.pca(df=Vm,center=TRUE,scale=TRUE,scan=FALSE)
# Eboulis des valeurs propres
screeplot(Vm.acp,type ="l", main = 'choix du nombre d axes')
nbaxes =3;
Vm.acp <-dudi.pca(df=Vm,center=T,scale=T,nf=nbaxes,scannf=FALSE)
s.corcircle(Vm.acp$co,xax=1,yax=3,clabel=0.7,sub="C. de corr des variables",possub="topright")
summary(Vm.acp)
New.vm<-Vm.acp$li
2. CAH sur les facteurs de l'ACP

# CAH sur les facteurs de l'ACP
preclus<-dist(New.vm)**2 #  La CAH op?re sur les distances
cah.vm<-hclust(d=preclus,method="ward")
# Visualisation de la CAH
par(mfrow=c(1,1))
nbclasses.fi =3
groups<-cutree(cah.vm,k=nbclasses.fi)
plot(cah.vm)
rect.hclust(cah.vm,k=nbclasses.fi,border="blue")
Les résultats de la classif

# Rattachement des groupes pour chacun des iris
> discrim.vm<-cbind(Vm,cluster=as.factor(groups))
> table(discrim.vm$cluster)

  1   2   3 
180 483 250 
 
3. Analyse discriminante pour affectation dans les classes
 
> index <- sample(nrow(discrim.vm), nrow(discrim.vm)*.70)
> #Echantillon d'apprentissage
> appren <- discrim.vm[index, ]
> #Echantillon de test
> test <- discrim.vm[-index, ]
> vm.disc= f.lda(appren[,-29],groups=appren$cluster)
> # Matrice de confusion
> pred.vm <-predict(vm.disc,newdata=test)
> # Taux  d'erreur
> Tx_err <- function(y,ypred){
+ mc <- table(y,ypred)
+ error <- 100*(mc[1,2]+mc[2,1])/sum(mc)
+ print(mc)
+ print(paste(round(error,2),"%",sep =""))
+ }
> Tx_err(test$cluster,pred.vm$class)
   ypred
y     1   2   3
  1  47   4   1
  2   3 142   2
  3   0  11  64
[1] "2.55%"
 
Une erreur de prévision de 2% avec une analyse discriminante
 
4. Avec une régression logistique multinomiale
 
# Avec une régression logistique multinomiale:
> # Avec une régression logistique multinomiale:
> library(nnet)
> logi=multinom(cluster~., data = appren)
# weights:  90 (58 variable)
initial  value 702.013252 
iter  10 value 113.648468
iter  20 value 75.506364
iter  30 value 69.704088
iter  40 value 65.256664
iter  50 value 62.358589
iter  60 value 61.152375
iter  70 value 60.246843
iter  80 value 59.528363
iter  90 value 58.957304
iter 100 value 58.452446
final  value 58.452446 
stopped after 100 iterations
> pp = predict(logi,newdata=test)
> Tx_err(test$cluster,pp)
   ypred
y     1   2   3
  1  51   0   1
  2   1 140   6
  3   1  12  62
[1] "0.36%"
 
 
Verdict : La reg. log. fait beaucoup mieux que l'analyse discriminante.
 
Certains auteurs  avaient déjà commencé à traiter la question. Ils voyaient l'analyse discriminante linéaire comme un cas particulier de la régression logistique (ce avec quoi je ne suis pas totalement d'accord, car aucune de ces deux méthodes ne s'affranchit de l'hypothèse forte de normalité et utilise les moments d'ordre deux à chaque fois)
Quelques conclusions néanmoins  intéressantes :

1. Lorsque le nombre de paramètre à estimer est important, le temps que met une reg log peut être 1,5 fois plus important que l'ADL du fait de l'algo itératif
2. La rég. log. peut être plus précise sur de petits échantillons, car du fait des modalités de référence, le nombre de paramètres à estimer est plus faible que si l'on utilise une ADL
 

GT
 

Jul 7, 2013

ggmap : Interesting toolbox for spatial analysis

ggmap is a new tool which enables such visualization by combining the spatial information of static maps from Google Maps, OpenStreetMap, Stamen Maps or CloudMade Maps with the layered grammar of graphics implementation of ggplot2

The library is developped by David Kahle and Hadley Wickham and in the latest R/Journal (Volume 5/1, June 2013), there is a whitepaper, very interesting that should to be read.

Let's use it to see where we can always buy tobacco at night in Paris.

Data are from data.ratp.fr, and give all authorized dealers by ratp. There's long and lat in this file, but this geocoding use lambert. So to show how ggmap work, we can extract adress, run geocoding using api/google and then ouput thematic map.

> require(ggmap)
> setwd("C:\\Users\\guibertt\\Desktop\\Geocodage - A faire")
> fic<-"commerces.csv"
> commerces<-read.csv(fic,header=T,sep=";", dec=",")
> head(commerces[,c(1,3:6)])
  DEA_CODE    ADRESSE_LIVRAISON DEA_CODE_POSTAL_LIVRAISON DEA_COMMUNE_LIVRAISON INSEE_LIVRAISON
1   510001       11   R. Mozart                     92230         Gennevilliers           92036
2   510002     81   Bd Voltaire                     92600    Asnières-sur-Seine           92004
3   510003 60   Av. Jean Moulin                     92390 Villeneuve-la-Garenne           92078
4   510004     2   Av. Michelet                     93400            Saint-Ouen           93070
5   510006      31   R. d'Anjou                     92600    Asnières-sur-Seine           92004
6   510007 45   R. Jules Larose                     92230         Gennevilliers           92036

We use the api google (limited at 2500 requests by day/ip adress) to geocode our places
>ad<-as.vector(tabac$adresse2)
>system.time(gc <- geocode(ad,output='latlona',messaging=FALSE))
> head(gc)
       lon      lat                                        address
1 2.267043 48.85066         108 avenue mozart, 75116 paris, france
2 2.297897 48.84552 56 rue de la croix nivert, 75015 paris, france
3 2.276529 48.84276             20 rue cauchy, 75015 paris, france
4 2.293242 48.83890  157 rue de la convention, 75015 paris, france
5 2.285362 48.83406       37 boulevard victor, 75015 paris, franc
>tabacp<-cbind(tabac,gc)
But, don't forget that ggmap is just a ggplot and we can get this, for example using OpenStreetMap
So we can now plot our places
png("paris4.png", width=800,height=600)
map <- get_map(location = 'paris',zoom=13,maptype="roadmap",color="color",source="google")
mymap = ggmap(map, darken = c(.3,'white'))
mymap+
  stat_bin2d(
aes(x = lon, y = lat, colour = typecommerce, fill = typecommerce),
size = .5, bins = 30, alpha = 1/2,
data = tabacp)
dev.off()



And at the end, with wrap, to split
png("paris5.png", width=800,height=600)
map <- get_map(location = 'paris',zoom=13,maptype="roadmap",color="bw",source="google")
mymap = ggmap(map, darken = c(.8,'white'))
mymap +  stat_bin2d(
aes(x = lon, y = lat, colour = TCO_LIBELLE, fill = TCO_LIBELLE),
size = .5, bins = 30, alpha = 1/2,
data = tabacp)+facet_wrap(~ TCO_LIBELLE)
dev.off()

Jun 9, 2013

How to read quickly large dataset in R?


Here, or there, I read many techniques to import a large dataset in R.
The option read.table or read.csv doesn't work anyway because, as discusshere, R load in memory. And sometimes, when we try to load a big dataset, we got this message :

Warning messages: 
1: Reached total allocation of 8056Mb: see help(memory.size)
2: Reached total allocation of 8056Mb: see help(memory.size) 

Many techniques can be used to load a large dataset. I found some there, or there. But there is two techniques that I never think before. 
Suppose that we have a large dataset with 10 millions rows

Comparing the methods for loading in R. 
- Using read.table

read.csv() performs a lot of analysis of the data it is reading, to determine the data types. So we can help R, by reading the first rows, determine the data type of the columns, and then, read the big data and provide the type of each columns and/or squeeze some of them if it doesn't need for analysis anyway;
Example
First we try to read a big data file (10 millions rows)
> system.time(df <-read.table(file="bigdf.csv",sep =",",dec=".")) Timing stopped at: 160.85 0.75 161.97 

 I let this run for a long period but no answer.

With this new method, we load the first rows, determine the data type and then, run read.table with indications of datatype.
> system.time (ds <- read.table("bigdf.csv", nrows=100, dec=".",sep=",")) user system elapsed 0 0 0 > classes <-sapply(ds, class) > classes V1 V2 V3 V4 "integer" "factor" "factor" "factor"

system.time(ds<-read.table("bigdf.csv",dec=".",sep=","colClasses=classes))
user  system elapsed 
234     432    128
As we see, this technique is not very interesting. It's also longer.
- We can use the package sqldf.

> require(sqldf)
> f <- file("bigdf.csv")
> system.time(SQLf <- sqldf("select * from f", dbname = tempfile(),
+                           file.format = list(header = T, row.names = F)))
Le chargement a nécessité le package : tcltk
   user  system elapsed 
  53.64    4.17   58.20 

Less of 1 minute  to import 10 millions rows  of an object of
> print(object.size(SQLf), units="Mb")
267 Mb

-          We can aslo used package read.table
> require(data.table)
Le chargement a nécessité le package : data.table
data.table 1.8.8  For help type: help("data.table")
> system.time(DT <- fread("bigdf.csv"))
   user  system elapsed 
 133.11    0.56  133.93 

But DT is a data.table format and a bit of transformation is require for use the table as dataframe using ddply from plyr package.

So. The point is : the package Sqldf is very useful to read quickly a large dataset in R. 10 millions rows in Less of a minute.


Jun 1, 2013

How logistic regression work ?

Discussing with a non statistician colleague, it seems that the logistic regression is not intuitive; Some basics questions like :
 - Why don't use the linear model?
 - What's logistic function?
 - How can we compute by hand, step by step to listen what is dealing by the glm function?

This post aims to answer that questions and may be this helps.

Suppose that we have this data : http://www.info.univ-angers.fr/~gh/wstat/pg.dar

 ID TAILLE GROUPE
1 A01    130      0
2 A02    140      0
3 C01    162      0
4 C02    160      1
5 A03    136      0
6 C03    165      1
 
and we want to predic the group according to the height. The problematic can be the level of risk according to age, or the customer segment according the amounts of transaction, etc. 
Let's remind.
When we compute a linear model (let's assume just one predictor : simple linear model), we have : E(y) =Cste + a1x1. Linear regression like all regressions focuses on the conditional probability distribution of Y given X.

The first think generally do is to draw the groupe = f(taille), we got :
The idea of Generaliszed Model (logistic regression is a particular ) is to replace E(Y) by something else.
For our example, we are interested by the probability of a person to be in group 0 or 1.
So, Instead of E(y) =Cste + a1x1, we seek P(Groupe==1) = a0 +a1*Taille. But, to solve the roblem, which is exactly the same to the other hand, we have to transform left hand side using a bijection between the interval[0,1]. That means to seek a "link" function that can help us to work in R.

The most useful function in logistic regression is : logit(p) = log(p/1-p). But one can also use the inverse of normal distribution(probit), the log-log distribution, or poisson distribution.

The method used to perform logistic regression is the maximization of likelihod estimator (MLE)
Read this post.
We sum up :
- Suppose in a population from which we are sampling, each individual has the same probability p to be in groupe 1 or groupe 0
- The likelihood is the joint probability of the data L = Product(P ** {Gourpe = 1} *(1 - p)**{Groupe = 0})
** Means power
For instance, we use log-likelihood.



How to interpret the likelihood :?

When we try to assign the group for a new id, it's natural to assign the group which have the best probability according to height.


Apply the MLE and perform logistic regression is done by


 



> Test = fit.logis(y=don$GROUPE,x=don$TAILLE)
> Test
  coef.est std.err
a  -27.190   8.885
b    0.181   0.058

 
 We can get the same output using glm function with "binomial" option.
>viaglm

Call:  glm(formula = don$GROUPE ~ don$TAILLE, family = "binomial", data = don)

Coefficients:
(Intercept)   don$TAILLE  
   -27.2103       0.1812  

Degrees of Freedom: 29 Total (i.e. Null);  28 Residual
Null Deviance:     38.19 
Residual Deviance: 10.89  AIC: 14.89


So, we can see that our optimisation via optim function is quite equivalent to glm function.
Just have a look

May be this helps to understand how it works !

May 18, 2013

Mining the last French presidential debate

After reading this post (thanks to him), I think it could be interesting to replicate this with some specific up of french language and to see and we can perform rapid view of the debate between Sarkozy and Hollande of the last 2nd round of presidential election.

Key words : TextMining, Elections, France, Debate, 2nd Round

We use the packages qdap from (Tyler Rinker) and tm to perform textmining analysis and the classical package like ggplot or RColorBrewer make our  graphics look pretty.
For Hollande






Top words From hollande
For Sarkozy
Top words from Sarkozy
As we can see, I've a problem to manage french accent. If somebody have any idea... We can also perform a quick Gantt plot basing on Qdap package and get some information about who lead the debate No surprise about the winner.

May 12, 2013

A new package : Quandl

Quandl is a new database management tool which seeks to become the place to find datasets. That is, each unique indicator is considered an independent data set. This helps them to seem to have a ginormous quantity of data sets. Source : Blog Econometric Simulation.
To load or find the datasets, we have to authentify using the API like with Twitter First, we need to set an account to receive the Quandl.auth





setwd("D:/PERSO/R_Working/Tutoriels/Probability")
library(Quandl)
Quandl.auth("XXXXXXXXXX")  ### Replace with yours
For example, we load the database about pollution and gdp and try to find the link.

plot(date, pollution, col = "red", type = "o", lwd = 2, ylim = c(70, 150), ylab = "", 
    main = "Evolution du PIB vs Pollution entre 94 et 2006 au Japon")
lines(date, pib, lty = 2, col = "purple", type = "o")
legend("topright", legend = c("PIB", "Pollution"), col = c("purple", "red"), 
    pch = 15, bty = "n", pt.cex = 2, cex = 0.8, text.col = "black", horiz = TRUE, 
    inset = c(0.1, 0.1))


plot of chunk unnamed-chunk-3
We fit a model to know if the growth of GDP is responsible of pollution's growth in japan between 1994 and 2006.

plot of chunk unnamed-chunk-4

Ecologists are not totaly wrong ! ! !

The entire code to run this post
Share it ! ! !