Ошибка при использовании интеграла (строка 85): A и B должны быть скалярами с плавающей запятой.

У меня есть следующий код, в котором я хочу оценить параметры настраиваемого распределения, используя функцию MATLAB mle(). Подробнее о распространении.

Основной код:

x           =  [0   0   0   0   0.000649967501624919    0.00569971501424929 0.0251487425628719  0.0693465326733663  0.155342232888356   0.284835758212089   0.458277086145693   0.658567071646418   0.908404579771011   1.17284135793210    1.43977801109945    1.71951402429879    1.98925053747313    2.27553622318884    2.57147142642868    2.80390980450977    3.03829808509575    3.26583670816459    3.45642717864107    3.65106744662767    3.81950902454877    3.98275086245688    4.11259437028149    4.24683765811709    4.35043247837608    4.43832808359582    4.58427078646068    4.62286885655717    4.68361581920904    4.75686215689216    4.80245987700615    4.84005799710015    4.86280685965702    4.91675416229189    4.92725363731813    4.90890455477226    4.96570171491425    4.92315384230789    4.95355232238388    4.92790360481976    4.93135343232838    4.90310484475776    4.90885455727214    4.86765661716914    4.85490725463727    4.81940902954852    4.81450927453627    4.78621068946553    4.74206289685516    4.71791410429479    4.69961501924904    4.65706714664267    4.63611819409030    4.60176991150443    4.57512124393780    4.53507324633768    4.48252587370631    4.47062646867657    4.43127843607820    4.39963001849908    4.37598120093995    4.29548522573871    4.31033448327584    4.21708914554272    4.21913904304785    4.18669066546673    4.16719164041798    4.09774511274436    4.07989600519974    4.02869856507175    3.98485075746213    3.95785210739463    3.93945302734863    3.90240487975601    3.87025648717564    3.81185940702965    3.78461076946153    3.74091295435228    3.71666416679166    3.67276636168192    3.65846707664617    3.61361931903405    3.58712064396780    3.55452227388631    3.53082345882706    3.49197540122994    3.48582570871456    3.46512674366282    3.41227938603070    3.36278186090695    3.35528223588821    3.31238438078096    3.27213639318034    3.23863806809660    3.24173791310434    3.19339033048348    3.20118994050298    3.16489175541223    3.10739463026849    3.09484525773711    3.08094595270237    3.02129893505325    3.02309884505775    2.99375031248438    2.95765211739413    2.93230338483076    2.89560521973901    2.87805609719514    2.85440727963602    2.82285885705715    2.80175991200440    2.79091045447728    2.73901304934753    2.72701364931753    2.73441327933603    2.71646417679116    2.68236588170592    2.65551722413879    2.63356832158392    2.60361981900905    2.58147092645368    2.57697115144243    2.54287285635718    2.53502324883756    2.47702614869257    2.50387480625969    2.46487675616219    2.45722713864307    2.42707864606770    2.41762911854407    2.39823008849558    2.38708064596770    2.34058297085146    2.35613219339033    2.32123393830309    2.30503474826259    2.27613619319034    2.27248637568122    2.25113744312784    2.24908754562272    2.22703864806760    2.20583970801460    2.17244137793110    2.15709214539273    2.16469176541173    2.12139393030348    2.12809359532023    2.11389430528474    2.09774511274436    2.07629618519074    2.07459627018649    2.05394730263487    2.04724763761812    2.01684915754212    2.01684915754212    2.00409979501025    1.98955052247388    1.96540172991350    1.95890205489726    1.93035348232588    1.92295385230738    1.90605469726514    1.89785510724464    1.87070646467677    1.88000599970002    1.86295685215739    1.84420778961052    1.82510874456277    1.80480975951202    1.80785960701965    1.80870956452177    1.77581120943953    1.76771161441928    1.77131143442828    1.76636168191590    1.75081245937703    1.73156342182891    1.69876506174691    1.70836458177091    1.70376481175941    1.67196640167992    1.68101594920254    1.66586670666467    1.66061696915154    1.64296785160742    1.63291835408230    1.62506874656267    1.62516874156292    1.60556972151392    1.59007049647518    1.59187040647968    1.57947102644868    1.57577121143943    1.54527273636318    1.57237138143093    1.54637268136593    1.54802259887006    1.50492475376231    1.52077396130193    1.50417479126044    1.50162491875406    1.50062496875156    1.48957552122394    1.47997600119994    1.47027648617569    1.44452777361132    1.45407729613519    1.44272786360682    1.43247837608120    1.41657917104145    1.40787960601970    1.39323033848308    1.40282985850707    1.39403029848508    1.38233088345583    1.37888105594720    1.37943102844858    1.36183190840458    1.34808259587021    1.34503274836258    1.33703314834258    1.33308334583271    1.32253387330633    1.32698365081746    1.29963501824909    1.30758462076896    1.29103544822759    1.29473526323684    1.27413629318534    1.26858657067147    1.27888605569722    1.26063696815159    1.27863606819659    1.25168741562922    1.23913804309785    1.24788760561972    1.22308884555772    1.24198790060497    1.22133893305335    1.20678966051697    1.20098995050247    1.20343982800860    1.18779061046948    1.19024048797560    1.17194140292985    1.17369131543423    1.16869156542173    1.15814209289536    1.15429228538573    1.15904204789761    1.12774361281936    1.15344232788361    1.13744312784361    1.12909354532273    1.12479376031198    1.11099445027749    1.11469426528674    1.11064446777661    1.10464476776161    1.10309484525774    1.10689465526724    1.07654617269137    1.07884605769712    1.07359632018399    1.06864656767162    1.07544622768862    1.06689665516724    1.04884755762212    1.06164691765412    1.04979751012449    1.04529773511324    1.02839858007100    1.03634818259087    1.01709914504275    1.02089895505225    1.01024948752562    1.01549922503875    1.01319934003300    1.01404929753512    1.00839958002100    0.995400229988501   0.989850507474626   0.978801059947003   0.977551122443878   0.980450977451127   0.975451227438628   0.969201539923004   0.964151792410380   0.964601769911504   0.958802059897005   0.955702214889256   0.948602569871506   0.960751962401880   0.941352932353382   0.928653567321634   0.949002549872506   0.937053147342633   0.913854307284636   0.916204189790510   0.915454227288636   0.902604869756512   0.909454527273636   0.895505224738763   0.898355082245888   0.894455277236138   0.902454877256137   0.883705814709265   0.888405579721014   0.876356182190891   0.881555922203890   0.878156092195390   0.868456577171141   0.870406479676016   0.863906804659767   0.862456877156142   0.858757062146893   0.851307434628269   0.851107444627769   0.833908304584771   0.843507824608770   0.831708414579271   0.836858157092145   0.829058547072646   0.828508574571272   0.822908854557272   0.820508974551273   0.815559222038898   0.819709014549273   0.809609519524024   0.813409329533523   0.800759962001900   0.806609669516524   0.806959652017399   0.792260386980651   0.787660616969152   0.783810809459527   0.794960251987401   0.771061446927654   0.788910554472276   0.789510524473776   0.763061846907655   0.776761161941903   0.767561621918904   0.773611319434028   0.750262486875656   0.765811709414529   0.765911704414779   0.748012599370032   0.741612919354032   0.757312134393280   0.752612369381531   0.741362931853407   0.742212889355532   0.741912904354782   0.743162841857907   0.732963351832408   0.732813359332033   0.733363331833408   0.721913904304785   0.716664166791661   0.726713664316784   0.709764511774411   0.700064996750163   0.710764461776911   0.717664116794160   0.707314634268287   0.707114644267787   0.705614719264037   0.709164541772911   0.696665166741663   0.680765961701915   0.686715664216789   0.694465276736163   0.683015849207540   0.681715914204290   0.694465276736163   0.688615569221539   0.680665966701665   0.672316384180791   0.672866356682166   0.656517174141293   0.665316734163292   0.671566421678916   0.666266686665667   0.652917354132293   0.663366831658417   0.651917404129794   0.663816809159542   0.661366931653417   0.647017649117544   0.655167241637918   0.647867606619669   0.636918154092295   0.645467726613669   0.633118344082796   0.640217989100545   0.634668266586671   0.618669066546673   0.635068246587671   0.632568371581421   0.623118844057797   0.623868806559672   0.623718814059297   0.621368931553422   0.623768811559422   0.608419579021049   0.616019199040048   0.609869506524674   0.606569671516424   0.614019299035048   0.610269486525674   0.596520173991300   0.595570221488926   0.593270336483176   0.596670166491675   0.598470076496175   0.597770111494425   0.593720313984301   0.592770361481926   0.585420728963552   0.580870956452177   0.584120793960302   0.580270986450677   0.577971101444928   0.579021048947553   0.572821358932053   0.585970701464927   0.572921353932303   0.567071646417679   0.569971501424929   0.571271436428179   0.568421578921054   0.567421628918554   0.569521523923804   0.563721813909305   0.558772061396930   0.562171891405430   0.557872106394680   0.549072546372681   0.558722063896805   0.536973151342433   0.561021948902555   0.544172791360432   0.552122393880306   0.553072346382681   0.546222688865557   0.551472426378681   0.540772961351932   0.541122943852807   0.542772861356932   0.530323483825809   0.526023698815059   0.529273536323184   0.524573771311435   0.525923703814809   0.524923753812309   0.516474176291185   0.527273636318184   0.527723613819309   0.518424078796060   0.517874106294685   0.516074196290186   0.517924103794810   0.523173841307935   0.514474276286186   0.513174341282936   0.498875056247188   0.518024098795060   0.507924603769812   0.505524723763812   0.507174641267937   0.502874856257187   0.502624868756562   0.500624968751562   0.510824458777061   0.490925453727314   0.492675366231688   0.489925503724814   0.478126093695315   0.485775711214439   0.491775411229439   0.489925503724814   0.491325433728314   0.487225638718064   0.485725713714314   0.485675716214189   0.477676116194190   0.483875806209690   0.478026098695065   0.470176491175441   0.471926403679816   0.483625818709065   0.469376531173441   0.474026298685066   0.467826608669567   0.462426878656067];

Censored    =  ones(1,size(x,2));% 
custpdf     =  @eval_custpdf;
custcdf     =  @eval_custcdf;
phat        =  mle(x,'pdf', custpdf,'cdf', custcdf,'start',[1 0.1 0.3 0.1 0.01 -0.3],...
               'lowerbound',[0 0 0 0 0 -inf],'upperbound',[inf inf inf inf inf inf],'Censoring',Censored);

% Cheking how close the estimated PDF and CDF match with those from the data x
t   =  0.001:0.001:0.5;
figure();
plot(t,x);hold on
plot(t,custpdf(t, phat(1), phat(2), phat(3), phat(4), phat(5), phat(6))) 

figure();
plot(t,cumsum(x)./sum(x));hold on
plot(t,custcdf(t, phat(1), phat(2), phat(3), phat(4), phat(5), phat(6)))

Функции:

function out = eval_custpdf(x,myalpha,mybeta,mytheta,a,b,c)
  first_integral      =  integral(@(x) eval_K(x,a,b,c),0,1).^-1;
  theta_t_ratio       =  (mytheta./x);
  incomplete_gamma    =  igamma(myalpha,theta_t_ratio.^mybeta);
  n_gamma             =  gamma(myalpha);
  exponent_term       =  exp(-theta_t_ratio.^mybeta-(c.*(incomplete_gamma./n_gamma)));


  numerator   =  first_integral.* mybeta.*incomplete_gamma.^(a-1).*...
               theta_t_ratio.^(myalpha*mybeta+1).*exponent_term;
  denominator =  mytheta.* n_gamma.^(a+b-1).* (n_gamma-incomplete_gamma.^mybeta).^(1-b);

  out         =  numerator./denominator;
end

function out = eval_custcdf(x,myalpha,mybeta,mytheta,a,b,c)
  first_integral      =  integral(@(x) eval_K(x,a,b,c),0,1).^-1;
  theta_t_ratio       =  (mytheta./x);
  incomplete_gamma    =  igamma(myalpha,theta_t_ratio.^mybeta);
  n_gamma             =  gamma(myalpha);
  second_integral     =  integral(@(x) eval_K(x,a,b,c),0, incomplete_gamma.^mybeta./n_gamma);
                                                         % |<----- PROBLEMATIC LINE ----->|

  out     =  first_integral*second_integral;

end    

function out = eval_K(x,a,b,c)

  out =  x.^(a-1).*(1-x).^(b-1).*exp(-c.*x);

end

Интеграл, вызывающий проблему, является вторым интегралом в функции eval_custcdf(), поскольку его верхний предел представляет собой массив (обозначенный ПРОБЛЕМАТИЧЕСКАЯ СТРОКА).

Есть ли способ взять одно значение из массива x, чтобы верхний предел оставался скалярным? А затем вычислить cdf таким образом, чтобы на выходе cdf был массив? Может быть, с forloop? Но я не могу понять, как это реализовать?

Как я могу обойти эту проблему?

Любая помощь будет оценена.

Заранее спасибо.


debugging matlab integral mle
person nashynash    schedule 11.06.2019    source источник
comment
Ничего общего с ответом. Но почему вы создаете анонимную функцию, которая вызывает ту же самую функцию? Также вы должны разделить вашу функцию на несколько строк, ваш код довольно сложно читать.   -  person obchardon    schedule 11.06.2019
comment
@obchardon Я попытался разбить функции на несколько строк, чтобы улучшить читаемость. Однако я не уверен, как не «создать анонимную функцию, вызывающую точно такую ​​же функцию». Не могли бы вы рассказать подробнее?   -  person nashynash    schedule 11.06.2019
comment
@obchardon Кажется, я исправил. Я избавился от лишнего звонка сейчас?   -  person nashynash    schedule 11.06.2019
comment
Если вы можете попробовать свой PDF-файл, тогда cdf = cumsum(pdf).   -  person Cris Luengo    schedule 11.06.2019

Ответы (1)


arrow_upward
2
arrow_downward

  • eval_custcdf - функция, которая должна возвращать 1D array of length n для заданного n ввода данных.
  • Я использую for loop для вычисления вывода для данного ввода, а затем возвращаю весь массив как вывод eval_custcdf

Я передал элементы входного массива по одному

Вот так eval_custcdf может выглядеть

function out = eval_custcdf(x,myalpha,mybeta,mytheta,a,b,c)
    out = zeros(size(x));
    for i = 1: length(x)
          first_integral      =  integral(@(w) eval_K(w,a,b,c),0,1).^-1;
          theta_t_ratio       =  (mytheta./x(i));
          incomplete_gamma    =  igamma(myalpha,theta_t_ratio.^mybeta);
          n_gamma             =  gamma(myalpha);
          second_integral     =  integral(@(w) eval_K(w,a,b,c),0, incomplete_gamma.^mybeta./n_gamma);

          out(i) = first_integral*second_integral;
    end
end
person Adam    schedule 11.06.2019
comment
Еще раз большое спасибо, Адам. Это хорошо работает. Теперь мне нужно найти способ найти хорошее предположение для начальных точек начала параметров. - person nashynash; 12.06.2019