Instantly share code, notes, and snippets.

# unaoya/beta.R

Created December 29, 2018 12:42
Show Gist options
• Save unaoya/109c1839a6327e106ad75a2a1148e3ea to your computer and use it in GitHub Desktop.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 import numpy as np import matplotlib.pyplot as plt print(np.random.rand(10)) def ord_sampling(N,n,p): sample = np.random.rand(N,n) sample_ord = np.sort(sample, axis=1)[:,p-1] fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.hist(sample_ord, bins=30) fig.show() N = int(1e+5) n = 10 p = 3 ord_sampling(N,n,p) def compare(N, n, p, fun): sample = np.random.rand(N,n) sample_ord = np.sort(sample, axis=1)[:,p-1] fig = plt.figure() ax = fig.add_subplot(1,1,1) x = np.linspace(0,1,100) ax.plot(x, fun(x)) ax.hist(sample_ord, bins=30, density=True) fig.show() n = 2 p = 1 fun = lambda x : 2 - 2 * x compare(N, n, p, fun)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 library(dplyr) library(ggplot2) set.seed(123) runif(10) runif(10) ord_sampling <- function(N,n,p){ x <- c() for(i in 1:N){ x[i] <- sort(runif(n))[p] } data.frame(x=x) %>% ggplot(aes(x=x)) + geom_histogram(binwidth = 0.05, boundary = 0, fill = "royal blue") + labs(title = paste("n =",as.character(n),", p = ",as.character(p))) } N <- 1e+5 n <- 10 p <- 3 ord_sampling(N,n,p) n <- 10 p <- 1 ord_sampling(N,n,p) n <- 10 p <- 2 ord_sampling(N,n,p) n <- 10 p <- 4 ord_sampling(N,n,p) n <- 10 p <- 5 ord_sampling(N,n,p) compare <- function(N,n,p,fun){ x <- c() for(i in 1:N){ x[i] <- sort(runif(n))[p] } data.frame(x=x) %>% ggplot(aes(x=x)) + geom_histogram(aes(y=..density..), fill = 'royal blue', binwidth = 0.05, boundary = 0) + stat_function(fun=fun) + labs(title = paste("n =",as.character(n),", p = ",as.character(p))) } N <- 1e+5 n <- 2 p <- 1 fun <- function(x) 2*(1-x) compare(N,n,p,fun)