abikoushi

library(dplyr)
library(cowplot)
library(deSolve)
set.seed(1234)
SIRsim <- function(beta,gamma){
mat <- matrix(0,12,12)
mat[sample(2:11,1),sample(2:11,1)] <- 1
out <- vector("list",100)
out[[1]] <- mat[2:11,2:11]
for(t in 2:100){

abikoushi

StatMA <- ggproto("StatMA", Stat, 
                  required_aes = c("x", "y"),
                  
                  compute_group = function(data, scales,windowsize) {
                    grid <- data.frame(x = data$x)
                    grid$y <- stats::filter(data$y, rep(1,windowsize)) / windowsize
                    grid
                  }
)

abikoushi

StatBinomCI <- ggproto("StatBinomCi", Stat, 
                       required_aes = c("x", "numerator", "denominator"),
                       
                       compute_group = function(data, scales, conf.level) {
                         CI <-t(mapply(function(x,n)binom.test(x=x, n=n, conf.level = conf.level)$conf.int,
                                       data$numerator,
                                       data$denominator))
                         grid <-as.data.frame(CI)
                         colnames(grid) <- c("ymin","ymax")
                         grid$x <- data$x

abikoushi

GeomLasagna <- ggproto("GeomLasagna", GeomRect,
                       extra_params = c("na.rm", "or_more"),
                       required_aes = c("x", "y","fill"),
                       setup_data = function(data, params) {
                         if (!is.na(params$or_more)) {
                           y2 <- with(data, ifelse(data$y >= params$or_more, paste0(params$or_more,"~"),data$y))
                           lev <-sort(unique(data$y[data$y < params$or_more]))
                           lev <-c(lev, paste0(params$or_more,"~"))
                           y2 <- factor(y2, levels = lev)
                           out=aggregate(fill~x+y2+PANEL+group,sum,data=data)

abikoushi

GeomDrilldownLabel <- ggproto("GeomDrilldownLabel", GeomLabel,
                             extra_params = c("na.rm"),
                             required_aes = c("x", "y","label"),
                             setup_data = function(data, params) {
                               as.data.frame(dplyr::mutate(dplyr::group_by(data,x,PANEL),y=cumsum(y)-y/2))
                             }
                             
)
geom_drilldown_label <- function(mapping = NULL, data = NULL, stat = "identity", position = "identity",
                                ..., na.rm = TRUE, show.legend = NA, inherit.aes = TRUE){

abikoushi

pathanal <-function(mcf_gadata,dimension="basicChannelGroupingPath"){
  mcf_gadata <-mcf_gadata %>% 
    mutate(totalConversions=as.numeric(totalConversions))

  path <- with(mcf_gadata,sapply(eval(parse(text=dimension)),function(x)strsplit(x," > ")))
  
  node <-unname(unlist(path))
  node <-gsub("CLICK|:|NA","",node)
  pathlength <- sapply(path,length)
  n <- length(path)

abikoushi

JacksonPollock <- function(n=5,r1=5,r2=0.5,drip=100,col="black"){
  old<-par(no.readonly = TRUE)
  par(mai=c(0,0,0,0))
  plot(c(0,1),c(0,1),type="n")
  pos<-locator(n)
  points(pos$x,pos$y,cex=rexp(n,1/r1),pch=16,col=col)
  d <-rpois(1,drip)
  u <-runif(d)
  for(i in 2:n){
  points(u*pos$x[i-1]+(1-u)*pos$x[i],u*pos$y[i-1]+(1-u)*pos$y[i],

abikoushi

NMF <- function(X,M=2,a_W=1,b_W=1,a_H=1,b_H=1,tol=1e-2,maxit=10000,
                seed=1234){
  set.seed(seed)
  D <-nrow(X)
  N <-ncol(X)
  W<-matrix(rgamma(D*M,a_W,b_W),D,M)
  H<-matrix(rgamma(M*N,a_H,b_H),M,N)
  
  logW0 <-log(W)
  logH0 <-log(H)

abikoushi

make_q_df <- function(sf,p=c(0.1,0.25,0.5,0.75,0.9)){
  sf_q <-quantile(sf,probs = p)
  if(is.null(sf$strata)){
    q_df<-sf_q$quantile %>% 
      as_data_frame() %>% 
      set_names(paste0("q",p))
    q_df[is.na(q_df)]<-Inf
  }else{
    q_df<-sf_q$quantile %>% 
      as_data_frame() %>% 

abikoushi

make_df4hist <-function(sf,bw){
  nb <- ceiling(max(sf$time)/bw)
  if(is.null(sf$strata)){
    sfdf <- data.frame(time=sf$time,survival=sf$surv,
                       stringsAsFactors = FALSE)
    df4hist <- sfdf %>% 
      mutate(time=cut(time,0:nb*bw)) %>% 
      group_by(time) %>% 
      summarise(TP=diff(range(survival))) %>% 
      mutate(density=TP/bw) %>%