Nothing
### DPelicit.R
### Function to Prior elicitation of the precision parameter of a
### Dirichlet process prior.
###
### Copyright: Alejandro Jara, 2006-2012.
###
### Last modification: 31-08-2006.
###
### This program is free software; you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation; either version 2 of the License, or (at
### your option) any later version.
###
### This program is distributed in the hope that it will be useful, but
### WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
### General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software
### Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
###
### The author's contact information:
###
### Alejandro Jara
### Department of Statistics
### Facultad de Matematicas
### Pontificia Universidad Catolica de Chile
### Casilla 306, Correo 22
### Santiago
### Chile
### Voice: +56-2-3544506 URL : http://www.mat.puc.cl/~ajara
### Fax : +56-2-3547729 Email: atjara@uc.cl
###
DPelicit<-function(n,method='JGL',a0=NULL,b0=NULL,mean=NULL,std=NULL)
UseMethod("DPelicit")
DPelicit.default<-function(n,method='JGL',a0=NULL,b0=NULL,mean=NULL,std=NULL)
{
cl <- match.call()
if(method=='JGL')
{
if(is.null(a0))
{
foo<-newton(c(2,2),n,mean,std**2,fn1)
inp<-matrix(c(foo[1],foo[2]),ncol=2)
out<-matrix(c(mean,std),ncol=2)
colnames(inp)<-c("a0","b0")
colnames(out)<-c("Mean","Std.")
rownames(inp)<-c(" ")
rownames(out)<-c(" ")
}
else
{
mean<-(a0/b0)*(digamma((a0+n*b0)/b0)-digamma(a0/b0))
var<-mean+(a0**2/b0**2)*(trigamma((a0+n*b0)/b0)-trigamma(a0/b0))+
(((a0/b0)*(trigamma((a0+n*b0)/b0)-trigamma(a0/b0))+digamma((a0+n*b0)/b0)-digamma(a0/b0))**2)*(a0/b0**2)
inp<-matrix(c(a0,b0),ncol=2)
out<-matrix(c(mean,sqrt(var)),ncol=2)
colnames(inp)<-c("a0","b0")
colnames(out)<-c("Mean","Std.")
rownames(inp)<-c(" ")
rownames(out)<-c(" ")
}
}
if(method=='KMQ')
{
if(is.null(a0))
{
foo<-newton(c(2,2),n,mean,std**2,fn2)
inp<-matrix(c(foo[1],foo[2]),ncol=2)
out<-matrix(c(mean,std),ncol=2)
colnames(inp)<-c("a0","b0")
colnames(out)<-c("Mean","Std.")
rownames(inp)<-c(" ")
rownames(out)<-c(" ")
}
else
{
mean<-(a0/b0)*log(1+n*b0/a0)
var<-(a0/b0)*log(1+n*b0/a0)-a0/b0+((log(1+n*b0/a0)-n*b0/
(a0+n*b0))**2)*(a0/b0**2)
inp<-matrix(c(a0,b0),ncol=2)
out<-matrix(c(mean,sqrt(var)),ncol=2)
colnames(inp)<-c("a0","b0")
colnames(out)<-c("Mean","Std.")
rownames(inp)<-c(" ")
rownames(out)<-c(" ")
}
}
output<-list(inp=inp,out=out,call=cl)
class(output)<-"DPelicit"
return(output)
}
fn1<-function(y,n,mean,var)
# evaluate the function at a0 and b0
# JGL
{
x<-exp(y)
fvec<-rep(0,2)
a0<-x[1]
b0<-x[2]
fvec[1]<-mean-(a0/b0)*(digamma((a0+n*b0)/b0)-digamma(a0/b0))
fvec[2]<-var-((a0/b0)*(digamma((a0+n*b0)/b0)-digamma(a0/b0))+(a0**2/b0**2)*(trigamma((a0+n*b0)/b0)-trigamma(a0/b0))+
(((a0/b0)*(trigamma((a0+n*b0)/b0)-trigamma(a0/b0))+digamma((a0+n*b0)/b0)-digamma(a0/b0))**2)*(a0/b0**2))
return(fvec)
}
fn2<-function(y,n,mean,var)
# evaluate the function at a0 and b0
# KMQ
{
x<-exp(y)
fvec<-rep(0,2)
fvec[1]<-mean-(x[1]/x[2])*log(1+n*x[2]/x[1])
fvec[2]<-var-((x[1]/x[2])*log(1+n*x[2]/x[1])-x[1]/x[2]+((log(1+n*x[2]/x[1])-n*x[2]/(x[1]+n*x[2]))**2)*(x[1]/x[2]**2))
return(fvec)
}
fjacb<-function(x,fvec,n,mean,var,fn)
# This function computes forward-difference approximation to Jacobian
# of the function fn.
# On input x is the point at which the Jacobian es to be evaluated,
# fvec is the vector of function values at the point.
# eps is the approximate square root of the machine precision
# A.J.V., 2006.
{
eps<-1e-08
fjac<-matrix(0,ncol=2,nrow=2)
for(j in 1:2)
{
temp<-x[j]
h<-eps*abs(temp)
if(h==0)h<-eps
x[j]<-temp+h
h<-x[j]-temp
f<-fn(x,n,mean,var)
x[j]<-temp
for(i in 1:2)
{
fjac[i,j]<-(f[i]-fvec[i])/h
}
}
return(fjac)
}
newton<-function(x0,n,mean,var,fn,ntrial=200,tolx=0.00001,tolf=0.00001)
# Given a initial guess x0 for a root, take ntrial Newton-Raphson steps
# to improve the root. Stop if the root converges in either summed absolute
# variable increments tolx or summed absolute function values tolf.
# A.J.V., 2006.
{
x<-log(x0)
p<-rep(0,2)
for(k in 1:ntrial)
{
fvec<-fn(x,n,mean,var)
fjac<-fjacb(x,fvec,n,mean,var,fn)
errf<-0
for(i in 1:2)
{
errf<-errf+abs(fvec[i])
}
if(errf<=tolf)break
for(i in 1:2)
{
p[i]<--fvec[i]
}
p<-solve(fjac, p)
errx<-0
for(i in 1:2)
{
errx<-errx+abs(p[i])
x[i]<-x[i]+p[i]
}
if(errx<=tolx)break
}
return(exp(x))
}
print.DPelicit<-function(x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n","Prior elicitation for the precision parameter","\n\nCall:\n", deparse(x$call), "\n\n", sep = "")
cat("Prior information for alpha in DP(alpha G0):\n")
print.default(format(x$inp, digits = digits), print.gap = 2,
quote = FALSE)
cat("\nPrior information for the number of clusters:\n")
print.default(format(x$out, digits = digits), print.gap = 2,
quote = FALSE)
cat("\n")
invisible(x)
}
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.