Раз уж разбираем задачу, доводим её до логического абсурда. Предположим, что у нас распределение описывается некой вероятностной функцией, являющейся суммой двух вероятностных функций с весами, сумма которых равна единице (условие нормировки). Если вспомнить про интегралы, что интеграл суммы равен сумме интегралов, а постоянный член выносится за интеграл, то функции можно описать, например, так.
Код:
> pMyDist<-function(q,lambda1,lambda2,w=0.5) w*ppois(q,lambda1)+(1-w)*ppois(q,lambda2)
> dMyDist<-function(x,lambda1,lambda2,w=0.5) w*dpois(x,lambda1)+(1-w)*dpois(x,lambda2)
> qMyDist<-function(p,lambda1,lambda2,w=0.5) w*qpois(p,lambda1)+(1-w)*qpois(p,lambda2)
Где w принимает значение [0,1] и является весовым коэффициентом.
Если имеем p(k)=f(k,lambda), где f распределение Пуассона, то результирующая функция описывается как p(k)=w*f(k,lambda1)+(1-w)*f(k,lambda2)
Графически результат выглядит вот так (в сравнении с обычным распределением Пуассона) для функций
p
1(k)=f(k,7,2)
p
2(k)=0.8*f(k,7,2)+0.2*f(k,0,8)
Код:
> plot(1:25,dpois(1:25,7.2),xlab="x",ylab="f(x)")
> lines(dMyDist(1:25,7.2,4,0.8),type="p",col="red",pch=2)
> legend(x=15,y=0.1,c("Распределение Пуассона","Сложное распределение"),pch=c(1,2),col=c("black","red"))
Счастье заключается в том, что теперь эту функцию можно попробовать подогнать под наше распределение.
Код:
> XX<-fitdist(LT,"MyDist",start=list(lambda1=9 ,lambda2=5 ), method="mle")
> summary(XX)
Fitting of the distribution ' MyDist ' by maximum likelihood
Parameters :
estimate Std. Error
lambda1 7.197402 0.05980069
lambda2 7.202613 0.05980405
Loglikelihood: -123149.4 AIC: 246302.8 BIC: 246320.5
Correlation matrix:
lambda1 lambda2
lambda1 1.0000000 -0.9227544
lambda2 -0.9227544 1.0000000
Как видно, принципиально лучше не стало. Отрицательный результат -- тоже результат. Но метод может пригодиться. Я его как-то удачно использовал.