In dariamed/gjamed: What the package does (short line)
rm(list=ls())
#Packages for gamma function estimation
library(pracma)
library(Brobdingnag)
library(copula)
library(plotly)
library(distrEx)
$V_{kn}$ functions and Prior cluster distribtion.
############### V for PY / Dir / NGG########################
v_py<- function(kval, sigma,theta,npoints){
c_v<-1:(kval-1)
v_nk<- (theta +sigma*c_v)
Vup<- prod(v_nk)
n_vec<- 1:(npoints-1)
Vlow<- prod(theta +n_vec)
V_t_nk<-Vup/Vlow
return(V_t_nk)
}
############### V for NGG########################
v_ng<- function(beta, sigma, kval, npoints){
sum<-0
coef_low<-as.brob(gamma(npoints))
coef_high<-as.brob(exp(beta)* sigma^(kval-1))
coef<- coef_high/coef_low
incv<- as.brob(vector(length=npoints))
for(i in (0:(npoints-1))){
gn<- as.brob(gammainc(kval - i/sigma,beta)[2])
#gn<- gamma_inc(kval - i/sigma,beta)
ckn<- as.brob(choose(npoints-1,i))
sum<- sum + ((-1)^i)*(beta^(i/sigma))*ckn*gn
incv[i+1]<- gn
}
sumf<- sum/coef
sumn<- as.numeric(sumf)
return(list(sum=sumf, incg= incv))
}
v_ng2<- function(beta, sigma, kval, npoints){
sum<-0
incv<- as.brob(vector(length=npoints))
for(i in (0:(npoints-1))){
coef_low<-as.brob(gamma(npoints))
coef_high<-as.brob(exp(beta)* sigma^(kval-1))
coef<- coef_high/coef_low
gn<- as.brob(gammainc(kval - i/sigma,beta)[2])
ckn<- as.brob(choose(npoints-1,i))
sum<- sum + ((-1)^i)*(beta^(i/sigma))*ckn*gn*coef
incv[i+1]<- gn
}
sumf<- sum
sumn<- as.numeric(sumf)
return(list(sum=sumn, incg= incv))
}
###############Generalized coefficient########################
gen_fac_coef<-function(kval,sigma,npoints){
sum<-0
kfac<-factorial(kval)
for(i in (0:kval)){
n_vec<- 0:(npoints-1)
sn<- prod(-i*sigma +n_vec)
ckn<- choose(kval,i)
#print((-1)^i*ckn*sn)
sum<- sum + ((-1)^i)*ckn*sn
}
sumf<- sum/kfac
return(sumf)
}
########### density for NGG #############################
prob_ng<- function(kg, sigma, npoints, beta){
pb_v_all<- v_ng2(beta, sigma, kg, npoints)
pb_v<-as.brob(pb_v_all$sum)
pb_gen<- as.brob(gen_fac_coef(kg,sigma, npoints))
prob<- (pb_v*pb_gen)/(as.brob(sigma^kg))
prob_num<- as.numeric(prob)
return(prob_num)
}
########### density for PY#############################
prob_py<- function(kg, npoints, sigma, theta){
pb_v<- v_py(kg,sigma, theta,npoints)
pb_gen<- gen_fac_coef(kg,sigma, npoints)
prob<- (pb_v*pb_gen)/(sigma^kg)
return(prob)
}
########### density for Dirichlet#############################
prob_dir<- function(k, npoints, theta){
n_vec<- 0:(npoints-1)
theta_n<- prod(theta +n_vec)
prob<- ((theta^k) *(abs(Stirling1(npoints,k))))/theta_n
return(prob)
}
prob_dir_large_dim<- function(k, npoints, theta){
n_vec<-as.brob( 0:(npoints-1))
theta_n<- prod(theta +n_vec)
stir<- as.brob(abs(Stirling1(npoints,k)))
powerk<- as.brob((theta^k))
prob_brob<- powerk*(stir/theta_n)
prob<- as.numeric(prob_brob)
return(prob)
}
#####################################################################################################################################
#prib_ng(kg, npoints, sigma, beta)
k_vec<-seq(1,50,by=5)
sigma_vec<-seq(0.2,0.8, by=0.1)
z<- outer(k_vec,sigma_vec,Vectorize(prob_ng),npoints=50, beta=1)
p<- plot_ly(showscale = TRUE) %>%
add_surface(x=k_vec, y=sigma_vec,z =z, cmin = min(z), cmax = max(z),colorbar=list(title='PY'), colorscale = list(c(0,1),c("rgb(255,112,184)","rgb(128,0,64)")),opacity = 0.98) %>%
layout(title="Prior distribution", scene = list(xaxis= list(title="K"),yaxis= list(title="sigma"),zaxis= list(title="N",range = c(min(z),max(z)))))
p
Data generation
Gaussian mixture:
$$Y \sim \frac{1}{2}N(1,0.2) +\frac{1}{2}N(5,0.2)$$
gen_data<- function(mu_vec=c(1,5,3,5), sigma_vec=rep(0.1, 4), ns=50, prob_vec=c(0.25,0.25,0.25,0.25) ){
l<- length(mu_vec)
components <- sample(1:l,prob=prob_vec,size=ns,replace=TRUE)
mus <- mu_vec
sds <- sqrt(sigma_vec)
samples <- rnorm(n=ns,mean=mus[components],sd=sds[components])
return(samples)
}
set.seed(123)
sample<- gen_data(mu_vec=c(1,10),sigma_vec=rep(0.2, 2),ns=10000,prob_vec=c(0.5,0.5))
hist(sample,probability=TRUE,breaks=20,col=grey(.9))
Model
\begin{equation}
\begin{aligned}
(X_i \mid (\mu_i,\tau_i)) &\stackrel{iid}{\sim} N(\mu_i,\tau_i^{-1})\
((\mu_i, \tau_i) \mid G) &\stackrel{iid}{\sim} G \
G &\sim \mathcal{P} .
\end{aligned}
\end{equation}
Firstly, assume $\tau_i=1$,then we have $X_i \sim N(\mu_i, 1)$, base measure $G_0 = N(0,1)$, $Y_i= (\mu_i)$. For sampling we use Polya Urn characterization
\begin{equation}
p(Y_i | {Y_{-i}}) = \frac{\alpha}{\alpha+n-1} G_0 + \frac{1}{\alpha + n -1} \sum_{j=1}^{n-1} \delta_{Y_j} \qquad \forall i=1,\ldots ,n.
\end{equation}
Gibbs sampler:
\begin{equation}
\begin{aligned}
G &\mid \theta, X \sim DP(\alpha +n,\frac{\alpha}{(\alpha + n)} G_0 + \frac{1}{\alpha + n)}\sum_j \delta_{\theta_j)}\
\theta_i &\mid \theta_{-i},X, G \sim \sum_{j\neq i}q_{ij} \delta_{\theta_j} + r_i G \
q_{ij} &= b N(x_,\theta_j)\
r_i &= b\alpha \int f(x_i, \theta)dG_0(\theta), \text{Here: } N(x_i, \theta,1)N(\theta, 0,1)\propto N(x_i/2, 1/2).\
\end{aligned}
\end{equation}
Prior
Gbase<- function(num, mean=0, sigma=1){
return( rnorm(n=1, mean, sigma))
}
PUrn_sample<- function(n, alpha, Gbase){
#theta.1 sample from base measure
theta<- Gbase(1)
group_vec <- c(1)
for (i in 2:n){
#probability to sample from base measure
p1<- alpha/ (alpha + i-1)
#probability to sample from existing atoms
p2<- 1/ (alpha + i-1)
pvec<- c(p1, p2*group_vec)
# prob <- c(alpha/(length(count.vector)+alpha), count.vector/(length(count.vector)+M))
group_num<- sample(0:length(group_vec), prob=pvec, size=1)
if (group_num == 0){
#add new atom
theta <- c(theta, Gbase(1))
group_vec <- c(group_vec, 1)
}else{
#update existing group for group_vec
group_vec[group_num] = group_vec[group_num] + 1
}
}
return(list(theta=theta, group= group_vec))
}
alpha_0<-1
n<- 500
x<- PUrn_sample(n, alpha_0,Gbase)
k_groups<- x$group
weights<- c(k_groups/(alpha_0+n), alpha_0/(alpha_0+n))
df_weights <- data.frame(matrix(NA, nrow =length(weights), ncol =1))
df_weights$pw<-weights
df_weights$tr<-1:length(weights)
pl_weigths<- ggplot(df_weights, aes(x=tr, y=pw)) +
geom_segment( aes(x=tr,xend=tr,y=0,yend=pw)) +
geom_point( size=0.5, color="red", fill=alpha("blue", 0.3), alpha=0.4, shape=21, stroke=2)+ labs(title=paste0("Prior weights"))+
theme_bw() + theme(axis.text.x = element_text(angle = 0, hjust = 1,size = 10), strip.text = element_text(size = 15),legend.position = "top", plot.title = element_text(hjust = 0.5))
pl_weigths
df_weights <- data.frame(matrix(NA, nrow =length(weights), ncol =1))
df_weights$pw<-weights
df_weights$tr<-c(x$theta,20)
pl_weigths<- ggplot(df_weights, aes(x=tr, y=pw)) +
geom_segment( aes(x=tr,xend=tr,y=0,yend=pw)) +
geom_point( size=0.5, color="red", fill=alpha("blue", 0.3), alpha=0.4, shape=21, stroke=2)+ labs(title=paste0("Prior weights"))+
theme_bw() + theme(axis.text.x = element_text(angle = 0, hjust = 1,size = 10), strip.text = element_text(size = 15),legend.position = "top", plot.title = element_text(hjust = 0.5))
pl_weigths
PolyaUrnGenerator <- function(num, M, G0.generator) {
X.vector <- G0.generator(n = 1)
count.vector <- 1
for (i in 2:num){
prob <- c(M/(length(count.vector)+M), count.vector/(length(count.vector)+M))
category <- which(rmultinom(n=1, size=1,prob=prob)==1) - 1
if (category == 0){
X.vector <- c(X.vector, G0.generator(n = 1))
count.vector <- c(count.vector, 1)
}else{
count.vector[category] = count.vector[category] + 1
}
}
return(list(X.vector, count.vector))
}
## Define a Polya Urn sampler
PolyaUrnSampler <- function(n, M, cutoff, num = 10000, G0.generator=rnorm){
output <- matrix(NA, nrow = n, ncol = length(cutoff) + 1)
cutoff <- c(-Inf, cutoff, Inf)
for(i in 1:n){
PU.sample <- PolyaUrnGenerator(num=num, M=M, G0.generator=G0.generator)
for(j in 1:ncol(output)){
output[i,j] <- sum(PU.sample[[2]][intersect(which(PU.sample[[1]] > cutoff[j]), which(PU.sample[[1]] <= cutoff[j+1]))])/num
}
}
return(output)
}
##
x<- PolyaUrnGenerator(1,1,rnorm)
y<- PUrn_sample(1, 1, Gbase)
Posterior_P_Urn_exact<- function(X,theta,i,alpha){
q_weights<- dnorm(X[i], theta)
r_weight<- alpha*dnorm(X[i], mean=0, sd=sqrt(2))
#sum should be one
num_norm<- sum(q_weights)+r_weight
q_weights_norm<- q_weights/num_norm
r_weight_norm<- r_weight/num_norm
label <- sample(0:length(theta),prob=c(r_weight_norm,q_weights_norm),size=1,replace=TRUE)
if (label==0){
theta_new<- rnorm(1, mean=X[i], sd=sqrt(1/2))
}
else{
theta_new <- theta[label]
}
}
sample<- gen_data(mu_vec=c(1,10),sigma_vec=rep(0.2, 2),ns=100,prob_vec=c(0.5,0.5))
X<- sample
theta.matrix <- matrix(NA, nrow = 10000, ncol = length(X))
alpha<-1
#weights.matrix <- matrix(NA, nrow = 100, ncol = length(X)+1)
theta.init<- runif(length(X),-3,3)
theta.matrix[1,] <- theta.init
for (i in 2:nrow(theta.matrix)){
theta_last<- theta.matrix[i-1,]
N <- length(X)
temp_theta <- c(theta_last, rep(NA, N))
for (j in 1:N){
temp_theta[N+j] <- Posterior_P_Urn_exact(X, temp_theta[(j+1):(j+N-1)],j,alpha)
}
theta_S<- temp_theta[(N+1):(N*2)]
theta.matrix[i,] <-theta_S
}
#burnin period
theta.final <- theta.matrix[5000:10000,]
#xgrid <- seq(from= -5, to=15,length=100)
#fgrid <- NULL
fbar.H <- function(xgrid,wh)
{ ## return a draw F ~ p(F | ...) (approx)
fx <- rep(0,length(xgrid))
for(h in 1:length(wh))
fx <- fx + dnorm(xgrid,m=wh[h],sd=1)
fx<- fx + M/(n+M)*dnorm(xgrid,m=m0,sd=sqrt(B0+sig))
return(fx/length(wh))
}
fbar.H_base <- function(xgrid,wh)
{ ## return a draw F ~ p(F | ...) (approx)
fx <- rep(0,length(xgrid))
for(h in 1:length(wh)){
fx <- fx + dnorm(xgrid,m=wh[h],sd=1)
}
fx_m<- fx/ (length(wh)+ alpha)
fx_m<- fx_m + alpha/(length(wh)+alpha)*dnorm(xgrid,m=0,sd=sqrt(2))
return(fx_m)
}
ecfbar.H_base <- function(xgrid,wh)
{ ## return a draw F ~ p(F | ...) (approx)
fx <- rep(0,length(xgrid))
for(h in 1:length(wh)){
fx <- fx + pnorm(xgrid,m=wh[h],sd=1)
}
fx_m<- fx/ (length(wh)+ alpha)
fx_m<- fx_m + alpha/(length(wh)+alpha)*pnorm(xgrid,m=0,sd=sqrt(2))
return(fx_m)
}
gibbs.H <- function(n.iter=10000){
xgrid <- seq(from= -5, to=12,length=50)
cgrid<- seq(from= -5, to=15,length=50)
fgrid <- NULL
ecfgrid <- NULL
plot(density(X),xlab="X",ylab="Y",bty="l",type="l",xlim=c(-5, 15),ylim=c(0,1), main="")
for(iter in 1:nrow(theta.final)){
## record draw F ~ p(F | th,sig,y) (approx)
f <- fbar.H_base(xgrid,theta.final[iter,])
lines(xgrid,f,col=iter,lty=3)
fgrid <- rbind(fgrid,f)
}
## add overall average (= posterior mean) to the plot
fbar <- apply(fgrid,2,mean)
lines(xgrid,fbar,lwd=3,col=2)
plot(ecdf(X),xlab="X",ylab="Y",bty="l", main="")
for(iter in 1:nrow(theta.final)){
## record draw F ~ p(F | th,sig,y) (approx)
ecf<- ecfbar.H_base(cgrid,theta.final[iter,])
lines(cgrid,ecf,col=iter,lty=3)
ecfgrid <- rbind(ecfgrid,ecf)
}
## add overall average (= posterior mean) to the plot
fbar <- apply(fgrid,2,mean)
ecfbar<- apply(ecfgrid,2,mean)
lines(cgrid,ecfbar,lwd=3,col=2)
return(list(f_ap=fbar, ecf=ecfbar))
}
f_approx<- gibbs.H()
ks.test(f_approx$ecf, ecdf(X))
Posterior_P_Urn_exact<- function(X,theta,i,alpha){
q_weights<- dnorm(X[i], theta)
r_weight<- alpha*dnorm(X[i], mean=0, sd=sqrt(2))
#sum should be one
num_norm<- sum(q_weights)+r_weight
q_weights_norm<- q_weights/num_norm
r_weight_norm<- r_weight/num_norm
label <- sample(0:length(theta),prob=c(r_weight_norm,q_weights_norm),size=1,replace=TRUE)
if (label==0){
theta_new<- rnorm(1, mean=X[i], sd=sqrt(1/2))
}
else{
theta_new <- theta[label]
}
}
sample<- gen_data(mu_vec=c(1,5),sigma_vec=rep(0.2, 2),ns=50,prob_vec=c(0.5,0.5))
X<- sample
theta.matrix <- matrix(NA, nrow = 10000, ncol = length(X))
alpha<-1
#weights.matrix <- matrix(NA, nrow = 100, ncol = length(X)+1)
theta.init<- runif(length(X),-3,3)
theta.matrix[1,] <- theta.init
for (i in 2:nrow(theta.matrix)){
theta_last<- theta.matrix[i-1,]
N <- length(X)
temp_theta <- c(theta_last, rep(NA, N))
for (j in 1:N){
temp_theta[N+j] <- Posterior_P_Urn_exact(X, temp_theta[(j+1):(j+N-1)],j,alpha)
}
theta_S<- temp_theta[(N+1):(N*2)]
theta.matrix[i,] <-theta_S
}
#burnin period
theta.final <- theta.matrix[5000:10000,]
xgrid <- seq(from= -10, to=2,length=50)
fgrid <- NULL
fbar.H <- function(xgrid,wh)
{ ## return a draw F ~ p(F | ...) (approx)
fx <- rep(0,length(xgrid))
for(h in 1:length(wh))
fx <- fx + dnorm(xgrid,m=wh[h],sd=1)
return(fx/length(wh))
}
gibbs.H <- function(n.iter=500){
xgrid <- seq(from= -10, to=10,length=50)
fgrid <- NULL
plot(density(X),xlab="X",ylab="Y",bty="l",type="l",xlim=c(-10, 10),ylim=c(0,1), main="")
for(iter in 1:nrow(theta.matrix)){
## record draw F ~ p(F | th,sig,y) (approx)
f <- fbar.H(xgrid,theta.matrix[iter,])
lines(xgrid,f,col=iter,lty=3)
fgrid <- rbind(fgrid,f)
}
## add overall average (= posterior mean) to the plot
fbar <- apply(fgrid,2,mean)
lines(xgrid,fbar,lwd=3,col=2)
}
Posterior_G<- function(n, alpha, Gbase, theta){
}
Approximation 1st order
$$\theta_{n} \mid \theta_{n-1},..\theta_1 \sim \frac{k_n \alpha}{n} H + \frac{1}{n} \sum_i^{k_n} (n_i - \alpha)\delta_{\theta^*_j}$$
Approximation 2nd order
$$\theta_{n} \mid \theta_{n-1},..\theta_1 \sim \frac{k_n \alpha + \beta_n}{n + \beta_n} H + \frac{1}{n + \beta_n} \sum_i^{k_n} (n_i - \alpha)\delta_{\theta^*_j}$$
For PY process $\beta_n = \theta$
polya_urn_model <- function(base_measure, N_ball, alpha) {
balls <- c()
for (i in 1:N_ball) {
if (runif(1) < alpha / (alpha + length(balls))) {
#Add a new ball color.
new_color <- base_measure()
balls <- c(balls, new_color)
} else {
#Pick out a ball from the urn, and add back a # ball of the same color
ball <- balls[sample(1:length(balls), 1)]
balls <- c(balls, ball)
} }
balls
}
alpha_vect <- c(1, 10, 100, 1000)
N_urns <- 3
sigma2 <- c(1,.4,.2)
N_draws <- 100
N_xaxis <- 200
x_axis <- seq(-3, 3, length = N_xaxis)
result <- NULL
for (alpha in alpha_vect) {
PU <- polya_urn_model(function() rnorm(1), N_draws, alpha)
for (u in 1:N_urns) {
res <- mapply(function(mean) dnorm(x_axis, rep(mean, N_xaxis),
rep(sigma2[u], N_xaxis)), PU)
res <- apply(res, 1, mean)
new_draw <- cbind(res, x_axis, alpha, sigma2[u])
result <- rbind(result, new_draw)
} }
result <- as.data.frame(result)
names(result) <- c("density", "x", "alpha", "sigma2")
DP_mixt <- qplot(data = result, y = density, x = x,
geom = c("line", "area")) +
facet_grid(alpha ~ sigma2, labeller =
label_bquote(rows = alpha == .(alpha),
cols = sigma^2 == .(sigma2))) +
aes(color = as.factor(alpha)) +
theme(legend.position = "none")
DP_mixt
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
dariamed/gjamed documentation built on Sept. 27, 2019, 5:31 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
rm(list=ls()) #Packages for gamma function estimation library(pracma) library(Brobdingnag) library(copula) library(plotly) library(distrEx)
$V_{kn}$ functions and Prior cluster distribtion.
############### V for PY / Dir / NGG######################## v_py<- function(kval, sigma,theta,npoints){ c_v<-1:(kval-1) v_nk<- (theta +sigma*c_v) Vup<- prod(v_nk) n_vec<- 1:(npoints-1) Vlow<- prod(theta +n_vec) V_t_nk<-Vup/Vlow return(V_t_nk) } ############### V for NGG######################## v_ng<- function(beta, sigma, kval, npoints){ sum<-0 coef_low<-as.brob(gamma(npoints)) coef_high<-as.brob(exp(beta)* sigma^(kval-1)) coef<- coef_high/coef_low incv<- as.brob(vector(length=npoints)) for(i in (0:(npoints-1))){ gn<- as.brob(gammainc(kval - i/sigma,beta)[2]) #gn<- gamma_inc(kval - i/sigma,beta) ckn<- as.brob(choose(npoints-1,i)) sum<- sum + ((-1)^i)*(beta^(i/sigma))*ckn*gn incv[i+1]<- gn } sumf<- sum/coef sumn<- as.numeric(sumf) return(list(sum=sumf, incg= incv)) } v_ng2<- function(beta, sigma, kval, npoints){ sum<-0 incv<- as.brob(vector(length=npoints)) for(i in (0:(npoints-1))){ coef_low<-as.brob(gamma(npoints)) coef_high<-as.brob(exp(beta)* sigma^(kval-1)) coef<- coef_high/coef_low gn<- as.brob(gammainc(kval - i/sigma,beta)[2]) ckn<- as.brob(choose(npoints-1,i)) sum<- sum + ((-1)^i)*(beta^(i/sigma))*ckn*gn*coef incv[i+1]<- gn } sumf<- sum sumn<- as.numeric(sumf) return(list(sum=sumn, incg= incv)) } ###############Generalized coefficient######################## gen_fac_coef<-function(kval,sigma,npoints){ sum<-0 kfac<-factorial(kval) for(i in (0:kval)){ n_vec<- 0:(npoints-1) sn<- prod(-i*sigma +n_vec) ckn<- choose(kval,i) #print((-1)^i*ckn*sn) sum<- sum + ((-1)^i)*ckn*sn } sumf<- sum/kfac return(sumf) } ########### density for NGG ############################# prob_ng<- function(kg, sigma, npoints, beta){ pb_v_all<- v_ng2(beta, sigma, kg, npoints) pb_v<-as.brob(pb_v_all$sum) pb_gen<- as.brob(gen_fac_coef(kg,sigma, npoints)) prob<- (pb_v*pb_gen)/(as.brob(sigma^kg)) prob_num<- as.numeric(prob) return(prob_num) } ########### density for PY############################# prob_py<- function(kg, npoints, sigma, theta){ pb_v<- v_py(kg,sigma, theta,npoints) pb_gen<- gen_fac_coef(kg,sigma, npoints) prob<- (pb_v*pb_gen)/(sigma^kg) return(prob) } ########### density for Dirichlet############################# prob_dir<- function(k, npoints, theta){ n_vec<- 0:(npoints-1) theta_n<- prod(theta +n_vec) prob<- ((theta^k) *(abs(Stirling1(npoints,k))))/theta_n return(prob) } prob_dir_large_dim<- function(k, npoints, theta){ n_vec<-as.brob( 0:(npoints-1)) theta_n<- prod(theta +n_vec) stir<- as.brob(abs(Stirling1(npoints,k))) powerk<- as.brob((theta^k)) prob_brob<- powerk*(stir/theta_n) prob<- as.numeric(prob_brob) return(prob) } ##################################################################################################################################### #prib_ng(kg, npoints, sigma, beta) k_vec<-seq(1,50,by=5) sigma_vec<-seq(0.2,0.8, by=0.1) z<- outer(k_vec,sigma_vec,Vectorize(prob_ng),npoints=50, beta=1) p<- plot_ly(showscale = TRUE) %>% add_surface(x=k_vec, y=sigma_vec,z =z, cmin = min(z), cmax = max(z),colorbar=list(title='PY'), colorscale = list(c(0,1),c("rgb(255,112,184)","rgb(128,0,64)")),opacity = 0.98) %>% layout(title="Prior distribution", scene = list(xaxis= list(title="K"),yaxis= list(title="sigma"),zaxis= list(title="N",range = c(min(z),max(z))))) p
Data generation
Gaussian mixture:
$$Y \sim \frac{1}{2}N(1,0.2) +\frac{1}{2}N(5,0.2)$$
gen_data<- function(mu_vec=c(1,5,3,5), sigma_vec=rep(0.1, 4), ns=50, prob_vec=c(0.25,0.25,0.25,0.25) ){ l<- length(mu_vec) components <- sample(1:l,prob=prob_vec,size=ns,replace=TRUE) mus <- mu_vec sds <- sqrt(sigma_vec) samples <- rnorm(n=ns,mean=mus[components],sd=sds[components]) return(samples) } set.seed(123) sample<- gen_data(mu_vec=c(1,10),sigma_vec=rep(0.2, 2),ns=10000,prob_vec=c(0.5,0.5)) hist(sample,probability=TRUE,breaks=20,col=grey(.9))
Model
\begin{equation} \begin{aligned} (X_i \mid (\mu_i,\tau_i)) &\stackrel{iid}{\sim} N(\mu_i,\tau_i^{-1})\ ((\mu_i, \tau_i) \mid G) &\stackrel{iid}{\sim} G \ G &\sim \mathcal{P} . \end{aligned} \end{equation}
Firstly, assume $\tau_i=1$,then we have $X_i \sim N(\mu_i, 1)$, base measure $G_0 = N(0,1)$, $Y_i= (\mu_i)$. For sampling we use Polya Urn characterization
\begin{equation} p(Y_i | {Y_{-i}}) = \frac{\alpha}{\alpha+n-1} G_0 + \frac{1}{\alpha + n -1} \sum_{j=1}^{n-1} \delta_{Y_j} \qquad \forall i=1,\ldots ,n. \end{equation}
Gibbs sampler: \begin{equation} \begin{aligned} G &\mid \theta, X \sim DP(\alpha +n,\frac{\alpha}{(\alpha + n)} G_0 + \frac{1}{\alpha + n)}\sum_j \delta_{\theta_j)}\ \theta_i &\mid \theta_{-i},X, G \sim \sum_{j\neq i}q_{ij} \delta_{\theta_j} + r_i G \ q_{ij} &= b N(x_,\theta_j)\ r_i &= b\alpha \int f(x_i, \theta)dG_0(\theta), \text{Here: } N(x_i, \theta,1)N(\theta, 0,1)\propto N(x_i/2, 1/2).\
\end{aligned} \end{equation}
Prior
Gbase<- function(num, mean=0, sigma=1){ return( rnorm(n=1, mean, sigma)) } PUrn_sample<- function(n, alpha, Gbase){ #theta.1 sample from base measure theta<- Gbase(1) group_vec <- c(1) for (i in 2:n){ #probability to sample from base measure p1<- alpha/ (alpha + i-1) #probability to sample from existing atoms p2<- 1/ (alpha + i-1) pvec<- c(p1, p2*group_vec) # prob <- c(alpha/(length(count.vector)+alpha), count.vector/(length(count.vector)+M)) group_num<- sample(0:length(group_vec), prob=pvec, size=1) if (group_num == 0){ #add new atom theta <- c(theta, Gbase(1)) group_vec <- c(group_vec, 1) }else{ #update existing group for group_vec group_vec[group_num] = group_vec[group_num] + 1 } } return(list(theta=theta, group= group_vec)) } alpha_0<-1 n<- 500 x<- PUrn_sample(n, alpha_0,Gbase) k_groups<- x$group weights<- c(k_groups/(alpha_0+n), alpha_0/(alpha_0+n)) df_weights <- data.frame(matrix(NA, nrow =length(weights), ncol =1)) df_weights$pw<-weights df_weights$tr<-1:length(weights) pl_weigths<- ggplot(df_weights, aes(x=tr, y=pw)) + geom_segment( aes(x=tr,xend=tr,y=0,yend=pw)) + geom_point( size=0.5, color="red", fill=alpha("blue", 0.3), alpha=0.4, shape=21, stroke=2)+ labs(title=paste0("Prior weights"))+ theme_bw() + theme(axis.text.x = element_text(angle = 0, hjust = 1,size = 10), strip.text = element_text(size = 15),legend.position = "top", plot.title = element_text(hjust = 0.5)) pl_weigths df_weights <- data.frame(matrix(NA, nrow =length(weights), ncol =1)) df_weights$pw<-weights df_weights$tr<-c(x$theta,20) pl_weigths<- ggplot(df_weights, aes(x=tr, y=pw)) + geom_segment( aes(x=tr,xend=tr,y=0,yend=pw)) + geom_point( size=0.5, color="red", fill=alpha("blue", 0.3), alpha=0.4, shape=21, stroke=2)+ labs(title=paste0("Prior weights"))+ theme_bw() + theme(axis.text.x = element_text(angle = 0, hjust = 1,size = 10), strip.text = element_text(size = 15),legend.position = "top", plot.title = element_text(hjust = 0.5)) pl_weigths PolyaUrnGenerator <- function(num, M, G0.generator) { X.vector <- G0.generator(n = 1) count.vector <- 1 for (i in 2:num){ prob <- c(M/(length(count.vector)+M), count.vector/(length(count.vector)+M)) category <- which(rmultinom(n=1, size=1,prob=prob)==1) - 1 if (category == 0){ X.vector <- c(X.vector, G0.generator(n = 1)) count.vector <- c(count.vector, 1) }else{ count.vector[category] = count.vector[category] + 1 } } return(list(X.vector, count.vector)) } ## Define a Polya Urn sampler PolyaUrnSampler <- function(n, M, cutoff, num = 10000, G0.generator=rnorm){ output <- matrix(NA, nrow = n, ncol = length(cutoff) + 1) cutoff <- c(-Inf, cutoff, Inf) for(i in 1:n){ PU.sample <- PolyaUrnGenerator(num=num, M=M, G0.generator=G0.generator) for(j in 1:ncol(output)){ output[i,j] <- sum(PU.sample[[2]][intersect(which(PU.sample[[1]] > cutoff[j]), which(PU.sample[[1]] <= cutoff[j+1]))])/num } } return(output) } ## x<- PolyaUrnGenerator(1,1,rnorm) y<- PUrn_sample(1, 1, Gbase)
Posterior_P_Urn_exact<- function(X,theta,i,alpha){ q_weights<- dnorm(X[i], theta) r_weight<- alpha*dnorm(X[i], mean=0, sd=sqrt(2)) #sum should be one num_norm<- sum(q_weights)+r_weight q_weights_norm<- q_weights/num_norm r_weight_norm<- r_weight/num_norm label <- sample(0:length(theta),prob=c(r_weight_norm,q_weights_norm),size=1,replace=TRUE) if (label==0){ theta_new<- rnorm(1, mean=X[i], sd=sqrt(1/2)) } else{ theta_new <- theta[label] } } sample<- gen_data(mu_vec=c(1,10),sigma_vec=rep(0.2, 2),ns=100,prob_vec=c(0.5,0.5)) X<- sample theta.matrix <- matrix(NA, nrow = 10000, ncol = length(X)) alpha<-1 #weights.matrix <- matrix(NA, nrow = 100, ncol = length(X)+1) theta.init<- runif(length(X),-3,3) theta.matrix[1,] <- theta.init for (i in 2:nrow(theta.matrix)){ theta_last<- theta.matrix[i-1,] N <- length(X) temp_theta <- c(theta_last, rep(NA, N)) for (j in 1:N){ temp_theta[N+j] <- Posterior_P_Urn_exact(X, temp_theta[(j+1):(j+N-1)],j,alpha) } theta_S<- temp_theta[(N+1):(N*2)] theta.matrix[i,] <-theta_S } #burnin period theta.final <- theta.matrix[5000:10000,] #xgrid <- seq(from= -5, to=15,length=100) #fgrid <- NULL fbar.H <- function(xgrid,wh) { ## return a draw F ~ p(F | ...) (approx) fx <- rep(0,length(xgrid)) for(h in 1:length(wh)) fx <- fx + dnorm(xgrid,m=wh[h],sd=1) fx<- fx + M/(n+M)*dnorm(xgrid,m=m0,sd=sqrt(B0+sig)) return(fx/length(wh)) } fbar.H_base <- function(xgrid,wh) { ## return a draw F ~ p(F | ...) (approx) fx <- rep(0,length(xgrid)) for(h in 1:length(wh)){ fx <- fx + dnorm(xgrid,m=wh[h],sd=1) } fx_m<- fx/ (length(wh)+ alpha) fx_m<- fx_m + alpha/(length(wh)+alpha)*dnorm(xgrid,m=0,sd=sqrt(2)) return(fx_m) } ecfbar.H_base <- function(xgrid,wh) { ## return a draw F ~ p(F | ...) (approx) fx <- rep(0,length(xgrid)) for(h in 1:length(wh)){ fx <- fx + pnorm(xgrid,m=wh[h],sd=1) } fx_m<- fx/ (length(wh)+ alpha) fx_m<- fx_m + alpha/(length(wh)+alpha)*pnorm(xgrid,m=0,sd=sqrt(2)) return(fx_m) } gibbs.H <- function(n.iter=10000){ xgrid <- seq(from= -5, to=12,length=50) cgrid<- seq(from= -5, to=15,length=50) fgrid <- NULL ecfgrid <- NULL plot(density(X),xlab="X",ylab="Y",bty="l",type="l",xlim=c(-5, 15),ylim=c(0,1), main="") for(iter in 1:nrow(theta.final)){ ## record draw F ~ p(F | th,sig,y) (approx) f <- fbar.H_base(xgrid,theta.final[iter,]) lines(xgrid,f,col=iter,lty=3) fgrid <- rbind(fgrid,f) } ## add overall average (= posterior mean) to the plot fbar <- apply(fgrid,2,mean) lines(xgrid,fbar,lwd=3,col=2) plot(ecdf(X),xlab="X",ylab="Y",bty="l", main="") for(iter in 1:nrow(theta.final)){ ## record draw F ~ p(F | th,sig,y) (approx) ecf<- ecfbar.H_base(cgrid,theta.final[iter,]) lines(cgrid,ecf,col=iter,lty=3) ecfgrid <- rbind(ecfgrid,ecf) } ## add overall average (= posterior mean) to the plot fbar <- apply(fgrid,2,mean) ecfbar<- apply(ecfgrid,2,mean) lines(cgrid,ecfbar,lwd=3,col=2) return(list(f_ap=fbar, ecf=ecfbar)) } f_approx<- gibbs.H() ks.test(f_approx$ecf, ecdf(X))
Posterior_P_Urn_exact<- function(X,theta,i,alpha){ q_weights<- dnorm(X[i], theta) r_weight<- alpha*dnorm(X[i], mean=0, sd=sqrt(2)) #sum should be one num_norm<- sum(q_weights)+r_weight q_weights_norm<- q_weights/num_norm r_weight_norm<- r_weight/num_norm label <- sample(0:length(theta),prob=c(r_weight_norm,q_weights_norm),size=1,replace=TRUE) if (label==0){ theta_new<- rnorm(1, mean=X[i], sd=sqrt(1/2)) } else{ theta_new <- theta[label] } } sample<- gen_data(mu_vec=c(1,5),sigma_vec=rep(0.2, 2),ns=50,prob_vec=c(0.5,0.5)) X<- sample theta.matrix <- matrix(NA, nrow = 10000, ncol = length(X)) alpha<-1 #weights.matrix <- matrix(NA, nrow = 100, ncol = length(X)+1) theta.init<- runif(length(X),-3,3) theta.matrix[1,] <- theta.init for (i in 2:nrow(theta.matrix)){ theta_last<- theta.matrix[i-1,] N <- length(X) temp_theta <- c(theta_last, rep(NA, N)) for (j in 1:N){ temp_theta[N+j] <- Posterior_P_Urn_exact(X, temp_theta[(j+1):(j+N-1)],j,alpha) } theta_S<- temp_theta[(N+1):(N*2)] theta.matrix[i,] <-theta_S } #burnin period theta.final <- theta.matrix[5000:10000,] xgrid <- seq(from= -10, to=2,length=50) fgrid <- NULL fbar.H <- function(xgrid,wh) { ## return a draw F ~ p(F | ...) (approx) fx <- rep(0,length(xgrid)) for(h in 1:length(wh)) fx <- fx + dnorm(xgrid,m=wh[h],sd=1) return(fx/length(wh)) } gibbs.H <- function(n.iter=500){ xgrid <- seq(from= -10, to=10,length=50) fgrid <- NULL plot(density(X),xlab="X",ylab="Y",bty="l",type="l",xlim=c(-10, 10),ylim=c(0,1), main="") for(iter in 1:nrow(theta.matrix)){ ## record draw F ~ p(F | th,sig,y) (approx) f <- fbar.H(xgrid,theta.matrix[iter,]) lines(xgrid,f,col=iter,lty=3) fgrid <- rbind(fgrid,f) } ## add overall average (= posterior mean) to the plot fbar <- apply(fgrid,2,mean) lines(xgrid,fbar,lwd=3,col=2) }
Posterior_G<- function(n, alpha, Gbase, theta){ }
Approximation 1st order
$$\theta_{n} \mid \theta_{n-1},..\theta_1 \sim \frac{k_n \alpha}{n} H + \frac{1}{n} \sum_i^{k_n} (n_i - \alpha)\delta_{\theta^*_j}$$
Approximation 2nd order
$$\theta_{n} \mid \theta_{n-1},..\theta_1 \sim \frac{k_n \alpha + \beta_n}{n + \beta_n} H + \frac{1}{n + \beta_n} \sum_i^{k_n} (n_i - \alpha)\delta_{\theta^*_j}$$
For PY process $\beta_n = \theta$
polya_urn_model <- function(base_measure, N_ball, alpha) { balls <- c() for (i in 1:N_ball) { if (runif(1) < alpha / (alpha + length(balls))) { #Add a new ball color. new_color <- base_measure() balls <- c(balls, new_color) } else { #Pick out a ball from the urn, and add back a # ball of the same color ball <- balls[sample(1:length(balls), 1)] balls <- c(balls, ball) } } balls } alpha_vect <- c(1, 10, 100, 1000) N_urns <- 3 sigma2 <- c(1,.4,.2) N_draws <- 100 N_xaxis <- 200 x_axis <- seq(-3, 3, length = N_xaxis) result <- NULL for (alpha in alpha_vect) { PU <- polya_urn_model(function() rnorm(1), N_draws, alpha) for (u in 1:N_urns) { res <- mapply(function(mean) dnorm(x_axis, rep(mean, N_xaxis), rep(sigma2[u], N_xaxis)), PU) res <- apply(res, 1, mean) new_draw <- cbind(res, x_axis, alpha, sigma2[u]) result <- rbind(result, new_draw) } } result <- as.data.frame(result) names(result) <- c("density", "x", "alpha", "sigma2") DP_mixt <- qplot(data = result, y = density, x = x, geom = c("line", "area")) + facet_grid(alpha ~ sigma2, labeller = label_bquote(rows = alpha == .(alpha), cols = sigma^2 == .(sigma2))) + aes(color = as.factor(alpha)) + theme(legend.position = "none") DP_mixt
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.