плотность вероятности

[type2b]

блин, ну какая там реализация!

$$ f(u) = \frac{d}{du}\int_{\{x: \, H(x) < u\}}\frac{1}{L}\,dx = \frac{d}{du}\int_0^L \theta(u-H(x))\frac{dx}{L} = \int_0^L \delta(u-H(x)) \frac{dx}{L} = \frac{1}{L}\sum_{x: H(x)=u} \frac{1}{H'(x)} $$

Здесь $\theta$ -- функция Хевисайда, которая равна 1 при положительном аргументе и 0 при отрицательном; мы воспользовались тем, что ее производная есть дельта-функция. Затем мы воспользовались тем, что $\int \delta(f(x))dx = \sum_{x: f(x)=0} \frac{1}{f'(x)}$, что очевидно из разложения $f$ в окрестности точек, где $f=0$.


Последнюю формулу в цепочке вы уже вывели из определения плотности вероятности. Только там надо было не забыть поделить на $\Delta h$, в соответствии с определением плотности вероятности.

[AndreyL]

Решение, конечно, компактное и понятное, но, к сожалению, это интеграл не берется аналитически. Численное его взятие оказывается не сильно устойчивым - почему-то функция распределения получается разрывной, хотя она точно непрерывная (если непрерывна $H(x)$). Решение с суммой (с разбиением на области монотонности) дает, хоть и громоздкое, но аналитическое решение, и функция распределения полностью аналогична моделированию Монте-Карлой.

[type2b]

(поправка: везде вместо $1/H'(x)$ должно быть $1/|H'(x)|$).

так плотность вероятности и есть разрывная. У нее особенности как минимум $~1/\sqrt{|h-h_0|}$ там, где $h_0$ совпадает с локальным экстремумом

[AndreyL]

В том то и дело, что не плотность, а функция распределения

$$ F(u) = \int_0^L \theta(u-H(x))\frac{dx}{L}$$

получается разрывной, что мне не совсем понятно. Скорее, это связано с вычислительными проблемами - я считал в лоб на Математике.

fu = Sin[x] + 1;
L = 10 Pi;
Plot[fu, {x, 0, L}]
{minh, maxh} = {MinValue[{fu, 0 <= x <= L}, x], MaxValue[{fu, 0 <= x <= L}, x]}

F[u_] := NIntegrate[HeavisideTheta[u - fu]/L, {x, 0, L}];
myTable = If[$KernelCount > 0, ParallelTable, Table];
dat = myTable[{u, F[u]}, {u, minh, maxh, L/10000}];
ListPlot[dat]

[type2b]

ага, это баг интегрирования. Видимо, из-за того, что $\theta$ почти везде постоянна, и стандартный алгоритм математики не может правильно разделить интервал интегрирования. Если вот так:

F[u_] := NIntegrate[HeavisideTheta[u - fu]/L, {x, 0, L}, 
   Method -> {"RiemannRule", "Type" -> "Left", "Points" -> 1000}, 
   MaxRecursion -> 0];

то выглядит получше.

[AndreyL]

Да, Вы правы

F[u_] := NIntegrate[HeavisideTheta[u - fu]/L, {x, 0, L}, 
   Method -> {"RiemannRule", "Type" -> "Left", "Points" -> 1000}, MaxRecursion -> 0];
Off[NIntegrate::ncvb]
gr1 = Plot[F[u], {u, minh, maxh}, PlotPoints -> 100, PlotStyle -> Red];
On[NIntegrate::ncvb]

myTable = If[$KernelCount > 0, ParallelTable, Table];
h = Sort[myTable[fu // N, {x, 0, L, L/10000}]];
n = Length[h]
Y = Table[i/n - 1/(2 n), {i, n}];
gr2 = ListPlot[Transpose[{h, Y}]];
Show[gr2, gr1]

Правда считает не быстро, но если это конечная задача (если не придется считать такие интегралы в цикле), то приемлемо. Спасибо!

[AndreyL]

Из спортивного интереса написал интегратор

MyInregral[func_, limits_] := 
 Block[{bonds =Sort[Union[Join[{limits[[2]], limits[[3]]}, 
       limits[[1]] /. Solve[{func == 0, limits[[2]] <= limits[[1]] <= limits[[3]]}, 
         limits[[1]]]]]], co}, co = Length[bonds];
  Total[Table[(bonds[[i + 1]] - bonds[[i]]) HeavisideTheta[ func /. {limits[[1]] -> (bonds[[i]] + bonds[[i + 1]])/2}], {i, co - 1}]]]
gr3 = Plot[MyInregral[u - fu, {x, 0, L}]/L, {u, minh, maxh}, PlotStyle -> Green];

Show[gr2, gr3]