Discussion :

[Artturc]

Bonjour,

Je suis en train de développer un algorithme d'optimisation d'un modèle de choix à classes latentes. Pour vérifier que mon algo fonctionne, je génère des données à partir de coefficients arbitraires, coefficients que j'espère retrouver avec mon algo.

Pour l'instant je ne retrouve pas les bons coefficients, et il est possible que ce soit ma génération de données elle-même qui fonctionne mal. Pour vérifier que ma génération de données fonctionne, je veux faire tourner des fonctions déjà implémentées dans R, notamment mlogit() (package mlogit). Seulement, j'obtiens l'erreur suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
> reg<-mlogit(choice~X1+X2+X3+X4+X5,data.shape)
Erreur dans drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) : 
  sous-programme Lapack dgesv : le système est exactement singulier

alors que mes données ne sont absolument pas colinéaires, puisque générées aléatoirement de façon indépendante les unes des autres.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
> head(data.shape)
            X1         X2         X3           X4        X5 choice id.var chid.var alt.var
1.1 -3.6705006  4.0189001  0.3673018 -0.001927193 -1.127631   TRUE      1        1       1
1.2  0.8894797 -5.7973798 -0.2573684 -4.532502341 -5.280530  FALSE      1        1       2
2.3 -3.1811889  2.0087362 -2.7873105 -1.358131276  5.201941   TRUE      1        2       3
2.4 -0.8415004 -0.7204216  3.6883280  0.726791957  2.776803  FALSE      1        2       4
3.5  0.2129098  2.0541943  1.3524543  3.112635710 -4.315230   TRUE      1        3       5
3.6  0.5743260 -2.8436047  2.6553177 -4.844616086 -2.785944  FALSE      1        3       6

avec id.var la variable identifiant un individu (j'ai 20 individus simulés), chid.var la variable identifiant un choix et alt.var celle identifiant une option, au sein d'un même choix. Les données se lisent ainsi : l'individu 1 a eu à choisir au choix 1 entre les options 1 et 2 (options caractérisées par des valeurs de X1-X5), et il a finalement choisi l'option 1 (choice=TRUE).

X=X1-X5 est bien de rang 5 de même que la matrice XX'. Par contre cbind(X,choice) n'est pas de rang 6 (mais je ne sais pas si ça a quoi que ce soit à faire dans le problème) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
> qr(X)$rank
[1] 5
> qr(X%*%t(X))$rank
[1] 5
> qr(cbind(X,choice))$rank
[1] 5

Les variables X1-X5 sont générées ainsi :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
> matrix(rnorm(5*K,mean=0,sd=2),nrow=K,ncol=5)
            [,1]        [,2]       [,3]       [,4]        [,5]
 [1,] -0.1388375  0.08751987  1.7515582  2.1411558 -2.07855080
 [2,]  5.6999010  1.68412965  2.6812098  1.5934430  3.53517888
 [3,] -1.6832822  1.32272432  1.2920483 -1.1747077 -3.24965556
 [4,] -1.0567598  0.28361254  0.4804106  0.1926974  0.89616776
 [5,] -1.6485392  3.03781734  2.6835301  1.0521919 -0.72972912
 [6,] -4.0830007  0.74012525 -0.8748868  0.8204626 -2.36926277
 [7,] -2.0079312 -1.69787389  1.6984261  2.0790443  3.08570832
 [8,] -1.8974952 -1.39385358  3.4356878  1.3855985 -1.33423274
 [9,]  1.5038651 -1.69747451  0.5177587 -2.3425068  0.08612213
[10,]  2.1027893  0.07387516 -0.5510225  1.2878800 -0.93083105
[11,] ..............

Pour chaque ligne, la valeur S=X1*b1 + X2*b2 + X3*b3 + X4*b4 + X5*b5 est calculée à partir de valeurs de b1 ... b5 fixées au préalable (ce sont les fameux coefficients à estimer). La variable 'choice' est TRUE pour le maximum de S au sein de chaque choix (variable chid.var).

Mon problème principal est que je ne sais pas à quel moment de mlogit() la fonction Lapack dgesv est appellée ni sur quelles matrices... L'appel à cette fonction est extérieur au code de la fonction mlogit, que je ne comprends qu'à moitié.

Merci d'avance pour votre aide.

[pitipoisson]

Bonjour,

Je ne connais pas cette fonction, mais je soupçonne que le problème puisse venir des valeurs de départ des coefficients utilisées par l'algorithme d'optimisation. Tu devrais peut-être essayer d'en spécifier d'autres avec l'argument start (par défaut

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

rep(0, K)

je crois).

Sinon, pour connaître la succession des appels qui ont menés à une erreur, tu peux utiliser la fonction

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

traceback()

après l'erreur en question.

[Artturc]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
>traceback()
9: .Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")
8: drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base"))
7: solve.default(H, g[!fixed])
6: solve(H, g[!fixed])
5: as.vector(solve(H, g[!fixed]))
4: mlogit.optim(method = "nr", print.level = 0, start = c(0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0), logLik = lnl.mlogits, X = X, y = y, 
       opposite = TRUE)
3: eval(expr, envir, enclos)
2: eval(opt, sys.frame(which = nframe))
1: mlogit(choice ~ X1 + X2 + X3 + X4 + X5, data.shape)

... et changer start ne sert à rien. J'ai essayé (1,1,1, ... 1,1), (1,2,1,2, ... 1,2) et (0.01,0.02, ... 2.99,3) sans succès. Pour information, il y a 300 zéros, ce qui correspond au nombre de lignes de data.shape.

[Artturc]

L'erreur a lieu à la ligne 47 de mlogit.optim() :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
> mlogit.optim
function (logLik, start, method = c("bfgs", "nr", "bhhh"), iterlim = 2000, 
    tol = 1e-06, ftol = 1e-08, steptol = 1e-10, print.level = 0, 
    constPar = NULL, ...) 
{
    method <- match.arg(method)
    param <- start
    callT <- match.call(expand.dots = TRUE)
    optimoptions <- c("iterlim", "tol", "method", "print.level", 
        "constPar", "ftol", "steptol")
    chi2 <- 1e+10
    i <- 0
    K <- length(param)
    d <- rep(0, K)
    fixed <- rep(FALSE, K)
    names(fixed) <- names(start)
    if (!is.null(constPar)) 
        fixed[constPar] <- TRUE
    f <- callT
    m <- match(optimoptions, names(callT), 0L)
    if (sum(m)) 
        f <- f[-m]
    f[[1]] <- as.name(f[[2]])
    f$gradient <- TRUE
    f$steptol <- steptol
    if (method == "nr") 
        f$hessian <- TRUE
    else f$hessian <- FALSE
    f[[2]] <- NULL
    names(f)[2] <- "param"
    x <- eval(f, parent.frame())
    if (print.level > 0) 
        cat(paste("Initial value of the function :", as.numeric(x), 
            "\n"))
    g <- attr(x, "gradient")
    if (method == "nr") 
        H <- attr(x, "hessian")[!fixed, !fixed]
    if (method == "bhhh") 
        H <- crossprod(attr(x, "gradi")[, !fixed])
    if (method == "bfgs") 
        Hm1 <- solve(crossprod(attr(x, "gradi")[, !fixed]))
    repeat {
        oldx <- x
        oldg <- g
        if (method == "bfgs") 
            d[!fixed] <- -as.vector(Hm1 %*% g[!fixed])
        else d[!fixed] <- -as.vector(solve(H, g[!fixed]))
        i <- i + 1
        if (i > iterlim) {
            code <- 4
            break
        }
        f$param <- param
        f$direction <- d
        f$initial.value <- oldx
        x <- eval(f, parent.frame())
        if (is.null(x)) {
            code = 3
            break
        }
        if (abs(x - oldx) < ftol) {
            code = 2
            break
        }
        g <- attr(x, "gradient")
        step <- attr(x, "step")
        param[!fixed] <- param[!fixed] + step * d[!fixed]
        if (method == "nr") 
            H <- attr(x, "hessian")[!fixed, !fixed]
        if (method == "bhhh") 
            H <- crossprod(attr(x, "gradi")[, !fixed])
        if (method == "bfgs") {
            incr <- step * d
            y <- g - oldg
            Hm1 <- Hm1 + outer(incr[!fixed], incr[!fixed]) * 
                (sum(y[!fixed] * incr[!fixed]) + as.vector(t(y[!fixed]) %*% 
                  Hm1 %*% y[!fixed]))/sum(incr[!fixed] * y[!fixed])^2 - 
                (Hm1 %*% outer(y[!fixed], incr[!fixed]) + outer(incr[!fixed], 
                  y[!fixed]) %*% Hm1)/sum(incr[!fixed] * y[!fixed])
        }
        chi2 <- -crossprod(d[!fixed], oldg[!fixed])
        if (print.level > 0) {
            chaine <- paste("iteration ", i, ", step = ", step, 
                ", lnL = ", round(x, 8), ", chi2 = ", round(chi2, 
                  8), "\n", sep = "")
            cat(chaine)
        }
        if (print.level > 1) {
            resdet <- rbind(param = param, gradient = g)
            print(round(resdet, 3))
            cat("--------------------------------------------\n")
        }
        if (abs(chi2) < tol) {
            code = 1
            break
        }
    }
    if (code == 3) 
        x <- oldx
    names(attr(x, "gradient")) <- colnames(attr(x, "gradi")) <- names(param)
    attr(x, "fixed") <- fixed
    est.stat = structure(list(elaps.time = NULL, nb.iter = i, 
        eps = chi2, method = method, code = code), class = "est.stat")
    result <- list(optimum = x, coefficients = param, est.stat = est.stat)
    result
}