search for: loglikelihood

Dear R experts, Attached you can find the expression of a loglikelihood function which I would like to maximize in R. So far, I have done maximization with the combined use of the mathematical programming language AMPL (www.ampl.com) and the solver SNOPT (http://www.sbsi-sol-optimize.com/manuals/SNOPT%20Manual.pdf). With these tools, maximization is carried out in a fe...

It has always been my understanding that deviance for GLMs is defined by; D = -2(loglikelihood(model) - loglikelihood(saturated model)) and this can be calculated by (or at least usually is); D = -2(loglikelihood(model)) As is done so in the code for 'polr' by Brian Ripley (in the package 'MASS') where the -loglikehood is minimised using optim; res <- optim...

Hi Rers, I have a silly question. I don't know how to express the loglikelihood function of 1/(x!) where x=x1,x2,....xn in R. Could anyone give me a hint? Thank you in advance. Chunhao Tu

Dear R-help, The documentation for systemfit shows that logLik() can be used to obtain loglikelihood values from linear systems estimated by systemfit(). It seems to me that logLik() cannot be used for nlsystemfit(). Does anyone know of any other packages that might let me obtain the loglikelihood of a model estimated with nlsystemfit()? Kind regards, Alex Olssen

I'm trying to estimate the parameters for GARCH(1,1) process. Here's my code: loglikelihood <-function(theta) { h=((r[1]-theta[1])^2) p=0 for (t in 2:length(r)) { h=c(h,theta[2]+theta[3]*((r[t-1]-theta[1])^2)+theta[4]*h[t-1]) p=c(p,dnorm(r[t],theta[1],sqrt(h[t]),log=TRUE)) } -sum(p) } Then I use nlminb to minimize the function loglikelihood: nlminb( c(0.001,0.001,0.001,0.001),...

Dear R users, I am estimating Poisson mixed models using glmmPQL (MASS) and lmer (lme4). We know that glmmPQL do not provide the correct loglikelihood for such models (it gives the loglike of a 'pseudo' or working linear mixed model). I would like to know how the loglike is calculated by lmer. A minor question is: why do glmmPQL and lmer give different degrees-of-freedom for the same estimated model? Does glmmPQL consider the scale parame...

...nction Chi-square 918.852 Degrees of freedom 36 P-value 0.000 Full model versus baseline model: Comparative Fit Index (CFI) 0.931 Tucker-Lewis Index (TLI) 0.896 Loglikelihood and Information Criteria: Loglikelihood user model (H0) -3737.745 Loglikelihood unrestricted model (H1) -3695.092 Number of free parameters 21 Akaike (AIC) 7517.490 Bayesian (BIC) 7595....

Hello, I am trying to understand a small quirk I came across in R. The following code results in an error: k <- c(2, 1, 1, 5, 5) f <- c(1, 1, 1, 3, 2) loglikelihood <- function(theta,k,f){ if( theta<1 && theta>0 ) return(-1*sum(log(choose(k,f))+f*log(theta)+(k-f)*log(1-theta))) return(NA) } nlm(loglikelihood ,0.5, k, f ) Running this code results in: Error in nlm(logliklihood, 0.5, k, f) : invalid function value in 'nlm&...

...thods (Nelder-Mead and BFGS) will need to evaluate the integrals, my second question is: why this difference? Below, I am providing the code that I use. Any help will be greatly appreciated. Deepankar ************ CODE START ******************* ############################# # COMPUTING THE LOGLIKELIHOOD # USING NUMERICAL INTEGRALS ############################# LLK <- function(alpha, y) { n <- length(y) lglik <- numeric(n) # TO BE SUMMED LATER TO GET THE LOGLIKELIHOOD lambda <- numeric(n-1) # GENERATING *lstar* for (i in 1:(n-1)) { # TO USE IN THE...

Hi,  I tried fitting loglinear model using the glm(catspec). The data used is FHtab. . An independence model was fitted. Here summary() and fitmacro( ) give different values for AIC.   I understand that fitmacro( ) takes the likelilhood ratio L2(deviance) to calculate AIC and uses the formula AIC= L2- d.f(deviance)*2 and this AIC is used for comparison of nested models. (Am I right?)   The value

Hi, I tried to compute maximum likelihood under gamma distribution, using nlm function. The code is this: __BEGIN__ vsamples<- c(103.9, 88.5, 242.9, 206.6, 175.7, 164.4) mlogl <- function(alpha, x) { if (length(alpha) > 1) stop("alpha must be scalar") if (alpha <= 0) stop("alpha must be positive") return(- sum(dgamma(x, shape = alpha, log = TRUE)))

Hi, I am trying to learn lognormal mixture models with EM. I was wondering how does one compute the log likelihood. The current implementation I have is as follows, which perform really bad in learning the mixture models. __BEGIN__ # compute probably density of lognormal. dens <- function(lambda, theta, k){ temp<-NULL meanl=theta[1:k] sdl=theta[(k+1):(2*k)]

Hi - I have a likelihood function that involves sums of two possions: L = a*dpois(Xi,theta1)*dpois(Yi,theta2)+b*(1-c)*a*dpois(Xi,theta1+theta3)*dpois(Yi,theta2) where a,b,c,theta1,theta2,theta3 are parameters to be estimated. (Xi,Yi) are observations. However, Xi and Yi are usually big (> 20000). This causes dpois to returns 0 depending on values of theta1, theta2 and theta3. My first

...lem without VGAM, using dnbinom in the wrapper instead. library(VGAM) phi <- 0.09 size <- 0.7 munb <- 72 x <- rzinb(200, phi=phi, size=size, munb=munb) #VGAM can fit using vglm... but I'd like to try optim > fit <- vglm(x~1, zinegbinomial,trace=TRUE) VGLM linear loop 1 : loglikelihood = -964.1915 VGLM linear loop 2 : loglikelihood = -964.1392 VGLM linear loop 3 : loglikelihood = -964.1392 > Coef(fit) phi munb k 0.1200317 62.8882874 0.8179183 > #build my wrapper function for dzinb > fncZiNB <- function(phi, size, munb){ + -...

Hi there, I'm trying to solve a ML problem where the likelihood function is a function of two numerical procedures and I'm having some problems figuring out how to do this. The log-likelihood function is of the form L(c,psi) = 1/T sum [log (f(c, psi)) - log(g(c,psi))], where c is a 2xT matrix of data and psi is the parameter vector. f(c, psi) is the transition density which can be

A question about the inner workings of glm() and dpois(): Suppose I call glm(y ~ x, family=poisson, weights = w) where y contains NON-INTEGER (but still nonnegative) values. (a) Does glm() still correctly maximise the weighted Poisson loglikelihood ? (i.e. the function given by the same formal expression as the weighted loglikelihood of independent Poisson variables Y_i except that the y_i are now nonnegative real numbers) (b) is this likely to change in future versions of R? In R version 1.3.0, glm() gives a warning that the responses...

...led each one of the outputs by using the assign and paste statements, i.e. assign(paste("group","subgroup",i),arima(...)). Works great, but here's what I need... I want to create a tables by group, and refer to particular information from each arima fit, such as AIC value, loglikelihood value, etc. My best guess is, assign(paste("group","subgroup"),get(paste("group",subgroup",i))$aic This won't work, due to the restrictions of my the assign function. If you can't understand this from the code, here it is verbally: I've got a bunch...

Hello, I am trying to use nlm to estimate the parameters that minimize the following function: Predict<-function(M,c,z){ + v = c*M^z + return(v) + } M is a variable and c and z are parameters to be estimated. I then write the negative loglikelihood function assuming normal errors: nll<-function(M,V,c,z,s){ n<-length(Mean) logl<- -.5*n*log(2*pi) -.5*n*log(s) - (1/(2*s))*sum((V-Predict(Mean,c,z))^2) return(-logl) } When I put the Mean and Variance (variables with 136 observations) into this function, and estimates for c,z, and s, it...

Hi Could anybody give me a bit of advice on some code I'm having trouble with? I've been trying to calculate the loglikelihood of a function iterated over a data of time values and I seem to be experiencing difficulty when I use the function expm(). Here's an example of what I am trying to do y<-c(5,10) #vector of 2 survival times p<-Matrix(c(1,0),1,2) #1x2 matrix Q<-Matrix(c(1,2,3,4),2,2) #2x2 ma...

...for my two treatment groups (p1 and p2), and then see whether the treatment affects p (i.e. whether p1 and p2 differ). # A subsample of the data, where 11 represent no response to 10 stimulations Nstim=c(1,1,1,1,1,1,2,2,3,4,5,6,8,11,11,11,11) # create the geometric likelihood function geometric.loglikelihood <- function(pi, y, n) { loglikelihood <- n*log(pi)-n*log(1-pi)+log(1-pi)*sum(y) return(loglikelihood) } # optimise the likelihood function test<-optim(c(0.5),geometric.loglikelihood,method="BFGS",hessian=TRUE,control=list(fnscale=-1),y=Nstim,n=length(Nstim)) # print par...