R help - Jun 2009 - implementing Maximum Likelihood with distrMod when only the PDF is known

Dear R users and Dear authors of the distr package and sequels

I am trying to use the (very nice) package distrMod as I want to 
implement maximum likelihood (ML) fit of some univariate data for which 
I have derived a theoretical continuous density (pdf). As it is a 
parametric density, I guess that I should implement myself a new 
distribution of class AbscontDistributions (as stated in the pdf on 
"creating new distributions in distr"), and then use MLEstimator()
from
the distrMod package. Is that correct or is there a simpler way to go? 
Note that I want to use the distr package because it allows me to 
implement simply the convolution of my theoretical pdf with some noise 
distribution that I want to model in the data, this is more difficult 
with fitdistr or mle.
It proved rather difficult for me to implement the new class following 
all the steps provided in the example, so I am asking if someone has an 
example of code he wrote to implement a parametric distribution from its 
pdf alone and then used it with MLEstimator().

I am sorry for the post is a bit long but it is a complicate question to 
me and I am not at all skillful in the handling of such notions as "S4 - 
class", etc.. so I am a bit lost here..

As a simple example, suppose my theoretical pdf is the skew normal 
distribution (available in package sn):

#skew normal pdf (default values = the standard normal N(0,1)

fsn<-function(x,mu=0,sd=1,d=0) {u = (x-mu)/sd;  f = dnorm(u)*pnorm(d*u); 
return(f/sd)}

# d = shape parameter (any real), mu = location (any real), sd = scale 
(positive real)

#to see what it looks like try
x<-seq(-1,4,length=200);plot(fsn(x,d=3),type="l")

#Now I tried to create the classes "SkewNorm" and
"SkewNormParameter"
copying the example for the binomial
##Class:parameters
setClass("SkewNormParameter",
representation=representation(mu="numeric",sd="numeric",d="numeric"),
prototype=prototype(mu=0,sd=1,d=0,name=gettext("Parameter of the Skew 
Normal distribution")),
contains="Parameter"
)

##Class: distribution (created using the pdf of the skew normal defined 
above)
setClass("SkewNorm",prototype = prototype(
     d = function(x, log = FALSE){fsn(x, mu=0, sd=1,d=0)},
     param = new("SkewNormParameter"),
     .logExact = TRUE,.lowerExact = TRUE),
contains = "AbscontDistribution"
)

#so far so good but then with
setMethod("mu", "SkewNormParameter", function(object) object
at mu)

#I get the following error message:

 > Error in setMethod("mu", "SkewNormParameter",
function(object)
object at mu) : no existing definition for function "mu"

I don't understand because to me mu is a parameter not a function... 
maybe that is too complex programming for me and I should switch to 
implementing my likelihood by hand with numerical convolutions and 
optim() etc., but I would like to know how to use distr, so if there is 
anyone who had the same problem and solved it, I would be very grateful 
for the hint !

All the best,
Guillaume



-- 
Guillaume MARTIN
Institut des Sciences de l'Evolution de Montpellier
ISEM UMR 5554, B?t. 22 Universit? Montpellier II
34 090 Montpellier, France

tel: (+33) 4 67 14 32 50

Dear Guillaume,

thanks for your interest in the distrMod package.

Regarding your question I took up your example and put a file under:

http://www.stamats.de/distrModExample.R

Hope that helps ...

Don't hesitate to contact me if you have further questions!

Best,
Matthias

guillaume.martin schrieb:
> Dear R users and Dear authors of the distr package and sequels
>
> I am trying to use the (very nice) package distrMod as I want to 
> implement maximum likelihood (ML) fit of some univariate data for 
> which I have derived a theoretical continuous density (pdf). As it is 
> a parametric density, I guess that I should implement myself a new 
> distribution of class AbscontDistributions (as stated in the pdf on 
> "creating new distributions in distr"), and then use
MLEstimator()
> from the distrMod package. Is that correct or is there a simpler way 
> to go? Note that I want to use the distr package because it allows me 
> to implement simply the convolution of my theoretical pdf with some 
> noise distribution that I want to model in the data, this is more 
> difficult with fitdistr or mle.
> It proved rather difficult for me to implement the new class following 
> all the steps provided in the example, so I am asking if someone has 
> an example of code he wrote to implement a parametric distribution 
> from its pdf alone and then used it with MLEstimator().
>
> I am sorry for the post is a bit long but it is a complicate question 
> to me and I am not at all skillful in the handling of such notions as 
> "S4 - class", etc.. so I am a bit lost here..
>
> As a simple example, suppose my theoretical pdf is the skew normal 
> distribution (available in package sn):
>
> #skew normal pdf (default values = the standard normal N(0,1)
>
> fsn<-function(x,mu=0,sd=1,d=0) {u = (x-mu)/sd;  f = 
> dnorm(u)*pnorm(d*u); return(f/sd)}
>
> # d = shape parameter (any real), mu = location (any real), sd = scale 
> (positive real)
>
> #to see what it looks like try
> x<-seq(-1,4,length=200);plot(fsn(x,d=3),type="l")
>
> #Now I tried to create the classes "SkewNorm" and
"SkewNormParameter"
> copying the example for the binomial
> ##Class:parameters
> setClass("SkewNormParameter",
>
representation=representation(mu="numeric",sd="numeric",d="numeric"),
> prototype=prototype(mu=0,sd=1,d=0,name=gettext("Parameter of the Skew 
> Normal distribution")),
> contains="Parameter"
> )
>
> ##Class: distribution (created using the pdf of the skew normal 
> defined above)
> setClass("SkewNorm",prototype = prototype(
>     d = function(x, log = FALSE){fsn(x, mu=0, sd=1,d=0)},
>     param = new("SkewNormParameter"),
>     .logExact = TRUE,.lowerExact = TRUE),
> contains = "AbscontDistribution"
> )
>
> #so far so good but then with
> setMethod("mu", "SkewNormParameter", function(object)
object at mu)
>
> #I get the following error message:
>
> > Error in setMethod("mu", "SkewNormParameter",
function(object)
> object at mu) : no existing definition for function "mu"
>
> I don't understand because to me mu is a parameter not a function... 
> maybe that is too complex programming for me and I should switch to 
> implementing my likelihood by hand with numerical convolutions and 
> optim() etc., but I would like to know how to use distr, so if there 
> is anyone who had the same problem and solved it, I would be very 
> grateful for the hint !
>
> All the best,
> Guillaume
>
>
>
-- 
Dr. Matthias Kohl
www.stamats.de