Tag Archives: algorithm

Classification from scratch, SVM 7/8

Seventh post of our series on classification from scratch. The latest one was on the neural nets, and today, we will discuss SVM, support vector machines.

A formal introduction

Here y takes values in \{-1,+1\}. Our model will be m(\mathbf{x})=\text{sign}[\mathbf{\omega}^T\mathbf{x}+b] Thus, the space is divided by a (linear) border\Delta:\lbrace\mathbf{x}\in\mathbb{R}^p:\mathbf{\omega}^T\mathbf{x}+b=0\rbrace

The distance from point \mathbf{x}_i to \Delta is d(\mathbf{x}_i,\Delta)=\frac{\mathbf{\omega}^T\mathbf{x}_i+b}{\|\mathbf{\omega}\|}If the space is linearly separable, the problem is ill posed (there is an infinite number of solutions). So consider

The strategy is to maximize the margin. One can prove that we want to solve \max_{\mathbf{\omega},m}\left\lbrace\frac{m}{\|\mathbf{\omega}\|}\right\rbrace
subject to y_i\cdot(\mathbf{\omega}^T\mathbf{x}_i)=m, \forall i=1,\cdots,n. Again, the problem is ill posed (non identifiable), and we can consider m=1: \max_{\mathbf{\omega}}\left\lbrace\frac{1}{\|\mathbf{\omega}\|}\right\rbrace
subject to y_i\cdot(\mathbf{\omega}^T\mathbf{x}_i)=1, \forall i=1,\cdots,n. The optimization objective can be written\min_{\mathbf{\omega}}\left\lbrace\|\mathbf{\omega}\|^2\right\rbrace

The primal problem

In the separable case, consider the following primal problem,\min_{\mathbf{w}\in\mathbb{R}^d,b\in\mathbb{R}}\left\lbrace\frac{1}{2}\|\mathbf{\omega}\|^2\right\rbracesubject to y_i\cdot (\mathbf{\omega}^T\mathbf{x}_i+b)\geq 1, \forall i=1,\cdots,n.

In the non-separable case, introduce slack (error) variables \mathbf{\xi} : if y_i\cdot (\mathbf{\omega}^T\mathbf{x}_i+b)\geq 1, there is no error \xi_i=0.

Let C denote the cost of misclassification. The optimization problem becomes\min_{\mathbf{w}\in\mathbb{R}^d,b\in\mathbb{R},{\color{red}{\mathbf{\xi}}}\in\mathbb{R}^n}\left\lbrace\frac{1}{2}\|\mathbf{\omega}\|^2 + C\sum_{i=1}^n\xi_i\right\rbracesubject to y_i\cdot (\mathbf{\omega}^T\mathbf{x}_i+b)\geq 1-{\color{red}{\xi_i}}, with {\color{red}{\xi_i}}\geq 0, \forall i=1,\cdots,n.

Let us try to code this optimization problem. The dataset is here

n = length(myocarde[,"PRONO"])
myocarde0 = myocarde
myocarde0$PRONO = myocarde$PRONO*2-1
C = .5

and we have to set a value for the cost C. In the (linearly) constrained optimization function in R, we need to provide the objective function f(\mathbf{\theta}) and the gradient \nabla f(\mathbf{\theta}).

f = function(param){
  w  = param[1:7]
  b  = param[8]
  xi = param[8+1:nrow(myocarde)]
  .5*sum(w^2) + C*sum(xi)}
grad_f = function(param){
  w  = param[1:7]
  b  = param[8]
  xi = param[8+1:nrow(myocarde)]

and (linear) constraints are written as \mathbf{U}\mathbf{\theta}-\mathbf{c}\geq \mathbf{0}

U = rbind(cbind(myocarde0[,"PRONO"]*as.matrix(myocarde[,1:7]),diag(n),myocarde0[,"PRONO"]),
C = c(rep(1,n),rep(0,n))

Then we use

constrOptim(theta=p_init, f, grad_f, ui = U,ci = C)

Observe that something is missing here: we need a starting point for the algorithm, \mathbf{\theta}_0. Unfortunately, I could not think of a simple technique to get a valid starting point (that satisfies those linear constraints).

Let us try something else. Because those functions are quite simple: either linear or quadratic. Actually, one can recognize in the separable case, but also in the non-separable case, a classic quadratic program\min_{\mathbf{z}\in\mathbb{R}^d}\left\lbrace\frac{1}{2}\mathbf{z}^T\mathbf{D}\mathbf{z}-\mathbf{d}\mathbf{z}\right\rbracesubject to \mathbf{A}\mathbf{z}\geq\mathbf{b}.

eps = 5e-4
y = myocarde[,"PRONO"]*2-1
X = as.matrix(cbind(1,myocarde[,1:7]))
n = length(y)
D = diag(n+7+1)
diag(D)[8+0:n] = 0 
d = matrix(c(rep(0,7),0,rep(C,n)), nrow=n+7+1)
A = Ui
b = Ci
sol = solve.QP(D+eps*diag(n+7+1), d, t(A), b, meq=1, factorized=FALSE)
qpsol = sol$solution
(omega = qpsol[1:7])
[1] -0.106642005446 -0.002026198103 -0.022513312261 -0.018958578746 -0.023105767847 -0.018958578746 -1.080638988521
(b     = qpsol[n+7+1])
[1] 997.6289927

Given an observation \mathbf{x}, the prediction is

y_pred = 2*((as.matrix(myocarde0[,1:7])%*%omega+b)>0)-1

Observe that here, we do have a classifier, depending if the point lies on the left or on the right (above or below, etc) the separating line (or hyperplane). We do not have a probability, because there is no probabilistic model here. So far.

The dual problem

The Lagrangian of the separable problem could be written introducing Lagrange multipliers \mathbf{\alpha}\in\mathbb{R}^n, \mathbf{\alpha}\geq \mathbf{0} as\mathcal{L}(\mathbf{\omega},b,\mathbf{\alpha})=\frac{1}{2}\|\mathbf{\omega}\|^2-\sum_{i=1}^n \alpha_i\big(y_i(\mathbf{\omega}^T\mathbf{x}_i+b)-1\big)Somehow, \alpha_i represents the influence of the observation (y_i,\mathbf{x}_i).

Consider the Dual Problem, with \mathbf{G}=[G_{ij}] and G_{ij}=y_iy_j\mathbf{x}_j^T\mathbf{x}_i
subject to \mathbf{y}^T\mathbf{\alpha}=\mathbf{0} and \mathbf{\alpha}\geq\mathbf{0}.

The Lagrangian of the non-separable problem could be written introducing Lagrange multipliers \mathbf{\alpha},{\color{red}{\mathbf{\beta}}}\in\mathbb{R}^n, \mathbf{\alpha},{\color{red}{\mathbf{\beta}}}\geq \mathbf{0}, and define the Lagrangian \mathcal{L}(\mathbf{\omega},b,{\color{red}{\mathbf{\xi}}},\mathbf{\alpha},{\color{red}{\mathbf{\beta}}}) as\frac{1}{2}\|\mathbf{\omega}\|^2+{\color{blue}{C}}\sum_{i=1}^n{\color{red}{\xi_i}}-\sum_{i=1}^n \alpha_i\big(y_i(\mathbf{\omega}^T\mathbf{x}_i+b)-1+{\color{red}{\xi_i}}\big)-\sum_{i=1}^n{\color{red}{\beta_i}}{\color{red}{\xi_i}}
Somehow, \alpha_i represents the influence of the observation (y_i,\mathbf{x}_i).

The Dual Problem become with \mathbf{G}=[G_{ij}] and G_{ij}=y_iy_j\mathbf{x}_j^T\mathbf{x}_i\min_{\mathbf{\alpha}\in\mathbb{R}^n}\left\lbrace\frac{1}{2}\mathbf{\alpha}^T\mathbf{G}\mathbf{\alpha}-\mathbf{1}^T\mathbf{\alpha}\right\rbrace
subject to \mathbf{y}^T\mathbf{\alpha}=\mathbf{0}, \mathbf{\alpha}\geq\mathbf{0} and \mathbf{\alpha}\leq {\color{blue}{C}}.
As previsouly, one can also use quadratic programming

eps = 5e-4
y = myocarde[,"PRONO"]*2-1
X = as.matrix(cbind(1,myocarde[,1:7]))
n = length(y)
Q = sapply(1:n, function(i) y[i]*t(X)[,i])
D = t(Q)%*%Q
d = matrix(1, nrow=n)
A = rbind(y,diag(n),-diag(n))
C = .5
b = c(0,rep(0,n),rep(-C,n))
sol = solve.QP(D+eps*diag(n), d, t(A), b, meq=1, factorized=FALSE)
qpsol = sol$solution

The two problems are connected in the sense that for all \mathbf{x}\mathbf{\omega}^T\mathbf{x}+b = \sum_{i=1}^n \alpha_i y_i (\mathbf{x}^T\mathbf{x}_i)+b

To recover the solution of the primal problem,\mathbf{\omega}=\sum_{i=1}^n \alpha_iy_i \mathbf{x}_ithus

omega = apply(qpsol*y*X,2,sum)
                           1                        FRCAR                        INCAR                        INSYS 
 0.0000000000000002439074265  0.0550138658687635215271960 -0.0920163239049630876653652  0.3609571899422952534486342 
                       PRDIA                        PAPUL                        PVENT                        REPUL 
-0.1094017965288692356695677 -0.0485213403643276475207813 -0.0660058643191372279579454  0.0010093656567606212794835

while b=y-\mathbf{\omega}^T\mathbf{x} (but actually, one can add the constant vector in the matrix of explanatory variables).

More generally, consider the following function (to make sure that D is a definite-positive matrix, we use the nearPD function).

svm.fit = function(X, y, C=NULL) {
 n.samples = nrow(X)
 n.features = ncol(X)
 K = matrix(rep(0, n.samples*n.samples), nrow=n.samples)
 for (i in 1:n.samples){
  for (j in 1:n.samples){
   K[i,j] = X[i,] %*% X[j,] }}
 Dmat = outer(y,y) * K
 Dmat = as.matrix(nearPD(Dmat)$mat) 
 dvec = rep(1, n.samples)
 Amat = rbind(y, diag(n.samples), -1*diag(n.samples))
 bvec = c(0, rep(0, n.samples), rep(-C, n.samples))
 res = solve.QP(Dmat,dvec,t(Amat),bvec=bvec, meq=1)
 a = res$solution 
 bomega = apply(a*y*X,2,sum)

On our dataset, we obtain

M = as.matrix(myocarde[,1:7])
center = function(z) (z-mean(z))/sd(z)
for(j in 1:7) M[,j] = center(M[,j])
bomega = svm.fit(cbind(1,M),myocarde$PRONO*2-1,C=.5)
y_pred = 2*((cbind(1,M)%*%bomega)>0)-1
obs  -1  1
  -1 27  2
  1   9 33

i.e. 11 misclassification, out of 71 points (which is also what we got with the logistic regression).

Kernel Based Approach

In some cases, it might be difficult to “separate” by a linear separators the two sets of points, like below,

It might be difficult, here, because which want to find a straight line in the two dimensional space (x_1,x_2). But maybe, we can distort the space, possible by adding another dimension

That’s heuristically the idea. Because on the case above, in dimension 3, the set of points is now linearly separable. And the trick to do so is to use a kernel. The difficult task is to find the good one (if any).

A positive kernel on \mathcal{X} is a function K:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R} symmetric, and such that for any n, \forall\alpha_1,\cdots,\alpha_n and \forall\mathbf{x}_1,\cdots,\mathbf{x}_n,\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_j k(\mathbf{x}_i,\mathbf{x}_j)\geq 0.
For example, the linear kernel is k(\mathbf{x}_i,\mathbf{x}_j)=\mathbf{x}_i^T\mathbf{x}_j. That’s what we’ve been using here, so far. One can also define the product kernel k(\mathbf{x}_i,\mathbf{x}_j)=\kappa(\mathbf{x}_i)\cdot\kappa(\mathbf{x}_j) where \kappa is some function \mathcal{X}\rightarrow\mathbb{R}.

Finally, the Gaussian kernel is k(\mathbf{x}_i,\mathbf{x}_j)=\exp[-\|\mathbf{x}_i-\mathbf{x}_j\|^2].

Since it is a function of \|\mathbf{x}_i-\mathbf{x}_j\|, it is also called a radial kernel.

linear.kernel = function(x1, x2) {
 return (x1%*%x2)
svm.fit = function(X, y, FUN=linear.kernel, C=NULL) {
 n.samples = nrow(X)
 n.features = ncol(X)
 K = matrix(rep(0, n.samples*n.samples), nrow=n.samples)
 for (i in 1:n.samples){
  for (j in 1:n.samples){
   K[i,j] = FUN(X[i,], X[j,])
 Dmat = outer(y,y) * K
 Dmat = as.matrix(nearPD(Dmat)$mat) 
 dvec = rep(1, n.samples)
 Amat = rbind(y, diag(n.samples), -1*diag(n.samples))
 bvec = c(0, rep(0, n.samples), rep(-C, n.samples))
 res = solve.QP(Dmat,dvec,t(Amat),bvec=bvec, meq=1)
 a = res$solution 
 bomega = apply(a*y*X,2,sum)

Link to the regression

To relate this duality optimization problem to OLS, recall that y=\mathbf{x}^T\mathbf{\omega}+\varepsilon, so that \widehat{y}=\mathbf{x}^T\widehat{\mathbf{\omega}}, where \widehat{\mathbf{\omega}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}
But one can also write y=\mathbf{x}^T\widehat{\mathbf{\omega}}=\sum_{i=1}^n \widehat{\alpha}_i\cdot \mathbf{x}^T\mathbf{x}_i
where \widehat{\mathbf{\alpha}}=\mathbf{X}[\mathbf{X}^T\mathbf{X}]^{-1}\widehat{\mathbf{\omega}}, or conversely, \widehat{\mathbf{\omega}}=\mathbf{X}^T\widehat{\mathbf{\alpha}}.

Application (on our small dataset)

One can actually use a dedicated R package to run a SVM. To get the linear kernel, use

df0 = df
df0$y = 2*(df$y=="1")-1
SVM1 = ksvm(y ~ x1 + x2, data = df0, C=.5, kernel = "vanilladot" , type="C-svc")

Since the dataset is not linearly separable, there will be some mistakes here

     -1 1
  -1  2 2
  1   1 5

The problem with that function is that it cannot be used to get a prediction for other points than those in the sample (and I could neither extract \omega nor b from the 24 slots of that objet). But it’s possible by adding a small option in the function

SVM2 = ksvm(y ~ x1 + x2, data = df0, C=.5, kernel = "vanilladot" , prob.model=TRUE, type="C-svc")

With that function, we convert the distance as some sort of probability. Someday, I will try to replicate the probabilistic version of SVM, I promise, but today, the goal is just to understand what is done when running the SVM algorithm. To visualize the prediction, use

pred_SVM2 = function(x,y){
return(predict(SVM2,newdata=data.frame(x1=x,x2=y), type="probabilities")[,2])}
vu = seq(-.1,1.1,length=251)
vv = outer(vu,vu,function(x,y) pred_SVM2(x,y))
contour(vu,vu,vv,add=TRUE,lwd=2,nlevels = .5,col="red")

Here the cost is C=.5, but of course, we can change it

SVM2 = ksvm(y ~ x1 + x2, data = df0, C=2, kernel = "vanilladot" , prob.model=TRUE, type="C-svc")
pred_SVM2 = function(x,y){
return(predict(SVM2,newdata=data.frame(x1=x,x2=y), type="probabilities")[,2])}
vu = seq(-.1,1.1,length=251)
vv = outer(vu,vu,function(x,y) pred_SVM2(x,y))
contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5,col="red")

As expected, we have a linear separator. But slightly different. Now, let us consider the “Radial Basis Gaussian kernel”

SVM3 = ksvm(y ~ x1 + x2, data = df0, C=2, kernel = "rbfdot" , prob.model=TRUE, type="C-svc")

Observe that here, we’ve been able to separare the white and the black points

     -1 1
  -1  4 0
  1   0 6
vu = seq(-.1,1.1,length=251)
vv = outer(vu,vu,function(x,y) pred_SVM3(x,y))
contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5,col="red")

Now, to be completely honest, if I understand the theory of the algorithm used to compute \omega and b with linear kernel (using quadratic programming), I do not feel confortable with this R function. Especially if you run it several times… you can get (with exactly the same set of parameters)


(to be continued…)

Machines, procedures and avoiding responsibility

Some people are trying to make us believe that artificial intelligence is a “revolution”. What if it wasn’t? Can we not simply see the logic of a process that goes back at least fifty years ago? Bureaucracy has pushed us to put in place simple procedures in all areas of everyday life, allowing everyone to avoid any responsibility, to no longer have to think and to be smart. Algorithms are scary, we wonder where the “human” is in these decision-making procedures… What if he had already disappeared long ago?

Insurance and bureaucratic procedures

As Crozier (1963) had shown, bureaucracy appears as a rational mode of organization, in which citizens are protected from clientelism and arbitrariness by the establishment of objective rules. Because bureaucracy is not only a hierarchical (and state) administration, it is above all a set of norms, procedures and formalities that encompass all human activities. Insurers, for example, quickly understood the value of these procedures. In the context of road accidents, the settlement of material accidents is made according to the “IRSA Convention” (Agreement on Direct Indemnification of the Insured and Recourse between Motor Insurance Companies) when two vehicles are involved (succeeding the so-called IDA Convention, created in 1968). After a damage assessment by an expert, the insurer establishes the liability of his insured and directly compensates him for the material damage and prejudice suffered. He then turns against the opposing insurer(s) according to the terms of recourse established by the agreement. This convention is based on the use of scales, corresponding to 13 cases of classic accidents: the case 51 for example relates to a collision between two vehicles, one of which was backing up or making a half-turn, and in this case, the responsibility is normally full (100%) for the driver of the vehicle reversing (or turning half-turn). In the case of personal injury accidents, the procedure is less transparent, but the so-called Badinter law proposed to set up a simple mechanism (for the victim but also, finally, for the driver’s insurer), with the aim of “improving the situation of road accident victims and speeding up compensation procedures”. To do this, several scales are used. For example, Partial Permanent Disability (P.P.I.) measures the “reduction in physical, psycho-sensory or intellectual potential that a victim remains with” and translates into a percentage of permanent disability, on a scale from 1 to 100, corresponding to an “indicative scale of functional deficits after sequellaires in common law”. If this scale is officially “indicative”, following it avoids having to justify oneself. The whole assessment of compensation will be based on “medical certificates”, which have become the basic administrative documents, themselves increasingly standardised, and meeting specific forms and standards. In a very bureaucratic logic, these procedures aim to avoid empathy being encouraged. The development of procedural standards also has this objective: to distance the victims, so that neutrality can be exercised, objectivity can be (supposedly) ensured and justice can be effective. Isn’t this form of indifference the basis of everyday insurance? Don’t claims managers want to be as far away from the caller as possible because the roof of his house has caught on fire, or because his wife is in hospital following a traffic accident?

The bureaucratic production of indifference

Bauman (2002) shows that this “production of indifference” often emerges from an extremely banal configuration, characteristic of all modern societies, namely the coupling of the “functional division of labour” and the “substitution of technical responsibility for moral responsibility”. In describing the bureaucracy in Pakistan, Hull (2012) explains that “these procedures are developed not because of a rationalization logic, but because public servants protect themselves by deploying them vigorously and widely. It is not surprising that Adolf Eichmann chose this line of defence at his trial in April 1961. In fact, when Arendt (1966), present during the trial, is reread about Adolf Eichmann’s responsibility in the implementation of the “final solution” (then sent by The New Yorker to Jerusalem to testify), she presents us with a “robotized” character, disubstantialized, having only obeyed orders. Without going back over the historical reality of Eichmann’s role and the staging of the trial, Hannah Arendt’s analysis of the bureaucracy is interesting. The subtitle of the book, the “banality of evil”, is characterized by the inability to be affected by what one does and the refusal to judge, like any bureaucrat analyzing a situation through the prism of a form, becoming a robot, whose responsibility would be called into question by mechanical obedience to orders. The latter would ultimately be “representative of a bureaucratic system, in which each individual is merely a blind cog… mechanically executing orders from respected authority”. Arendt’s thesis (in Cesarani’s words (2010)) is that “Eichmann was telling the truth when he presented himself as a civil servant without passion, as a tiny cog in the vast exterminating machine, and when he claimed that he could very easily have been replaced by someone else. Obeying orders, following a procedure then makes Eichmann guilty, but the question of responsibility remains open (in particular by distinguishing individual responsibility, and collective responsibility). Returning precisely to the question of responsibility a few years later, Arendt (2005) wonders, “how to judge without clinging to standards, preconceived norms and general rules under which to subsume particular cases? ». Bureaucracy can reassure by its rationality, but frighten by the indifference it generates.

Because bureaucracy is nothing new. The ancient scribes were the first bureaucrats, as Wilford (2001) states. More recently, if Karl Marx has looked at the bureaucracy of industry to study the domination of the bourgeoisie and capitalism, Max Weber has shown that the bureaucracy accommodates all forms of power. He writes that “real domination[is exercised] in the maintenance of daily administration”, “bureaucracy is characterized by a much greater impossibility than one has to escape it”. The company is the privileged place for the development of a bureaucracy, “the requirement of calculability and predictability as rigorous as possible promotes the growth of a special layer of administrators and imposes a certain type of structuring on it” (quoted by Claude Lefort). For Max Weber, bureaucracy is not a parasite (as Marx thought) but a fundamental component of capitalism. Subcontracting, outsourcing, just-in-time work are only possible thanks to practices based on greater bureaucracy. All information should be codified as accurately as possible, no approximation should be made in decisions, and the division of tasks should be certified and standardised. Partitioning protects employees from a sense of responsibility. The “silos” have become the “experts’ paradise”. The segmented and sequential form of the work then protects the members of the organization as Dupuy (2011) shows.

The endless search for standards

This production of indifference was studied at length by Michael Herzfeld, who showed the specificity of state bureaucracy as a public power, and concretizing the denial of any difference: “it offers the effective and generalized capacity to reject those who do not fit into pre-established categories and considered normal” as Hibou (2012) notes. Sainati & Schalchi (2007) also notes this by studying the importance of bureaucracy and standards in the security drift of recent years, “each individual citizen is apprehended according to the category of offender to which he or she is supposed to belong. Everyone is necessarily a suspect in having committed, wanting or being able to commit an offence. This policy of social intolerance will complete the transformation of justice (especially criminal justice) into a total bureaucratic system…”.

Normative inflation, which is observed in the world of justice, has also been observed in the world of finance, with the decline of state authorities (central banks, financial market authorities, various regulators) which no longer intervene directly, but through the imposition of increasingly strict administrative rules. These include Basel II-type management rules, but also the different partitions between the various activities (deposit and merchant banks, advisory activities and market activities, for example). Risk management is fundamentally bureaucratic, involving standards, grids and codes that will generate automatic reactions. This reporting gives a very simplified vision of the activity, but this synthesis of information makes it possible to take decisions more quickly. This phase seems essential given the specificities of each branch, almost preventing a global vision of the bank’s activity. This set of rules and procedures also provides “protection”: in the face of (judicial) uncertainty, the best way to defend yourself is to respect procedures and rules. The goal is not to avoid bankruptcy, but to protect oneself in the event of an accusation. Respect for the rules then becomes more important than their purpose. These rules seem to be a political response to the various crises experienced by the banking world. This bureaucratic inflation brings a form of tranquillity and comfort, while creating a form of dilution of responsibilities. As Hibou (2012) notes, “in the name of individual responsibility, everyone must respect the standards, but compliance with the standards is worth failing responsibility in the event of a problem” (recalling the Kerviel affair in passing). This is noted, for example, by Brunson & Jacobsson (2002) when they state that the audit culture is certainly a culture of responsibility, but of individual responsibility. Collective responsibility is all the more diluted as governments delegate their regulatory power to private actors in a vague manner if the standard is adopted, traceability techniques will make it possible to trace back to the individual responsible for the act causing a failure; whereas if the standard is not, no one will be responsible.

And this search for standards is endless. Thus, many opponents of industrial agriculture opposed the standardization of food products, and in response, developed organic consumption standards. The latter were then questioned by local networks, which wanted to support the concept of “eating local”, and to be recognized, adopted new standards. The response to proceedings is an escalation of proceedings. Max Weber said it in 1920, “when those who are subject to bureaucratic control seek to escape the influence of existing bureaucratic apparatus, normally this is only possible by creating a proper organization that will also be subject to bureaucratization.

From man to machine or seeing machines as human beings…

MacDonald has long been committed to ensuring a consistent product throughout the world. Following the logic of the Taylorism method of organization, the company has written a procedure guide explaining the correct steps to take when salting French fries, filling a glass of soda, etc. In call centers, procedures are followed very scrupulously by employees, with call queue management, scenarios displayed on the screen, and the employee only has to unfold it. Simone Weil spoke of a “scientific organisation of work”, where work was dehumanised, reduced to a state of mechanical energy. Taylorism was the expression in factories of this fascination for science, seeing the human being as a machine. But there is nothing new, since in the middle of the 17th century Thomas Hobbes wrote “since life is nothing more than a movement of limbs, the beginning of which is in some way an inner main part, why could we not say that all automata (machines that move themselves, by springs and wheels, like a watch) have an artificial life? ».

Computers were invented to relieve a number of repetitive tasks by unwinding an algorithm implemented by a human being. Hanoi tricks are a very simple, very repetitive game… but solved by an incredibly simple algorithm. This is probably why this “game” is still taught in all computer science and algorithmics courses, because the “game” itself is actually very boring very quickly. Many training courses no longer teach creativity, but sets of procedures to follow. The Box & Jenkins method of forecasting is a long procedure that we simply follow to the letter, mechanically: we park the series, we model it by an autoregressive process, then we validate the hypotheses. And if it doesn’t work, we do it again.

The actuarial profession, as a science, is based on a set of simple procedures. For example, to construct a rate, we start by building a base, using underwriting information, we will look at the number of claims for each policy, and we will use a model to describe this counting variable (typically a Poisson regression). If we wish to do more advanced things, we will use part of the data to build the model, and another part to test the predictions of our model. We will do the same thing with claims costs. The approach is simple: we collect data, we estimate a model, we test the model, possibly we retain the best if we have the choice between several. It’s so simple that a computer can almost do it itself…

Machines everywhere!?

For decades, organizations have tried to put in place procedures to avoid arbitrary (humane?) decision-making. At the same time, engineers have developed increasingly powerful machines to repeat routine tasks over and over again. From this double observation, we cannot be surprised to see the machines more and more present, everywhere. At least this is the thesis defended by Susskind & Susskind (2015) which anticipates many transformations of the majority of professions, going beyond a simple robotization of routine tasks (including for legal professions, as Remus & Levy sees it (2015)).

But we must not be mistaken, if machines replace men, they are not men for all that. For example, responsibility always lies with one person: the designer of the machine, its operator – who has some technical knowledge – or the user – who often has none. Abiteboul (2016) goes even further by asking the underlying ethical question, namely “can a computer system be assigned a responsibility” (and joining Hannah Arendt’s questions). But beyond responsibility, one can argue that the machines have no intention. They sometimes perform complex tasks, but because they have been programmed to do so. That’s the beauty of engineering: a phone (as smart as it is) is just a set of electronic components that can perform increasingly complex calculations (sometimes it can locate it on a map and offer me a quick route to a restaurant), but the phone wants nothing.

An automatic translator makes it possible to translate a text in a few seconds, but only because it has been programmed for that, whereas if a child learns a language, it is because he understands that it is an essential step to communicate with his parents. If I type a sentence in an unknown language into a translator, I am impressed to get an answer that makes sense. But if I master both languages (even a little), I am on the contrary often disappointed, probably hoping better. These machines are often very predictable, which is both a quality and a defect. Isn’t that what you ask every engineer? The machine must be reliable, obey the finger and the eye. The first difference between man and machine is that man can disobey. And this is his greatest wealth!

A word that comes up again and again when we practice data science is the word “learning” (we speak of “machine learning”). We are told that the machine tries to reproduce the human walk of the human being when he learns. The child learns to recognize letters, then learns to put them end to end to write words, in a fairly conscious way since he is usually guided by his teacher. He also learns to recognize faces, often unconsciously this time. The machine also learns to recognize handwriting or faces. The machine will make mistakes, but it will also “learn from its mistakes”, and it will improve. Like a human being? It’s not as simple as that. The child will get discouraged several times before being able to read correctly, and he will persevere. We find the “conatus” dear to Spinoza, the effort that everything that actually exists makes, the perseverance that every living being shows. And often, this learning is done in pain. If a machine makes a mistake, recognizing a “5” when it was a “3”, it makes a mistake, and for it it stops there. For the child who hoped to have 5 candies and has only 3, he will really feel the mistake. What about this person who thought he had an important appointment at 3:00, and who arrives two hours late?

To pretend that a telephone can be “smart” is probably indicative of what we think intelligence is. But machines are far from being able to replace human beings. They don’t want to, and fortunately, as long as they’re not programmed to, that’s not going to happen.


Abiteboul, Serge 2016. Informaticiens, tous coupables ? Le Monde, 1er septembre 2016

Arendt, Hannah 1966. Eichmann à Jérusalem : Rapport sur la banalité du mal. Folio Histoire.

Arendt, Hannah 2005. Responsabilité et jugement. 2005.

Bauman, Zygmunt 2002. Modernité et Holocauste. La Fabrique.

Brunson, Niels & Jacobsson, Bengt. . 2002. A World of Standards. Oxford University Press.

Cesarani. David. 2010. Adolf Eichmann .Tallandier.

Crozier, Michel. 1963. Le Phénomène bureaucratique. Seuil.

Dupuy, François. 2011. Lost in Management.: La vie quotidienne des entreprises au XXIe siècle. Seuil.

Herzfeld, Michael. 1992. The Social Production of Indifference. The University of Chicago Press.

Hibou, Isabelle. 2012. La bureaucratisation du monde à l’ère néolibérale, La Découverte.

Hobbes, Thomas. 1651. Leviathan. (version Folio Essais)

Hull, Matthew. 2012. Government of Paper The Materiality of Bureaucracy in Urban Pakistan. University of California Press.

Lefort, Claude. 1979. Éléments d’une critique de la bureaucratie, Paris, Gallimard.

Remus, Dana & Levy, Frank Can Robots Be Lawyers? Computers, Lawyers, and the Practice of Law. Working Paper

Sainati, Gilles & Schalchli, Ulrich. 2007. La décadence sécuritaire. La Fabrique.

Susskind Richard & Susskind Dana. 2015. The Future of the Professions: How Technology Will Transform the Work of Human Experts. Oxford University Press.

Weber, Max. 1921. La domination légale à direction administrative bureaucratique. Paru dans Économie et Société.

Weil, Simone 1999. L’Enracinement, Œuvres, Quarto Gallimard.

Wilford, John 2011. Greek Tablet May Shed Light on Early Bureaucratic Practices. New York Times, 4 Avril 2011

i One may wonder whether the situationists were right when they criticized the Fordist (bureaucratic) society, in which the certainty of not dying of hunger had been exchanged for the risk of dying of boredom.


Regression tree using Gini’s index

In order to illustrate the construction of regression tree (using the CART methodology), consider the following simulated dataset,

> set.seed(1)
> n=200
> X1=runif(n)
> X2=runif(n)
> P=.8*(X1<.3)*(X2<.5)+
+   .2*(X1<.3)*(X2>.5)+
+   .8*(X1>.3)*(X1<.85)*(X2<.3)+
+   .2*(X1>.3)*(X1<.85)*(X2>.3)+
+   .8*(X1>.85)*(X2<.7)+
+   .2*(X1>.85)*(X2>.7) 
> Y=rbinom(n,size=1,P)  
> B=data.frame(Y,X1,X2)

with one dichotomos varible (the variable of interest, ), and two continuous ones (the explanatory ones  and ).

> tail(B)
    Y        X1        X2
195 0 0.2832325 0.1548510
196 0 0.5905732 0.3483021
197 0 0.1103606 0.6598210
198 0 0.8405070 0.3117724
199 0 0.3179637 0.3515734
200 1 0.7828513 0.1478457

The theoretical partition is the following

Here, the sample can be plotted below (be careful, the first variate is on the y-axis above, and the x-axis below) with blue dots when  equals one, and red dots when  is null,

> plot(X1,X2,col="white")
> points(X1[Y=="1"],X2[Y=="1"],col="blue",pch=19)
> points(X1[Y=="0"],X2[Y=="0"],col="red",pch=19)

In order to construct the tree, we need a partition critera. The most standard one is probably Gini’s index, which can be writen, when ‘s are splited in two classes, denoted here 

L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/arbre-comp-04.png” ne peut être affichée car elle contient des erreurs.

or when ‘s are splited in three classes, denoted 

etc. Here,  are just counts of observations that belong to partition  such that  takes value . But it is possible to consider other criteria, such as the chi-square distance,


where, classically

when we consider two classes (one knot) or, in the case of three classes (two knots)

Here again, the idea is to maximize that distance: the idea is to discriminate, so we want samples as not independent as possible. To compute Gini’s index consider

> GINI=function(y,i){
+ T=table(y,i)
+ nx=apply(T,2,sum)
+ pxy=T/matrix(rep(nx,each=2),2,ncol(T))
+ vxy=pxy*(1-pxy)
+ zx=apply(vxy,2,sum)
+ n=sum(T)
+ -sum(nx/n*zx)
+ }

We simply construct the contingency table, and then, compute the quantity given above. Assume, first, that there is only one explanatory variable. We split the sample in two, with all possible spliting values , i.e.

Then, we compute Gini’s index, for all those values. The knot is the value that maximizes Gini’s index. Once we have our first knot, we keep it (call it, from now on ). And we reiterate, by seeking the best second choice: given one knot, consider the value that splits the sample in three, and give the highest Gini’s index, Thus, we consider either the following partition

or this one

I.e. we cut either below, or above the previous knot. And we iterate. The code can be something like that,

> X=X2
> u=(sort(X)[2:n]+sort(X)[1:(n-1)])/2
> knot=NULL
> for(s in 1:4){
+ vgini=rep(NA,length(u))
+ for(i in 1:length(u)){
+ kn=c(knot,u[i])
+ F=function(x){sum(x<=kn)}
+ I=Vectorize(F)(X)
+ vgini[i]=GINI(Y,I)
+ }
+ plot(u,vgini)
+ k=which.max(vgini)
+ cat("knot",k,u[k],"\n")
+ knot=c(knot,u[k])
+ u=u[-k]
+ }
knot 69 0.3025479 
knot 133 0.5846202 
knot 72 0.3148172 
knot 111 0.4811517

At the first step, the value of Gini’s index was the following,

which was maximal around 0.3. Then, this value is considered as fixed. And we try to construct a partition in three parts (spliting either below or above 0.3). We get the following plot for Gini’s index (as a function of this second knot)

 which is maximum when the split the sample around 0.6 (which becomes our second knot). Etc. Now, let us compare our code with the standard R function,

> tree(Y~X2,method="gini")
node), split, n, deviance, yval
      * denotes terminal node

 1) root 200 49.8800 0.4750  
   2) X2 < 0.302548 69 12.8100 0.7536 *
   3) X2 > 0.302548 131 28.8900 0.3282  
     6) X2 < 0.58462 65 16.1500 0.4615  
      12) X2 < 0.324591 7  0.8571 0.1429 *
      13) X2 > 0.324591 58 14.5000 0.5000 *
     7) X2 > 0.58462 66 10.4400 0.1970 *

We do obtain similar knots: the first one is 0.302 and the second one 0.584. So, constructing tree is not that difficult…

Now, what if we consider our two explanatory variables? The story remains the same, except that the partition is now a bit more complex to write. To find the first knot, we consider all values on the two components, and again, keep the one that maximizes Gini’s index,

> n=nrow(B)
> u1=(sort(X1)[2:n]+sort(X1)[1:(n-1)])/2
> u2=(sort(X2)[2:n]+sort(X2)[1:(n-1)])/2
> gini=matrix(NA,nrow(B)-1,2)
> for(i in 1:length(u1)){
+ I=(X1<u1[i])
+ gini[i,1]=GINI(Y,I)
+ I=(X2<u2[i])
+ gini[i,2]=GINI(Y,I)
+ }
> mg=max(gini)
> i=1+sum(mg==max(gini[,2]))
> par(mfrow = c(1, 2))
> plot(u1,gini[,1],ylim=range(gini),col="green",type="b",xlab="X1",ylab="Gini index")
> abline(h=mg,lty=2,col="red")
> if(i==1){points(u1[which.max(gini[,1])],mg,pch=19,col="red")
+          segments(u1[which.max(gini[,1])],mg,u1[which.max(gini[,1])],-100000)}
> plot(u2,gini[,2],ylim=range(gini),col="green",type="b",xlab="X2",ylab="Gini index")
> abline(h=mg,lty=2,col="red")
> if(i==2){points(u2[which.max(gini[,2])],mg,pch=19,col="red")
+          segments(u2[which.max(gini[,2])],mg,u2[which.max(gini[,2])],-100000)}
> u2[which.max(gini[,2])]
[1] 0.3025479

The graphs are the following: either we split on the first component (and we obtain the partition on the right, below),

or we split on the second one (and we get the following partition),

Here, it is optimal to split on the second variate, first. And actually, we get back to the one-dimensional case discussed previously: as expected, it is optimal to split around 0.3. This is confirmed with the code below,

> library(tree)
> arbre=tree(Y~X1+X2,data=B,method="gini")
> arbre$frame[1:4,]
     var   n       dev      yval splits.cutleft splits.cutright
1     X2 200 49.875000 0.4750000      <0.302548       >0.302548
2     X1  69 12.811594 0.7536232      <0.800113       >0.800113
4 <leaf>  57  8.877193 0.8070175                               
5 <leaf>  12  3.000000 0.5000000

For the second knot, four cases should be considered: spliting on the second variable (again), either above, or below the previous knot (see below on the left) or spliting on the first one. Then whe have wither a partition below or above the previous knot (see below on the right),

Etc. To visualize the tree, the code is the following

> plot(arbre)
> text(arbre)
> partition.tree(arbre)


Note that we can also visualize the partition. Nice, isn’t it?

To go further, the book Classification and Regression Trees by Leo Breiman (and co-authors) is awesome. Note that there are also interesting sections in the bible Elements of Statistical Learning: Data Mining, Inference, and Prediction by Trevor Hastie, Robert Tibshirani and Jerome Friedman (which can be downloaded from http://www.stanford.edu/~hastie/…)