How to use cdf function using beta function with more than the 2 standard parameter (shape parameters)?

My question is how to get cdf beta with more than the two standard shape parameter so I can define the probability of the observed data from the theoretical distribution of random permutated data?

I have 1,0000 random permutation of ~44,000 gene allele frequencies to generate Fst, which creates a 1000x44,000 matrix. Each column vector is a random set of Fst values across 44,000 genes from one permutation, Each row vector are Fst values for 1,000 different random permutations for the same gene. The districution of column is best described (bess fit) to a beta distribution with 4 paramters, not the two standard shape parater (a, b). The extra two parameter (lets call them p & q) define the minimum and maximum. Or location is p, and the scale is q-p. Thus, the best fit is beta distrigution with a, b, p, q.

These two extra (non-standard) parameter provide a more robust fit to the data.

HOWEVER, there seems to be no option in MatLab for a beta distribution with these extra scalar parameters. If I could use them, I could define the most likely beta distribution and thus the p-value of the observed data.

In summary can I fit a beta distribution that include the two extra scalar parameters to more precise determine the probility of the observed data?