不幸的是，我现在不能花时间来发布所有的细节，但我可以确认，解决3个参数（a，b和尺度参数）的3个方程的一般策略确实有效;我能够恢复OP显示的解决方案。
(1)三方程组：当x_1 = 100，x_2 = 250，x_3 = 500，分位数_1 = 1/10，分位数_2 = 1/2，分位数_3 = 9/10时，求出cdf（x_k）=分位数_k。
一般β分布CDF

beta_incomplete_regularized(a, b, (x - x0)/L)

其中x 0是偏移量（上述33.3），L是比例。（Scipy符号可能不同。）
(2)雅可比矩阵。对cdf（x_k）关于a、b和L求微分，给出3个梯度。构造一个矩阵，以这些梯度为行。最大值表示导数为

[beta_incomplete_regularized(a,b,(x-x0)/L)
  *(log((x-x0)/L)+psi[0](b+a)-psi[0](a))
  -(gamma(a)*hypergeometric_regularized([a,a,1-b],[a+1,a+1],(x-x0)/L)
            *gamma(b+a)*(x-x0)^a)
   /(L^a*gamma(b)),
 (gamma(b)*gamma(b+a)
          *hypergeometric_regularized([b,b,1-a],[b+1,b+1],1-(x-x0)/L)
          *(1-(x-x0)/L)^b)
  /gamma(a)
  +beta_incomplete_regularized(b,a,1-(x-x0)/L)
   *((-log(1-(x-x0)/L))-psi[0](b+a)+psi[0](b)),
 -(L^((-a)-1)*(1-(x-x0)/L)^(b-1)*(x-x0)^a)/beta(a,b)]

正则化超几何函数是3F 2函数除以gamma（n[k]）的乘积，其中n是第二参数列表。
在网上快速搜索一下，Scipy似乎没有实现3F 2，所以这可能是一个障碍。
(3)为多维牛顿算法提供方程和雅可比矩阵。
很抱歉，我不能给予一个完整的解决方案，或者用Python，现在。希望这对所有的帮助一样。
编辑：这里有一个完整的解决方案，为Maxima。

display2d: false;
linel: 65;

/* construct cdf for stretched beta */
assume (x > 0, L > 0, x0 > 0);
assume (x > x0, x < x0 + L);
assume (a > 0, b > 0);
load (distrib);
cdf: cdf_beta ((x - x0)/L, a, b);

/* construct gradient wrt free parameters a, b, L */
grad_e: [diff (cdf, a), diff (cdf, b), diff (cdf, L)];
expand_hypergeometric_regularized: lambda ([aa, bb, xx], hypergeometric (aa, bb, xx) / product (gamma (bb[k]), k, 1, length (bb)));
grad_e1: subst (hypergeometric_regularized = expand_hypergeometric_regularized, grad_e);

/* construct equations given user-specified data */
x0_assumed: 333/10;
F90: 100;
Median: 250;
F10: 500;
eqs_with_x0: [subst (x = F90, cdf) = 1/10, subst (x = Median, cdf) = 1/2, subst (x = F10, cdf) = 9/10];
eqs: subst (x0 = x0_assumed, eqs_with_x0);

/* construct Jacobian matrix by evaluating gradient at data values */
DF_with_x0: apply ('matrix, makelist (subst (x = xx, grad_e1), xx, [F90, Median, F10]));
DF: subst (x0 = x0_assumed, DF_with_x0);

/* solve equations via multidimensional Newton algorithm */
load (mnewton);
newtondebug: true;
/* need initial guess; assume cdf is skewed to the right, so a < b
 * and guess L is somewhat bigger than interquantile range
 */
solution: mnewton (eqs, [a, b, L], [2, 3, 600], DF);

/* test solution in equations */
float (subst (solution[1], eqs));

/* display pdf of beta distribution */
pdf: subst (x0 = x0_assumed, subst (solution[1], pdf_beta ((x - x0)/L, a, b)));
plot2d (pdf, [x, x0_assumed, x0_assumed + assoc (L, solution[1])]);

/* mean and sd */
printf (true, "Mean ~,1f", subst (solution[1], mean_beta (a, b) * L + x0_assumed));
printf (true, "Std ~,1f", subst (solution[1], std_beta (a, b) * L));

/* table of quantiles */
quantile: subst (solution[1], lambda ([q], quantile_beta (q, a, b) * L + x0_assumed));
for p in [100, 95, 90, 75, 50, 25, 10, 5, 0]
    do printf (true, "F~d ~,1f~%", p, quantile(1 - p/100));

下面是我得到的解决方案：

[a = 1.782967949218066,b = 7.660029561011594,
         L = 1301.704394701864]

和输出：

Mean 279.1
Std 157.6
F100 33.3
F95 76.6
F90 100.0
F75 157.2
F50 250.0
F25 371.8
F10 500.0
F5 581.7
F0 1335.0