Since $P(within 1~4)$ is the ultimate parameter of concern, that is why I came up with the wish to test for this parameter. I did go for the $\chi^2$ test as well. Though, with either $\chi^2$-test or simulation method ($p=0.283$), so far I did not see any statistically significant difference. Meanwhile, as you mentioned, $HCp$ is systematically larger than $LCp$, this might be the only thing I can argue on then. All in all, that is a fair solution above. Thanks.

@Roland Thanks for your response. I agree that bootstrapping can be a compromise in my case. Should have thought of that. However, I also feel a bit "uncomfortable" about the gap between samples during bootstrapping. Even though the true distribution is unknown, the variables should be well approximated by the criteria above (#1-5). So, if possible, it seems preferable to simulate a distribution over doing bootstrapping. This is not a must here, though (just trying to learn more, like when I learnt to simulate correlated multivariates).

I also tried some exponential function but I have not been able to solve them...

Actually, I am finding an equation, not necessarily $f(x)=ax^3+bx^2+cx+d$, which I can map $(p,q,r)$ on and has to satisfy (1)-(4), $f''(x)$ positive (yes, the slopes should increase with x), AND something that I forgot to mention: $f'''(x)$ positive. I tried $f(x)=ax^2+bx+c$ before, but this does not generate a fitting that I want, because I want the increase in slope to be mild at first and more vigorous later (i.e. $f''(x)$ will increase with x). So I came to this cubic equation, which did not give me what I want because now $f''(x)$ is not always positive.

@JamesPhillips Thank you very much! Let me look into your solution first

@James sure no problems. I have made some edits at the end of the post.