B站视频BV1Us4y1s7E1以及BV1wY4y1d7pB所用到的代码

6stars.r

library(ggplot2)
library(dplyr)

prob_per_trial = function(x){
  p_ori = 0.02
  p_add = 0.02
  if (x==0){
    return(1)
  }else if (x<=50){
    return(p_ori)
  }else if (x<=99){
    return(p_ori+(x-50)*p_add)
  }else{
    return(prob_per_trial(x %% 99))
  }
}
prob_per_trial2 = function(x){
  p_ori = 0.006
  p_add = 0.06
  if (x<76){
    return(p_ori)
  }else if (x<90){
    return(p_ori+(x-75)*p_add)
  }else if (x==90 | x==0){
    return(1)
  }else{
    return(prob_per_trial2(x %% 90))
  }
}

pdf_x = function(x0){
  if(x0>1){
    p_agg_mul = 1
    for (i in 1:(x0-1)){
      p_agg_mul = p_agg_mul*(1-prob_per_trial(i))
    }
    p_agg_mul = p_agg_mul*prob_per_trial(x0)
    return(p_agg_mul)
  }else if(x0==1){
    return(prob_per_trial(x0))
  }
  else{
    return(NULL)
  }
}

pdf_x2 = function(x0){
  if(x0>1){
    p_agg_mul = 1
    for (i in 1:(x0-1)){
      p_agg_mul = p_agg_mul*(1-prob_per_trial2(i))
    }
    p_agg_mul = p_agg_mul*prob_per_trial2(x0)
    return(p_agg_mul)
  }else if(x0==1){
    return(prob_per_trial2(x0))
  }
  else{
    return(NULL)
  }
}


x = 1:300
p_val = sapply(x,pdf_x)
p_val2 = sapply(x,pdf_x2)
df = data.frame(x=x,p1=p_val,p2=p_val2) %>% mutate(cdf1=cumsum(p1),cdf2=cumsum(p2))

df %>% ggplot(aes(x=x))+geom_line(aes(y=p1),color='red')+geom_line(aes(y=cdf1),color='red')+geom_line(aes(y=p2),color='blue')+geom_line(aes(y=cdf2),color='blue')+ylim(0,1)+xlim(1,90)+ylab('probability')+xlab('trials before first success')

expectation = sum(x*p_val)

monte_carlo.r

library(ggplot2)
library(dplyr)

prob_per_trial = function(x){
  p_ori = 0.02
  p_add = 0.02
  if (x==0){
    return(1)
  }else if (x<=50){
    return(p_ori)
  }else if (x<=99){
    return(p_ori+(x-50)*p_add)
  }else{
    return(prob_per_trial(x %% 99))
  }
}

gacha_trial = function(trial_number){
  if(trial_number %% 1000==1){
    print(paste('trial number',trial_number))
  }
  combo = 0
  for(t in 1:300){
    
    # draw a card
    combo = combo+1
    p_6star = prob_per_trial(combo)
    rst=sample(c(0,1),size=1,replace=TRUE,prob=c(1-p_6star,p_6star))
    
    if(rst==1){
      # get a 6 star, reset combo  
      combo=0
      final_rst = sample(c(0,1),size=1,replace=TRUE,prob=c(1-0.35,0.35))
      if(final_rst==1){
        # get you b**** !
        # stop this trial and return total draws
        return(t)
      }
    }
  }
  # if after 300 draws, still no luck
  return(t)
}

wai_trial = function(trial_number){
  p_prev = 0.03*(15/68)
  if(trial_number %% 1000==1){
    print(paste('trial number',trial_number))
  }
  combo = 0
  for(t in 1:300){
    
    # draw a card
    combo = combo+1
    p_6star = prob_per_trial(combo)
    rst=sample(c(0,1),size=1,replace=TRUE,prob=c(1-p_6star,p_6star))
    
    if(rst==1){
      # get a 6 star, reset combo  
      combo=0
      final_rst = sample(c(0,1),size=1,replace=TRUE,prob=c(1-0.35,0.35))
      if(final_rst==1){
        # NOPE
        return(-1)
      }else{
        wai_rst = sample(c(0,1),size=1,replace=TRUE,prob=c(1-p_prev,p_prev))
        if(wai_rst==1){
          return(t)
        }
      }
    }
  }
  # if after 300 draws, still no luck
  return(t)
}

draws = sapply(1:100000,gacha_trial)
draws2 = sapply(1:100000,wai_trial)
df = data.frame(1:100000,draw1=draws,draw2=draws2)
draws2_filtered = (df %>% filter(draws2>0))$draw2
density_1 = density(draws)
density_2 = density(draws2_filtered)

# histogram comparison
ggplot(data=NULL)+geom_histogram(aes(x=draws),bins=100,fill='gray')+geom_histogram(aes(x=draws2),bins=100)+xlim(1,300)
# density plot comparison
ggplot()+geom_density(aes(x=draws),fill='gray')+geom_density(aes(x=draws2_filtered),fill='black')