[R] LogLikelihood of a Distribution Given Fixed Parameters

Rolf Turner r.turner at auckland.ac.nz
Tue Apr 29 02:07:40 CEST 2014

Indeed it should be lower=FALSE.  Duhhhh!  Sorry 'bout that!



On 29/04/14 02:51, Jacob Warren (RIT Student) wrote:
> Thanks Rolf. That took care of it. It should be lower=FALSE through
> right? I want the upper tail because my values are right censored?
> Regards,
> Jake
> On Fri, Apr 25, 2014 at 12:50 AM, Rolf Turner <r.turner at auckland.ac.nz
> <mailto:r.turner at auckland.ac.nz>> wrote:
>     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 <http://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