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.
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]
FindFit
gives several different error messages for Mathematica 10.4.1 Windows 10 and does not give a result.) $\endgroup$data
in your question doesn't match the figure. There are 11 peaks rather than 7 peaks in the supplied data. $\endgroup$