Probability Density Function (PDF) di una funzione data

[nickwing]

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

[nickwing]

"cenzo":Nel grafico però non vedo il modello logistico classico (Dose ad covariate) che ha un AIC=119.07
(la funzione link è riportata su Bentzen).

Inoltre ho provato anche un altro modello - complementary log log (*) - con funzione link (espressa rispetto ai tuoi parametri):

\( \displaystyle p=1- 0.5^{ e^{ \frac{-2 \gamma_{50}}{\log2} \left( 1-\frac{\text{EUD}}{\text{TCD}_{50}} \right) } } \)

Con questo risulta TCD50=40.012, γ50=2.408, log(L)=−55.287 e AIC=114.57 (il più basso finora)

(*) La sua espressione originaria è \( \log(-\log(1-p))=a+bx \)

Ho messo solo Log-Dose (Niemierko), la curva nera
Tocca che faccio una prova pure col complementary log come dicevi tu. Nel frattempo mi sono attrezzato per il calcolo dell'intervallo di confidenza di tali valori individuati ($gamma_50$ e $TCD50$). Ho seguito il metodo descritto da Zhu e riportato in dettaglio nelle appendici di questo articolo di Dawson. Ossia utilizzando il cosiddetto "profile likelihood method" sottraggo al valore di MLE 3.84/2 come spiegato qui:

The uncertainty in the fitted parameters was assessed by the profile likelihood method. We varied each parameter separately around its optimum value and determined the lower and upper bound value of the parameter at which the natural log of the likelihood function dropped by 3.84/2. The parameter value was varied while the likelihood function was maximized with the two other parameters. The value of 3.84 corresponds to a confidence level of 95% with 1 degree of freedom for a $chi^2$ distribution

.
Prima di implementare il tutto nel programma principale ho scritto un programmino "ad hoc" per calcolare la MLE dei due parametri famosi leggendo la nostra ormai arcinota serie di dati simulati: nei due screenshot che ti posto ho eseguito il calcolo dei parametri $TCD50$ e $gamma_50$ per il modello di Lyman (quello che deriva dalla error function ). I risultati sono i due intervalli di confidenza calcolati per i due parametri riportando l'altro parametro al suo valore "ottimale". Per calcolarli ho provato prima con la soluzione numerica di un'equazione non lineare ma con gli algoritmi in mio possesso ottenevo una serie di errori di programma, alla fine mi sono deciso a scrivere un algoritmo iterativo molto semplice che cerca i valori estremi dell'intervallo di confidenza muovendo il valore della variabile da analizzare in alto e in basso e ponendo l'uguaglianza del risultato con quello calcolato come $MLE - 3.84/2$ che trovi anche riportato negli screenshot:





Secondo te il metodo gira in questi termini ?

[cenzo1]

"nickwing":[quote="cenzo"]Nel grafico però non vedo il modello logistico classico (Dose ad covariate) che ha un AIC=119.07
(la funzione link è riportata su Bentzen).

Sarebbe la curva nera (Niemierko)
Tocca che faccio una prova pure col complementary log come dicevi tu.[/quote]
La curva nera Niemierko dovrebbe essere il logit in cui si assume come covariate il logaritmo della dose (la funzione link è quella della tabella 1 di Bentzen, riga logit, colonna logarithm of dose as covariate, in pratica la prima formula che hai riportato nel topic). AIC = 126.01. Giusto ?
[EDIT]

Il logit classico invece, stessa tabella Bentzen, stessa riga Logit ma prima colonna (Dose as covariate). Trovi la funzione link:
\( p=\frac{1}{1+e^{4 \gamma_{50} \left( 1- \frac{EUD}{TCD_{50}} \right) } } \) con la quale hai AIC = 119.07.

Sarebbe molto interessante verificare se la cloglog (AIC = 114.57) si comporta bene anche con altri set di dati.
Finora non mi sembra di averla vista utilizzata in altri articoli...

Sto anche lavorando ad un altro modello ancora (cauchit) che sembra ancora migliore avendo AIC = 110.29.
Appena riesco a ricostruire la funzione link parametrizzata su TCD50 e gamma50 la posto.

Per gli intervalli di confidenza ho sempre lasciato fare tutto al programma R (di cui mi fido al 99.999%).
Leggo con calma il tuo intervento e i riferimenti, poi ti faccio sapere.

[nickwing]

Hai documentazione su questo modello Cauchit? Da quel poco che sono riuscito a trovare ho visto che potrebbe essere molto utile in caso si abbiano a disposizione eventi posizionati agli estremi dell'intervallo osservato (per noi, della dose EUD), cosa non infrequente (molte osservazioni e pochi eventi positivi all'estremo più alto della casistica) nella pratica clinica reale. Tra l'altro non è stato mai implementato, da quel che ho visto, in radiobiologia, se lo utilizzassimo per analizzare qualche casistica reale ne potremmo ricavare un bel lavoro pilota da pubblicare

[cenzo1]

Con ordine:

1) Sorry, sulla Cauchit non ti so dare riferimenti. L'ho scoperta grazie ad R. Tra le funzioni link di default, oltre alla logit e probit, figurano la cloglog e la cauchit. Ho pensato allora che fossero in un certo senso importanti e ho provato per curiosità sui nostri dati e ho visto che avevano un AIC più basso.
La cauchit si basa sulla CDF della distribuzione di Cauchy (come la probit sulla CDF normale).
La funzione link parametrizzata su tcd50 e gamma50 è questa:

\( \displaystyle p=\frac{1}{2} + \frac{1}{\pi} \arctan \left[ \pi \ \gamma_{50} \left( \frac {\text{EUD}}{\text{TCD}_{50}}-1 \right) \right] \)

Mi risulta: tcd50=39.564 ; gamma50=7.002 ; logL=-53.144 ; AIC=110.29
Prova a disegnare il grafico: è molto diverso dagli altri, molto più ripido per i valori centrali di EUD (infatti ha un gamma50 molto alto).
Che ne pensi ?
Il grafico della cloglog (AIC=114.57) invece è in linea con gli altri.

2) Nell'immagine che hai riportato più sopra leggo, per il modello probit (Lyman), una gamma50=2.898.
Nick, non mi torna. Io ottengo gamma50=2.049.
Ho provato sia con la formula di Zhu (CDF normale) sia con la formula di Bentzen (uguale alla (17) di Cancers, erf function). Sempre 2.049.
Credo che l'errore è su Cancers quando riporta $m=1/(\gamma_{50} sqrt(\pi))$ oppure $\gamma_{50}=1/(m sqrt(\pi))$
Sono entrambe errate. A me torna invece:
$m=1/(\gamma_{50} sqrt(2*\pi))$ oppure $\gamma_{50}=1/(m sqrt(2*\pi))$

3) Per gli intervalli di confidenza col metodo del profile likelihood ti consiglio di usare, invece di 3.84 un valore con più cifre significative:
\( \chi_{0.95,1}^2= 3.841459\)
Forse è questo che giustifica la piccola differenza che ottengo sul CI di tcd50: 36.769 - 40.747 (contro il tuo 36.783 - 40.725)
Il metodo mi sembra giusto, anzi: sei stato bravissimo a scriverti un algoritmo iterativo ad hoc. Sempre in Delphi ?

Confesso che ho problemi a calcolare il CI. Per il momento ho scritto una macro in Excel che lancia il solver ...
Volevo sfruttare i risultati del CI di R. Il problema è che in R la parametrizzazione dei modelli è diversa e non so come riportare il CI dei parametri che usa R (diciamo a e b di un modello lineare a+bx) in termini di gamma50 e tcd50.
Ho la relazione che lega (a,b) con (tcd50,gamma50) ma al momento non so come convertire i CI...

PS: ho provato anche il modello Weibull (vedi Cancers): mi risulta AIC=117.51 (sui "nostri" dati)

[cenzo1]

Dimenticavo: sono curioso di dare uno sguardo all'articolo di Dawson (e all'appendice), ma non ho accesso a Science Direct...

[nickwing]

Ho lasciato i sorgenti modificati in ospedale. Domani e Ferragosto e tornerò lì mercoledì perchè me ne vò tre giorni al mare
Ho problemi a riscrivere l'algoritmo iterativo per il modello di Lyman, becco sempre un errore nel calcolo della funzione... Sarà pure che ho avuto una settimana lunga...
Ti posto il codice dell'algoritmo per Niemierko, ossia il Logit Questo funziona per il calcolo degli estremi del C.I.

procedure TForm1.Button3Click(Sender: TObject);
var
  m, n, p: Integer; // contatori
  CIup, CIdown: Float;			// estremi dell'intervallo di confidenza
  TempUP, TempDOWN: Float;	       // termini temporanei per il calcolo
  TempMLE: Float;						// valore temporaneo di MLE
begin
	LblMLE.Caption:= FloatToStr(MLE_N);
  LblPrMLE.Caption:= FloatToStr(MLE_N - 1.9207295);
	// calcolo intervallo di confidenza nel modello di Niemierko
  case Form1.RadioGroup1.ItemIndex of
  	0: begin  // gamma50
      m:= 0; p:= 1;
      TempUP:= gamma50;
      // calcolo valore basso di CI di gamma50
      while m < 14 do begin
      	TempDOWN:= TempUP - Power(10, -m)*p;
        TempMLE:= 0;
        for n:= 1 to 200 do begin TempMLE:= TempMLE + DataList[n].Outcome * Ln(1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))) +
                (1 - DataList[n].Outcome) * (Ln(1 - 1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))));
        end;
        if TempMLE < MLE_N - 1.9207295 then begin
          inc(m);
          p:= 1;
        end else begin
          inc(p);
          TempUP:= TempDOWN;
        end;
      end;
      CIdown:= RoundN((TempUP + TempDOWN)/2, 13);
      // caloolo valore alto di CI gamma50
      m:= 0; p:= 1;
      TempUP:= gamma50;
      while m < 14 do begin
      	TempDOWN:= TempUP + Power(10, -m)*p;
        TempMLE:= 0;
        for n:= 1 to 200 do begin TempMLE:= TempMLE + DataList[n].Outcome * Ln(1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))) +
                (1 - DataList[n].Outcome) * (Ln(1 - 1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))));
        end;
        if TempMLE < MLE_N - 1.9207295 then begin
          inc(m);
          p:= 1;
        end else begin
          inc(p);
          TempUP:= TempDOWN;
        end;
      end;
      CIup:= RoundN((TempUP + TempDOWN)/2, 13);
      // aggiunge gli elementi alla string sul CI
      Lbl95CI.Caption:= FloatToStr(CIdown) + ' - ' + FloatToStr(CIup);
    end;

    1: begin  // TCD50
      m:= 0; p:= 1;
      TempUP:= TCD50;
      // calcolo valore basso di CI di TCD50
        TempDOWN:= TempUP - Power(10, -m)*p;
        TempMLE:= 0;
        for n:= 1 to 200 do begin TempMLE:= TempMLE + DataList[n].Outcome * Ln(1/(1+ Power(TempDown/DataList[n].EUD, 4*gamma50))) +
                (1 - DataList[n].Outcome) * (Ln(1 - 1/(1+ Power(TempDOWN/DataList[n].EUD, 4*gamma50))));
        end;
        if TempMLE < MLE_N - 1.9207295 then begin
          inc(m);
          p:= 1;
        end else begin
          inc(p);
          TempUP:= TempDOWN;
        end;
      end;
      CIdown:= RoundN((TempUP + TempDOWN)/2, 13);
      // caloolo valore alto di CI gamma50
      m:= 0; p:= 1;
      TempUP:= TCD50;
      // calcolo valore basso di CI di TCD50
      while m < 14 do begin
      	TempDOWN:= TempUP + Power(10, -m)*p;
        TempMLE:= 0;
        for n:= 1 to 200 do begin TempMLE:= TempMLE + DataList[n].Outcome * Ln(1/(1+ Power(TempDown/DataList[n].EUD, 4*gamma50))) +
                (1 - DataList[n].Outcome) * (Ln(1 - 1/(1+ Power(TempDown/DataList[n].EUD, 4*gamma50))));
        end;
        if TempMLE < MLE_N - 1.9207295 then begin
          inc(m);
          p:= 1;
        end else begin
          inc(p);
          TempUP:= TempDOWN;
        end;
      end;
      CIup:= RoundN((TempUP + TempDOWN)/2, 13);
      // aggiunge gli elementi alla string sul CI
      Lbl95CI.Caption:= FloatToStr(CIdown) + ' - ' + FloatToStr(CIup);
    end;
  end;
end;

non è commentato in ogni dettaglio ma credo tu ti possa orientare

[nickwing]

Allora ci ho lavorato un pò su ancora e con il valore identificato ci siamo, con l'intervallo di confidenza ancora no... boh?!? Mi sa che devo dare una ritoccata all'algoritmo

[nickwing]

l'articolo di Dawson te lo puoi scaricare da qui...

[nickwing]

ahhh, alla fine ho provato anche Cauchit, effettivamente ha una forma... un pò strana ma trovo i tuoi stessi valori

[cenzo1]

Ho scritto una tabella riassuntiva con le formule relative ai 10 modelli analizzati (la puoi scaricare qui).

Questi sono:
1) Logit (Suit, Bentzen)
2) log Logit (Niemierko)
3) Probit (Lyman)
4) log Probit
5) Poisson (Munro Gilbert, Kallman...)
6) log Poisson
7) Cloglog
8) Weibull o log Cloglog (Cancers)
9) Cauchit
10) log Cauchit

Ogni modello dei 5 ha un suo "gemello" in cui compare il logaritmo della dose come variabile esplicativa (come fa Bentzen).
Su Bentzen ci sono 1,2,3,4,5,6. Su Cancers compare anche il Weibull (equivalente a un log Cloglog).
La colonna "link function" secondo me fa cogliere molto bene le relazioni tra i gemelli e tra i vari modelli in generale, molto più di quanto si capisca con la parametrizzazione D50, gamma50.

Il Cloglog e i due Cauchit non li ho visti negli articoli che mi hai passato; potenzialmente quindi sarebbero una novità.. tra l'altro sono i tre che presentano il migliore adattamento (AIC) sulla base dati inventata..

Possibili critiche ai due Cauchit (che hanno grafici quasi sovrapposti):
- gamma50 molto alto: presentano un brusco salto nella probabilità per un intervallo ristretto di EUD (questo però può anche essere un pregio);
- convergono lentamente a p=0 per EUD bassa e a p=1 per EUD alta (rispetto agli altri).

[cenzo1]

"nickwing":l'articolo di Dawson te lo puoi scaricare da ...

Grazie!
Speravo di trovare qualche indicazione su come "trasferire" i CI da (a,b) ai parametri D50,gamma50.
Purtroppo non ho trovato quello che speravo.

Come ti dicevo prima ho trovato i CI per D50 e gamma50 con una macro Excel...
Non mi fido moltissimo dei risultati. La terza cifra decimale può essere non accurata di una unità.
In qualche raro caso (non mi ricordo più in quale modello in particolare, forse il Cauchit) anche la seconda.
(Dico questo perchè ho calcolato con Excel anche i CI sui param a,b per confrontare con i risultati di R di cui mi fido).

I risultati sono in questa immagine:



Una osservazione sui modelli: tra ogni coppia dei "gemelli" il modello in cui figura il logaritmo della dose si comporta sempre peggio del rispettivo che dipende dalla dose (dose as covariate). Almeno su questa base dati.

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...

[cenzo1]

"nickwing":Ti posto il codice dell'algoritmo per Niemierko, ossia il Logit Questo funziona per il calcolo degli estremi del C.I.

      m:= 0; p:= 1;
      TempUP:= gamma50;
      // calcolo valore basso di CI di gamma50
      while m < 14 do begin
        TempDOWN:= TempUP - Power(10, -m)*p;
        TempMLE:= 0;
        for n:= 1 to 200 do begin TempMLE:= TempMLE + DataList[n].Outcome * Ln(1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))) +
                (1 - DataList[n].Outcome) * (Ln(1 - 1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))));
        end;
        if TempMLE < MLE_N - 1.9207295 then begin
          inc(m);
          p:= 1;
        end else begin
          inc(p);
          TempUP:= TempDOWN;
        end;
      end;
      CIdown:= RoundN((TempUP + TempDOWN)/2, 13);

Il codice non mi convince. Secondo me non implementa il metodo del profile Likelihood, per come l'ho capito.
La cosa essenziale è che la logL1 che calcoli tiene sempre fisso il parametro D50; invece dovresti massimizzare la logL1 variando solo D50 e tenendo fisso gamma50=TempDOWN nello step corrente.
Al termine del tuo algoritmo, è vero che è soddisfatta la relazione MLE - logL1- 1.92 = 0, però non sei in condizioni da avere massimizzato la logL1 (D50 deve cambiare).

Se ti può aiutare, la macro che ho fatto in Excel lavora con questa logica:
(stiamo cercando l'estremo inferiore del CI per gamma 50)

1) Fisso i valori iniziali: D50i=D50 (quello ottimo da MLE); gamma50i=un valore un po' inferiore al gamma50 stimato con MLE
2) massimizzo la logL1 tenendo fisso gamma50i e variando solo D50i
3) do while ( abs (MLE - logL1- 1.92..) > err) ''' err = 0.000000001
4) aggiorno il valore di gamma50i tale che (MLE - logL1- 1.92 = 0) lasciando fisso D50i (sempre col solver)
5) massimizzo la logL1 tenendo fisso gamma50i e variando solo D50i
6) loop

In pratica per come l'ho impostato si lavora su due livelli: massimizzazione della logL1 funzione solo di D50i e trovare un gamma50i che soddisfa la relazione MLE - logL1- 1.92 = 0 a parità di D50i.
Ovviamente le cose sono legate: componendone una, si scompone l'altra. Perciò il loop. L'algoritmo mi risulta sempre convergere.
Esco con gamma50lower=gamma50i ultimo calcolato.

Ho avuto qualche problema con alcuni modelli per via che, con probabilità molto piccole o molto alte, viene p=0 oppure 1-p=0 e giustamente Excel mi da un errore per calcolare la logL1 (c'è il logaritmo di 0 ).
Non entro nei dettagli per ovviare al problema ...

Sarebbe molto utile se trovassi una libreria come per la massimizzazione!

[nickwing]

"cenzo":[quote="nickwing"]Ti posto il codice dell'algoritmo per Niemierko, ossia il Logit Questo funziona per il calcolo degli estremi del C.I.

      m:= 0; p:= 1;
      TempUP:= gamma50;
      // calcolo valore basso di CI di gamma50
      while m < 14 do begin
        TempDOWN:= TempUP - Power(10, -m)*p;
        TempMLE:= 0;
        for n:= 1 to 200 do begin TempMLE:= TempMLE + DataList[n].Outcome * Ln(1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))) +
                (1 - DataList[n].Outcome) * (Ln(1 - 1/(1+ Power(TCD50/DataList[n].EUD, 4*TempDOWN))));
        end;
        if TempMLE < MLE_N - 1.9207295 then begin
          inc(m);
          p:= 1;
        end else begin
          inc(p);
          TempUP:= TempDOWN;
        end;
      end;
      CIdown:= RoundN((TempUP + TempDOWN)/2, 13);

Il codice non mi convince. Secondo me non implementa il metodo del profile Likelihood, per come l'ho capito.
La cosa essenziale è che la logL1 che calcoli tiene sempre fisso il parametro D50; invece dovresti massimizzare la logL1 variando solo D50 e tenendo fisso gamma50=TempDOWN nello step corrente.
Al termine del tuo algoritmo, è vero che è soddisfatta la relazione MLE - logL1- 1.92 = 0, però non sei in condizioni da avere massimizzato la logL1 (D50 deve cambiare).[/quote]

Ho capito. In effetti io pensavo che si trattasse solamente di trovare l'intersezione fra una retta orizzontale in cui $y=MLE - 1.92$, e per fare questo, non essendo riuscito a trovare una soluzione algebrica all'equazione ho fatto la ricerca in maniera iterativa, ossia creando degli intervalli e "stringendo" i valori, prima di gamma50 e poi di D50, intorno al valore dato, mantenendo però fisso il valore dell'altra costante nel modello.

"cenzo":Se ti può aiutare, la macro che ho fatto in Excel lavora con questa logica:
(stiamo cercando l'estremo inferiore del CI per gamma 50)

1) Fisso i valori iniziali: D50i=D50 (quello ottimo da MLE); gamma50i=un valore un po' inferiore al gamma50 stimato con MLE
2) massimizzo la logL1 tenendo fisso gamma50i e variando solo D50i
3) do while ( abs (MLE - logL1- 1.92..) > err) ''' err = 0.000000001
4) aggiorno il valore di gamma50i tale che (MLE - logL1- 1.92 = 0) lasciando fisso D50i (sempre col solver)
5) massimizzo la logL1 tenendo fisso gamma50i e variando solo D50i
6) loop

In pratica per come l'ho impostato si lavora su due livelli: massimizzazione della logL1 funzione solo di D50i e trovare un gamma50i che soddisfa la relazione MLE - logL1- 1.92 = 0 a parità di D50i.
Ovviamente le cose sono legate: componendone una, si scompone l'altra. Perciò il loop. L'algoritmo mi risulta sempre convergere.
Esco con gamma50lower=gamma50i ultimo calcolato.

Posso provare ad implementarlo in Delphi . Mi serve un pò di lavoro però

"cenzo":Sarebbe molto utile se trovassi una libreria come per la massimizzazione!

La libreria ce l'ho, è quella che ho utilizzato per calcolare la MLE, in realtà funziona cercando il minimo della funzione con il metodo Broyden-Fletcher-Goldfarb-Shanno, altrimenti detto Quasi-Newton. Lo trovi descritto a pagina 56, paragrafo 7.2.4 di questo file, che è la documentazione della libreria che ho usato. Ovviamente nel mio programma faccio calcolare il minimo dell'inverso della funzione MLE (a questo ci arrivo con le mie conoscenze matematiche ).

[nickwing]

un'altra domanda, da ignorante quale sono: ma come hai fatto a tradurre la forma "canonica" dei vari modelli nella forma che contiene i parametri $D_50$ e $gamma_50$? Nel senso, che equazioni poni per trasformare i termini $a$ e $b$ della forma classica dei diversi modelli nella forma contenete solo le costanti che ci interessano?

[cenzo1]

"nickwing":come hai fatto a tradurre la forma "canonica" dei vari modelli nella forma che contiene i parametri $D_50$ e $gamma_50$?

Ho utilizzato il significato geometrico di $D_50$ e $gamma_50$:
$p(x=D_50)=0.50$
$(dp)/(dx) (x=D_50)=\gamma_{50} / D_{50}$

Prendiamo ad esempio il modello log Logit (Niemierko): $"logit " p= ln(p/(1-p))=a+blnx->p/(1-p)=e^(a+blnx)->p=e^(a+blnx)/(1+e^(a+blnx))=1/(1+e^(-(a+blnx)))$
Prima equazione: $ln(0.5/(1-0.5))=a+bln D_50 ->a+bln D_50 =0->ln D_50=-a/b->D_50=e^(-a/b)$

Seconda equazione: $(dp)/(dx) =(-e^(-(a+blnx))*(-b/x))/(1+e^(-(a+blnx)))^2->(dp)/(dx) (x=D_50)=(-e^(-(a+blnD_50))*(-b/D_50))/(1+e^(-(a+blnD_50)))^2=(-e^(-0)*(-b/D_50))/(1+e^(-0))^2=b/(4D_50)$
Dovendo essere appunto per la seconda equazione: $(dp)/(dx) (x=D_50)=b/(4D_50)=\gamma_{50} / D_{50}->b=4gamma_50->gamma_50=b/4$

Riassumendo, ho ottenuto le seguenti equazioni:
$D_50=e^(-a/b)$
$gamma_50=b/4$

Ne calcolo le inverse:
$b=4gamma_50$
$a=-blnD_50=-4gamma_50*lnD_50$

Infine sostituisco queste espressioni nella funzione $p=1/(1+e^(-(a+blnx)))$:
$p=1/(1+e^(-(-4gamma_50*lnD_50+4gamma_50 lnx)))=1/(1+e^(-(-4gamma_50(lnD_50-lnx))))=1/(1+e^(4gamma_50 ln(D_50/x)))=1/( 1+e^( ln (D_50/x)^(4 gamma_50)) )=1/( 1+(D_50/x)^(4 gamma_50) )$
che è proprio l'espressione a noi familiare del modello Niemierko

In modo del tutto analogo si procede per gli altri modelli.

[nickwing]

elementare... non ci avevo pensato

[cenzo1]

"nickwing":[quote="cenzo"]Sarebbe molto utile se trovassi una libreria come per la massimizzazione!

La libreria ce l'ho, è quella che ho utilizzato per calcolare la MLE, in realtà funziona cercando il minimo della funzione con il metodo Broyden-Fletcher-Goldfarb-Shanno, altrimenti detto Quasi-Newton. Lo trovi descritto a pagina 56, paragrafo 7.2.4 di questo file, che è la documentazione della libreria che ho usato.[/quote]
Mi riferivo ad una libreria che implementa direttamente il metodo del profile likelihood. In effetti penso sia difficile esista qualcosa di specifico per Delphi. Ma sono convinto che in qualche modo si possa risolvere con la libreria che hai.

Per caso hai modo di scaricare (e condividere ) questo articolo ?

"nickwing":Ovviamente nel mio programma faccio calcolare il minimo dell'inverso della funzione MLE

Volevi dire l'opposto -MLE ?

[cenzo1]

"cenzo":Per caso hai modo di scaricare (e condividere ) questo articolo ?

Come non detto. Sono riuscito ad averlo.. però l'approccio mi sembra troppo complicato..

Invece il metodo suggerito in questo breve articolo mi sembra più facilmente implementabile. Dagli uno sguardo.

[nickwing]

"cenzo":[quote="nickwing"]Ovviamente nel mio programma faccio calcolare il minimo dell'inverso della funzione MLE

Volevi dire l'opposto -MLE ? [/quote]

ehhmmm... si... "le parole sono importanti" (cit. N. Moretti)

[nickwing]

A proposito: mandami la mail tramite PM, così se abbiamo altri articoli e file da scambiare evitiamo di metterli su un server pubblico tipo zshare...