0
$\begingroup$

I wrote a code in Mathematica, which tries to match the data by using a weighted sum of 4-parameter beta distribution. I have real-world frequency data, that look like the below image. In the latter, I manually sketched some “sub-distributions”, total 7 in general. I thought using a highly flexible 4-parameter beta distribution is a way to go. However, my FindFit command gives an error message, and you can see from the plot the fit is far from perfect. I suppose the problem arises due to the poor guess values, but I have no clue how to rectify the issue.

Can somebody please help? Please note I must have analytic expression in the end, that is why I used a weighted sum of known densities. The mathematica code is given below.enter image description here Thank you. 

data = {{0.`, 0.14943079650402`}, {0.1`, 
0.00021827389719`}, {0.1111111`, 0.00032284128247`}, {0.125`, 
0.00105345470365`}, {0.1428571`, 0.00189604749903`}, {0.1666667`, 
0.00496507203206`}, {0.2`, 0.0128972902894`}, {0.2222222`, 
0.00094202003675`}, {0.25`, 0.03224607929587`}, {0.2857143`, 
0.00478106550872`}, {0.3`, 0.00093960686354`}, {0.3333333`, 
0.08392683416605`}, {0.375`, 0.00326516595669`}, {0.4`, 
0.02893583476543`}, {0.4285714`, 0.0072281844914`}, {0.4444444`, 
0.00200031138957`}, {0.5`, 0.22689934074879`}, {0.5555556`, 
0.00208450458013`}, {0.5714286`, 0.00720308022574`}, {0.6`, 
0.02667851746082`}, {0.625`, 0.00338809774257`}, {0.6666667`, 
0.10289531946182`}, {0.7`, 0.00141532090493`}, {0.7142857`, 
0.00444651674479`}, {0.75`, 0.03364257141948`}, {0.7777778`, 
0.00085208012024`}, {0.8`, 0.01232228334993`}, {0.8333333`, 
0.00435024732724`}, {0.8571429`, 0.00152650044765`}, {0.875`, 
0.00095871940721`}, {0.8888889`, 0.00024728401331`}, {0.9`, 
0.00045259558829`}, {1.`, 0.23558811843395`}};

   f[x_, a_, b_, p_, q_] := 
  Piecewise[{{(Gamma[
         p + q]*((x - a)^(p - 1))*((b - x)^(q - 1)))/(Gamma[p]*
        Gamma[q]*((b - a)^(p + q - 1))), a <= x <= b}, {0, 
     x > b && x < a}}];

model = w1*f[x, a1, b1, p1, q1] + w2*f[x, a2, b2, p2, q2] + 
   w3*f[x, a3, b3, p3, q3] + w4*f[x, a4, b4, p4, q4] + 
   w5*f[x, a5, b5, p5, q5] + w6*f[x, a6, b6, p6, q6] + 
   w7*f[x, a7, b7, p7, q7];

result = FindFit[
  data, {model, 
   w1 > 0 && w2 > 0 && w3 > 0 && w4 > 0 && w5 > 0 && w6 > 0 && 
    w7 > 0 && w1 + w2 + w3 + w4 + w5 + w6 + w7 == 1 && p1 > 0 && 
    p2 > 0 && p3 > 0 && p4 > 0 && p5 > 0 && p6 > 0 && p7 > 0 && 
    q1 > 0 && q2 > 0 && q3 > 0 && q4 > 0 && q5 > 0 && q6 > 0 && 
    q7 > 0 && b1 > a1 && b2 > a2 && b3 > a3 && b4 > a4 && b5 > a5 && 
    b6 > a6 && b7 > a7 && a1 >= 0 && a2 >= 0 && a3 >= 0 && a4 >= 0 && 
    a5 >= 0 && a6 >= 0 && a7 >= 0}, {{w1, 0.10}, {a1, 0.01}, {b1, 
    0.22}, {p1, 3}, {q1, 2}, {w2, 0.20}, {a2, 0.22}, {b2, 0.30}, {p2, 
    2.9}, {q2, 3.12}, {w3, 0.21}, {a3, 0.30}, {b3, 0.38}, {p3, 
    4}, {q3, 3}, {w4, 0.33}, {a4, 0.38}, {b4, 0.45}, {p4, 4.5}, {q4, 
    5.3}, {w5, 0.08}, {a5, 0.45}, {b5, 0.65}, {p5, 40}, {q5, 11}, {w6,
     0.07}, {a6, 0.65}, {b6, 0.78}, {p6, 51}, {q6, 2.1}, {w7, 
    0.01}, {a7, 0.78}, {b7, 0.99}, {p7, 14}, {q7, 10}}, x]

 w1est = w1 /. result[[1]]; a1est = a1 /. result[[2]]; b1est = 
 b1 /. result[[3]]; p1est = p1 /. result[[4]]; q1est = 
 q1 /. result[[5]]; w2est = w2 /. result[[6]]; a2est = 
 a2 /. result[[7]]; b2est = b2 /. result[[8]]; p2est = 
 p2 /. result[[9]]; q2est = q2 /. result[[10]]; w3est = 
 w3 /. result[[11]]; a3est = a3 /. result[[12]]; b3est = 
 b3 /. result[[13]]; p3est = p3 /. result[[14]]; q3est = 
 q3 /. result[[15]]; w4est = w4 /. result[[16]]; a4est = 
 a4 /. result[[17]]; b4est = b4 /. result[[18]]; p4est = 
 p4 /. result[[19]]; q4est = q4 /. result[[20]]; w5est = 
 w5 /. result[[21]]; a5est = a5 /. result[[22]]; b5est = 
 b5 /. result[[23]]; p5est = p5 /. result[[24]]; q5est = 
 q5 /. result[[25]]; w6est = w6 /. result[[26]]; a6est = 
 a6 /. result[[27]]; b6est = b6 /. result[[28]]; p6est = 
 p6 /. result[[29]]; q6est = q6 /. result[[30]]; w7est = 
 w7 /. result[[31]]; a7est = a7 /. result[[32]]; b7est = 
 b7 /. result[[33]]; p7est = p7 /. result[[34]]; q7est = 
 q7 /. result[[35]];

modelest[x_] = 
  w1est*f[x, a1est, b1est, p1est, q1est] + 
   w2est*f[x, a2est, b2est, p2est, q2est] + 
   w3est*f[x, a3est, b3est, p3est, q3est] + 
   w4est*f[x, a4est, b4est, p4est, q4est] + 
   w5est*f[x, a5est, b5est, p5est, q5est] + 
   w6est*f[x, a6est, b6est, p6est, q6est] + 
   w7est*f[x, a7est, b7est, p7est, q7est];

Plot[modelest[x], {x, 0, 1}]~Show~ListPlot[data]

ListPlot[data, Filling -> Axis]
Improve this question
$\endgroup$
6
  • $\begingroup$ I think you're getting a pretty good fit given that you have 33 data points and 35 parameters to be estimated. You might want to ask the more general question as "What should I do with this data? Fit a probability distribution or perform a regression? Or something else?" at CrossValidated and then come back to this site as to how to implement the solution. (And the line with FindFit gives several different error messages for Mathematica 10.4.1 Windows 10 and does not give a result.) $\endgroup$
    – JimB
    Jun 17, 2017 at 17:45
  • $\begingroup$ Dear Jim, thank you for your reply. My message might have been confusing, sorry about that. I am getting a pretty bad fit actually. That image above (the red lines), are hand-drawn by me. If you run my Mathematica code, it will show it. $\endgroup$
    – Alex
    Jun 17, 2017 at 21:54
  • $\begingroup$ Your question is not confusing. However it is impossible to do what you think you want to do. You have 33 data points and 35 parameters. Doesn't that statement make alarm bells sound? Also the beta distributions provide you with probability "densities" and not relative frequencies for values with unequal distances between apparently binned values (or the values are simply a finite set of values and a continuous mixture distribution will likely not make sense). And as I said before, you're mixing up regression and fitting a probability distribution. Get thee to a statistician! $\endgroup$
    – JimB
    Jun 17, 2017 at 22:41
  • $\begingroup$ Dear Jim, thank you for your message. Now I see your point. I definitely need to fit a probability distribution. Do I need to use MixtureDistribution? I need a closed form expression for a fitted density, and obviously the integral under the curve must add up to 1. $\endgroup$
    – Alex
    Jun 18, 2017 at 2:46
  • $\begingroup$ The data in your question doesn't match the figure. There are 11 peaks rather than 7 peaks in the supplied data. $\endgroup$
    – Jack LaVigne
    Jun 18, 2017 at 3:03

1 Answer 1

Reset to default
1
$\begingroup$

This is an extended comment rather than an answer.

I'm convinced that the values you show must be ratios of a small number of relatively small integers meaning that attempting to fit a continuous mixture distribution is not warranted. Here's why I think that. Suppose we have the following combinations of numerators and denominators:

Table of ratios

The resulting ratios match your data perfectly and explains the varying distances between values: these are the only values possible which means there was no grouping of the data into unequally-sized bins. You have a discrete distribution with 33 unique values.

I understand that likely someone of authority over you has probably insisted that you construct an analytic expression. But that makes no sense in this case. The meaning of the numerators and denominators needs to be described to make headway in summarizing what must be a complete census (as you've described in other posts at Wolfram Community). Short of that your best bet is to provide the complete table of 33 ratios along with their relative frequencies.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Dear Jim, Thank you for your message. Your are right, it is a discrete distribution. But the task still remains. I would appreciate any feedback, regards. $\endgroup$
    – Alex
    Jun 18, 2017 at 12:49
  • $\begingroup$ You're certainly not making things easy for us. You supply the wrong graph for the data presented (as @JackLaVigne has pointed out) and it appears that you have a mixture of binomial distributions all with a proportion around 0.5 (with the number of trials taking on values 6, 7, 8, 9, and 10) but with an overabundance of zeros and ones. The solution is simply to present the table of values (and describe what the values mean). Forget the fit. Tell your boss he's crazy to ask for an analytic description. (Well, that's easy for me to say: it's not my boss.) $\endgroup$
    – JimB
    Jun 18, 2017 at 17:40
  • $\begingroup$ Dear Jim, thank you for your message. Regarding the graph, I am not sure why it came up wrong as you say. I use a standard Mathematica command to generate that graph (ListPlot[data, Filling -> Axis]) and I am not sure what I should have to make it correctly looking. I also appreciate your feedback regarding forgetting the fit. If you take it as a general problem, with a bunch of values of of a random variables and associated frequencies, how would you best fit a density to it? $\endgroup$
    – Alex
    Jun 18, 2017 at 17:54
  • $\begingroup$ I know you keep insisting it is wrong, but technically, for ANY given values of X and associated probabilities, one can try to best fit even if one does not know what X represents and where the numbers came from. Maybe the fit would not be the best, but one can at least try. $\endgroup$
    – Alex
    Jun 18, 2017 at 17:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.