Matrice singulière dans mlogit

Version imprimable

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:

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:

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:

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:

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.

09/11/2011, 09h45
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:

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:

traceback()

après l'erreur en question.

Code:

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.

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

Code:

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
}