-
|
| .
|
|
|
|
In altro topic, vedi Weibull vs Rayleigh , è spuntata la funzione Gamma di Eulero per la determinazione di uno dei parametri della distribuzione di Weibull.
Approfondendo quella questione ho visto che la funzione diretta Gamma() di Eulero non è implementata tra le funzioni Built-in di Excel che invece ha la funzione Ln.Gamma().
Utilizzando quest'ultima è comunque possibile risalire al valore restituito dalla funzione Gamma().
Adesso vorrei poter scrivere una funzione UDF che implementi direttamente la funzione Gamma di Eulero.
E' lo scopo di questo topic.Attached Image
|
|
|
| .
|
-
|
| .
|
|
|
|
i cosiddetti "workaround" senza scrivere una funzione dedicata sono i seguenti:
Γ(x) = EXP(GAMMALN(x))
Γ(x) = EXP(-1) / GAMMADIST(1, x, 1, FALSE)
il primo fa uso della built-in function GAMMALN(x), il secondo fa uso della distribuzione GAMMA.
non ho volutamente cercato su internet la UDF specifica perché ritengo di notevole interesse l'esposizione che afazio intende fare in questo post.
|
|
|
| .
|
-
|
| .
|
|
|
|
CITAZIONE (reversi @ 18/3/2016, 11:35) i cosiddetti "workaround" senza scrivere una funzione dedicata sono i seguenti:
Γ(x) = EXP(GAMMALN(x))
Γ(x) = EXP(-1) / GAMMADIST(1, x, 1, FALSE)
il primo fa uso della built-in function GAMMALN(x), il secondo fa uso della distribuzione GAMMA.
non ho volutamente cercato su internet la UDF specifica perché ritengo di notevole interesse l'esposizione che afazio intende fare in questo post.
Ringraziandoti per l'interesse manifestato, spero di non deludere le tue aspettative. Ti anticipo però che io non sono un esperto e tranne la funzione che è di mia scrittura (e che proporrò alla fine) quanto riferirò lo si può facilmente trovare su internet.
Per quanto riguarda la UDF in VBA per excel, le mie ricerche non hanno dato esito positivo tranne il fatto di aver trovato la sua implementazione nel pacchetto di Leo Volpi XNumber.dll.
Essendo una DLL non ho potuto leggere il codice e poi quel pacchetto implementa oggetti matematici con aritmetica a grande precisione.
|
|
|
| .
|
-
|
| .
|
|
|
|
Un po' di storia.
La prima volta che mi sono imbattuto, con interesse, su Christian Goldbach è stata in occasione della lettura del libro "L'enigma dei numeri primi". Prima forse ne avevo letto il nome ma non si era destato il mio interesse tanto da poterlo ricordare.
Famosa è la sua congettura sui numeri primi che afferma che ogni naturale maggiore di 2 è sempre esprimibile come somma tra due o più numeri primi.
Per esempio 28 lo si può esprimere come 23+3+2. Provate con qualsiasi numero abbordabile e troverete sempre una somma tra numeri primi eventualmente anche duplicati che da come risultato il numero dato.
La congettura non è stata ancora dimostrata e per chi riesce a farlo c'è un premio ancora valido.
Vi starete chiedendo cosa c'entra Goldbach con la questione della funzione gamma di Eulero.
Goldbach, che nasce a Königsberg (vi dice niente questa cittadina?) ebbe intenso scambio epistolare con il nostro Eulero e proprio da uno dei quesiti proposti da Goldbach ad Eulero che nasce il famoso problema dei sette ponti di Königsberg.
Ma oltre alla sua congettura sui numeri primi ed alla questione del percorso unico sui ponti della sua città natale, Goldbach sollecitò Eulero sulla questione dei fattoriali.
Da qui nasce la funzione Gamma di Eulero.
In sostanza la funzione Gamma altro non è che il fattoriale di un numero reale. Quindi Eulero ha esteso il concetto di fattoriale oltre che agli interi anche ai numeri reali.
3! = 1*2*3 = 6
3.25! =Γ(3.25+1)= 2.54925696671849...
Edited by afazio - 19/3/2016, 23:44
|
|
|
| .
|
-
|
| .
|
|
|
|
Non riporterò le trattazioni presenti nei numerosi siti di matematica che spiegano la funzione Gamma poiché non è questo lo scopo del topic e poichè potrete tranquillamente leggerle da voi, voglio invece applicare usando excel e la sua funzione EXP(LN.GAMMA(x)] il concetto di fattoriale di un numero reale.
Costruite, pertanto, una tabella con la prima colonna dedicata al numero reale x di cui vogliamo il GAMMA(x) e la seconda dedicata a contenere la funzione =EXP(LN.GAMMA(x)]
Proviamo quindi a calcolare GAMMA(3.25).
Il risultato, come vi ho anticipato, è 2.54925696671849...
Per capire come spunta il fattoriale, arricchiamo la tabella inserendo i reali precedenti e successivi che differiscono di 1 partendo da 3.25,
quindi la sequenza 0.25, 1.25; 2.25; 3.25; 4.25; 5.25... basta fermiamoci qui.
Otteniamo la seguente tabella:
Questa ancora non ci fa notare nulla. Aggiungete una terza colonna col prodotto tra il numero x ed il suo GAMMA(x)
Otteniamo:
Da qui notiamo subito che
0.25*Γ(0.25) =Γ(1.25)
Completando la colonna otteniamo la seguente matrice:
Vediamo confermata la relazione:
Γ(x+1)=x*Γ(x)
Ma se ci fate caso, il valore di Γ(5.25) lo possiamo ottenere mediante il prodotto di tutti i reali che precedono il numero e che differiscono di 1 a partire da esso, per il Gamma del più piccolo di essi, cioè:
Γ(5.25) = 4.25*3.25*2.25*1.25*0.25*Γ(0.25)
Basta quindi conoscere il Gamma della parte decimale del numero dato per determinare mediante prodotto delle sequenza ottenuta partendo dal decimale sommando sempre 1 col gamma della parte decimale , il gamma del reale dato.
E' più facile comporre la tabella che non a scrivere il concetto.
Provate.
|
|
|
| .
|
-
|
| .
|
|
|
|
Che succede se invece inseriamo degli interi?
Dovremmo vedere duplicata la funzione fattoriale che conosciamo. Quindi per n=4 dovremmo trovare 1*2*3*4=24
Invece otteniamo il fattoriale del precedente cioe 3!=6:
Questo semplicemente perche la relazione è:
Γ(n)=(n-1)!
oppure
Γ(n+1)=n!
|
|
|
| .
|
-
|
| .
|
|
|
|
E se invece tentiamo di calcolare la funzione Gamma per un reale negativo?
Nulla. Excel non riesce a calcolarli e questo semplicemente per il fatto che noi siamo dovuti passare attraverso la funzione LN.GAMMA() in mancanza della funzione diretta Gamma():
Eppure guardando il grafico della funzione postata nel primo messaggio, vediamo che i valori della funzione esistono pure per reali negativi esclusi gli interi negativi.
Edited by afazio - 19/3/2016, 23:47
|
|
|
| .
|
-
|
| .
|
|
|
|
Vediamo adesso di iniziare a ragionare con un po' di codice.
Anche se non è questo lo scopo, penso sia utile iniziare a vedere come poter affrontare la cosa con la funzione gamma applicata agli interi.
Questo non dovrebbe offrire particolari difficoltà potendo procedere in maniera ricorsiva come per la funziona fattoriale.
In questo caso la ricorsione si basa sulle seguenti relazioni:
La funzione ricorsiva, cioè che chiama se stessa fino a trovare il punto di uscita, è la seguente:
CODICE Private Function Fatt_ricorsivo(n As Long) As Long
If n <= 1 Then
Fatt_ricorsivo = 1
Else
Fatt_ricorsivo = n * Fatt_ricorsivo(n - 1)
End If
End Function
Una funzione di una estrema semplicità che poi è quella utilizzata nei corsi specifici per spiegare il concetto della "ricorsione".
Come per il fattoriale, anche per la funzione Gamma() applicata agli interi possiamo individuare un procedimento ricorsivo basato sulle seguenti relazioni:
con queste è immediato il codice che segue:
CODICE Public Function Gamma_ricorsivo(n As Long) As Long
If n <= 1 Then
Gamma_ricorsivo = 1
Else
Gamma_ricorsivo = (n - 1) * Gamma_ricorsivo(n - 1)
End If
End Function
Notate che il codice è sostanzialmente uguale a quello del fattoriale ricorsivo tranne l'adattamento dei fattori.
Ma a noi serve la funzione gamma() applicata ai reali e questa è poco utile allo scopo.
Edited by afazio - 19/3/2016, 23:47
|
|
|
| .
|
-
|
| .
|
|
|
|
Volendo scrivere un codice ricorsivo per la funzione gamma() applicata ad un reale dovremmo riferirci ad una formula che contiene se stessa anche al secondo membro. Ma gia la definizione di gamma, contiene gamma dato che come abbiamo visto risulta:
estendendo coi prodotti dei termini successivi abbiamo:
Ci manca però il punto di uscita della funzione ricorsiva, cioè un valore noto della stessa funzione che lo si dovrebbe determinare attraverso l'integrale che definisce la funzione stessa:
Ma questo integrale, salvo casi particolari, non è di facile soluzione.
Uno dei casi particolari risolti da Eulero è, per esempio, il valore di Γ(1/2)
In questo caso l'integrale si scrive come:
Il valore di Γ(1/2)=radq(PI) può essere quindi assunto come punto di uscita per il calcolo di tutti i reali esprimibili come:
x= n+1/2
con n= 1, 2 , 3....
quindi per i numeri reali: 1.5; 2.5; 3.5; 4.5; ....
In particolare avremo:
Γ(4.5)= 3.5*2.5*1.5*0.5*Radq(PI) = 11.63173...
che corrisponde al valore ottenuto con la funzione Excel=EXP(LN.GAMM(4.5))
|
|
|
| .
|
-
|
| .
|
|
|
|
vedo che afazio è partito con la creazione della udf.
mi inserisco nel topic per postare - con il consenso dello stesso afazio, perché non volevo interferire con il suo sviluppo - una funzione VB che calcola il gamma di eulero:
CODICE Public Function Gamma(ByVal Dou1 As Double) As Double
Dim I As Integer
Dim Dou2 As Double
Dim Dou3 As Double
If Fix(Dou1) = Dou1 Then
If Dou1 = 0 Or Dou1 = 1 Then
Gamma = 1
Exit Function
Else
Dou2 = 1
For I = 1 To Dou1 - 1
Dou2 = Dou2 * I
Next
Gamma = Dou2
Exit Function
End If
ElseIf Dou1 > 3 Then
Dou2 = 1
For I = 1 To (Fix(Dou1) - 1)
Dou2 = Dou2 * (Dou1 - 1)
Dou1 = Dou1 - 1
Next
Dou3 = Exp(-0.57721566 * Dou1) / Dou1
For I = 1 To 9999
Dou3 = Dou3 * ((1 + (Dou1 / I)) ^ (-1)) * Exp(Dou1 / I)
Next
Dou3 = Dou2 * Dou3
Gamma = Dou3
Else
Dou3 = Exp(-0.577215664901533 * Dou1) / Dou1
For I = 1 To 9999
Dou3 = Dou3 * ((1 + (Dou1 / I)) ^ (-1)) * Exp(Dou1 / I)
Next
Gamma = Dou3
End If
End Function
afazio nota che questo codice fornisce un risultato che differisce (oltre un certo numero di decimali) da quello fornito dalla funzione:
Γ(x) = EXP(GAMMALN(x))
è vero. ma non nel senso che afazio mi ha manifestato.
cioè, non è Γ(x) = EXP(GAMMALN(x)) a fornire un risultato più preciso, ma è vero l'opposto: questa funzione così costruita fornisce risultati la cui precisione non supera le 9, massimo 10 cifre significative. pare inoltre che abbia problemi anche nell'intervallo 1-2.
è comunque bene precisare che non è possibile ottenere alcun risultato esatto con qualunque implementazione di gamma per la stessa ragione per cui non è possibile determinare il valore esatto di pi greco.
|
|
|
| .
|
-
|
| .
|
|
|
|
CITAZIONE (reversi @ 18/3/2016, 21:16) è comunque bene precisare che non è possibile ottenere alcun risultato esatto con qualunque implementazione di gamma per la stessa ragione per cui non è possibile determinare il valore esatto di pi greco.
Giusto. Come per pi.greco o anche e (numero di Nepero), non è possibile ottenere un valore esatto in assoluto per la funzione gamma(x), tuttavia per pigreco sono noti un certo numero di decimali esatti (qualche centinaio o forse migliaio non ricordo) lo stesso dicasi per il numero di Nepero.
Quindi la corretta dicitura sarebbe "esatto fino alla tot cifra decimale" e possiamo solo parlare di "approssimazione alla tot cifra decimale"
E' chiaro che non riuscendo a risolvere l'integrale in maniera esatta si deve fare ricorso a qualche metodo numerico o a qualche formula che approssima.
Edited by afazio - 20/3/2016, 12:37
|
|
|
| .
|
-
|
| .
|
|
|
|
aggiunta
il numero 0.577215664901533 che compare nel codice che ho postato è la cosiddetta "costante di eulero-mascheroni" legata alla funzione gamma.
lorenzo mascheroni è stato un matematico che seppe dare un contributo originale ed innovativo anche nel campo della scienza delle costruzioni indicando, per primo, il reale meccanismo di collasso degli archi in muratura, mediante la formazione di cerniere plastiche. meccanismo confermato poi da navier e definitivamente sistematizzato - poco più di 30 anni fa - da heyman.
CITAZIONE (afazio @ 18/3/2016, 21:24) ... tuttavia per pigreco sono noti un certo numero di decimali esatti (qualche centinaio o forse migliaio non ricordo) lo stesso dicasi per il numero di Nepero.
centinaia di migliaia.
una volta mi chiesi: a che pro? la risposta fu: testano la potenza di calcolo dei computer.
|
|
|
| .
|
-
|
| .
|
|
|
|
Le approssimazioni di gamma(x)
La funzione gamma, salvo casi particolari, non è determinabile con esattezza e questo per il fatto che non possediamo la soluzione dell'integrale che la definisce.
Per calcolare gamma(x) con x reale dobbiamo quindi fare ricorso ad un qualche metodo numerico che ci permette di valutare l'integrale, oppure in alternativa a qualche formula che lo approssima con una fissata precisione.
Molti autori hanno proposto delle formule approssimate. La più semplice è quella proposta da Stirling:
Ho messo insieme su tre diverse colonne i valori restituiti da EXP(LN.GAMMA(x)), quello restituito dalla funzione postata da Reversi e quello ottenuto con l'approssimazione di Stirling:
Notiamo che:
- l'approssimazione di Stirling è troppo grossolana. Tuttavia essa viene tranquillamente applicata come valida approssimazione di Gamma nelle applicazioni statistiche;
- l'approssimazione di Stirling non calcola il gamma per reali inferiori ad 1
- la funzione proposta da reversi si scosta da quella di Excel dopo qualche decimale.
In particolare anche per Γ(1/2) per il quale abbiamo la soluzione esatta [ Γ(1/2)=radq(PI)] i risultati differiscono dalla quinta cifra decimale in poi.
funzione di reversi: Γ(1/2)=1.77243169463232...
funzione EXP(LN.Gamma(1/2))=1.77245385102051
In effetti abbiamo:
RadQ(PI.GRECO) = 1.77245385090552..
Sono costretto a smentire Reversi infatti confrontando il primo con il terzo ottengo:
excel 1.77245385102051
esatt 1.77245385090552..
confrontando il secondo col terzo ottengo:
Rever 1.77243169463232
esatto 1.77245385090552
Il numero di cifre decimali uguali nel primo confronto è 8, il numero di cifre decimali uguali nel secondo confronto è 4.
Si deduce che la funzione excel da un risultato "più preciso" di quello fornito dalla funzione.
Ho anche voluto fare il calcolo "a mano" dei gamma dei reali proposti ed ecco la tabella riassuntiva:
Con calcolo a mano intendo l'applicazione della formula:
Γ(4.5)=3.5*2.5*1.5*0.5*radq(PI.Greco())
------
Purtroppo la precisione di qualsiasi calcolo viene limitata in excel alle prime 15 cifre complessive (parte intera e parte decimale), dopo di che è inutile cercare altro.
Edited by afazio - 20/3/2016, 12:44
|
|
|
| .
|
-
|
| .
|
|
|
|
Per curiosità ho voluto calcolare gamma(0.5) in WXmaxima con la funzione
Function: bffac (expr, n)
Bigfloat version of the factorial (shifted gamma) function. The second argument is how many digits to retain and return, it's a good idea to request a couple of extra.
Categories: Gamma and factorial functions Numerical evaluation
bffac(0.5,50); ho praticamente imposto 50 cifre significative e secondo quanto riportato nell'help il risultato dovrebbe essere fedele fino alla 48 esima cifra.
Il risultato è: 8.8622692545276280387098406224793732052316680157007b-1 che si traduce in:
Wxmaxima: 0.88622692545276280387098406224793732052316680157007
Noto che questo però coincide con il gamma(1.5) che abbiamo finora determinato. Ecco spiegato lo "(shifted gamma)" presente nella definizione.
vediamone l'implacabile confronto con quanto ci restituisce Excel:
EXP(LN.GAMM(1.5)) = 0.886226925395262
UDF Reversi..............= 0.88612723053155
a mano:...................= 0.886226925452758
Ho evidenziato le cifre coincidenti.
Avendo compreso solo dopo il significato della shifted gamma, ho fatto determinare il valore bffac(-0.5,50)
ottenendo: 1.7724538509055064468543668422381998109837106913866b0
che corrisponderebbe a gamma(0.5)
Il confronto:
EXP(LN.GAMM(0.5)) = 1.77245385102051
UDF Reversi..............= 1.77243169463232
a mano:...................= 1.77245385090552
Per poter ottenere risultati con grande numero di cifre superando i limiti imposti da Excel dovremmo utilizzare la libreria XNumber.DLL di Leo volpi che permette di impostare il numero di cifre significative a piacimento (o quasi). Non ho fatto la prova ma se qualcuno ne ha voglia puo' farlo e riportare qui i risultati.
Edited by afazio - 18/3/2016, 23:10
|
|
|
| .
|
-
|
| .
|
|
|
|
L'approssimazione Quasi-iperbolica
Questo metodo di approssimazione è basato sul fatto che la parte di curva nell'intervallo [0 , 1] è quasi una iperbole.
Partendo quindi dall'equazione dell'iperbole gli autori di questo metodo hanno determinatioi coefficienti correttivi da applicare per ottenere la migliore interpolazione di quella parte della curva mediante un polinomio correttivo di settimo grado dell'equazione dell'iperbole 1/x.
Si sfrutterà il valore determinato con questa approssimazione come punto di uscita di una funzione ricorsiva.
Per esempio, volendo calcolare Gamma(8.372) possiamo scrivere il seguente prodotto:
gamm(8.372)=7.372*6.372*5.372*4.372*3.372*2.372*1.372*0.372*gamma(0.372).
Ecco lo stralcio del documento con riportata la descrizione ed i coefficienti.
Di questo metodo non ho scritto nulla e forse cercherò di scrivere una funzione che si basa su esso alla fine.
Edited by afazio - 20/3/2016, 12:57Attached Image
|
|
|
| .
|
33 replies since 18/3/2016, 10:16 1595 views
.
.
.
