Hi everyone!
I already posted
http://r.789695.n4.nabble.com/Declaring-a-density-function-with-for-loop-td4635699.html
a question on finding density values of a new Binomial like distribution
which has the following pmf:
http://r.789695.n4.nabble.com/file/n4636134/kb.png
Thank fully
http://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=user_nodes&user=124474
Berend Hasselman and
http://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=user_nodes&user=125777
David Carlson cleared out my problems in a timely manner.
Now I have a question on estimating the parameters of this distribution
when a data set is provided.
For example, I have a data from a journal paper:
*> values<- c(0,1,2,3,4,5,6,7) #Binomial variable values
> frequency<-c(47,54,43,40,40,41,39,95) #frequency counts of the above
> binomial values, totally I have 399 samples
*
Now to estimate the parameters, I first declared the negative log
likelihood function of this new distribution as below(I used two methods,
one with a loop and one without the loop)
#without any loops
*> Loglik.newdis1<- function(x,a,b,fre,n){
KBLL1<-sum(fre*(log(a*b)+lchoose(n,x)+log(sum((-1)**(0:1000) *
choose(b-1, 0:1000) * beta(x+a+a*0:1000, n-x+1)))))
return(-KBLL1)
}*
#here a and b are the parameters to be estimated and x=binomial
values, fre=frequencies and n=binomial trials
#and since the inner series is a convergent infinite series, a large
number like 1000 is defined as the maximum
#with for loop
*>Loglik.newdis2<-function(x,a,b,fre,n,imax) {
term<-0
for (i in 0:imax) {
term=term+(((-1)**i)*(choose(b-1,i))*(beta(x+a+a*i,n-x+1)))
}
dens=a*b*choose(n,x)*term
KBLL2<-sum(fre*log(dens))
return(-KBLL2)
} *
#Now to estimate the parameters, I used
http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=bbmle:mle2 mle2(from the
package bbmle) , which is a wrapper around the optim function
*>estimates1=mle2(Loglik.newdis1, start=list(a=1,b=1),
data=list(x=values,fre=frequency, n=7, imax=2000)) *
Error in optim(par = c(1, 1), fn = function (p) :
non-finite finite-difference value [1]
In addition: There were 50 or more warnings (use warnings() to see the first
50)
/#this gives the above error message/
>estimates2=mle2(Loglik.newdis2, start=list(a=1,b=1),
data=list(x=values,fre=frequency, n=7, imax=2000))
> estimates2
Call:
mle2(minuslogl = Loglik.newdis2, start = list(a = 1, b = 1),
data = list(x = values, fre = frequency, n = 7, imax = 2000))
Coefficients:
a b
0.7574478 0.6399205
Log-likelihood: -815.48
/#here the values are not equal to the vallues provided in the journal
paper
#also thee inner series doesn't converge even for very large number
#values given in the paper are
a= 0.70 and b= 0.59 and -loglikelihood = 813.6939/
*
Please point out my mistakes and help me on this r-code?
Thank you very much!!*
--
View this message in context:
http://r.789695.n4.nabble.com/declaring-negative-log-likelihood-of-a-distribution-tp4636134.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
