In apeterson91/fhmm:

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(fhmm)
library(tidyverse)

set.seed(3431)
R <- 5
d <- seq(from=0,to=R,by=0.01)
f_1 <- function(x) 3*((1/2)*dbeta(x/R,1,8) + (1/2)*dbeta(x/R,6,1))
f_2 <- function(x) 3*((1/5)*dbeta(x/R,3,2) + (2/3)*dbeta(x/R,3,1) + (2/15)*dbeta(x/R,1,1))
f_3 <- function(x) 3*((1/2)*dbeta(x/R,8,2) + (1/2)*dbeta(x/R,30,50))
pltdf <- tibble(Distance = d,
       Intensity = f_1(d),
       Cluster = 1) %>% 
    rbind(.,tibble(Distance = d,
                   Intensity = f_2(d),
                   Cluster = 2)) %>% 
    rbind(.,tibble(Distance = d,
                   Intensity = f_3(d),
                   Cluster = 3))
pltdf %>% ggplot(aes(x=Distance,y=Intensity)) + geom_line() + 
    facet_wrap(~Cluster) + theme_bw() + 
    theme(strip.background = element_blank()) + 
    ggtitle("Intensity Functions")

num_schools <- 50
schools_1 <- BEFcluster::rnhpp(nsim = num_schools,lambda = f_1,
                   interval = c(0,R),seed = 3431,
                   max =max(f_1(d)))
schools_2 <- BEFcluster::rnhpp(nsim = num_schools,lambda = f_2,interval = c(0,R),seed = 3431,max = max(f_2(d)))
schools_3 <- BEFcluster::rnhpp(nsim = num_schools,lambda = f_3 ,interval = c(0,R),seed = 3431,max = max(f_3(d)))
school_df <- as_tibble(schools_1)
school_df <- rbind(school_df,as_tibble(schools_2) %>% mutate(sim_id = sim_id + num_schools)) %>% 
    mutate(density = ifelse(sim_id<=num_schools,1,2) )
school_df <- rbind(school_df,schools_3 %>% mutate(sim_id = sim_id + 2*num_schools, density=3))
school_df %>% ggplot(aes(x=event_times)) + geom_density() + facet_wrap(~density) + 
    theme_bw() + theme(strip.background = element_blank()) + xlab("Distance")

r <- school_df %>% arrange(sim_id) %>% 
    select(event_times) %>% pull()

n_j <- school_df %>% group_by(sim_id) %>% count() %>% 
    ungroup() %>% mutate(start = (cumsum(n) ) ) %>% 
    mutate(start_ = replace_na(dplyr::lag(start),0) ) %>% select(-start) %>% 
    rename(start=start_,go =n) %>% 
    select(start,go) %>% as.matrix()

fit <- fhmm(r = r, n_j = n_j,
            L = 4, K = 3,
            mu_0 = 0, kappa_0 = 1, ## Normal Mean Base Measure Hyperparameters
            sigma_0 = 1, nu_0 = 1, ## Inverse Chi square prior Hyperparameters
            iter_max = 1.5E4,
            warm_up = 1.3E4,
            thin = 1,
            seed = 34143)