diff --git a/CHANGELOG_UNRELEASED.md b/CHANGELOG_UNRELEASED.md
index eee3614a7..e48778c6e 100644
--- a/CHANGELOG_UNRELEASED.md
+++ b/CHANGELOG_UNRELEASED.md
@@ -58,6 +58,50 @@
   + `measurable_function.v`
   + `measure.v`
 
+- in `pseudometric_normed_Zmodule.v`:
+  + lemma `continuous_comp_cvg`
+
+- in `derive.v`:
+  + lemma `derive1_onem`
+
+- in `ftc.v`:
+  + lemmas `integration_by_substitution_onem`, `Rintegration_by_substitution_onem`
+
+- in `probability.v`:
+  + lemmas `continuous_onemXn`, `onemXn_derivable`,
+    `derivable_oo_continuous_bnd_onemXnMr`, `derive_onemXn`,
+    `Rintegral_onemXn`
+  + definition `XMonemX`
+  + lemmas `XMonemX_ge0`, `XMonemX_le1`, `XMonemX0n`, `XMonemXn0`,
+    `XMonemX00`, `XMonemXC`, `XMonemXM`
+  + lemmas `continuous_XMonemX`, `within_continuous_XMonemX`,
+    `measurable_XMonemX`, `bounded_XMonemX`, `integrable_XMonemX`,
+    `integrable_XMonemX_restrict`, `integral_XMonemX_restrict`
+  + definition `beta_fun`
+  + lemmas `EFin_beta_fun`, `beta_fun_sym`, `beta_fun0n`, `beta_fun00`,
+    `beta_fun1S`, `beta_fun11`, `beta_funSSS`, `beta_funSS`, `beta_fun_fact`
+  + lemmas `beta_funE`, `beta_fun_gt0`, `beta_fun_ge0`
+  + definition `beta_pdf`
+  + lemmas `measurable_beta_pdf`, `beta_pdf_ge0`, `beta_pdf_le_beta_funV`,
+    `integrable_beta_pdf`, `bounded_beta_pdf_01`
+  + lemma `invr_nonneg_proof`, definition `invr_nonneg`
+  + definition `beta_prob`
+  + lemmas `integral_beta_pdf`, `beta_prob01`, `beta_prob_fin_num`,
+    `beta_prob_dom`, `beta_prob_uniform`,
+    `integral_beta_prob_bernoulli_prob_lty`,
+    `integral_beta_prob_bernoulli_prob_onemX_lty`,
+    `integral_beta_prob_bernoulli_prob_onem_lty`, `beta_prob_integrable`,
+    `beta_prob_integrable_onem`, `beta_prob_integrable_dirac`,
+    `beta_prob_integrable_onem_dirac`, `integral_beta_prob`
+  + definition `div_beta_fun`
+  + lemmas `div_beta_fun_ge0`, `div_beta_fun_le1`
+  + definition `beta_prob_bernoulli_prob`
+  + lemma `beta_prob_bernoulli_probE`
+
+
+- in `unstable.v`:
+  + lemmas `leq_prod2`, `leq_fact2`, `normr_onem`
+
 ### Changed
 
 - in `lebesgue_stieltjes_measure.v` specialized from `numFieldType` to `realFieldType`:
@@ -305,6 +349,14 @@
 
 ### Renamed
 
+
+- in `probability.v`:
+  + `bernoulli` -> `bernoulli_prob`
+  + `bernoulli_probE` -> `bernoulliE`
+  + `integral_bernoulli_prob` -> `integral_bernoulli_prob`
+  + `measurable_bernoulli` -> `measurable_bernoulli_prob`
+  + `measurable_bernoulli2` -> `measurable_bernoulli_prob2`
+
 ### Generalized
 
 - in `pseudometric_normed_Zmodule.v`:
diff --git a/classical/unstable.v b/classical/unstable.v
index 9b57974e4..088c1f90c 100644
--- a/classical/unstable.v
+++ b/classical/unstable.v
@@ -76,6 +76,35 @@ elim: r => [|a l]; first by rewrite !big_nil exp1n.
 by rewrite !big_cons; case: ifPn => // Pa <-; rewrite expnMn.
 Qed.
 
+Lemma leq_prod2 (x y n m : nat) : (n <= x)%N -> (m <= y)%N ->
+  (\prod_(m <= i < y) i * \prod_(n <= i < x) i <= \prod_(n + m <= i < x + y) i)%N.
+Proof.
+move=> nx my; rewrite big_addn -addnBA//.
+rewrite [in leqRHS]/index_iota -addnBAC// iotaD big_cat/=.
+rewrite mulnC leq_mul//.
+  by apply: leq_prod; move=> i _; rewrite leq_addr.
+rewrite subnKC//.
+rewrite -[in leqLHS](add0n m) big_addn.
+rewrite [in leqRHS](_ : y - m = ((y - m + x) - x))%N; last first.
+  by rewrite -addnBA// subnn addn0.
+rewrite -[X in iota X _](add0n x) big_addn -addnBA// subnn addn0.
+by apply: leq_prod => i _; rewrite leq_add2r leq_addr.
+Qed.
+
+Lemma leq_fact2 (x y n m : nat) : (n <= x) %N -> (m <= y)%N ->
+  (x`! * y`! * ((n + m).+1)`! <= n`! * m`! * ((x + y).+1)`!)%N.
+Proof.
+move=> nx my.
+rewrite (fact_split nx) -!mulnA leq_mul2l; apply/orP; right.
+rewrite (fact_split my) mulnCA -!mulnA leq_mul2l; apply/orP; right.
+rewrite [leqRHS](_ : _ =
+    (n + m).+1`! * \prod_((n + m).+2 <= i < (x + y).+2) i)%N; last first.
+  by rewrite -fact_split// ltnS leq_add.
+rewrite mulnA mulnC leq_mul2l; apply/orP; right.
+do 2 rewrite -addSn -addnS.
+exact: leq_prod2.
+Qed.
+
 Section max_min.
 Variable R : realFieldType.
 
@@ -280,6 +309,14 @@ Qed.
 End onem.
 Notation "`1- r" := (onem r) : ring_scope.
 
+Lemma normr_onem {R : realDomainType} (x : R) :
+  (0 <= x <= 1 -> `| `1-x | <= 1)%R.
+Proof.
+move=> /andP[x0 x1]; rewrite ler_norml; apply/andP; split.
+  by rewrite lerBrDl lerBlDr (le_trans x1)// lerDl.
+by rewrite lerBlDr lerDl.
+Qed.
+
 Lemma onemV (F : numFieldType) (x : F) : x != 0 -> `1-(x^-1) = (x - 1) / x.
 Proof. by move=> ?; rewrite mulrDl divff// mulN1r. Qed.
 
diff --git a/theories/derive.v b/theories/derive.v
index 3e9e0608f..5b2a0d643 100644
--- a/theories/derive.v
+++ b/theories/derive.v
@@ -1267,6 +1267,12 @@ Proof. by have /derivableP := @derivable_id x v; rewrite derive_val. Qed.
 
 End derive_id.
 
+Lemma derive1_onem {R : numFieldType} :
+  (fun x => 1 - x : R^o)^`()%classic = cst (-1).
+Proof.
+by apply/funext => x; rewrite derive1E deriveB// derive_id derive_cst sub0r.
+Qed.
+
 Section Derive_lemmasVR.
 Variables (R : numFieldType) (V : normedModType R).
 Implicit Types (f g : V -> R) (x v : V).
@@ -1332,6 +1338,15 @@ rewrite scalerA -scalerDl mulrCA -[f x * _]exprS.
 by rewrite [in LHS]mulr_natl exprfctE -mulrSr mulr_natl.
 Qed.
 
+(*Global Instance is_deriveX f n x v (df : R) :
+  is_derive x v f df -> is_derive x v (f ^+ n.+1) ((n.+1%:R * f x ^+n) *: df).
+Proof.
+move=> dfx; elim: n => [|n ihn]; first by rewrite expr1 expr0 mulr1 scale1r.
+rewrite exprS; apply: is_derive_eq.
+rewrite scalerA -scalerDl mulrCA -[f x * _]exprS.
+by rewrite [in LHS]mulr_natl exprfctE -mulrSr mulr_natl.
+Qed.*)
+
 Lemma derivableX f n x v : derivable f x v -> derivable (f ^+ n) x v.
 Proof. by case: n => [_|n /derivableP]; [rewrite expr0|]. Qed.
 
@@ -1339,6 +1354,10 @@ Lemma deriveX f n x v : derivable f x v ->
   'D_v (f ^+ n) x = (n%:R * f x ^+ n.-1) *: 'D_v f x.
 Proof. by move=> /derivableP df; rewrite derive_val. Qed.
 
+(*Lemma deriveX f n x v : derivable f x v ->
+  'D_v (f ^+ n.+1) x = (n.+1%:R * f x ^+ n) *: 'D_v f x.
+Proof. by move=> /derivableP df; rewrite derive_val. Qed.*)
+
 Fact der_inv f x v : f x != 0 -> derivable f x v ->
   (fun h => h^-1 *: (((fun y => (f y)^-1) \o shift x) (h *: v) - (f x)^-1)) @
   0^' --> - (f x) ^-2 *: 'D_v f x.
diff --git a/theories/ftc.v b/theories/ftc.v
index 157a707d5..1e5b940be 100644
--- a/theories/ftc.v
+++ b/theories/ftc.v
@@ -1793,6 +1793,46 @@ Qed.
 
 End integration_by_substitution.
 
+Section integration_by_substitution_onem.
+Context {R : realType}.
+Let mu := (@lebesgue_measure R).
+Local Open Scope ereal_scope.
+
+Lemma integration_by_substitution_onem (G : R -> R) (r : R) :
+  (0 < r <= 1)%R ->
+  {within `[0%R, r], continuous G} ->
+  (\int[mu]_(x in `[0%R, r]) (G x)%:E =
+  \int[mu]_(x in `[(1 - r)%R, 1%R]) (G (1 - x))%:E).
+Proof.
+move=> r01 cG.
+have := @integration_by_substitution_decreasing R (fun x => 1 - x)%R G (1 - r) 1.
+rewrite subKr subrr => -> //.
+- by apply: eq_integral => x xr; rewrite !fctE derive1_onem opprK mulr1.
+- by rewrite ltrBlDl ltrDr; case/andP : r01.
+- by move=> x y _ _ xy; rewrite ler_ltB.
+- by rewrite derive1_onem; move=> ? ?; exact: cvg_cst.
+- by rewrite derive1_onem; exact: is_cvg_cst.
+- by rewrite derive1_onem; exact: is_cvg_cst.
+- split => /=.
+  + by move=> x xr1; exact: derivableB.
+  + apply: cvg_at_right_filter; rewrite subKr.
+    apply: (@continuous_comp_cvg _ R^o R^o _ (fun x => 1 - x)%R)=> //=.
+      by move=> x; apply: (@continuousB _ R^o)  => //; exact: cvg_cst.
+    by under eq_fun do rewrite subKr; exact: cvg_id.
+  + by apply: cvg_at_left_filter; apply: (@cvgB _ R^o) => //; exact: cvg_cst.
+Qed.
+
+Lemma Rintegration_by_substitution_onem (G : R -> R) (r : R) :
+  (0 < r <= 1)%R ->
+  {within `[0%R, r], continuous G} ->
+  (\int[mu]_(x in `[0%R, r]) (G x) =
+  \int[mu]_(x in `[(1 - r)%R, 1%R]) (G (1 - x)))%R.
+Proof.
+by move=> r01 cG; rewrite [in LHS]/Rintegral integration_by_substitution_onem.
+Qed.
+
+End integration_by_substitution_onem.
+
 Section ge0_integralT_even.
 Context {R : realType}.
 Let mu := @lebesgue_measure R.
diff --git a/theories/normedtype_theory/pseudometric_normed_Zmodule.v b/theories/normedtype_theory/pseudometric_normed_Zmodule.v
index 07329928f..dcbdec1bc 100644
--- a/theories/normedtype_theory/pseudometric_normed_Zmodule.v
+++ b/theories/normedtype_theory/pseudometric_normed_Zmodule.v
@@ -817,6 +817,20 @@ Arguments cvgr_neq0 {R V T F FF f}.
 #[global] Hint Extern 0 (ProperFilter _^'+) =>
   (apply: at_right_proper_filter) : typeclass_instances.
 
+Lemma continuous_comp_cvg {R : numFieldType} (U V : pseudoMetricNormedZmodType R)
+  (f : U -> V) (g : R -> U) (r : R) (l : V) : continuous f ->
+  (f \o g) x @[x --> r] --> l -> f x @[x --> g r] --> l.
+Proof.
+move=> cf fgl; apply/(@cvgrPdist_le _ V) => /= e e0.
+have e20 : 0 < e / 2 by rewrite divr_gt0.
+move/(@cvgrPdist_le _ V) : fgl => /(_ _ e20) fgl.
+have /(@cvgrPdist_le _ V) /(_ _ e20) fgf := cf (g r).
+rewrite !near_simpl/=; near=> t.
+rewrite -(@subrK _ (f (g r)) l) -(addrA (_ + _)) (le_trans (ler_normD _ _))//.
+rewrite (splitr e) lerD//; last by near: t.
+by case: fgl => d /= d0; apply; rewrite /ball_/= subrr normr0.
+Unshelve. all: by end_near. Qed.
+
 Section closure_left_right_open.
 Variable R : realFieldType.
 Implicit Types z : R.
diff --git a/theories/probability.v b/theories/probability.v
index 6927b7e6e..f86d5a8de 100644
--- a/theories/probability.v
+++ b/theories/probability.v
@@ -8,7 +8,7 @@ From mathcomp Require Import exp numfun lebesgue_measure lebesgue_integral.
 From mathcomp Require Import reals interval_inference ereal topology normedtype.
 From mathcomp Require Import sequences derive esum measure exp trigo realfun.
 From mathcomp Require Import numfun lebesgue_measure lebesgue_integral kernel.
-From mathcomp Require Import ftc gauss_integral hoelder.
+From mathcomp Require Import charge ftc gauss_integral hoelder.
 
 (**md**************************************************************************)
 (* # Probability                                                              *)
@@ -47,7 +47,7 @@ From mathcomp Require Import ftc gauss_integral hoelder.
 (*                                                                            *)
 (* ```                                                                        *)
 (*      bernoulli_pmf p == Bernoulli pmf with parameter p : R                 *)
-(*          bernoulli p == Bernoulli probability measure when 0 <= p <= 1     *)
+(*     bernoulli_prob p == Bernoulli probability measure when 0 <= p <= 1     *)
 (*                         and \d_false otherwise                             *)
 (*     binomial_pmf n p == binomial pmf with parameters n : nat and p : R     *)
 (*    binomial_prob n p == binomial probability measure when 0 <= p <= 1      *)
@@ -68,6 +68,11 @@ From mathcomp Require Import ftc gauss_integral hoelder.
 (*   exponential_prob r == exponential probability measure                    *)
 (*      poisson_pmf r k == pmf of the Poisson distribution with parameter r   *)
 (*       poisson_prob r == Poisson probability measure                        *)
+(*          XMonemX a b := x ^+ a * `1-x ^+ b                                 *)
+(*         beta_fun a b := \int[mu]_x (XMonemX a.-1 b.-1 \_`[0,1] x)          *)
+(*             beta_pdf == probability density function for beta              *)
+(*            beta_prob == beta probability measure                           *)
+(* div_beta_fun a b c d := beta_fun (a + c) (b + d) / beta_fun a b            *)
 (* ```                                                                        *)
 (*                                                                            *)
 (******************************************************************************)
@@ -1037,13 +1042,16 @@ Proof.
 by apply/measurable_funTS/measurable_fun_if => //=; exact: measurable_funB.
 Qed.
 
-Definition bernoulli {R : realType} (p : R) : set bool -> \bar R := fun A =>
-  if (0 <= p <= 1)%R then \sum_(b \in A) (bernoulli_pmf p b)%:E else \d_false A.
+Definition bernoulli_prob {R : realType} (p : R) : set bool -> \bar R :=
+  fun A => if (0 <= p <= 1)%R then
+    \sum_(b \in A) (bernoulli_pmf p b)%:E
+  else
+    \d_false A.
 
 Section bernoulli.
 Context {R : realType} (p : R).
 
-Local Notation bernoulli := (bernoulli p).
+Local Notation bernoulli := (bernoulli_prob p).
 
 Let bernoulli0 : bernoulli set0 = 0.
 Proof.
@@ -1109,10 +1117,10 @@ Section bernoulli_measure.
 Context {R : realType}.
 Variables (p : R) (p0 : (0 <= p)%R) (p1 : ((NngNum p0)%:num <= 1)%R).
 
-Lemma bernoulli_dirac : bernoulli p = measure_add
+Lemma bernoulli_dirac : bernoulli_prob p = measure_add
   (mscale (NngNum p0) \d_true) (mscale (1 - (Itv01 p0 p1)%:num)%:nng \d_false).
 Proof.
-apply/funext => U; rewrite /bernoulli; case: ifPn => [p01|]; last first.
+apply/funext => U; rewrite /bernoulli_prob; case: ifPn => [p01|]; last first.
   by rewrite p0/= p1.
 rewrite measure_addE/= /mscale/=.
 have := @subsetT _ U; rewrite setT_bool => UT.
@@ -1130,10 +1138,10 @@ have [->|->|->|->] /= := subset_set2 UT.
 Qed.
 
 End bernoulli_measure.
-Arguments bernoulli {R}.
+Arguments bernoulli_prob {R}.
 
 Lemma eq_bernoulliV2 {R : realType} (P : probability bool R) :
-  P [set true] = P [set false] -> P =1 bernoulli 2^-1.
+  P [set true] = P [set false] -> P =1 bernoulli_prob 2^-1.
 Proof.
 move=> Ptrue_eq_false; apply/eq_bernoulli.
 have : P [set: bool] = 1%E := probability_setT P.
@@ -1148,11 +1156,12 @@ Context {R : realType}.
 Variables (p : R) (p01 : (0 <= p <= 1)%R).
 Local Open Scope ereal_scope.
 
-Lemma bernoulliE A : bernoulli p A = p%:E * \d_true A + (`1-p)%:E * \d_false A.
+Lemma bernoulli_probE A :
+  bernoulli_prob p A = p%:E * \d_true A + (`1-p)%:E * \d_false A.
 Proof. by case/andP : p01 => p0 p1; rewrite bernoulli_dirac// measure_addE. Qed.
 
-Lemma integral_bernoulli (f : bool -> \bar R) : (forall x, 0 <= f x) ->
-  \int[bernoulli p]_y (f y) = p%:E * f true + (`1-p)%:E * f false.
+Lemma integral_bernoulli_prob (f : bool -> \bar R) : (forall x, 0 <= f x) ->
+  \int[bernoulli_prob p]_y (f y) = p%:E * f true + (`1-p)%:E * f false.
 Proof.
 move=> f0; case/andP : p01 => p0 p1; rewrite bernoulli_dirac/=.
 rewrite ge0_integral_measure_sum// 2!big_ord_recl/= big_ord0 adde0/=.
@@ -1166,8 +1175,8 @@ Local Open Scope ring_scope.
 Variable R : realType.
 Implicit Type p : R.
 
-Lemma measurable_bernoulli :
-  measurable_fun setT (bernoulli : R -> pprobability bool R).
+Lemma measurable_bernoulli_prob :
+  measurable_fun setT (bernoulli_prob : R -> pprobability bool R).
 Proof.
 apply: (measurability (@pset _ _ _ : set (set (pprobability _ R)))) => //.
 move=> _ -[_ [r r01] [Ys mYs <-]] <-; apply: emeasurable_fun_infty_o => //=.
@@ -1181,14 +1190,15 @@ apply: emeasurable_sum => n; move=> k Ysk; apply/measurableT_comp => //.
 exact: measurable_bernoulli_pmf.
 Qed.
 
-Lemma measurable_bernoulli2 U : measurable U ->
-  measurable_fun setT (bernoulli ^~ U : R -> \bar R).
+Lemma measurable_bernoulli_prob2 U : measurable U ->
+  measurable_fun setT (bernoulli_prob ^~ U : R -> \bar R).
 Proof.
-by move=> ?; exact: (measurable_kernel (kprobability measurable_bernoulli)).
+move=> mU.
+exact: (measurable_kernel (kprobability measurable_bernoulli_prob)).
 Qed.
 
 End measurable_bernoulli.
-Arguments measurable_bernoulli {R}.
+Arguments measurable_bernoulli_prob {R}.
 
 Section binomial_pmf.
 Local Open Scope ring_scope.
@@ -1317,9 +1327,9 @@ End binomial_probability.
 
 Lemma integral_binomial_prob (R : realType) n p U : (0 <= p <= 1)%R ->
   (\int[binomial_prob n p]_y \d_(0 < y)%N U =
-  bernoulli (1 - `1-p ^+ n) U :> \bar R)%E.
+  bernoulli_prob (1 - `1-p ^+ n) U :> \bar R)%E.
 Proof.
-move=> /andP[p0 p1]; rewrite bernoulliE//=; last first.
+move=> /andP[p0 p1]; rewrite bernoulli_probE//=; last first.
   rewrite subr_ge0 exprn_ile1//=; [|exact/onem_ge0|exact/onem_le1].
   by rewrite lerBlDr addrC -lerBlDr subrr; exact/exprn_ge0/onem_ge0.
 rewrite (@integral_binomial _ n p _ _ (fun y => \d_(1 <= y)%N U))//.
@@ -2032,3 +2042,768 @@ apply: ge0_emeasurable_sum.
   by move=> k x/= [_ x01] _; rewrite lee_fin poisson_pmf_ge0.
 by move=> k Ysk; apply/measurableT_comp => //; exact: measurable_poisson_pmf.
 Qed.
+
+(**md lemmas about the function $x \mapsto (1 - x)^n$ *)
+(* TODO: move? *)
+Section about_onemXn.
+Context {R : realType}.
+Implicit Types x y : R.
+
+Lemma continuous_onemXn n x : {for x, continuous (fun y => `1-y ^+ n)}.
+Proof.
+apply: (@continuous_comp _ _ _ (@onem R) (fun x => x ^+ n)).
+  by apply: (@cvgB _ R^o); [exact: cvg_cst|exact: cvg_id].
+exact: exprn_continuous.
+Qed.
+
+Lemma onemXn_derivable n x : derivable (fun y => `1-y ^+ n) x 1.
+Proof.
+have := @derivableX _ _ (@onem R) n x 1.
+by rewrite fctE; apply; exact: derivableB.
+Qed.
+
+Lemma derivable_oo_continuous_bnd_onemXnMr n x :
+  derivable_oo_continuous_bnd (fun y => `1-y ^+ n * x) 0 1.
+Proof.
+split.
+- by move=> y y01; apply: derivableM => //=; exact: onemXn_derivable.
+- apply: cvgM; last exact: cvg_cst.
+  apply: cvg_at_right_filter.
+  apply: (@cvg_comp _ _ _ (fun x => `1-x) (fun x => x ^+ n)).
+    by apply: cvgB; [exact: cvg_cst|exact: cvg_id].
+  exact: exprn_continuous.
+- apply: cvg_at_left_filter.
+  apply: cvgM; last exact: cvg_cst.
+  apply: (@cvg_comp _ _ _ (fun x => `1-x) (fun x => x ^+ n)).
+    by apply: cvgB; [exact: cvg_cst|exact: cvg_id].
+  exact: exprn_continuous.
+Qed.
+
+Lemma derive_onemXn n x :
+  (fun y => `1-y ^+ n)^`()%classic x = - n%:R * `1-x ^+ n.-1.
+Proof.
+rewrite (@derive1_comp _ (@onem _) (fun x => x ^+ n))//; last first.
+  exact: exprn_derivable.
+rewrite derive1E exp_derive// derive1E deriveB// -derive1E.
+by rewrite derive1_cst derive_id sub0r mulrN1 [in RHS]mulNr scaler1.
+Qed.
+
+Lemma Rintegral_onemXn n :
+  \int[lebesgue_measure]_(x in `[0, 1]) (`1-x ^+ n) = n.+1%:R^-1 :> R.
+Proof.
+rewrite /Rintegral.
+rewrite (@continuous_FTC2 _ _ (fun x => `1-x ^+ n.+1 / - n.+1%:R))//=.
+- rewrite onem1 expr0n/= mul0r onem0 expr1n mul1r sub0r.
+  by rewrite -invrN -2!mulNrn opprK.
+- by apply: continuous_in_subspaceT => x x01; exact: continuous_onemXn.
+- exact: derivable_oo_continuous_bnd_onemXnMr.
+- move=> x x01; rewrite derive1Mr//; last exact: onemXn_derivable.
+  by rewrite derive_onemXn mulrAC divff// mul1r.
+Qed.
+
+End about_onemXn.
+
+(**md about the function $x \mapsto x^a  (1 - x)^b$ *)
+Section XMonemX.
+Context {R : numDomainType}.
+Implicit Types (x : R) (a b : nat).
+
+Definition XMonemX a b := fun x => x ^+ a * `1-x ^+ b.
+
+Lemma XMonemX_ge0 a b x : x \in `[0, 1] -> 0 <= XMonemX a b x.
+Proof.
+by rewrite in_itv=> /andP[? ?]; rewrite mulr_ge0 ?exprn_ge0 ?subr_ge0.
+Qed.
+
+Lemma XMonemX_le1 a b x : x \in `[0, 1] -> XMonemX a b x <= 1.
+Proof.
+rewrite in_itv/= => /andP[t0 t1].
+by rewrite mulr_ile1// ?(exprn_ge0,onem_ge0,exprn_ile1,onem_le1).
+Qed.
+
+Lemma XMonemX0n n x : XMonemX 0 n x = `1-x ^+ n.
+Proof. by rewrite /XMonemX/= expr0 mul1r. Qed.
+
+Lemma XMonemXn0 n x : XMonemX n 0 x = x ^+ n.
+Proof. by rewrite /XMonemX/= expr0 mulr1. Qed.
+
+Lemma XMonemX00 x : XMonemX 0 0 x = 1.
+Proof. by rewrite XMonemX0n expr0. Qed.
+
+Lemma XMonemXC a b x : XMonemX a b (1 - x) = XMonemX b a x.
+Proof. by rewrite /XMonemX [in LHS]/onem opprB addrCA subrr addr0 mulrC. Qed.
+
+Lemma XMonemXM a b a' b' x :
+  XMonemX a' b' x * XMonemX a b x = XMonemX (a + a') (b + b') x.
+Proof. by rewrite mulrCA -mulrA -exprD mulrA -exprD (addnC b'). Qed.
+
+End XMonemX.
+
+Section XMonemX_realType.
+Context {R : realType}.
+Local Notation XMonemX := (@XMonemX R).
+
+Lemma continuous_XMonemX a b : continuous (XMonemX a b).
+Proof.
+by move=> x; apply: cvgM; [exact: exprn_continuous|exact: continuous_onemXn].
+Qed.
+
+Lemma within_continuous_XMonemX A a b : {within A, continuous (XMonemX a b)}.
+Proof. by apply: continuous_in_subspaceT => x _; exact: continuous_XMonemX. Qed.
+
+Lemma measurable_XMonemX A a b : measurable_fun A (XMonemX a b).
+Proof.
+apply/measurable_funM => //; apply/measurable_funX => //.
+exact: measurable_funB.
+Qed.
+
+Lemma bounded_XMonemX a b : [bounded XMonemX a b x | x in `[0, 1]%classic].
+Proof.
+exists 1; split; [by rewrite num_real|move=> x x1 /= y y01].
+rewrite ger0_norm//; last by rewrite XMonemX_ge0.
+move: y01; rewrite in_itv/= => /andP[? ?].
+rewrite (le_trans _ (ltW x1))// mulr_ile1 ?exprn_ge0//.
+- by rewrite subr_ge0.
+- by rewrite exprn_ile1.
+- by rewrite exprn_ile1 ?subr_ge0// lerBlDl addrC -lerBlDl subrr.
+Qed.
+
+Local Notation mu := lebesgue_measure.
+
+Lemma integrable_XMonemX a b : mu.-integrable `[0, 1] (EFin \o XMonemX a b).
+Proof.
+apply: continuous_compact_integrable => //; first exact: segment_compact.
+by apply: continuous_in_subspaceT => x _; exact: continuous_XMonemX.
+Qed.
+
+Lemma integrable_XMonemX_restrict a b :
+  mu.-integrable [set: R] (EFin \o XMonemX a.-1 b.-1 \_`[0,1]).
+Proof.
+rewrite -restrict_EFin; apply/integrable_restrict => //=.
+by rewrite setTI; exact: integrable_XMonemX.
+Qed.
+
+Lemma integral_XMonemX_restrict U a b :
+  (\int[mu]_(x in U) (XMonemX a b \_ `[0, 1] x)%:E =
+   \int[mu]_(x in U `&` `[0%R, 1%R]) (XMonemX a b x)%:E)%E.
+Proof.
+rewrite [RHS]integral_mkcondr /=; apply: eq_integral => x xU /=.
+by rewrite restrict_EFin.
+Qed.
+
+End XMonemX_realType.
+
+Section beta_fun.
+Context {R : realType}.
+Notation mu := (@lebesgue_measure _).
+Local Open Scope ring_scope.
+Local Notation XMonemX := (@XMonemX R).
+
+Definition beta_fun a b : R := \int[mu]_x (XMonemX a.-1 b.-1 \_`[0,1]) x.
+
+Lemma EFin_beta_fun a b :
+  ((beta_fun a b)%:E = \int[mu]_x (XMonemX a.-1 b.-1 \_`[0,1] x)%:E)%E.
+Proof.
+rewrite fineK//; apply: integrable_fin_num => //=.
+under eq_fun.
+  move=> x.
+  rewrite /= -/((EFin \o ((XMonemX a.-1 b.-1) \_ _)) x) -restrict_EFin.
+  over.
+by apply/integrable_restrict => //=; rewrite setTI; exact: integrable_XMonemX.
+Qed.
+
+Lemma beta_fun_sym a b : beta_fun a b = beta_fun b a.
+Proof.
+rewrite -[LHS]Rintegral_mkcond Rintegration_by_substitution_onem//=.
+- rewrite subrr -[RHS]Rintegral_mkcond; apply: eq_Rintegral => x x01.
+  by rewrite XMonemXC.
+- by rewrite ltr01 lexx.
+- exact: within_continuous_XMonemX.
+Qed.
+
+Lemma beta_fun0n b : (0 < b)%N -> beta_fun 0 b = b%:R^-1.
+Proof.
+move=> b0; rewrite -[LHS]Rintegral_mkcond.
+under eq_Rintegral do rewrite XMonemX0n.
+by rewrite Rintegral_onemXn// prednK.
+Qed.
+
+Lemma beta_fun00 : beta_fun 0 0 = 1%R.
+Proof.
+rewrite -[LHS]Rintegral_mkcond.
+under eq_Rintegral do rewrite XMonemX00.
+rewrite Rintegral_cst//= mul1r lebesgue_measure_itv/= lte_fin ltr01.
+by rewrite oppr0 adde0.
+Qed.
+
+Lemma beta_fun1S b : beta_fun 1 b.+1 = b.+1%:R^-1.
+Proof.
+rewrite /beta_fun -Rintegral_mkcond.
+under eq_Rintegral do rewrite XMonemX0n.
+by rewrite Rintegral_onemXn.
+Qed.
+
+Lemma beta_fun11 : beta_fun 1 1 = 1.
+Proof. by rewrite (beta_fun1S O) invr1. Qed.
+
+Lemma beta_funSSS a b :
+  beta_fun a.+2 b.+1 = a.+1%:R / b.+1%:R * beta_fun a.+1 b.+2.
+Proof.
+rewrite -[LHS]Rintegral_mkcond.
+rewrite (@Rintegration_by_parts _ _ (fun x => `1-x ^+ b.+1 / - b.+1%:R)
+    (fun x => a.+1%:R * x ^+ a)); last 7 first.
+  exact: ltr01.
+  apply/continuous_subspaceT => x.
+  by apply: cvgM; [exact: cvg_cst|exact: exprn_continuous].
+  split.
+    by move=> x x01; exact: exprn_derivable.
+    by apply: cvg_at_right_filter; exact: exprn_continuous.
+    by apply: cvg_at_left_filter; exact: exprn_continuous.
+  by move=> x x01; rewrite derive1E exp_derive scaler1.
+  by apply/continuous_subspaceT => x x01; exact: continuous_onemXn.
+  exact: derivable_oo_continuous_bnd_onemXnMr.
+  move=> x x01; rewrite derive1Mr; last exact: onemXn_derivable.
+  by rewrite derive_onemXn mulrAC divff// mul1r.
+rewrite {1}/onem !(expr1n,mul1r,expr0n,subr0,subrr,mul0r,oppr0,sub0r)/=.
+transitivity (a.+1%:R / b.+1%:R * \int[mu]_(x in `[0, 1]) XMonemX a b.+1 x).
+  under [in LHS]eq_Rintegral.
+    move=> x x01.
+    rewrite mulrA mulrC mulrA (mulrA _ a.+1%:R) -(mulrA (_ * _)%R).
+    over.
+  rewrite /=.
+  rewrite RintegralZl//=; last exact: integrable_XMonemX.
+  by rewrite -mulNrn -2!mulNr -invrN -mulNrn opprK (mulrC _ a.+1%:R).
+by rewrite Rintegral_mkcond.
+Qed.
+
+Lemma beta_funSS a b : beta_fun a.+1 b.+1 =
+  a`!%:R / (\prod_(b.+1 <= i < (a + b).+1) i)%:R * beta_fun 1 (a + b).+1.
+Proof.
+elim: a b => [b|a ih b].
+  by rewrite fact0 mul1r add0n /index_iota subnn big_nil invr1 mul1r.
+rewrite beta_funSSS [in LHS]ih !mulrA; congr *%R; last by rewrite addSnnS.
+rewrite -mulrA mulrCA 2!mulrA.
+rewrite -natrM (mulnC a`!) -factS -mulrA -invfM; congr (_ / _).
+rewrite big_add1 [in RHS]big_nat_recl/=; last by rewrite addSn ltnS leq_addl.
+by rewrite -natrM addSnnS.
+Qed.
+
+Lemma beta_fun_fact a b :
+  beta_fun a.+1 b.+1 = (a`! * b`!)%:R / (a + b).+1`!%:R.
+Proof.
+rewrite beta_funSS beta_fun1S natrM -!mulrA; congr *%R.
+(* (b+1 b+2 ... b+1 b+a)^-1 / (a+b+1) = b! / (a+b+1)! *)
+rewrite factS [in RHS]mulnC natrM invfM mulrA; congr (_ / _).
+rewrite -(@invrK _ b`!%:R) -invfM; congr (_^-1).
+apply: (@mulfI _ b`!%:R).
+  by rewrite gt_eqF// ltr0n fact_gt0.
+rewrite mulrA divff// ?gt_eqF// ?ltr0n ?fact_gt0 ?mul1r//.
+rewrite [in RHS]fact_prod -natrM; congr (_%:R).
+rewrite fact_prod -big_cat/= /index_iota subn1 -iotaD.
+by rewrite subSS addnK subn1 addnC.
+Qed.
+
+Lemma beta_funE a b : beta_fun a b =
+  if (a == 0)%N && (0 < b)%N then
+    b%:R^-1
+  else if (b == 0)%N && (0 < a)%N then
+    a%:R^-1
+  else
+    a.-1`!%:R * b.-1`!%:R / (a + b).-1`!%:R.
+Proof.
+case: a => [|a].
+  rewrite eqxx/=; case: ifPn => [|].
+    by case: b => [|b _] //; rewrite beta_fun0n.
+  rewrite -leqNgt leqn0 => /eqP ->.
+  by rewrite beta_fun00 eqxx ltnn/= fact0 mul1r divr1.
+case: b => [|b].
+  by rewrite beta_fun_sym beta_fun0n// fact0 addn0/= mulr1 divff.
+by rewrite beta_fun_fact natrM// addnS.
+Qed.
+
+Lemma beta_fun_gt0 a b : 0 < beta_fun a b.
+Proof.
+rewrite beta_funE.
+case: ifPn => [/andP[_ b0]|]; first by rewrite invr_gt0 ltr0n.
+rewrite negb_and => /orP[a0|].
+  case: ifPn => [/andP[_]|]; first by rewrite invr_gt0// ltr0n.
+  rewrite negb_and => /orP[b0|].
+    by rewrite divr_gt0// ?mulr_gt0 ?ltr0n ?fact_gt0.
+  by rewrite -leqNgt leqn0 (negbTE a0).
+rewrite -leqNgt leqn0 => /eqP ->; rewrite eqxx/=.
+case: ifPn; first by rewrite invr_gt0 ltr0n.
+by rewrite -leqNgt leqn0 => /eqP ->; rewrite fact0 mul1r divr1.
+Qed.
+
+Lemma beta_fun_ge0 a b : 0 <= beta_fun a b.
+Proof. exact/ltW/beta_fun_gt0. Qed.
+
+End beta_fun.
+
+Section beta_pdf.
+Local Open Scope ring_scope.
+Context {R : realType}.
+Variables a b : nat.
+
+Local Notation XMonemX := (@XMonemX R).
+
+Definition beta_pdf t : R := XMonemX a.-1 b.-1 \_`[0,1] t / beta_fun a b.
+
+Lemma measurable_beta_pdf : measurable_fun [set: R] beta_pdf.
+Proof.
+apply: measurable_funM => //; apply/measurable_restrict => //.
+by rewrite setTI; exact: measurable_XMonemX.
+Qed.
+
+Lemma beta_pdf_ge0 t : 0 <= beta_pdf t.
+Proof.
+rewrite /beta_pdf divr_ge0 ?beta_fun_ge0//.
+rewrite patchE; case: ifPn => //=.
+by rewrite inE/= => ?; exact: XMonemX_ge0.
+Qed.
+
+Lemma beta_pdf_le_beta_funV x : beta_pdf x <= (beta_fun a b)^-1.
+Proof.
+rewrite /beta_pdf ler_pdivrMr ?beta_fun_gt0// mulVf ?gt_eqF ?beta_fun_gt0//.
+by rewrite patchE; case: ifPn => //; rewrite inE/= => ?; exact: XMonemX_le1.
+Qed.
+
+Local Notation mu := lebesgue_measure.
+
+Let int_beta_pdf01 : (\int[mu]_(x in `[0%R, 1%R]) (beta_pdf x)%:E =
+                      \int[mu]_x (beta_pdf x)%:E :> \bar R)%E.
+Proof.
+rewrite integral_mkcond/=; apply: eq_integral => /=x _.
+by rewrite /beta_pdf/= !patchE; case: ifPn => [->//|_]; rewrite mul0r.
+Qed.
+
+Lemma integrable_beta_pdf : mu.-integrable [set: _] (EFin \o beta_pdf).
+Proof.
+apply/integrableP; split.
+  by apply/measurable_EFinP; exact: measurable_beta_pdf.
+under eq_integral.
+  move=> /= x _.
+  rewrite ger0_norm//; last by rewrite beta_pdf_ge0.
+  over.
+simpl.
+rewrite -int_beta_pdf01.
+apply: (@le_lt_trans _ _ (\int[mu]_(x in `[0%R, 1%R]) (beta_fun a b)^-1%:E)%E).
+  apply: ge0_le_integral => //=.
+  - by move=> x _; rewrite lee_fin beta_pdf_ge0.
+  - apply/measurable_funTS/measurable_EFinP => /=.
+    exact: measurable_beta_pdf.
+  - by move=> x _; rewrite lee_fin invr_ge0// beta_fun_ge0.
+  - by move=> x _; rewrite lee_fin beta_pdf_le_beta_funV.
+rewrite integral_cst//= lebesgue_measure_itv//=.
+by rewrite lte01 oppr0 adde0 mule1 ltry.
+Qed.
+
+Lemma bounded_beta_pdf_01 : [bounded beta_pdf x | x in `[0%R, 1%R]%classic].
+Proof.
+exists (beta_fun a b)^-1; split; first by rewrite num_real.
+move=> // y y1.
+near=> M => /=.
+rewrite (le_trans _ (ltW y1))//.
+near: M => M /=; rewrite in_itv/= => /andP[M0 M1].
+rewrite ler_norml; apply/andP; split.
+  rewrite lerNl (@le_trans _ _ 0%R)// ?invr_ge0 ?beta_fun_ge0//.
+  by rewrite lerNl oppr0 beta_pdf_ge0.
+rewrite /beta_pdf ler_pdivrMr ?beta_fun_gt0//.
+rewrite mulVf ?gt_eqF ?beta_fun_gt0//.
+by rewrite patchE; case: ifPn => //; rewrite inE => ?; exact: XMonemX_le1.
+Unshelve. all: by end_near. Qed.
+
+End beta_pdf.
+
+(* TODO: move? *)
+Lemma invr_nonneg_proof (R : numDomainType) (p : {nonneg R}) :
+  (0 <= (p%:num)^-1)%R.
+Proof. by rewrite invr_ge0. Qed.
+
+Definition invr_nonneg (R : numDomainType) (p : {nonneg R}) :=
+  NngNum (invr_nonneg_proof p).
+
+Section beta.
+Local Open Scope ring_scope.
+Context {R : realType}.
+Variables a b : nat.
+
+Local Notation mu := (@lebesgue_measure R).
+Local Notation XMonemX := (@XMonemX R).
+
+Let beta_num (U : set (measurableTypeR R)) : \bar R :=
+  \int[mu]_(x in U) (XMonemX a.-1 b.-1 \_`[0,1] x)%:E.
+
+Let beta_numT : beta_num [set: (measurableTypeR R)] = (beta_fun a b)%:E.
+Proof. by rewrite /beta_num/= EFin_beta_fun. Qed.
+
+Let beta_num_lty U : measurable U -> (beta_num U < +oo)%E.
+Proof.
+move=> mU.
+apply: (@le_lt_trans _ _ (\int[mu]_(x in U `&` `[0%R, 1%R]) 1)%E); last first.
+  rewrite integral_cst//= ?mul1e//.
+    rewrite (le_lt_trans (measureIr _ _ _))//= lebesgue_measure_itv//= lte01//.
+    by rewrite EFinN sube0 ltry.
+  exact: measurableI.
+rewrite /beta_num integral_XMonemX_restrict ge0_le_integral//=.
+- exact: measurableI.
+- by move=> x [_ ?]; rewrite lee_fin XMonemX_ge0.
+- by apply/measurable_funTS/measurableT_comp => //; exact: measurable_XMonemX.
+- by move=> x [_ ?]; rewrite lee_fin XMonemX_le1.
+Qed.
+
+Let beta_num0 : beta_num set0 = 0%:E.
+Proof. by rewrite /beta_num integral_set0. Qed.
+
+Let beta_num_ge0 U : (0 <= beta_num U)%E.
+Proof.
+rewrite /beta_num integral_ge0//= => x Ux; rewrite lee_fin.
+by rewrite patchE; case: ifPn => //; rewrite inE/= => x01; exact: XMonemX_ge0.
+Qed.
+
+Let beta_num_sigma_additive : semi_sigma_additive beta_num.
+Proof.
+move=> /= F mF tF mUF; rewrite /beta_num; apply: cvg_toP.
+  apply: ereal_nondecreasing_is_cvgn => m n mn.
+  apply: lee_sum_nneg_natr => // k _ _; apply: integral_ge0 => /= x Fkx.
+  rewrite patchE; case: ifPn => //; rewrite inE/= => x01.
+  by rewrite lee_fin XMonemX_ge0.
+rewrite ge0_integral_bigcup//=.
+- apply/measurable_funTS/measurableT_comp => //=.
+  by apply/measurable_restrict => //=; rewrite setTI; exact: measurable_XMonemX.
+- move=> x [? _ ?]; rewrite lee_fin.
+  by rewrite patchE; case: ifPn => //; rewrite inE/= => x0; exact: XMonemX_ge0.
+Qed.
+
+HB.instance Definition _ := isMeasure.Build _ _ _ beta_num
+  beta_num0 beta_num_ge0 beta_num_sigma_additive.
+
+Definition beta_prob :=
+  @mscale _ _ _ (invr_nonneg (NngNum (beta_fun_ge0 a b))) beta_num.
+
+HB.instance Definition _ := Measure.on beta_prob.
+
+Let beta_prob_setT : beta_prob setT = 1%:E.
+Proof.
+rewrite /beta_prob /= /mscale /= beta_numT.
+by rewrite -EFinM mulVf// gt_eqF// beta_fun_gt0.
+Qed.
+
+HB.instance Definition _ :=
+  @Measure_isProbability.Build _ _ _ beta_prob beta_prob_setT.
+
+Lemma integral_beta_pdf U : measurable U ->
+  (\int[mu]_(x in U) (beta_pdf a b x)%:E = beta_prob U :> \bar R)%E.
+Proof.
+move=> mU.
+rewrite /beta_pdf.
+under eq_integral do rewrite EFinM/=.
+rewrite ge0_integralZr//=.
+- by rewrite /beta_prob/= /mscale/= muleC.
+- apply/measurable_funTS/measurableT_comp => //.
+  by apply/measurable_restrict => //=; rewrite setTI; exact: measurable_XMonemX.
+- move=> x Ux; rewrite patchE; case: ifPn => //; rewrite inE/= => x01.
+  by rewrite lee_fin XMonemX_ge0.
+- by rewrite lee_fin invr_ge0// beta_fun_ge0.
+Qed.
+
+Lemma beta_prob01 : beta_prob `[0, 1] = 1%:E.
+Proof.
+rewrite -beta_prob_setT/= /beta_prob/= /mscale/= /beta_num/=.
+rewrite [in RHS]integral_XMonemX_restrict/= setTI.
+by rewrite [in LHS]integral_XMonemX_restrict/= setIid.
+Qed.
+
+Lemma beta_prob_fin_num U : measurable U -> beta_prob U \is a fin_num.
+Proof.
+move=> mU; rewrite ge0_fin_numE//.
+rewrite (@le_lt_trans _ _ (beta_prob setT))//=.
+  by rewrite le_measure// inE.
+by rewrite probability_setT ltry.
+Qed.
+
+Lemma beta_prob_dom : beta_prob `<< mu.
+Proof.
+move=> A mA muA0; rewrite /beta_prob /mscale/=.
+apply/eqP; rewrite mule_eq0 eqe invr_eq0 gt_eqF/= ?beta_fun_gt0//; apply/eqP.
+rewrite /beta_num integral_XMonemX_restrict.
+apply/eqP; rewrite eq_le; apply/andP; split; last first.
+  by apply: integral_ge0 => x [_ ?]; rewrite lee_fin XMonemX_ge0.
+apply: (@le_trans _ _ (\int[mu]_(x in A `&` `[0%R, 1%R]) 1)%E); last first.
+  rewrite integral_cst ?mul1e//=; last exact: measurableI.
+  by rewrite -[leRHS]muA0 measureIl.
+apply: ge0_le_integral => //=; first exact: measurableI.
+- by move=> x [_ x01]; rewrite lee_fin XMonemX_ge0.
+- by apply/measurable_funTS/measurableT_comp => //; exact: measurable_XMonemX.
+- by move=> x [_ ?]; rewrite lee_fin XMonemX_le1.
+Qed.
+
+End beta.
+Arguments beta_prob {R}.
+
+Lemma beta_prob_uniform {R : realType} :
+  beta_prob 1 1 = uniform_prob (@ltr01 R).
+Proof.
+apply/funext => U.
+rewrite /beta_prob /uniform_prob.
+rewrite /mscale/= beta_fun11 invr1 !mul1e.
+rewrite integral_XMonemX_restrict integral_uniform_pdf.
+apply: eq_integral => /= x.
+rewrite inE => -[Ux/=]; rewrite in_itv/= => x10.
+rewrite /XMonemX !expr0 mul1r.
+by rewrite /uniform_pdf x10 subr0 invr1.
+Qed.
+
+Lemma integral_beta_prob_bernoulli_prob_lty {R : realType} a b (f : R -> R) U :
+  measurable_fun setT f ->
+  (forall x, x \in `[0, 1]%R -> 0 <= f x <= 1)%R ->
+  (\int[beta_prob a b]_x `|bernoulli_prob (f x) U| < +oo :> \bar R)%E.
+Proof.
+move=> mf /= f01.
+apply: (@le_lt_trans _ _ (\int[beta_prob a b]_x cst 1 x))%E.
+  apply: ge0_le_integral => //=.
+    apply: measurableT_comp => //=.
+    by apply: (measurableT_comp (measurable_bernoulli_prob2 _)).
+  by move=> x _; rewrite gee0_abs// probability_le1.
+by rewrite integral_cst//= mul1e -ge0_fin_numE// beta_prob_fin_num.
+Qed.
+
+Lemma integral_beta_prob_bernoulli_prob_onemX_lty {R : realType} n a b U :
+  (\int[beta_prob a b]_x `|bernoulli_prob (`1-x ^+ n) U| < +oo :> \bar R)%E.
+Proof.
+apply: integral_beta_prob_bernoulli_prob_lty => //=.
+  by apply: measurable_funX => //; exact: measurable_funB.
+move=> x; rewrite in_itv/= => /andP[x0 x1].
+rewrite exprn_ge0 ?subr_ge0//= exprn_ile1// ?subr_ge0//.
+by rewrite lerBlDl -lerBlDr subrr.
+Qed.
+
+Lemma integral_beta_prob_bernoulli_prob_onem_lty {R : realType} n a b U :
+  (\int[beta_prob a b]_x `|bernoulli_prob (1 - `1-x ^+ n) U| < +oo :> \bar R)%E.
+Proof.
+apply: integral_beta_prob_bernoulli_prob_lty => //=.
+  apply: measurable_funB => //.
+  by apply: measurable_funX => //; exact: measurable_funB.
+move=> x; rewrite in_itv/= => /andP[x0 x1].
+rewrite -lerBlDr opprK add0r.
+rewrite andbC lerBlDl -lerBlDr subrr.
+rewrite exprn_ge0 ?subr_ge0//= exprn_ile1// ?subr_ge0//.
+by rewrite lerBlDl -lerBlDr subrr.
+Qed.
+
+Lemma beta_prob_integrable {R :realType} a b a' b' :
+  (beta_prob a b).-integrable `[0, 1]
+    (fun x : measurableTypeR R => (XMonemX a' b' x)%:E).
+Proof.
+apply/integrableP; split.
+  by apply/measurableT_comp => //; exact: measurable_XMonemX.
+apply: (@le_lt_trans _ _ (\int[beta_prob a b]_(x in `[0%R, 1%R]) 1)%E).
+  apply: ge0_le_integral => //=.
+    by do 2 apply/measurableT_comp => //; exact: measurable_XMonemX.
+  move=> x; rewrite in_itv/= => /andP[x0 x1].
+  rewrite lee_fin ger0_norm; last first.
+    by rewrite !mulr_ge0// exprn_ge0// onem_ge0.
+  rewrite mulr_ile1// ?exprn_ge0 ?onem_ge0// exprn_ile1// ?onem_ge0//.
+  exact: onem_le1.
+rewrite integral_cst//= mul1e.
+by rewrite -ge0_fin_numE// beta_prob_fin_num.
+Qed.
+
+Lemma beta_prob_integrable_onem {R : realType} a b a' b' :
+  (beta_prob a b).-integrable `[0, 1]
+    (fun x : measurableTypeR R => (`1-(XMonemX a' b' x))%:E).
+Proof.
+apply: (eq_integrable _ (cst 1 \- (fun x : measurableTypeR R =>
+  (XMonemX a' b' x)%:E))%E) => //.
+apply: integrableB => //=.
+  apply/integrableP; split => //.
+  rewrite (eq_integral (fun x => (\1_setT x)%:E))/=; last first.
+    by move=> x _; rewrite /= indicT normr1.
+  rewrite integral_indic//= setTI /beta_prob /mscale/= lte_mul_pinfty//.
+    by rewrite lee_fin invr_ge0 beta_fun_ge0.
+  rewrite (_ : integral _ _ _ = \int[lebesgue_measure]_x
+    (((@XMonemX R a.-1 b.-1) \_ `[0, 1]) x)%:E)%E; last first.
+    rewrite integral_mkcond/=; apply: eq_integral => /= x _.
+    by rewrite !patchE; case: ifPn => // ->.
+  have /integrableP[_] := @integrable_XMonemX_restrict R a b.
+  under eq_integral.
+    move=> x _.
+    rewrite gee0_abs//; last first.
+      rewrite lee_fin patchE; case: ifPn => //; rewrite inE/= => x01.
+      by rewrite XMonemX_ge0.
+    over.
+  by [].
+exact: beta_prob_integrable.
+Qed.
+
+Lemma beta_prob_integrable_dirac {R : realType} a b a' b' (c : bool) U :
+  (beta_prob a b).-integrable `[0, 1]
+    (fun x : measurableTypeR R => (XMonemX a' b' x)%:E * \d_c U)%E.
+Proof.
+apply: integrableMl => //=; last first.
+  exists 1; split => // x x1/= _ _; rewrite (le_trans _ (ltW x1))//.
+  by rewrite ger0_norm// indicE; case: (_ \in _).
+exact: beta_prob_integrable.
+Qed.
+
+Lemma beta_prob_integrable_onem_dirac {R : realType} a b a' b' (c : bool) U :
+  (beta_prob a b).-integrable `[0, 1]
+    (fun x : measurableTypeR R => (`1-(XMonemX a' b' x))%:E * \d_c U)%E.
+Proof.
+apply: integrableMl => //=; last first.
+  exists 1; split => // x x1/= _ _; rewrite (le_trans _ (ltW x1))//.
+  by rewrite ger0_norm// indicE; case: (_ \in _).
+exact: beta_prob_integrable_onem.
+Qed.
+
+Section integral_beta_prob.
+Context {R : realType}.
+Local Notation mu := (@lebesgue_measure R).
+
+Lemma integral_beta_prob a b f U : measurable U -> measurable_fun U f ->
+  (\int[beta_prob a b]_(x in U) `|f x| < +oo)%E ->
+  (\int[beta_prob a b]_(x in U) f x =
+   \int[mu]_(x in U) (f x * (beta_pdf a b x)%:E))%E.
+Proof.
+move=> mU mf finf.
+rewrite -(Radon_Nikodym_change_of_variables (beta_prob_dom a b))//=; last first.
+  by apply/integrableP; split.
+apply: ae_eq_integral => //.
+- apply: emeasurable_funM => //; apply: (measurable_int mu).
+  apply: (integrableS _ _ (@subsetT _ _)) => //=.
+  by apply: Radon_Nikodym_integrable; exact: beta_prob_dom.
+- apply: emeasurable_funM => //=; apply/measurableT_comp => //=.
+  by apply/measurable_funTS; exact: measurable_beta_pdf.
+- apply: ae_eqe_mul2l => /=.
+  rewrite Radon_NikodymE//=; first exact: beta_prob_dom.
+  move=> ?.
+  case: cid => /= h [h1 h2 h3].
+(* uniqueness of Radon-Nikodym derivative up to equal on non null sets of mu *)
+  apply: integral_ae_eq => //=.
+  + apply: integrableS h2 => //.
+    apply/measurable_funTS/measurableT_comp => //.
+    exact: measurable_beta_pdf.
+  + by move=> E E01 mE; rewrite -h3//= integral_beta_pdf.
+Qed.
+
+End integral_beta_prob.
+
+Section beta_prob_bernoulliE.
+Context {R : realType}.
+Local Notation mu := lebesgue_measure.
+Local Open Scope ring_scope.
+
+Definition div_beta_fun a b c d : R := beta_fun (a + c) (b + d) / beta_fun a b.
+
+Lemma div_beta_fun_ge0 a b c d : 0 <= div_beta_fun a b c d.
+Proof. by rewrite /div_beta_fun divr_ge0// beta_fun_ge0. Qed.
+
+Lemma div_beta_fun_le1 a b c d : (0 < a)%N -> (0 < b)%N ->
+  div_beta_fun a b c d <= 1.
+Proof.
+move=> a0 b0.
+rewrite /div_beta_fun ler_pdivrMr// ?mul1r ?beta_fun_gt0//.
+rewrite !beta_funE.
+rewrite addn_eq0 (gtn_eqF a0)/=.
+rewrite addn_eq0 (gtn_eqF b0)/=.
+rewrite ler_pdivrMr ?ltr0n ?fact_gt0//.
+rewrite mulrAC.
+rewrite ler_pdivlMr ?ltr0n ?fact_gt0//.
+rewrite -!natrM ler_nat.
+move: a a0 => [//|a _].
+rewrite addSn.
+move: b b0 => [//|b _].
+rewrite [(a + c).+1.-1]/=.
+rewrite [a.+1.-1]/=.
+rewrite [b.+1.-1]/=.
+rewrite addnS.
+rewrite [(_ + b).+1.-1]/=.
+rewrite (addSn b d).
+rewrite [(b + _).+1.-1]/=.
+rewrite (addSn (a + c)).
+rewrite [_.+1.-1]/=.
+rewrite addSn addnS.
+by rewrite leq_fact2// leq_addr.
+Qed.
+
+Definition beta_prob_bernoulli_prob a b c d U : \bar R :=
+  \int[beta_prob a b]_(y in `[0, 1])
+    bernoulli_prob ((@XMonemX R c d \_`[0, 1])%R y) U.
+
+Lemma beta_prob_bernoulli_probE a b c d U : (a > 0)%N -> (b > 0)%N ->
+  beta_prob_bernoulli_prob a b c d U = bernoulli_prob (div_beta_fun a b c d) U.
+Proof.
+move=> a0 b0.
+rewrite /beta_prob_bernoulli_prob.
+under eq_integral => x.
+  rewrite inE/= in_itv/= => x01.
+  rewrite bernoulli_probE/=; last first.
+  rewrite patchE; case: ifPn => x01'.
+    by rewrite XMonemX_ge0//= XMonemX_le1.
+  by rewrite lexx ler01.
+  over.
+rewrite /= [in RHS]bernoulli_probE/= ?div_beta_fun_ge0 ?div_beta_fun_le1//=.
+under eq_integral => x x01.
+  rewrite patchE x01/=.
+  over.
+rewrite /= integralD//=; last 2 first.
+  exact: beta_prob_integrable_dirac.
+  exact: beta_prob_integrable_onem_dirac.
+congr (_ + _).
+  rewrite integralZr//=; last exact: beta_prob_integrable.
+  congr (_ * _)%E.
+  rewrite integral_beta_prob//; last 2 first.
+    by apply/measurableT_comp => //; exact: measurable_XMonemX.
+    by have /integrableP[_] := @beta_prob_integrable R a b c d.
+  rewrite /beta_pdf.
+  under eq_integral do rewrite EFinM -muleA muleC -muleA.
+  rewrite /=.
+  transitivity ((beta_fun a b)^-1%:E * \int[mu]_(x in `[0%R, 1%R])
+      (@XMonemX R (a + c).-1 (b + d).-1 \_`[0,1] x)%:E)%E.
+    rewrite -integralZl//=; last first.
+      by apply: (integrableS measurableT) => //=; exact: integrable_XMonemX_restrict.
+    apply: eq_integral => x x01.
+    rewrite muleA muleC muleA -(EFinM (x ^+ c)) -/(XMonemX c d x) -EFinM mulrA.
+    rewrite !patchE x01 XMonemXM// -EFinM mulrC.
+    by move: a a0 b b0 => [//|a] _ [|b].
+  rewrite /div_beta_fun mulrC EFinM; congr (_ * _)%E.
+  rewrite /beta_fun integral_mkcond/= fineK; last first.
+    by apply: integrable_fin_num => //; exact: integrable_XMonemX_restrict.
+  by apply: eq_integral => /= x _; rewrite !patchE; case: ifPn => // ->.
+under eq_integral do rewrite muleC.
+rewrite /= integralZl//=; last exact: beta_prob_integrable_onem.
+rewrite muleC; congr (_ * _)%E.
+rewrite integral_beta_prob//=; last 2 first.
+  apply/measurableT_comp => //=.
+  by apply/measurable_funB => //; exact: measurable_XMonemX.
+  by have /integrableP[] := @beta_prob_integrable_onem R a b c d.
+rewrite /beta_pdf.
+under eq_integral do rewrite EFinM muleA.
+rewrite integralZr//=; last first.
+  apply: integrableMr => //=.
+  - by apply/measurable_funB => //=; exact: measurable_XMonemX.
+  - apply/ex_bound => //.
+    + apply: (@globally_properfilter _ _ 0%R) => //=.
+      by apply: inferP; rewrite in_itv/= lexx ler01.
+    + exists 1 => t.
+      rewrite /= in_itv/= => t01.
+      apply: normr_onem; apply/andP; split.
+        by rewrite mulr_ge0// exprn_ge0// ?onem_ge0//; case/andP: t01.
+      by rewrite mulr_ile1// ?exprn_ge0 ?exprn_ile1// ?onem_ge0 ?onem_le1//;
+        case/andP: t01.
+  - exact: integrableS (integrable_XMonemX_restrict _ _).
+transitivity ((\int[mu]_x ((@XMonemX R a.-1 b.-1 \_`[0,1] x)%:E -
+   (@XMonemX R (a + c).-1 (b + d).-1 \_`[0,1] x)%:E)) * (beta_fun a b)^-1%:E)%E.
+  congr (_ * _)%E; rewrite integral_mkcond/=; apply: eq_integral => x _.
+  rewrite !patchE; case: ifPn => [->|]; last by rewrite EFinN subee.
+  rewrite /onem -EFinM mulrBl mul1r EFinB EFinN; congr (_ - _)%E.
+  rewrite XMonemXM.
+  by move: a a0 b b0 => [|a]// _ [|b].
+rewrite integralB_EFin//=; last 2 first.
+  exact: integrableS (integrable_XMonemX_restrict _ _).
+  exact: integrableS (integrable_XMonemX_restrict _ _).
+rewrite EFinB muleBl//; last by rewrite -!EFin_beta_fun.
+by rewrite -!EFin_beta_fun -EFinM divff// gt_eqF// beta_fun_gt0.
+Qed.
+
+End beta_prob_bernoulliE.