[R] LogLikelihood of a Distribution Given Fixed Parameters
Rolf Turner
r.turner at auckland.ac.nz
Fri Apr 25 06:50:28 CEST 2014
As usual I am too lazy to fight my way through the rather convoluted
code presented, but it seems to me that you just want to calculate a log
likelihood. And that is bog-simple:
The log likelihood for i.i.d. data is just the sum of log f(y_i) where
the y_i are your observed values and f() is the density function of the
distribution that you have in mind.
Where there is (right) censoring you take the sum of log f(y_i) over all
the non-censored values and then add k*(1-F(cens.time)) where k is the
number of censored values and F() is the cumulative distribution
function corresponding to f().
In your case it would appear that f(y) = dlnorm(y,1.66,0.25) and
F(y) = plnorm(y,1.66,0.25). Note that instead of using 1-F(cens.time)
you can use plnorm(cens.time,1.66,0.25,lower=TRUE) and that instead of
taking logs explicitly you can set log=TRUE in the calls to dlnorm() and
plnorm().
cheers,
Rolf Turner
On 25/04/14 07:27, Jacob Warren (RIT Student) wrote:
> I'm trying to figure out if there is a way in R to get the loglikelihood of
> a distribution fit to a set of data where the parameter values are fixed.
> For example, I want to simulate data from a given alternate lognormal
> distribution and then I will fit it to a lognormal distribution with null
> parameter values to see what the likelihood of the null distribution is
> given random data from the alternate distribution.
>
> I have been using fitdistrplus for other purposes but I cannot use it to
> fix both parameter values.
>
> Here is an example of what I've been working with...
>
> nullmu<-1.66 #set null mu
> altmu<-1.58 #set alt mu
> sd.log<-0.25 #set common sigma
> cens.time<-6 #if simulated times are greater than this turn them into right
> censored times
>
> #simulating lognormal data (time) from altnative dist
> (sim<-rlnorm(n=samplesize, meanlog=altmu, sdlog=sd.log))
> #if the time was > cens.time replace time with cens.time
> (sim[which(sim>cens.time)]<-cens.time)
> sim
>
> #create a variable indicating censoring
> (cens<-sim)
> cens[which(sim==cens.time)]<-NA
> cens
>
> #create the data frame to be passed to fitdistcens and fitdist
> (x<-data.frame(left=sim,right=cens))
>
>
> #if there is censored data use fitdistcens else use fitdist
> ifelse(length(which(is.na(cens)))>0,
> simfit<-fitdistcens(censdata=x, distr="lnorm"),
> simfit<-fitdist(data=x[,1], distr="lnorm")
> )
>
> #Now I can get the loglikelihood of the MLE fitted distribution
> simfit$loglik
>
> #I want to get the loglikelihood of the distribution with the null
> parameterization
> #This is what I can't get to work
> #I can't seem to find any function that allows me to set both parameter
> values
> #so I can get to loglikelihood of the of the parameterization given the data
> nulldist<-fitdistcens(censdata=x, distr="lnorm", start=list(meanlog=nullmu,
> sdlog=sd.log)
>
> #Then I want to do a likelihood ratio test between the two distributions
> pchisq((-2*simfit$loglik--2*nulldist$loglik), df=2, lower.tail=FALSE)
More information about the R-help
mailing list