Nothing
## penalized likelihood
pml_ll_der = function(para,X,offset,Y,n_one,ytwo,fid,cumsumy,posind,posindy,nb,nind,k,sigmas)
{
betas = para[1:nb]
logw = para[(nb+1):(nb+k)]
temp = pml_ll_der_eigen(X,offset,Y,fid-1,as.double(cumsumy),posind-1,posindy,nb,nind,k,betas,logw,sigmas)
return(list("objective"=temp$fn,"gradient"=temp$gr))
}
pml_ll_der2 = function(para,X,offset,Y,n_one,ytwo,fid,cumsumy,posind,posindy,nb,nind,k,sigmas)
{
betas = para[1:nb]
logw = para[(nb+1):(nb+k)]
temp = pml_ll_der_eigen(X,offset,Y,fid-1,as.double(cumsumy),posind-1,posindy,nb,nind,k,betas,logw,sigmas)
# return(list("objective"=temp$fn,"gradient"=temp$gr))
objective = temp$fn
attr(objective, "gradient") = temp$gr
objective
}
pql_ll = function(para,X,offset,Y,n_one,ytwo,fid,cumsumy,posind,posindy,nb,nind,k,betas,reml,eps,ord,intcol)
{
exps = exp(para[1])
alpha = 1/(exps-1)
gamma = para[2]
lambda = 1/(sqrt(exps)*(exps-1))
betas[intcol] = betas[intcol] - para[1]/2
repml = opt_pml(X,offset,Y,fid-1,as.double(cumsumy),posind-1,posindy,nb,nind,k,betas,para,reml,eps,ord)
if(repml$second < (-1))
{stop()}
loglik_temp = ifelse(is.nan(repml$loglik),repml$loglikp,repml$loglik)
loglik = loglik_temp + nind*gamma*log(gamma) + k*alpha*log(lambda) - k*Lgamma(alpha)
loglik = loglik + (sum(Lgamma(ytwo+gamma))-(length(posindy)-sum(n_one))*Lgamma(gamma) + sum(n_one)*log(gamma)+ n_one[2]*log(gamma+1))
loglik = loglik - 0.5*repml$logdet + log(1+repml$second)
-loglik
}
pql_gamma_ll = function(para,X,offset,Y,n_one,ytwo,fid,cumsumy,posind,posindy,nb,nind,k,betas,reml,eps,gamma,ord,intcol)
{
exps = exp(para[1])
alpha = 1/(exps-1)
lambda = 1/(sqrt(exps)*(exps-1))
betas[intcol] = betas[intcol] - para[1]/2
repml = opt_pml(X,offset,Y,fid-1,as.double(cumsumy),posind-1,posindy,nb,nind,k,betas,c(para[1],gamma),reml,eps,ord)
if(repml$second < (-1))
{stop()}
loglik_temp = ifelse(is.nan(repml$loglik),repml$loglikp,repml$loglik)
loglik = loglik_temp + nind*gamma*log(gamma) + k*alpha*log(lambda) - k*Lgamma(alpha)
loglik = loglik + (sum(Lgamma(ytwo+gamma))-(length(posindy)-sum(n_one))*Lgamma(gamma) + sum(n_one)*log(gamma)+ n_one[2]*log(gamma+1))
loglik = loglik - 0.5*repml$logdet + log(1+repml$second)
-loglik
}
pql_nbm_ll = function(para,X,offset,Y,n_one,ytwo,fid,cumsumy,posind,posindy,nb,nind,k,betas,reml,eps,ord,intcol)
{
alpha = para[1]
gamma = para[2]
betas[intcol] = betas[intcol] - para[1]/2
repml = opt_pml_nbm(X,offset,Y,fid-1,as.double(cumsumy),posind-1,posindy,nb,nind,k,betas,para,reml,eps,1)
loglik = repml$loglik + nind*gamma*log(gamma) + sum(Lgamma(Y+gamma)) -length(Y)*Lgamma(gamma)
loglik = loglik - 0.5*repml$logdet
-loglik
}
pql_nbm_gamma_ll = function(para,X,offset,Y,n_one,ytwo,fid,cumsumy,posind,posindy,nb,nind,k,betas,reml,eps,gamma,ord,intcol)
{
alpha = para[1]
betas[intcol] = betas[intcol] - para[1]/2
repml = opt_pml_nbm(X,offset,Y,fid-1,as.double(cumsumy),posind-1,posindy,nb,nind,k,betas,c(para[1],gamma),reml,eps,1)
loglik = repml$loglik + nind*gamma*log(gamma) + sum(Lgamma(Y+gamma)) -length(Y)*Lgamma(gamma)
loglik = loglik - 0.5*repml$logdet
-loglik
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.