Лабораторная 4

Задание 1. Моделирование одномерной случайной величины

f <- function(x) ifelse(x <= 4, 2*x+3, -x +15 )
df <- data.frame(x=1:10, y=f(1:10))
ggplot(data=df, aes(x=x, y=y)) + geom_line() + geom_point()

norm_k <- 72
f_norm <- function(x) f(x)/norm_k
ggplot(data=df, aes(x=x, y=f_norm(x))) + geom_line() + geom_point()

STEP <- 0.01
x <- seq(1, 10, STEP)
y <- f_norm(x)
ss <- Reduce(function(a,b) c(a, tail(a,1) + b*STEP), y[2:length(y)], c(0))
df <- data.frame(x=x, y=y, ss=ss)

Сгенерируем 10000 значений с помощью нашего генератора и построим гистограмму выборочной плотности рампределения полученной выбоки.

f_gen = function(x) df[df$ss >= x, 'x'][[1]]
gen_dataset <- data.frame(x=unlist(Map(f_gen, runif(10000))))
ggplot(data=gen_dataset) + geom_histogram(aes(x=x,..density..), binwidth = 0.4, fill='green', col='red', alpha = .2) + geom_line(data=df, aes(x=x, y=f_norm(x))) + geom_point(data=df, aes(x=x, y=f_norm(x)))

Задание 2. Моделирование двумерной случайной величины

Функция плотности определена на прямоугольнике [1, 5] x [0, 2]. Прямоугольник разделен диагональю, проходящей через точки (1, 0) и (5, 2), на два треугольника. Вероятность попадания в левый треугольник равна 0.25, а в правый треугольник 0.75.

inverseGen <- function (N, fun, xlim){ 
  STEP <- 0.01
  x <- seq(xlim[1], xlim[2], STEP)
  y <- unlist(Map(fun, x))
  ss <- Reduce(function(a,b) c(a, tail(a,1) + b*STEP), y[2:length(y)], c(0))
  df <- data.frame(x=x, y=y, ss=ss)
  f_gen = function(x) tail(df[df$ss < x, 'x'],1)[[1]]
  data.frame(x=unlist(Map(f_gen, runif(N))))
}

f <- function(x) ifelse(x[2] >= 0.5*x[1] - 0.5, 0.25, 0.75 )
res <- adaptIntegrate(f, lowerLimit = c(1, 0), upperLimit = c(2, 5), tol = 1e-3)
k_norm = res$integral
f_norm <- function(x) f(x)/k_norm
takeIntergal = function(b1,b2) adaptIntegrate(f_norm, lowerLimit = c(1, 0), upperLimit = c(b1, b2), tol = 1e-3)$integral
s1<- seq(1,2,length.out = 20)
s2<- seq(0,5,length.out = 20)
CDF_Table <- adply(expand.grid(x=s1, y=s2), 1, function(row) data.frame(ss=takeIntergal(row$x, row$y)))
CDF_Table_sort <- CDF_Table[order(CDF_Table$ss), ]
f_gen = function(x) CDF_Table_sort[CDF_Table_sort$ss >= x, ][1, c('x', 'y')]

N_SAMPLES = 10000
#gen_dataset <- data.frame(matrix(unlist(Map(f_gen, runif(N_SAMPLES))), nrow=N_SAMPLES, byrow = TRUE))
#colnames(gen_dataset) <- c('x', 'y')
generator2 <- function(x){
  x <- inverseGen(1, function(x) x/16 + 1/16, c(1,5))[[1]]
  y <- inverseGen(1, function(p) ifelse(p >= 0.5*x - 0.5, 1/(x+1), 3/(x+1) ), c(0,2))[[1]]
  data.frame(x=x, y=y)
}
gen_dataset <- data.frame(matrix(unlist(Map(generator2, 1:10000)), ncol=2, byrow=TRUE))
colnames(gen_dataset) <- c('x', 'y')

Попробуем графически изобразить полученный результат. Построим гистограмму на плоскости. Насыщенностью цвета будем обозначить количество попаданий в карманы.