Можно ли восстановить распределение?

[Neytrall]

Здравствуйте,

У меня немного нестандартный вопрос. Допустим у меня есть вектор из распределения Гаусса, но проблема в том, что некоторые числа из вектора равны нулю (алгоритм обнулял все числа, которые были на каком-то расстояние от нуля). Гистограмма выглядит так:


image hosting websites

Мне известно математическое ожидание изначального распределения (допустим $\mu$), то что это распределение Гаусса и что числа не равные нулю относятся к хвостам этого распределения. Есть ли какой-нибудь алгоритм, способный "восстановить" это распределение? То бишь скорректировать эти нули, так чтобы они заполнили пространство между хвостами и получившийся вектор имел нормальное распределение с заданным математическим ожиданием и оценкой среднеквадратического отклонения.
Мне не нужно восстановить точно каждый ноль.

[Александрович]

Пусть всего $n$ значений с.в., а левее $x=-0,2$ попало $m$ значений, тогда $x=-0,2$ это $100 \cdot \frac{m}{n}$%-ная квантиль нормального распределения $N[\mu;\sigma]$. Тоже самое можно сказать о значении $x=0,2$.

[Neytrall]

Спасибо, если я правильно понял то:
если $X\sim N(\mu, \sigma)$, то $F(\frac{x-\mu}{\sigma}<neg)=\frac{m_{neg}}{N}$ и $F(\frac{x-\mu}{\sigma}>pos)=\frac{m_{pos}}{N}$, где $neg$ это максимальное значение на отрицательной стороне (примерно -0.2) и $pos$ это минимальное значение на положительной стороне(примерно 0.2), то мне надо решить систему уравнений:

$$\int_{-\infty}^{neg} f\left(\frac{x-\hat{\mu}}{\sigma}\right) dx=\frac{m_{neg}}{N}$$

и

$$\int_{pos}^{\infty} f\left(\frac{x-\hat{\mu}}{\sigma}\right) dx=\frac{m_{pos}}{N}$$

для $\sigma$. Так?
Есть ли в R функция способная быстро и безболезненно получить ответ?

[Neytrall]

Я попробовал multiroot():

tryfunc=function(sigma){
  funcd=function(x){dnorm(x,mean=b.mean, sd=sigma)}
  y1=integrate(funcd, lower=-Inf, upper=max.neg)$value
  y2=integrate(funcd, lower=min.pos, upper=Inf)$value
  y=y1-y2+mpos/N-mneg/N
 return(y) 
}
multiroot(start=b.sd,tryfunc)

Но код выдал ошибку:

 Error in integrate(funcd, lower = -Inf, upper = mneg) : 
  non-finite function value In addition: Warning message:
In dnorm(x, mean = b.mean, sd = s) : NaNs produced

Почему то y1 и y2 всегда равны нулю...

[realeugene]

Мне кажется, по квантилям очень плохо оценивать параметры в вашей задаче. Наиболее информативны для вас бины гистограммы вблизи порога отсечения, но именно там гистограмма искажается сильнее всего (в один бин попала часть отсчётов). Думаю, методы теориии оценивания могут вам помочь. Рассмотрите комплексную гипотезу, включающую два неизвестных параметра: дисперсию вашего нормального распределения и порог отсечения. Матожидание и, вероятно, число испытаний вы знаете. Выведите выражение для вероятности получить ту или иную гистограмму в зависимости от ваших известных и неизвестных параметров. Воспользуйтесь методом максимального правдоподобия для поиска наилучшей оценки неизвестных параметров исходя из вашей гистограммы.

Это, конечно, если у вас гистограмма получена экспериментально. Если же она посчитана аналитически для бесконечного числа испытаний, то любых двух несимметричных относительно матожидания бинов выше порога отсечения вам должно хватить для решения простого уравнения относительно дисперсии.

[Александрович]

Neytrall, система нужна для определения неизвестных мо и ско. А здесь вы получаете ско из каждого уравнения. Они должны быть одинаковыми.

[Neytrall]

Александрович в сообщении #1242326 писал(а):

Neytrall, система нужна для определения неизвестных мо и ско. А здесь вы получаете ско из каждого уравнения. Они должны быть одинаковыми.

А разве

funcd=function(x){dnorm(x,mean=b.mean, sd=sigma)}

не предполагает, что используется одно и то де ско?

[Neytrall]

Я разобрался! Просто запутался в функции. Она должна быть:

tryfunc=function(sigma){
  funcd=function(x){dnorm(x,mean=b.mean, sd=sigma)}
  y1=integrate(funcd, lower=-Inf, upper=mneg)$value
  y2=integrate(funcd, lower=mpos, upper=Inf)$value
  y=npos/n+nneg/n-y1-y2;
  c(y)
}
multiroot(start=b.sd,tryfunc)

[Александрович]

Ещё один вариант. Определить ско из максимума критерия Пирсона.

[Neytrall]

Как?

[Александрович]

Как при проверке гипотезы о виде распределения св.

[Евгений Машеров]

В целом это задача оценивания по цензурированным данным (только цензурирование не то, которое чаще всего встречается, с одной стороны, а "дырка", цензурируются данные в определённом интервале).
А в частности я бы $\chi^2$ попробовал бы. С ячейками, соответствующими бинам гистограммы (или объединённым, если наблюдений мало), а одна ячейка - все обнулённые. Поскольку по условию оценивается один параметр - можно даже простым перебором по сетке пройтись.

[Neytrall]

Я не уверен, что $\chi^2$ является оптимальным. Гистограмма это результат бутстрепа и если размер выборки достаточный, то можно предположить, что вектор должен иметь нормальное распределение. Поскольку в данном случае есть оба хвоста, то использование отсечении по $m^{-}=\max(x<0)$ и по $m^{+}=\min(x>0)$ и поиск корня по интегралам выглядит оптимальнее.
Я сейчас прогнал симуляцию, и она показала, что этот метод выдает лучший результат, чем построение доверительного интервала на основание сечения $2.5\%$ и $97.5\%$ (percentile CI).

[Евгений Машеров]

Ну, в данном частном случае, когда где-то с две трети выборки попали в "диапазон обнуления", и матожидание предполагается известным, самым простым, но при этом дающим достаточную точность может оказаться оценка, использующая только долю p попавших в "диапазон обнуления" [L,R]
$\Phi(\frac {R-\mu}\sigma)-\Phi(\frac {L-\mu}\sigma)=p$
В качестве неизвестной использовать обратную к стандартному отклонению, начальное приближение к стандартному отклонению взять, как одну шестую размаха, решать Ньютоном (впрочем, эти советы не обязательны)