Как сгенерировать данные для гауссовских распределений в этих двух сценариях в R?

В «Элементах статистического обучения» Тибширани при сравнении моделей наименьших квадратов / линейных моделей и знания эти 2 сценария указаны:

Сценарий 1. Обучающие данные в каждом классе были сгенерированы из двумерных распределений Гаусса с некоррелированными компонентами и разными средними значениями.

Сценарий 2: обучающие данные в каждом классе были получены из смеси 10 гауссовских распределений с низкой дисперсией, при этом отдельные средние сами были распределены по Гауссу.

Идея состоит в том, что первая лучше подходит для моделей наименьших квадратов / линейных моделей, а вторая - для моделей типа knn (тех, которые отличаются от того, что я понимаю, поскольку knn учитывает ближайшие точки, а не все точки).

Как в R смоделировать данные для обоих сценариев?

Конечная цель состоит в том, чтобы иметь возможность воспроизвести оба сценария, чтобы доказать, что эффективно 1-й сценарий лучше объясняется линейной моделью, чем 2-й.

Спасибо!


person scc    schedule 04.09.2015    source источник
comment
Хотя у этого вопроса есть статистический компонент - было бы интересно объяснить различие между двумя сценариями - его акцент на R заставил его собирать чисто программно-ориентированные ответы, что сделало его вне темы CV.   -  person whuber    schedule 04.09.2015
comment
На самом деле я хочу иметь возможность воспроизвести два сценария, чтобы лучше понять различия и то, как они применяются к линейным моделям.   -  person scc    schedule 04.09.2015


Ответы (3)


Это может быть сценарий 1

library(mvtnorm)

N1 = 50
N2 = 50
K = 2

mu1 = c(-1,3)
mu2 = c(2,0)

cov1 = 0
v11 = 2
v12 = 2
Sigma1 = matrix(c(v11,cov1,cov1,v12),nrow=2)

cov2 = 0
v21 = 2
v22 = 2
Sigma2 = matrix(c(v21,cov2,cov2,v22),nrow=2)

x1 = rmvnorm(N1,mu1,Sigma1)
x2 = rmvnorm(N2,mu2,Sigma2)

Это может быть кандидатом для моделирования из гауссовой смеси:

BartSimpson <- function(x,n = 100){ 
   means <- as.matrix(sort(rnorm(10)))
   dens <- .1*rowSums(apply(means,1,dnorm,x=x,sd=.1)) 
   rBartSimpson <- c(apply(means,1,rnorm,n=n/10,sd=.1))
   return(list("thedensity" = dens,"draws" = rBartSimpson))
}

x <- seq(-5,5,by=.01)

plot(x,BartSimpson(x)$thedensity,type="l",lwd=4,col="yellow2",xlim=c(-4,4),ylim=c(0,0.6))
person Christoph Hanck    schedule 04.09.2015
comment
Спасибо :) Почему сорт? - person scc; 04.09.2015
comment
вполне может быть лишним, я что-то поправлял и пробовал вот это, но не в этом суть, как оказалось - person Christoph Hanck; 04.09.2015

В приведенном ниже коде я сначала создаю 10 различных средств классов, а затем использую средства для извлечения случайных значений из этих средств. Код для этих двух сценариев идентичен, но вам придется отрегулировать дисперсию внутри и между классами, чтобы получить желаемые результаты.

Сценарий 1:

Здесь вы хотите сгенерировать 10 классов с разными средствами (я предполагаю, что средства соответствуют двумерному гауссовскому распределению). Разница между классами намного меньше, чем разница внутри классов.

library(MASS)
n <- 20
# subjects per class
classes <- 10
# number of classes
mean <- 100
# mean value for all classes
var.between <- 25
# variation between classes
var.within <- 225
# variation within classes
covmatrix1 <- matrix(c(var.between,0,0,var.between), nrow=2)
# covariance matrix for the classes
means <- mvrnorm(classes, c(100,100), Sigma=covmatrix1)
# creates the means for the two variables for each class using variance between classes
covmatrix2 <- matrix(c(var.within,0,0,var.within), nrow=2)
# creates a covariance matrix for the subjects
class <- NULL
values <- NULL
for (i in 1:10) {
  temp <- mvrnorm(n, c(means[i], means[i+classes]), Sigma=covmatrix2)
  class <- c(class, rep(i, n))
values <- c(values, temp)
}
# this loop uses generates data for each class based on the class means and variance within classes
valuematrix <- matrix(values, nrow=(n*classes))
data <- data.frame (class, valuematrix)
plot(data$X1, data$X2)

В качестве альтернативы, если вы не заботитесь об указании дисперсии между классами и не хотите никакой корреляции внутри классов, вы можете просто сделать это:

covmatrix <- matrix(c(225, 0, 0, 225), nrow=2)
# specifies that the variance in both groups is 225 and no covariance
values <- matrix(mvrnorm(200, c(100,100), Sigma=covmatrix), nrow=200)
# creates a matrix of 200 individuals with two values each.

Сценарий 2:

Здесь единственное отличие состоит в том, что различия между классами больше, чем различия внутри классов. Попробуйте поменять значение переменной var.bet между примерно 500 и переменной var.ithin до 25, и вы увидите четкую кластеризацию на диаграмме рассеивания:

n <- 20
# subjects per class
classes <- 10
# number of classes
mean <- 100
# mean value for all classes
var.between <- 500
# variation between classes
var.within <- 25
# variation within classes
covmatrix1 <- matrix(c(var.between,0,0,var.between), nrow=2)
# covariance matrix for the classes
means <- mvrnorm(classes, c(100,100), Sigma=covmatrix1)
# creates the means for the two variables for each class using variance between classes
covmatrix2 <- matrix(c(var.within,0,0,var.within), nrow=2)
# creates a covariance matrix for the subjects
class <- NULL
values <- NULL
for (i in 1:10) {
  temp <- mvrnorm(n, c(means[i], means[i+classes]), Sigma=covmatrix2)
  class <- c(class, rep(i, n))
values <- c(values, temp)
}
# this loop uses generates data for each class based on the class means and variance within classes
valuematrix <- matrix(values, nrow=(n*classes))
data <- data.frame (class, valuematrix)
plot(data$X1, data$X2)

График должен подтверждать кластеризацию данных.

Надеюсь это поможет!

person JonB    schedule 04.09.2015
comment
Я ожидал, что сценарий 1 будет лучше объяснен применением линейной модели (model = lm (X2 ~ X1, data = data)), чем сценарий 2, но это не похоже на правду, хотя я вижу, что в сценарии 2 данные сгруппированы. Даже изменение var.betweeen и var.within в сценарии 1 на низкие значения (например, 2), похоже, не всегда обеспечивает более высокие значения p для линейной модели в сценарии 1. Вы знаете, почему? - person scc; 04.09.2015
comment
Извините, но я мало что знаю о моделях Knn. В основном я пытался помочь с созданием данных. - person JonB; 07.09.2015

С помощью обоих ответов здесь я использовал это:

mixed_dists = function(n, n_means, var=0.2) {
    means = rnorm(n_means, mean=1, sd=2)
    values <- NULL
    class <- NULL
    for (i in 1:n_means) {
        temp <- rnorm(n/n_means, mean=means[i], sd=0.2)
        class <- c(class, rep(i, n/n_means))
        values <- c(values, temp)
    }
    return(list(values, class));
}

N = 100

#Scenario 1: The training data in each class were generated from bivariate Gaussian distributions 
#with uncorrelated components and different means.
scenario1 = function () { 
    var = 0.5
    n_groups = 2
    m = mixed_dists(N, n_groups, var=var)
    x = m[[1]]
    group = m[[2]]
    y = mixed_dists(N, n_groups, var=var)[[1]]
    data = matrix(c(x,y, group), nrow=N, ncol=3)
    colnames(data) = c("x", "y", "group")
    data = data.frame(data)
    plot(x=data$x,y=data$y, col=data$group)
    model = lm(y~x, data=data)
    summary(model) 
}



#Scenario 2: The training data in each class came from a mixture of 10 
#low-variance Gaussian distributions, with individual means themselves
#distributed as Gaussian.
scenario2 = function () {
    var = 0.2 # low variance
    n_groups = 10
    m = mixed_dists(N, n_groups, var=var)
    x = m[[1]]
    group = m[[2]]
    y = mixed_dists(N, n_groups, var=var)[[1]]
    data = matrix(c(x,y, group), nrow=N, ncol=3)
    colnames(data) = c("x", "y", "group")
    data = data.frame(data)
    plot(x=data$x,y=data$y, col=data$group)
    model = lm(y~x, data=data)
    summary(model)
}
# scenario1()
# scenario2()

Таким образом, в основном данные в сценарии 1 четко разделены на 2 класса, а данные в сценарии 2 имеют около 10 кластеров и не могут быть четко разделены с помощью прямой линии. Действительно, запустив линейную модель для обоих сценариев, можно увидеть, что в среднем она будет лучше применяться к сценарию 1, чем к сценарию 2.

person scc    schedule 06.09.2015