Подобрать плотность вероятности

[Genie]

Добрый день!
Пожалуйста, помогите разобраться с такой задачей:
Имеются некоторые экспериментальные данные, и общий вид функции плотности вероятности. Нужно подобрать параметры этой функции.
Функция плотности: $p(x)=ce^{-ax}x^b$, где $c$ вычисляется из условия $\int_0^\infty p(x)\,dx=1$
Используя Matlab, я расчитала, что $c=\frac{a^{1+b}}{\Gamma(1+b)}$

Далее, я хочу решить систему методом наименьших квадратов. Система выглядит так:

$$\int_{i_1}^{i_2} p(x)\,dx=\frac{N_1}{N_a}$$

$$\int_{i_2}^{i_3} p(x)\,dx=\frac{N_5}{N_a}$$

...

$$\int_{i_5}^{i_6} p(x)\,dx=\frac{N_5}{N_a}$$

Где $N_k$ - количество попаданий в интервал от $i_k$ до $i_k+1$, $N_a$ - общее количество наблюдений.

Поскольку неизвестных две ($a$,$b$), а уравнений 5, я хочу использовать метод наименьших квадратов. Но, когда я вычисляю $\int p(x)\,dx$, получаю выражение, содержащее функцию WhittakerM(1/2 b, 1/2 b + 1/2, a x). Я гуглила, но ничего толкового про эту функцию не нашла. Где можно про нее почитать? Можно ли взять производные от нее по $b$ и $a$? Из того, что я наша, я поняла, что это табличная функция – получается, производные взять нельзя?

Еще я пробовала решать другим способом – просто ввести имеющиеся данные в программу statistica и попробовать подобрать параметры с помощью нелинейной регрессии, с помощью функции, определяемой пользователем. Если функция - $y(x)= e^{-ax}x^b$, то параметры подбираются, а если написать нужную мне функцию

$$y(x)= \frac{e^{-ax}x^b a^{1+b}}{\Gamma(1+b)}$$

, то выдается ошибка «predictions are probably very redundant; estimates suspect”. Как избавиться от этой ошибки?

Буду благодарна за любую помощь.

[ASA]

Genie писал(а):

Подобрать плотность вероятности
Функция плотности: $p(x)=ce^{-ax}x^b$, где $c$ вычисляется из условия $\int_0^\infty p(x)\,dx=1$

$p(x)=ce^{-ax}x^b$ при $x\geqslant 0$ и $p(x)=0$ при $x<0$. Может в этом дело?

Добавлено спустя 8 минут 12 секунд:

Genie писал(а):

Но, когда я вычисляю $\int p(x)\,dx$, получаю выражение, содержащее функцию WhittakerM(1/2 b, 1/2 b + 1/2, a x)

В Матлабе? Численно или используя символьные вычисления?

[Genie]

ASA писал(а):

Genie писал(а):

Подобрать плотность вероятности
Функция плотности: $p(x)=ce^{-ax}x^b$, где $c$ вычисляется из условия $\int_0^\infty p(x)\,dx=1$

$p(x)=ce^{-ax}x^b$ при $x\geqslant 0$ и $p(x)=0$ при $x<0$. Может в этом дело?

Да, функция плотности такая, как Вы написали. Но я не очень понимаю, в чем тут дело? Ведь интеграл плотности я вычисляю от нуля до бесконечности.

ASA писал(а):

Genie писал(а):

Но, когда я вычисляю $\int p(x)\,dx$, получаю выражение, содержащее функцию WhittakerM(1/2 b, 1/2 b + 1/2, a x)

В Матлабе? Численно или используя символьные вычисления?

В Матлабе, используя символьные вычисления. Определяю $a,b,x$ как символы и считаю интеграл

$$\int\frac{e^{-ax}x^b a^{1+b}}{\Gamma(1+b)}$$

, результат:

$$\frac{a^b\,x^b\,(ax)^{-\frac{1}{2}b}\,e^{-\frac{1}{2}ax}WhittakerM(\frac{1}{2}b,\frac{1}{2}b+\frac{1}{2},ax)}{(1+b)\Gamma(1+b)}$$

Как в Матлабе взять производную отдельно по a, и отдельно по b?

Пока писала, сообразила что может быть я не правильно беру интеграл. Я вызываю команду int(y), где

$$y=\frac{e^{-ax}x^b a^{1+b}}{\Gamma(1+b)}$$

. По какой переменной "берется" этот интеграл? Пожалуйста, посоветуйте хорошую книгу по matlab - я в нем совсем чайник.

[ASA]

Genie в сообщении #191674 писал(а):

По какой переменной "берется" этот интеграл?

По умолчанию, вроде бы $x$, но лучше написать явно.

Genie в сообщении #191674 писал(а):

выражение, содержащее функцию WhittakerM

Это, видимо, функция Уиттекера или какая-то ее разновидность. Но это - такие дебри.

А если ограничиться целым неотрицательным $b$? Тогда $\Gamma(1+b)=b!$

[Genie]

ASA писал(а):

Это, видимо, функция Уиттекера или какая-то ее разновидность. Но это - такие дебри.

А если ограничиться целым неотрицательным $b$? Тогда $\Gamma(1+b)=b!$

Да, насколько я понимаю, это она. Но я не смогла найти толковой информации про нее - не подскажете, где можно про нее почитать?

Да, при целых неотрицательных $b$ $\Gamma(1+b)=b!$, но, к сожалению, от функции Уиттекера это не избавляет.

[GAA]

Genie писал(а):

Но я не смогла найти толковой информации про нее - не подскажете, где можно про нее почитать?

См. ссылки внизу страницы Whittaker Function.

Добавлено спустя 16 минут 37 секунд:

Книга Уиттекер Э.Т., Ватсон Дж.Н. Курс современного анализа. Часть 2. Трансцендентные функции (2-е изд.) [перевод издания 1927г.] свободно доступна в электронном виде, например, на EqWorld. [В частности, см. гл 16 Вырожденная гипергеометрическая функция.]

Добавлено спустя 56 минут 54 секунды:

Чтение постановки, изложенной в первом сообщение Genie темы «Подобрать плотность вероятности», не приводит к мысли о разумности применения м.н.к. Вполне возможно, написано не все. Но я бы задумался о применении метода минимума $\chi^2$ или об оценке методом максимального правдоподобия параметров полиномиального (по другому — мультиномиального) распределения. Это особенно осмысленно, если результат подгонки будет сопровождаться проверкой гипотезы о принадлежности выборки гипотетическому семейству при помощи критерия Пирсона.

[--mS--]

GAA писал(а):

Чтение постановки, изложенной в первом сообщение Genie темы «Подобрать плотность вероятности», не приводит к мысли о разумности применения м.н.к. Вполне возможно, написано не все. Но я бы задумался о применении метода минимума $\chi^2$ или об оценке методом максимального правдоподобия параметров полиномиального (по другому — мультиномиального) распределения.

Это почти то же самое, что делает автор. В любом случае нужно минимизировать по параметрам функционал вида $\sum (p_i(a,b)-N_i/N)^2$.

[Genie]

Спасибо за советы! Буду дальше разбираться с задачей.

[GAA]

Спасибо, --mS--. Очень я бестолково написал. Хотел сказать следующее.
1. Желательно уточнить, действительно ли надо искать оценку асимптотически эквивалентную оценке метода минимума $\chi^2$, поскольку её придется искать численно. Может можно найти оценки по исходной выборке (до группировки)? Например, методом моментов. Оценку метода минимума $\chi^2$ придется искать, если будет использоваться критерий Пирсона.
2. Если ничего по теме не преподавалось, то может быть можно поступить просто так: взять прямоугольную параметрическую область и вычислить значения $x^2(a, b)$ на «достаточно» плотной сетке. [Обратите внимание: интегралы Вам надо брать определенные, а не неопределенные, но сильно это дело не упрощает].
3. Если значение оценки следует получить с приличной точностью, то уточнить значение можно при помощи какого-нибудь метода основанного на производных.

Если мне не изменяет память, в пакете Statistica 5.5 выполнить это не удастся. Матлабом я не владею. Так как ничего больше не знаю, то попробовал бы вычисления выполнить в пакете Maple. Maple 12, например, умеет вычислять не только необходимые в этой задаче интегралы, но и производные от них.