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
| Option Compare Binary
Option Explicit
Option Base 1
#Const DEBUG_ = False 'Compilation conditionnelle
Private Const EPSILON As Double = 0.00000001 'Précision minimale
Private Const MAX_ITER As Long = 1000 'Nombre maximal d'itérations
Private Const EM As Double = 0.577215664901533 'La constante d'Euler-Mascheroni
'Cette fonction retourne la valeur de la fonction mathématique digamma (psi):
Private Function digamma(x As Double) As Double
Dim s As Double, n As Long
For n = 0 To MAX_ITER
s = s + 1 / (n + 1) - 1 / (n + x)
Next
digamma = s - EM
End Function
'Cette fonction retourne la valeur de la fonction mathématique trigamma (psi-prime):
Private Function trigamma(x As Double) As Double
Dim n As Long, s As Double
For n = 0 To MAX_ITER
s = s + (n + x) ^ -2
Next
trigamma = s
End Function
'Procédure principale:
Public Function calcEMV(ByRef r As Range) As Double()
Dim p(1 To 2) As Double 'p = (alpha, tau)
Dim g(1 To 2) As Double 'le gradient
Dim h(1 To 2, 1 To 2) As Double 'la matrice hessienne
Dim s(1 To 3) As Double
Dim c As Long, t As Variant, k As Variant
Dim n As Long, mu As Double, var As Double
With Application.WorksheetFunction
n = r.Cells.Count
s(1) = 0
s(2) = 0
s(3) = 0
For Each k In r.Cells
s(1) = s(1) + Log(k.Value)
s(2) = s(2) + Log(k.Value) ^ 2
s(3) = s(3) + Log(Log(1 / k.Value))
Next
mu = (-1 / n) * s(1)
var = (s(2) / n) - mu ^ 2
p(1) = mu ^ 2 / var 'alpha-MM
p(2) = mu / var 'tau-MM
c = 0
Do
#If DEBUG_ Then
'À compléter...
#End If
c = c + 1
g(1) = Log(p(2)) - digamma(p(1)) + s(3)
g(2) = p(1) / p(2) + s(1)
h(1, 1) = -n * trigamma(p(1))
h(1, 2) = n / p(2)
h(2, 1) = h(1, 2)
h(2, 2) = -n * (p(1) / p(2) ^ 2)
t = .Transpose(.MMult(.MInverse(h), .Transpose(g)))
p(1) = p(1) - t(1)
p(2) = p(2) - t(2)
Loop Until Sqr(.SumSq(g)) <= EPSILON Or c >= MAX_ITER
End With
calcEMV = p
Exit Function
GESTION_ERREUR:
'À compléter...
End Function
Sub test111()
Dim rng As Range
Set rng = Feuil9.Range("A2:A5")
calcEMV (rng)
End Sub |
Partager