samples.R (1451B)
1 # Three simple functions 2 DEFAULT_N = 1000000 3 normal <- function (mean, std, n=DEFAULT_N){ 4 return(rnorm(n, mean, std)) 5 } 6 7 lognormal <- function(meanlog, sdlog, n=DEFAULT_N){ 8 return(rlnorm(n, meanlog = meanlog, sdlog = sdlog)) 9 } 10 11 to <- function(low, high, n=DEFAULT_N){ 12 normal95confidencePoint = 1.6448536269514722 13 logLow = log(low) 14 logHigh = log(high) 15 meanlog = (logLow + logHigh)/2 16 sdlog = (logHigh - logLow) / (2.0 * normal95confidencePoint) 17 return(lognormal(meanlog, sdlog, n)) 18 } 19 20 mixture <- function(samples_list, weights_array, n=DEFAULT_N){ # note that this takes a list, not an array 21 normalized_weights = weights_array/sum(weights_array) 22 cummulative_sums = cumsum(normalized_weights) 23 helper_probs = runif(n) 24 results = vector(mode='numeric', length=n) 25 for(i in c(1:n)){ 26 helper_which_list = which(cummulative_sums > helper_probs[i]) 27 # helper_loc = ifelse(is.na(helper_which_list[1]), 1, helper_which_list[1]) 28 if(is.na(helper_which_list[1])){ 29 print("This should never happen") 30 } 31 helper_loc = helper_which_list[1] 32 target_samples = samples_list[[helper_loc]] 33 result = sample(target_samples, 1) 34 results[i] = result 35 } 36 return(results) 37 } 38 39 # Example 40 p_a = 0.8 41 p_b = 0.5 42 p_c = p_a * p_b 43 44 dists = list(c(0), c(1), to(1, 3), to(2, 10)) 45 # print(dists) 46 weights = c((1 - p_c), p_c/2, p_c/4, p_c/4) 47 # print(weights) 48 result = mixture(dists, weights) 49 mean_result = mean(result) 50 print(mean_result)