HiI want to construct a logliikelood function in RHere is the situationy=number
of particles emitted in 1 hr period~pois(30)p=probability of detection of
radiation particlesx=number of particles detected by a radiation
detector~pois(30p)where p~beta(a,1)I have to calculate the loglikehood for a
for the range a(2,50)I wish to simulate 100 random samples for each aHere is my
code:-m=481n=100x = c(15, 36, 29, 28, 37, 32, 25, 27, 31, 21, 25, 27, 28, 31,
28, 20, 34, 25, 20, 34, 15,21, 28, 24, 31, 19, 34, 29, 18, 25, 16, 19, 44, 26,
34, 31, 21, 28, 11, 31, 21, 34, 25, 25,30, 23, 21, 35, 36, 21, 27, 29, 30, 22,
25, 30, 24, 27, 28, 22, 36, 29, 33, 35, 30, 32, 27,26, 25, 27, 23, 21, 39, 33,
24, 21, 19, 34, 32, 28, 27, 28, 23, 20, 24, 29, 21, 22, 31, 28,27, 28, 29, 21,
30, 28, 31, 22, 29, 18)a=c(1:m)*0.1+1.9p<-rbeta(n,a,1)loglik<-numeric(m)for(i
in
1:n){p[i]<-rbeta(n,a[i],1)x[i]<-rpois(n,p[i]*30)loglik<-function(x,p)logLik(a=a(i))<-n*log(a)+(log(30)*sum(x))-sum(log(factorial(x)))-30*(!
sum(p[i]))+sum(x[i]*log(p[i]))+(a[j]-1)*sum(log(p[i]))}But somehow I am not
getting the answer please HELPP.S. the data for x is a random sample of number
of particles detected in 1 hr (given in the question)Really appreciate the
helpThanks
[[alternative HTML version deleted]]
______________________________________________
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.
