Probability Density Function (PDF) di una funzione data
Salve a tutti! Il mio primo post qui.
Inizio con un problema per me molto complesso che cerco di esplicitare in breve. Ho una curva "dose/risposta" che risponde ad una legge data dalla funzione:
$ f(x):1/(1+(theta/x)^gamma) $
Ho una serie di valori sperimentali funzione di $ x $ che dovrebbero "fittare" tale funzione. In base a questi valori sperimentali devo calcolare il valore di $ theta $ e $gamma $ che sono costanti della funzione.
Dopo essermi informato leggendo di tutto e di più ho trovato che la maniera per risolvere questo problema risiede nella determinazione della cosiddetta "stima di massima verosimiglianza" (Maximum Likelihood Extimation - MLE) che consente di ricavare i valori delle costanti più probabili a determinare la mia serie di dati.
Il primo passo del problema dovrebbe essere, se ho capito bene, ricavare la "funzione di verosimiglianza" (Likelihood function: $ L $) che a sua volta dipende da una una "probability density function" (PDF) per tale distribuzione di valori.
Quindi:
1) Come posso ricavare la PDF a partire dalla mia $ f(x) $ ?
Una volta determinata $ PDF = p(x; theta, gamma) $ posso ricavare la $ L $ in questo modo:
$ L=prod_(i = 1)^(n) p(x;theta,gamma) $
Mi piacerebbe sapere se fino ad ora ho scritto castronerie o sono andato bene. Dopodichè mi servirebbe ottenere la $ PDF $ per iniziare i calcoli...
Grazie
Inizio con un problema per me molto complesso che cerco di esplicitare in breve. Ho una curva "dose/risposta" che risponde ad una legge data dalla funzione:
$ f(x):1/(1+(theta/x)^gamma) $
Ho una serie di valori sperimentali funzione di $ x $ che dovrebbero "fittare" tale funzione. In base a questi valori sperimentali devo calcolare il valore di $ theta $ e $gamma $ che sono costanti della funzione.
Dopo essermi informato leggendo di tutto e di più ho trovato che la maniera per risolvere questo problema risiede nella determinazione della cosiddetta "stima di massima verosimiglianza" (Maximum Likelihood Extimation - MLE) che consente di ricavare i valori delle costanti più probabili a determinare la mia serie di dati.
Il primo passo del problema dovrebbe essere, se ho capito bene, ricavare la "funzione di verosimiglianza" (Likelihood function: $ L $) che a sua volta dipende da una una "probability density function" (PDF) per tale distribuzione di valori.
Quindi:
1) Come posso ricavare la PDF a partire dalla mia $ f(x) $ ?
Una volta determinata $ PDF = p(x; theta, gamma) $ posso ricavare la $ L $ in questo modo:
$ L=prod_(i = 1)^(n) p(x;theta,gamma) $
Mi piacerebbe sapere se fino ad ora ho scritto castronerie o sono andato bene. Dopodichè mi servirebbe ottenere la $ PDF $ per iniziare i calcoli...
Grazie
Risposte
"cenzo":
Nota importante: fidarsi poco di Excel!...
Per Excel $(\chi^2_(0.95,1))/2=$1.9207295747449 ( =INV.CHI(0,05;1)/2 )
Per R lo stesso valore vale invece 1.9207294103471 ( options(digits=14); qchisq(0.95,1)/2 )
Ergo: ho adottato il valore di R per il calcolo dei CI... ma.. chissà che altro combina Excel con logaritmi, distrib.norm.st e altro...
A proposito la funzione InvKhi2 della libreria TPmath che utilizzo mi dà il seguente risultato: 1,92072941034706, coerente con quello di R
"cenzo":
il metodo suggerito in questo breve articolo mi sembra più facilmente implementabile.
Ho provato a implementare l'idea di questo articolo.
Non sarà il masimo dell'efficienza ma mi sembra semplice ed efficace e penso che puoi sfruttare in gran parte la libreria citata.
Ipotesi: dobbiamo trovare il limite inferiore di $D_50$. Lo chiamo D50L.
Il valore iniziale per D50L è D50 (stimato con MLE).
Il valore iniziale di maxlogL(D50L) quindi è uguale proprio a MLE.
La funzione da annullare all'inizio vale f = maxlogL - MLE +1.92 = 1.92..
Dobbiamo esplorare la maxlogL per valori D50L più piccoli di D50.
Fissiamo uno step di -0.5.
Se D50=38.68, il primo valore da analizzare sarà D50L=38.18
In questo valore calcolo maxlogL, cioè massimizzo la logL con D50L=38.18 fisso al variare dell'altro parametro (gamma50).
Calcolo quindi la f(D50L).
L'obiettivo è di trovare un intervallo per D50L in cui f cambia segno.
Se f è sempre positiva, faccio un altro step. D50L=37.68. Ricalcolo f e così via.
Ad un certo punto per un certo valore D50L f diventa per la prima volta negativa.
Fermo l'eplorazione, dato che ho un intervallo per D50L in cui f cambia segno (gli ultimi due valori calcolati).
A questo punto lancio l'algoritmo del metodo delle secanti (oppure di bisezione), partendo dagli ultimi due valori trovati dalla griglia.
Li hai nella libreria a pag. 61 e 62.
La funzione che gli devi passare è la f = maxlogL(D50L) - MLE +1.92
in cui maxlogL massimizza la verosimiglianza per D50L fissato, al variare solo dell'altro parametro.
Ho fatto una prova con Excel (modello Probit):
Nel caso di D50 Upper lo step della griglia dovrebbe essere positivo, ad esempio +0.5
Per gamma50 adotterei uno step della griglia più piccolo, ad esempio 0.2.
Magari per velocizzare puoi anche provare con step più grandi, ma attenzione alla stabilità dell'algoritmo, soprattutto nel caso del metodo delle secanti (che dovrebbe essere più veloce del metodo di bisezione, ma la cui convergenza non è garantita).
Che ne pensi: si può fare ?
"cenzo":
Che ne pensi: si può fare ?
Ci sto già lavorando
seguo il tuo "algoritmo" visto che mi hai risparmiato un pò di spremitura di meningi in più 
"cenzo":
La funzione da annullare all'inizio vale f = maxlogL - MLE +1.92 = 1.92..
Scusami, volevi dire f = maxlogL - MLE + 1.92 = 0 ?
"cenzo":
Ho provato a implementare l'idea di questo articolo.
Non sarà il masimo dell'efficienza ma mi sembra semplice ed efficace e penso che puoi sfruttare in gran parte la libreria citata.
Ipotesi: dobbiamo trovare il limite inferiore di $D_50$. Lo chiamo D50L.
Il valore iniziale per D50L è D50 (stimato con MLE).
Il valore iniziale di maxlogL(D50L) quindi è uguale proprio a MLE.
La funzione da annullare all'inizio vale f = maxlogL - MLE +1.92 = 1.92..
[cut]
Quindi stai impementando questo secondo algoritmo ?
"nickwing":
[quote="cenzo"]La funzione da annullare all'inizio vale f = maxlogL - MLE +1.92 = 1.92..
Scusami, volevi dire f = maxlogL - MLE + 1.92 = 0 ?[/quote]
No, credo di avere scritto giusto. All'inizio parto con D50L=D50, quindi:
f = maxlogL(D50L)-MLE+1.92=maxlogL(D50)-MLE+1.92=MLE-MLE+1.92=1.92
In altre parole quando D50L=D50, la maxlogL è proprio la MLE !
Ovviamente alla fine dell'algoritmo avremo f=0. Ma D50L sarà l'estremo inferiore del CI, non sarà più uguale al D50 iniziale.
Nell'immagine allegata si dovrebbero seguire abbastanza bene i vari passaggi (spero!). Nel caso chiedi pure.
Anche se già vedo che ho dimenticato di aggiungere il "diviso 2" nell'intestazione della colonna dove compare $\chi^2_(0.95,1)$.
Si intende che debba essere $(\chi^2_(0.95,1))/2=1.92...$
Per chiarezza aggiungo che MLE è un valore costante, il massimo della likelihood che ha consentito di stimare entrambi i parametri D50 e gamma50.
maxlogL è una funzione che dipende da D50L. Abbasso un po' D50L rispetto al valore iniziale D50 e faccio una nuova massimizzazione della verosimiglianza tenendo fisso D50L, al variare del parametro gamma50.
Mi sta fumando il cervello ma non ne vengo a capo... Mi sa che è più complicato di quanto pensassi, nel senso che il procedimento che mi hai descritto sembra chiaro, tradurlo in un algoritmo in Pascal che sia fluido e, soprattutto, funzionante, un pò meno...
Ti riporto l'implementazione dell'ultimo algoritmo in R per il calcolo di d50Lower nel modello Probit.
Per la minimizzazione ho usato una funzione di R (nlminb)
Ho fatto uso del metodo delle secanti per l'ultima parte.
La formattazione non è perfetta perchè mi manda a capo dei commenti: se copi e incolli in un editor di testo, le cose dovrebbero tornare a posto (uso notepad++ che ha una formattazione predefinita anche per R, il codice risulta molto leggibile).
Può darsi che i problemi che hai incontrato siano riconducibili alla minimizzazione della logL con un parametro fisso.
Ho dovuto creare due funzioni diverse per adattarmi a come lavora la funzione nlminb che minimizza.
Come ho scritto alla fine lo stesso codice funziona anche per d50Upper (estremo sup. del CI) cambiando il passo della griglia con lo stesso valore ma positivo.
Per gli estremi di g50 invece il codice va rifatto (in parte), soprattutto introducendo un'altra funzione da minimizzare (la chiamerei probit.nLL.g50), che tenga conto dello scambio tra parametro fisso e parametro su cui minimizzare.
Vedi se ne cavi qualcosa..
Per la minimizzazione ho usato una funzione di R (nlminb)
Ho fatto uso del metodo delle secanti per l'ultima parte.
La formattazione non è perfetta perchè mi manda a capo dei commenti: se copi e incolli in un editor di testo, le cose dovrebbero tornare a posto (uso notepad++ che ha una formattazione predefinita anche per R, il codice risulta molto leggibile).
# x e y contengono rispettivamente la dose e la risposta della nostra base dati
# negative log likelihood funzione del vettore (d50,g50)
# funzione la cui minimizzione fornisce MLE e la stima d50, g50
probit.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- pnorm(g50*sqrt(2*pi)*(x/d50-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
# minimizzazione della funzione .nll
out <- nlminb(c(40,2.5),probit.nLL) # valori di partenza fissati in d50=40 e g50=2.5 (uguale per tutti i modelli)
MLE <- -out$objective # massima verosimiglianza (segno negativo)
d50 <- out$par[1] # stima d50 con MLE
g50 <- out$par[2] # stima g50 con MLE
# negative log likelihood funzione del solo vettore (g50), su cui minimizzo.
# d50 resterà costante nelle chiamate di questa funzione
# è stato necessario separare d50 e g50 perchè la funzione di minimizzazione lavora col primo parametro
# della funzione e io voglio minimizzare solo su g50, tenendo fisso d50
probit.nLL.d50 <- function (g50,d50) {
p <- pnorm(g50*sqrt(2*pi)*(x/d50-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
pd50 <- 0.5 # passo della griglia. negativo CI Lower, positivo per CI Upper
ca1 <- qchisq(0.95,1)/2 # 1.92072941034706
# Inizializzazioni
d50L1 <- d50 # valore iniziale per d50L1 uguale a d50 stimato prima da MLE
f1 <- ca1 # f1=1.92072941034706
d50L2 <- d50L1 + pd50 # d50L2 è ad un passo della griglia
out <- nlminb(2.5,d50=d50L2,probit.nLL.d50) # minimizzo .nll solo variando g50 (starting value sempre 2.5), con un fissato d50=d50L2
f2 <- -out$objective - MLE + ca1 # calcolo la f2, cioè la f nel punto d50L2
# Ricerca sulla griglia per un intervallo in cui f cambia segno
while(f2 > 0) { # esco quando f2 diventa negativa
d50L1 <- d50L2 # aggiorno d50L1 e f1
f1 <- f2
d50L2 <- d50L1 + pd50 # nuovo passo sulla griglia per d50L2
out <- nlminb(2.5,d50=d50L2,probit.nLL.d50) # minimizzo .nll con d50L2 fisso
f2 <- -out$objective - MLE + ca1 # f2, cioè f calcolata in d50L2
}
# Metodo delle secanti per individuare il punto di nullo della funzione f
while(abs(f2)>1e-14 & abs(d50L2-d50L1)>1e-12 ) { # esco appena è soddisfatta una delle due condizioni (puoi anche provare ad imporre entrambe con un OR)
tmp <- (f2*d50L1-f1*d50L2)/(f2-f1) # memorizzo il valore successivo per d50L2
d50L1 <- d50L2 # aggiorno d50L1, d50L2, f1
d50L2 <- tmp
f1 <- f2
out <- nlminb(2.5,d50=d50L2,probit.nLL.d50)
f2 <- -out$objective - MLE + ca1 # ricalcolo f2, cioè f in d50L2
}
d50Lower <- d50L2 # estremo inferiore del CI per d50
d50Lower
[1] 36.7690291136
d50
[1] 38.683526410575
# Stesso codice ma con pd50 <- 0.5 (positivo invece che negativo)
d50Upper
[1] 40.746565391597Può darsi che i problemi che hai incontrato siano riconducibili alla minimizzazione della logL con un parametro fisso.
Ho dovuto creare due funzioni diverse per adattarmi a come lavora la funzione nlminb che minimizza.
Come ho scritto alla fine lo stesso codice funziona anche per d50Upper (estremo sup. del CI) cambiando il passo della griglia con lo stesso valore ma positivo.
Per gli estremi di g50 invece il codice va rifatto (in parte), soprattutto introducendo un'altra funzione da minimizzare (la chiamerei probit.nLL.g50), che tenga conto dello scambio tra parametro fisso e parametro su cui minimizzare.
Vedi se ne cavi qualcosa..
Grazie, il problema più grosso era che nel realizzare l'operazione $ln(phi^(p))$ oppure $ln((1-phi)^(1-p))$ da parte dell'algoritmo che cerca il minimo (il goldsearch, per una sola incognita ossia quella da massimizzare durante le operazioni di ricerca della MLE con valori modificati di $TCD50$ e $gamma_50$) la ricerca "incocciava" in valori di dell'argomento del logaritmo che erano pari a zero (ad es. se se $phi=1$) e quindi ne ricavavo una marea di eccezioni. Ora dovrei essere alle porte della soluzione (spero ).
).
"nickwing":
$ln(phi^(p))$ oppure $ln((1-phi)^(1-p))$
Intendevi dire \( \ln(p^y) \) e \( \ln(1-p)^{1-y} \) ?
Mi confermi che stai provando il secondo algoritmo ?
Per il goldsearch, che valori X e Y gli passi in input ?
Curiosità: perchè non sei rimasto sul BFGS, facendolo lavorare su una funzione ad una sola variabile? Non funziona ?
"cenzo":
[quote="nickwing"]$ln(phi^(p))$ oppure $ln((1-phi)^(1-p))$
Intendevi dire \( \ln(p^y) \) e \( \ln(1-p)^{1-y} \) ?[/quote]
esattamente, laddove per $phi$ intendo la funzione che calcola la probabilità cumulativa
come sempre difetto nella forma... ahhh se passassi a studiare matematica sul serio ne dovrei fare di strada (all'indietro) per recuperare un pò di allenamento 
"cenzo":
Mi confermi che stai provando il secondo algoritmo ?
sto testando il "tuo".
"cenzo":
Per il goldsearch, che valori X e Y gli passi in input ?
valori "coerenti" con i parametri attesi (per $gamma$ valori tra 0.1 e 10, per TCD fra 10 e 90).
"cenzo":
Curiosità: perchè non sei rimasto sul BFGS, facendolo lavorare su una funzione ad una sola variabile? Non funziona ?
ci avevo pensato ma l'ho fatto solo per pigrizia, perchè BFGS (che lavora eccellentemente visto che mi da i tuoi stessi valori per MLE $gamma$ e $TCD50$ necessita di scrivere una ulteriore funzione "gradiente" e non mi andava di "affollare" troppo il listato, sennò con una finetrella così già comincio a perdermici dentro...
Ok 
Anch'io se ricordi avevo incontrato errori del tipo \( \ln(0) \) lavorando in Excel.
Questo però accadeva solo con il Cloglog e il log-Cloglog (aka Weibull). Gli altri modelli non mi hanno dato questo problema.
valori "coerenti" con i parametri attesi (per $gamma$ valori tra 0.1 e 10, per TCD fra 10 e 90).[/quote]
Il problema potrebbe (?) essere nei valori X e Y che passi alla Goldsearch.
Possibilmente dovrebbero essere valori abbastanza "stretti" intorno ad un intervallo che contiene la soluzione ottima.
Per valori che se ne discostano molto siamo molto lontani dal centro della curva dose-risposta, per cui le probabilità vanno a zero o ad uno, da cui il problema col logaritmo.
E' per questo che avevo pensato di partire dal D50 ottimo (MLE) e diminuirlo un po' alla volta (passo -0,5), controllando ogni volta quanto vale la f. Non appena f cambia segno, possiedo un intervallo per D50 abbastanza stretto in cui lanciare il metodo delle secanti (o di bisezione).
Se ho capito bene, il vantaggio di BFGS è che non gli dai un intervallo tra X e Y, ma un solo punto iniziale in cui si calcola il gradiente.
Ad esempio con la funzione di minimizzazione che ho usato in R, gli ho imposto sempre di partire con \( \gamma_{50}=2.5 \) (vedi il codice), perchè sappiamo che è un valore plausibile per quasi tutti i modelli.
EDIT: ora che ci penso, forse sarebbe stato ancora meglio impostare il \( \gamma_{50} \) iniziale uguale proprio al valore che esce dalla stima MLE iniziale per il modello!
Anch'io se ricordi avevo incontrato errori del tipo \( \ln(0) \) lavorando in Excel.
Questo però accadeva solo con il Cloglog e il log-Cloglog (aka Weibull). Gli altri modelli non mi hanno dato questo problema.
"nickwing":
[quote="cenzo"]Per il goldsearch, che valori X e Y gli passi in input ?
valori "coerenti" con i parametri attesi (per $gamma$ valori tra 0.1 e 10, per TCD fra 10 e 90).[/quote]
Il problema potrebbe (?) essere nei valori X e Y che passi alla Goldsearch.
Possibilmente dovrebbero essere valori abbastanza "stretti" intorno ad un intervallo che contiene la soluzione ottima.
Per valori che se ne discostano molto siamo molto lontani dal centro della curva dose-risposta, per cui le probabilità vanno a zero o ad uno, da cui il problema col logaritmo.
E' per questo che avevo pensato di partire dal D50 ottimo (MLE) e diminuirlo un po' alla volta (passo -0,5), controllando ogni volta quanto vale la f. Non appena f cambia segno, possiedo un intervallo per D50 abbastanza stretto in cui lanciare il metodo delle secanti (o di bisezione).
Se ho capito bene, il vantaggio di BFGS è che non gli dai un intervallo tra X e Y, ma un solo punto iniziale in cui si calcola il gradiente.
Ad esempio con la funzione di minimizzazione che ho usato in R, gli ho imposto sempre di partire con \( \gamma_{50}=2.5 \) (vedi il codice), perchè sappiamo che è un valore plausibile per quasi tutti i modelli.
EDIT: ora che ci penso, forse sarebbe stato ancora meglio impostare il \( \gamma_{50} \) iniziale uguale proprio al valore che esce dalla stima MLE iniziale per il modello!
sto lavorando in questa direzione... e la soluzione mi sembra vicina
Passo al BFGS, è più "stabile" e "converge" meglio (spero). Con GoldSearch ho valori inesatti e il CI si chiude prima del dovuto, non so perchè. Lo avevo testato in un'altra ricerca su un singolo parametro e sembrava funzionare...
Vabbè... riscriverò un pò di codice e poi vediamo
Vabbè... riscriverò un pò di codice e poi vediamo
OK
Nel frattempo ho apportato modifiche al mio algoritmo. Sono passato alla bisezione al posto del metodo delle secanti.
Sarà meno efficiente (peraltro non si nota nemmeno) ma ho sperimentato maggiore stabilità.
Ho adottato un'unica condizione di terminazione sulla differenza tra due soluzioni consecutive per l'estremo dell'intervallo di confidenza.
Ho sempre il problema di \( \ln(0) \) per il cloglog che sto cercando di aggirare..
Su probit e cauchit invece tutto OK. Poi vedrò anche sugli altri, in modo da confermare (o smentire) i risultati ottenuti con Excel.
PS: la settimana prossima sarò via per un breve viaggetto
Nel frattempo ho apportato modifiche al mio algoritmo. Sono passato alla bisezione al posto del metodo delle secanti.
Sarà meno efficiente (peraltro non si nota nemmeno) ma ho sperimentato maggiore stabilità.
Ho adottato un'unica condizione di terminazione sulla differenza tra due soluzioni consecutive per l'estremo dell'intervallo di confidenza.
Ho sempre il problema di \( \ln(0) \) per il cloglog che sto cercando di aggirare..
Su probit e cauchit invece tutto OK. Poi vedrò anche sugli altri, in modo da confermare (o smentire) i risultati ottenuti con Excel.
PS: la settimana prossima sarò via per un breve viaggetto
"cenzo":
OK
Nel frattempo ho apportato modifiche al mio algoritmo. Sono passato alla bisezione al posto del metodo delle secanti.
Sarà meno efficiente (peraltro non si nota nemmeno) ma ho sperimentato maggiore stabilità.
Ho adottato un'unica condizione di terminazione sulla differenza tra due soluzioni consecutive per l'estremo dell'intervallo di confidenza.
nella ricerca del valore non uso la soluzione numerica di una equazione non lineare (ossia il metodo delle bisezioni o quello delle secanti) per "semplicità" di programmazione. Dovrei "nidificare" una funzione che vado a massimizzare (magari con il BFGS, e quindi mi tocca creare una funzione "satellite" per il gradiente) in un'altra funzione che mi calcola una soluzione numerica. Magari potrebbe pure essere un'idea ma preferisco fare le cose in maniera più semplice... senza perdermi...
"cenzo":
PS: la settimana prossima sarò via per un breve viaggetto
E buon viaggio!!!
come non detto... sto provando le bisezioni... un centinaio di linee di codice ancora...
Siiiiiiiiii!!! Ce l'ho fatta
Ecco nello screenshot l'intervallo di confidenza per il $gamma_50$ del modello Lyman (probit):
ed ecco a seguire quello del TCD50
Ecco nello screenshot l'intervallo di confidenza per il $gamma_50$ del modello Lyman (probit):
ed ecco a seguire quello del TCD50
direi che corrispondono ai tuoi
Alla fine ho usato il GoldSearch per massimizzare la MLE temporanea con i valori di $gamma_50$ e $TCD_50$ alternativamente passati da un'altra funzione che ho provveduto a risolvere numericamente con la bisezione, come hai fatto tu. L'algoritmo impiega una manciata di secondi a trovare la soluzione, ho posto un limite di 10000 iterazioni per la funzione bisezione e ho provveduto a passare i parametri in maniera attenta per fare in modo che i limiti andassero rispettivamente in basso (vettore di ricerca del GoldSearch dal $gamma_50temp$ al $gamma_50$ impiegato per il limite inferiore del C.I. e dal $gamma_50$ al $gamma_50temp$ per il limite superiore, similmente ho fatto per il $TCD50$). Mi è costato due settimane di lavoro ma ce l'ho fatta. Bisogna ora capire se l'algoritmo è stabile in tutte le condizioni di dati posti in analisi... ma il più è fatto
Alla fine ho usato il GoldSearch per massimizzare la MLE temporanea con i valori di $gamma_50$ e $TCD_50$ alternativamente passati da un'altra funzione che ho provveduto a risolvere numericamente con la bisezione, come hai fatto tu. L'algoritmo impiega una manciata di secondi a trovare la soluzione, ho posto un limite di 10000 iterazioni per la funzione bisezione e ho provveduto a passare i parametri in maniera attenta per fare in modo che i limiti andassero rispettivamente in basso (vettore di ricerca del GoldSearch dal $gamma_50temp$ al $gamma_50$ impiegato per il limite inferiore del C.I. e dal $gamma_50$ al $gamma_50temp$ per il limite superiore, similmente ho fatto per il $TCD50$). Mi è costato due settimane di lavoro ma ce l'ho fatta. Bisogna ora capire se l'algoritmo è stabile in tutte le condizioni di dati posti in analisi... ma il più è fatto
"nickwing":
Siiiiiiiiii!!! Ce l'ho fatta
Very good!
Posso confermare i risultati ottenuti con Excel per tutti i modelli.
Per Cloglog e log-cloglog è stato necessario qualche accorgimento per gestire le situazioni "ln(0)", come risulta dal codice seguente.
Se vuoi verificare la stabilità del tuo codice, potresti testarlo con questi due modelli.
Sono riuscito ad utilizzare un'unica funzione CI() per tutti e 4 gli estremi.
# Leggo i dati (x dose; y risposta binaria)
ca1 <- qchisq(0.95,1)/2
logit.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- 1/(1+exp(4*g50*(1-x/d50)))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logit.nLL.d50 <- function (a,b) {
p <- 1/(1+exp(4*a*(1-x/b)))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logit.nLL.g50 <- function (a,b) {
p <- 1/(1+exp(4*b*(1-x/a)))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
loglogit.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- 1/(1+(d50/x)^(4*g50))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
loglogit.nLL.d50 <- function (a,b) {
p <- 1/(1+(b/x)^(4*a))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
loglogit.nLL.g50 <- function (a,b) {
p <- 1/(1+(a/x)^(4*b))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
probit.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- pnorm(g50*sqrt(2*pi)*(x/d50-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
probit.nLL.d50 <- function (a,b) {
p <- pnorm(a*sqrt(2*pi)*(x/b-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
probit.nLL.g50 <- function (a,b) {
p <- pnorm(b*sqrt(2*pi)*(x/a-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logprobit.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- pnorm(g50*sqrt(2*pi)*log(x/d50))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logprobit.nLL.d50 <- function (a,b) {
p <- pnorm(a*sqrt(2*pi)*log(x/b))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logprobit.nLL.g50 <- function (a,b) {
p <- pnorm(b*sqrt(2*pi)*log(x/a))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
poisson.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- 0.5^exp(2*g50/log(2)*(1-x/d50))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
poisson.nLL.d50 <- function (a,b) {
p <- 0.5^exp(2*a/log(2)*(1-x/b))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
poisson.nLL.g50 <- function (a,b) {
p <- 0.5^exp(2*b/log(2)*(1-x/a))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logpoisson.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- 0.5^((d50/x)^(2*g50/log(2)))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logpoisson.nLL.d50 <- function (a,b) {
p <- 0.5^((b/x)^(2*a/log(2)))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logpoisson.nLL.g50 <- function (a,b) {
p <- 0.5^((a/x)^(2*b/log(2)))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
# Problema log(0): calcolo basato su cp=1-p e scarto casi noti logL=0
cloglog.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
cp <- 0.5^exp(2*g50/log(2)*(x/d50-1))
ind <- !((y==0 & cp==1) | (y==1 & cp==0))
return(-sum(y[ind]*log(1-cp[ind])+(1-y[ind])*log(cp[ind])))
}
cloglog.nLL.d50 <- function (a,b) {
cp <- 0.5^exp(2*a/log(2)*(x/b-1))
ind <- !((y==0 & cp==1) | (y==1 & cp==0))
return(-sum(y[ind]*log(1-cp[ind])+(1-y[ind])*log(cp[ind])))
}
cloglog.nLL.g50 <- function (a,b) {
cp <- 0.5^exp(2*b/log(2)*(x/a-1))
ind <- !((y==0 & cp==1) | (y==1 & cp==0))
return(-sum(y[ind]*log(1-cp[ind])+(1-y[ind])*log(cp[ind])))
}
# Problema log(0): calcolo basato su cp=1-p e scarto casi noti logL=0
logcloglog.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
cp <- 0.5^((x/d50)^(2*g50/log(2)))
ind <- !((y==0 & cp==1) | (y==1 & cp==0))
return(-sum(y[ind]*log(1-cp[ind])+(1-y[ind])*log(cp[ind])))
}
logcloglog.nLL.d50 <- function (a,b) {
cp <- 0.5^((x/b)^(2*a/log(2)))
ind <- !((y==0 & cp==1) | (y==1 & cp==0))
return(-sum(y[ind]*log(1-cp[ind])+(1-y[ind])*log(cp[ind])))
}
logcloglog.nLL.g50 <- function (a,b) {
cp <- 0.5^((x/a)^(2*b/log(2)))
ind <- !((y==0 & cp==1) | (y==1 & cp==0))
return(-sum(y[ind]*log(1-cp[ind])+(1-y[ind])*log(cp[ind])))
}
cauchit.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- 0.5+1/pi*atan(pi*g50*(x/d50-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
cauchit.nLL.d50 <- function (a,b) {
p <- 0.5+1/pi*atan(pi*a*(x/b-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
cauchit.nLL.g50 <- function (a,b) {
p <- 0.5+1/pi*atan(pi*b*(x/a-1))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logcauchit.nLL <- function (z) {
d50 <- z[1]
g50 <- z[2]
p <- 0.5+1/pi*atan(pi*g50*log(x/d50))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logcauchit.nLL.d50 <- function (a,b) {
p <- 0.5+1/pi*atan(pi*a*log(x/b))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
logcauchit.nLL.g50 <- function (a,b) {
p <- 0.5+1/pi*atan(pi*b*log(x/a))
return(-sum(y*log(p)+(1-y)*log(1-p)))
}
# Costruisco estremo intervallo confidenza 95% (vedi ca1) - profile likelihood
CI <- function(nLL,MLE,parm1,parm2,arg) {
# Inizializzazioni
if (arg=="d50L") dx <- -0.5
else if (arg=="d50U") dx <- +0.5
else if (arg=="g50L") dx <- -0.2
else if (arg=="g50U") dx <- +0.2
x1 <- parm1
f1 <- ca1
x2 <- x1 + dx
out <- nlminb(parm2,nLL,b=x2)
f2 <- -out$objective - MLE + ca1
# Grid search (cerco intervallo in cui f cambia segno)
while(f2 > 0) {
x1 <- x2
f1 <- f2
x2 <- x1 + dx
out <- nlminb(parm2,nLL,b=x2)
f2 <- -out$objective - MLE + ca1
}
# Metodo Bisezione
while(abs(x2-x1)>1e-12) {
xc <- (x1+x2)/2
out <- nlminb(parm2,nLL,b=xc)
fc <- -out$objective - MLE + ca1
if (fc>0) {
x1 <- xc
f1 <- fc
}
else {
x2 <- xc
f2 <- fc
}
}
return((x1+x2)/2)
}
fit <- function(mod) {
if (mod=="mo.logit") {m1 <- logit.nLL; m2 <- logit.nLL.d50; m3 <- logit.nLL.g50}
else if (mod=="mo.loglogit") {m1 <- loglogit.nLL; m2 <- loglogit.nLL.d50; m3 <- loglogit.nLL.g50}
else if (mod=="mo.probit") {m1 <- probit.nLL; m2 <- probit.nLL.d50; m3 <- probit.nLL.g50}
else if (mod=="mo.logprobit") {m1 <- logprobit.nLL; m2 <- logprobit.nLL.d50; m3 <- logprobit.nLL.g50}
else if (mod=="mo.poisson") {m1 <- poisson.nLL; m2 <- poisson.nLL.d50; m3 <- poisson.nLL.g50}
else if (mod=="mo.logpoisson") {m1 <- logpoisson.nLL; m2 <- logpoisson.nLL.d50; m3 <- logpoisson.nLL.g50}
else if (mod=="mo.cloglog") {m1 <- cloglog.nLL; m2 <- cloglog.nLL.d50; m3 <- cloglog.nLL.g50}
else if (mod=="mo.logcloglog") {m1 <- logcloglog.nLL; m2 <- logcloglog.nLL.d50; m3 <- logcloglog.nLL.g50}
else if (mod=="mo.cauchit") {m1 <- cauchit.nLL; m2 <- cauchit.nLL.d50; m3 <- cauchit.nLL.g50}
else if (mod=="mo.logcauchit") {m1 <- logcauchit.nLL; m2 <- logcauchit.nLL.d50; m3 <- logcauchit.nLL.g50}
# MLE
out <- nlminb(c(40,2.5),m1)
MLE <- -out$objective
d50 <- out$par[1]
g50 <- out$par[2]
# CI (intervalli confidenza)
d50L <- CI(m2,MLE,d50,g50,"d50L")
d50U <- CI(m2,MLE,d50,g50,"d50U")
g50L <- CI(m3,MLE,g50,d50,"g50L")
g50U <- CI(m3,MLE,g50,d50,"g50U")
return(list(MLE=MLE,"D50 (95% CI)"=c(d50,d50L,d50U),"g50 (95% CI)"=c(g50,g50L,g50U)))
}
Questi i risultati, in linea con quelli ottenuti con (la macro) Excel:
> fit("mo.logit")
$MLE
[1] -57.53318
$`D50 (95% CI)`
[1] 38.95585 37.07019 40.95478
$`g50 (95% CI)`
[1] 2.395000 1.818529 3.111523
> fit("mo.loglogit")
$MLE
[1] -61.00301
$`D50 (95% CI)`
[1] 38.30490 36.38319 40.33705
$`g50 (95% CI)`
[1] 2.216040 1.690060 2.867823
> fit("mo.probit")
$MLE
[1] -58.56481
$`D50 (95% CI)`
[1] 38.68353 36.76903 40.74657
$`g50 (95% CI)`
[1] 2.048982 1.612506 2.548257
> fit("mo.logprobit")
$MLE
[1] -62.59998
$`D50 (95% CI)`
[1] 37.78386 35.84774 39.87780
$`g50 (95% CI)`
[1] 1.867548 1.479811 2.310868
> fit("mo.poisson")
$MLE
[1] -65.14399
$`D50 (95% CI)`
[1] 37.48554 35.39330 39.82364
$`g50 (95% CI)`
[1] 1.641092 1.331775 1.972332
> fit("mo.logpoisson")
$MLE
[1] -71.65757
$`D50 (95% CI)`
[1] 36.44288 34.28337 38.93023
$`g50 (95% CI)`
[1] 1.440549 1.156031 1.744810
> fit("mo.cloglog")
$MLE
[1] -55.28668
$`D50 (95% CI)`
[1] 40.01209 38.15666 41.90072
$`g50 (95% CI)`
[1] 2.407601 1.828375 3.083486
> fit("mo.logcloglog")
$MLE
[1] -56.75386
$`D50 (95% CI)`
[1] 39.27007 37.38875 41.20812
$`g50 (95% CI)`
[1] 2.258346 1.727629 2.890386
> fit("mo.cauchit")
$MLE
[1] -53.14368
$`D50 (95% CI)`
[1] 39.56398 38.18472 41.25986
$`g50 (95% CI)`
[1] 7.001754 4.055993 12.864001
> fit("mo.logcauchit")
$MLE
[1] -54.36018
$`D50 (95% CI)`
[1] 39.48045 38.12085 41.15144
$`g50 (95% CI)`
[1] 6.961346 4.015679 12.848006
E buon viaggio!!!
Grazie!
Ecco il codice delle due funzioni
e questo è il codice associato al bottone che lancia il calcolo:
Gira tutto qui. Ora devo testarlo con Niemierko e poi Poisson (sono i tre principali modelli che mi interessano perchè i più impiegati)
La precisione di calcolo che ho scelto è un pò alta ( 1.0E-14) ma volevo esplorare quanto fosse "dispendioso" in termini di tempi di calcolo l'algoritmo. Ovviamente una decina di cifre decimali in meno andrebbero ugualmente bene per i miei impieghi.
function MLECILyman(X: Float) : Float; // funzione che calcola i valori massimizzati di gamma50 e TCD50 per il CI
var n: Integer;
Temp, t: Float;
begin
Temp:= 0;
case Form1.RadioGroup1.ItemIndex of
0: begin // calcolo del CI per trovare il CI di gamma50 - trova il massimo temporaneo di TCD50
for n:= 1 to 200 do begin
t:= (DataList[n].EUD - X)/(X/(Sqrt2Pi * TempGamma));
Temp:= Temp + Ln(Power(FNorm(t), DataList[n].Outcome)) + Ln(Power(1 - FNorm(t), 1 - DataList[n].Outcome));
end;
TempMLE:= Temp;
MLECILyman:= - Temp;
end;
1: begin // calcolo del CI per trovare il CI di TCD50 - trova il massimo temporaneto di gamma50
for n:= 1 to 200 do begin
t:= (DataList[n].EUD - TempTCD50)/(TempTCD50*1/(Sqrt2Pi * X));
Temp:= Temp + Ln(Power(FNorm(t), DataList[n].Outcome)) + Ln(Power(1 - FNorm(t), 1 - DataList[n].Outcome));
end;
TempMLE:= Temp;
MLECILyman:= - Temp;
end;
end;
end;
function CalculateCILyman (X: Float) : Float; // funzione che da la soluzione numerica dell'intervallo di confidenza
var
BoundA, BoundB: Float;
Xmin, Ymin: Float;
begin
// limite inferiore del CI di gamma50
if (CI_Bound = low) and (Form1.RadioGroup1.ItemIndex = 0) then begin
TempGamma:= X; // la funzione da risolvere qui è funzione di gamma50
BoundA:= TCD50 - 2;
BoundB:= TCD50;
end;
if (CI_Bound = high) and (Form1.RadioGroup1.ItemIndex = 0) then begin
TempGamma:= X; // la funzione da risolvere qui è funzione di gamma50
BoundA:= TCD50;
BoundB:= TCD50 + 2;
end;
if (CI_Bound = low) and (Form1.RadioGroup1.ItemIndex = 1) then begin
TempTCD50:= X; // la funzione da risolvere qui è funzione di TCD50
BoundA:= gamma50 - 1;
BoundB:= gamma50;
end;
if (CI_Bound = high) and (Form1.RadioGroup1.ItemIndex = 1) then begin
TempTCD50:= X; // la funzione da risolvere qui è funzione
BoundA:= gamma50;
BoundB:= gamma50 + 1;
end;
GoldSearch(MLECILyman, BoundA, BoundB, 10000, 1.0E-14, Xmin, Ymin);
CalculateCILyman:= TempMLE - MLE_L + HalfChi2; // risultato della funzione da risolvere con bisezione
// il risultato di questa funzione viene posto uguale a 0
end;
e questo è il codice associato al bottone che lancia il calcolo:
procedure TForm1.btnLymanCIClick(Sender: TObject);
var
X, Y, F: Float;
begin
Screen.Cursor:= crHourGlass;
case RadioGroup1.ItemIndex of
0: begin
// calcolo del limite inferiore di CI di gamma50
X:= gamma50 - 1; Y:= gamma50;
Bisect(CalculateCILyman, X, Y, 10000, 1E-14, F);
MinCIGamma50:= X;
// calcolo del limite superiore di CI di gamma50
X:= gamma50; Y:= gamma50 + 1;
Bisect(CalculateCILyman, X, Y, 10000, 1E-14, F);
MaxCIGamma50:= X;
Lbl95CI.Caption:= FloatToStr(MinCIGamma50) + ' - ' + FloatToStr(MaxCIGamma50);
end;
1: begin
// calcolo del limite inferiore di CI di TCD50
X:= TCD50 - 2; Y:= TCD50;
Bisect(CalculateCILyman, X, Y, 10000, 1E-14, F);
MinCITCD50:= X;
// calcolo del limite inferiore di CI di TCD50
X:= TCD50; Y:= TCD50 + 2;
Bisect(CalculateCILyman, X, Y, 10000, 1E-14, F);
MaxCITCD50:= X;
Lbl95CI.Caption:= FloatToStr(MinCITCD50) + ' - ' + FloatToStr(MaxCITCD50);
end;
end;
Screen.Cursor:= crDefault;
end;
Gira
tutto qui. Ora devo testarlo con Niemierko e poi Poisson (sono i tre principali modelli che mi interessano perchè i più impiegati)La precisione di calcolo che ho scelto è un pò alta ( 1.0E-14) ma volevo esplorare quanto fosse "dispendioso" in termini di tempi di calcolo l'algoritmo. Ovviamente una decina di cifre decimali in meno andrebbero ugualmente bene per i miei impieghi.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo