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
| Function foo(y As Range, x As Range)
'compter les variables
Dim i As Long, j As Long, jj As Long
'dimensions des donn�es
Dim K As Long, N As Long
N = y.Rows.Count
K = x.Columns.Count
'initialisation des vecteurs coeff b et score bx avec la moyenne de y
Dim b() As Double, bx() As Double, ybar As Double
ReDim b(1 To K) ' les 6 coefficients du score
ReDim bx(1 To N) ' score en soi par emprunteur
ybar = Application.WorksheetFunction.Average(y)
b(1) = Log(ybar / (1 - ybar))
For i = 1 To N
bx(i) = b(1)
Next i
'Definir variables dans Newton
Dim sens As Double, maxiter As Integer, iter As Integer, delta As Double
Dim lambda() As Double, logV() As Double, dlogV() As Double, hesse() As Double, hinv(), hinvg()
ReDim lambda(1 To N)
'Parametrage de boucle
sens = 1 * 10 ^ (-11): maxiter = 50
ReDim logV(1 To maxiter)
delta = sens + 1: iter = 1: logV(1) = 0
'Boucle pour iterations de Newton sur iter
Do While Abs(delta) > sens And iter < maxiter
iter = iter + 1
'Mise � Zero deriv�e premi�re et seconde
Erase dlogV, hesse
ReDim dlogV(1 To K): ReDim hesse(1 To K, 1 To K)
'calculer prediction Lambda, gradient dlogV, Hessien hesse, et log vraisembl logV
For i = 1 To N
lambda(i) = 1 / (1 + Exp(-bx(i)))
'calcul de la d�riv�e premi�re
For j = 1 To K
dlogV(j) = dlogV(j) + (y(i) - lambda(i)) * x(i, j)
' calcul de la d�riv�e seconde
For jj = 1 To K
hesse(jj, j) = hesse(jj, j) - lambda(i) * (1 - lambda(i)) * x(i, jj) * x(i, j)
Next jj
Next j
'Calcul du Log de la vraisemblance pour cette iteration
logV(iter) = logV(iter) + y(i) * Log(1 / (1 + Exp(-bx(i)))) + (1 - y(i)) * Log(1 - 1 / (1 + Exp(-bx(i))))
Next i
'calcul de l'inverse du Hessien ( = hinv) et le multiplier avec dlogV
hinv = Application.WorksheetFunction.MInverse(hesse)
hinvg = Application.WorksheetFunction.MMult(dlogV, hinv)
' Si converge, sortie et garder le b avec le hessien estim�
delta = logV(iter) - logV(iter - 1)
If Abs(delta) <= sens Then Exit Do
'appliquer Newton schema pour mise a jour coeff b
For j = 1 To K
b(j) = b(j) - hinvg(j)
Next j
'Calcul nouveau score (bx)
For i = 1 To N
bx(i) = 0
For j = 1 To K
bx(i) = bx(i) + b(j) * x(i, j)
Next j
Next i
' Fin de la boucle While
Loop
foo = b
End Function |
