0

I have the following code that I wish to estimate the parameters of a custom distribution using MATLAB's function mle(). For more details on the distribution.

The main code is:

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)))

The functions are:

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

The integral that is causing the problem is the second intergral in the function eval_custcdf() as its upper limit is an array (denoted by PROBLEMATIC LINE).

Is there a way to take a single value from the array x such that the upper limit remains a scalar? And then calculate the cdf such that the output of the cdf is an array? Using a forloop, maybe? But I cannot seem to figure how to implement that?

How can I work around this problem?

Any help would be appreciated.

Thanks in advance.

nashynash
  • 375
  • 2
  • 19
  • 1
    Nothing to do with the answer. But why do you create an anonymous function that call call the exact same function ? Also you should split your function into several lines, it's pretty hard to read your code. – obchardon Jun 11 '19 at 16:54
  • @obchardon I tried to split the functions into several lines to improve readability. However, I am not sure how not to ` create an anonymous function that calls the exact same function`. Could you elaborate more? – nashynash Jun 11 '19 at 17:40
  • @obchardon I think I fixed it. Did I get rid of the redundant call now? – nashynash Jun 11 '19 at 17:54
  • If you can sample your PDF, then `cdf = cumsum(pdf)`. – Cris Luengo Jun 11 '19 at 18:20

1 Answers1

2
  • eval_custcdf is a function expected to return 1D array of length n for a given n data input.
  • I use a for loop to compute the output for a given input, then return the whole array as output of eval_custcdf

I passed the input array elements one at a time

This is how eval_custcdf may look like

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
Adam
  • 2,726
  • 1
  • 9
  • 22
  • 1
    Thank you very much once again, Adam. That works well. Now I need to figure out a way to find the good guess for the initial start points of the parameters. – nashynash Jun 12 '19 at 01:26