TransformedDistribution using integration with parameter doesn't work

Question

I would like to ask for kind help with a following problem using TransformedDistribution in Mathematica. First here is the beginning of the code, that works:

n = 5;
k = 4;
f1 = OrderDistribution[{UniformDistribution[{0, 1}], n}, k]
f2 = OrderDistribution[{f1, n}, k]
w[v_?NumberQ] := NIntegrate[(v - b2)*PDF[f2, b2], {b2, 0, v}]
b[v_?NumberQ] := v - w[v]

The above code works so far, and getting any evaluated number from function b[v] works as well, for example when I write b[0.5] it gives me a number 0.49981412080261206, which is ok. I will not explain the context, but what I want to do is get a transformed distribution:

f3 := TransformedDistribution[b[v], v \[Distributed] f1]

So I want to transform each point in f1 by the function b[v]. Unfortunately this doesn't work, since when I try to plot pdf of f3 I get an empty graph, also checking the domain of f3 gives me [-inf,+inf], while clearly the domain here should be [0,1], since the basic function f1 is on [0,1], and the transformation is not decreasing each point by more than it's value. In short, there must be something wrong with the way I am using TransformedDistribution in this context. I am guessing the problem arises, because of how the integral limit depends on the parameter itself. I would very much appreciate if I could get any help on how to fix this issue, and how to make it work, so that I can get the transformed distribution.

Answer 1

n = 5;
k = 4;

f1 = OrderDistribution[{UniformDistribution[{0, 1}], n}, k];

f2 = OrderDistribution[{f1, n}, k];

The definition of w[v] can be done exactly with Integrate rather than NIntegrate

w[v_] = Assuming[Element[v, Reals],
   Integrate[(v - b2)*PDF[f2, b2], {b2, 0, v}] // Simplify];

b[v_] = v - w[v] // Simplify

f3 = TransformedDistribution[b[v], v \[Distributed] f1];

The mean and standard deviation for f3 are

{m, s} = (#[f3] // N) & /@ {Mean, StandardDeviation}

(*  {0.628867, 0.140471}  *)

The PDF for f3 is zero outside of the interval

PDF[f3, x][[1, 1, -1]]

(*  0 < x < 51351184/66927861  *)

% // N

(*  0. < x < 0.767262  *)

Plot[Evaluate@PDF[f3, x], {x, 0, 51351184/66927861}, Frame -> True, 
 Mesh -> {Range[m - 4 s, m, s], None}, 
 MeshStyle -> {Red, AbsolutePointSize[4]}, Exclusions -> True]

NIntegrate[PDF[f3, x], {x, 0, 51351184/66927861 - 10^-75}, 
 WorkingPrecision -> 50, MaxRecursion -> 100]

(*  0.99999999999999999999999999999879953560281933889629  *)