Дискретное приближение к двумерному нормальному распределению

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

введите описание изображения здесь

Вот что я сделал до сих пор.

library(mvtnorm)
library(graphics)

euclide = function(x,y){sqrt(x^2+y^2)}
maxdist = 40
sigma = diag(2)
m = matrix(0,ncol=maxdist*2 + 1, nrow=maxdist*2 + 1)
for (row in -maxdist:maxdist){
    for (col in -maxdist:maxdist){
        if ( euclide(abs(row), abs(col)) < maxdist ){
            lower = c(row-0.5, col-0.5)
            upper = c(row+0.5, col+0.5)
            p = pmvnorm(lower = lower , upper = upper, mean = c(0,0), sigma = sigma)    
        } else {
            p = 0
        }
        m[row + maxdist + 1,col + maxdist + 1] = p
    }
}
m = m[rowSums(m)!=0,colSums(m)!=0]
contour(m, levels = exp(-20:0), xlim=c(0.3,0.7), ylim=c(0.3,0.7))

введите описание изображения здесь

Работает нормально. Однако это довольно медленно (для больших maxdist), и я надеюсь сократить время его вычислений. Но это не моя главная проблема ...

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


r statistics distribution normal-distribution
person Remi.b    schedule 04.07.2015    source источник
comment
Вы слышали о Tauchen? А генерирование суставов нормальных вариаций от N(0,1) розыгрышей? (например, слайд 269-272 здесь) Или вы работаете на этом упражнении для другой конкретной цели?   -  person MichaelChirico    schedule 04.07.2015
comment
Могу я спросить, с какой целью вы хотите использовать строящуюся матрицу?   -  person Robert Dodier    schedule 04.07.2015
comment
Это ядро ​​рассредоточения (вероятность рассредоточения к ячейке, находящейся на заданном расстоянии) с немного более сложным набором деталей, поскольку внутри ячеек есть подъячейки. С вычислительной точки зрения имитировать людей, находящихся в клетках и спаривающихся случайным образом, быстрее, чем в пространственном континууме. Вот почему я хотел это дискретное приближение.   -  person Remi.b    schedule 04.07.2015
comment
@ Remi.b Хорошо, спасибо за информацию. Мне интересно, будет ли поэтому важная операция сверткой (этой переходной матрицы с распределением по текущему состоянию, чтобы получить распределение по следующему состоянию). Если это так, рассмотрите возможность вычисления свертки с помощью быстрого преобразования Фурье для большей эффективности.   -  person Robert Dodier    schedule 05.07.2015

Ответы (2)


arrow_upward
2
arrow_downward

Вот простая реализация. Как @DanielJohnson говорит, что вы можете просто использовать одномерную нормальную форму cdf, но она должна быть такой же, как при использовании pmvnorm, как показано ниже. Версия, использующая pnorm, намного быстрее.

## Choose the matrix dimensions
yticks <- xticks <- seq(-3, 3, length=100)
side <- diff(yticks[1:2])  # side length of squares
sigma <- diag(2)               # standard devs. for f2
mu <- c(0,0)                # means

## Using pnorm
f <- Vectorize(function(x, y, side, mu1, mu2, s1, s2)
    diff(pnorm(x+c(-1,1)*side/2, mu1, s1)) * diff(pnorm(y+c(-1,1)*side/2, mu2, s2)),
    vec=c("x", "y"))

## Using pmvnorm
f2 <- Vectorize(function(x, y, side, mu, sigma)
    pmvnorm(lower=c(x,y)-side/2, upper=c(x,y)+side/2, mean=mu, sigma=sigma),
                vec=c("x", "y"))

## get prob. of squares, mu are means, s are standards devs.
mat <- outer(xticks, yticks, f, side=side, mu1=0, mu2=0, s1=1,s2=1)
mat2 <- outer(xticks, yticks, f2, side=side, mu=mu, sigma=sigma)

## test equality
all(abs(mat2-mat) < 1e-11)  # TRUE
all.equal(mat2, mat)        # TRUE

## See how it looks
library(lattice)
persp(mat, col="lightblue", theta=35, phi=35, shade=0.1)

введите описание изображения здесь

person Rorschach    schedule 04.07.2015

arrow_upward
3
arrow_downward

Я не человек R, но уверен, что для нормального распределения есть функция CDF. Если вам нужна буквально матрица вероятностей попадания в каждый квадрат, мы можем использовать эту функцию CDF, чтобы получить ответ. Поскольку двумерное нормальное распределение имеет независимые маргинальные распределения, вопрос здесь сводится к заданию двух вопросов для каждого квадрата, описываемого положениями оси [x_left, x_right] и [y_left, y_right]:

  1. Какова вероятность того, что нормальная одномерная случайная величина окажется в интервале [x_left, x_right]?
  2. Какова вероятность того, что другая независимая одномерная нормальная случайная величина окажется в интервале [y_left, y_right]?

Поскольку они независимы, полная вероятность квадрата равна:

P = (CDF(x_right) - CDF(x_left))*(CDF(y_right) - CDF(y_left))

Это точный ответ, поэтому время вычислений не должно быть проблемой!

РЕДАКТИРОВАТЬ: Я также должен сказать, что вы можете выбрать сетку с большим количеством отметок на каждой оси, близким к нулю, чтобы получить желаемое разрешение. Приведенная выше формула вероятности для каждого квадрата все еще верна.

person Daniel Johnson    schedule 04.07.2015