5
$\begingroup$

I am trying to use TransformedDistribution for a data analysis task. The distribution performs well with RandomVariate applied. It is slow for CDF, and unusable with PDF where it cannot return a single value given 5 minutes. I would be grateful for any advice as to how I can make this work.

expr = ArcTan[Cos[α] Cos[θ] - Cos[ϕ] Sin[α] Sin[θ], 
              Sqrt[(Cos[θ] Sin[α] + Cos[α] Cos[ϕ] Sin[θ])^2 + Sin[θ]^2 Sin[ϕ]^2]];

dist = TransformedDistribution[expr /. {θ -> 10. °, α -> 1. °},
                               ϕ \[Distributed] UniformDistribution[{-π, π}]];

(* no problem with RandomVariate *)
Timing[t1 = RandomVariate[dist, 100000]/Degree;]

(* {0.`,Null} *)

Histogram[t1, 50, "PDF", Frame -> True]

enter image description here

(* CDF takes a LONG time *)
Timing[CDF[dist, 10 Degree]]

(* {2.6676171`,0.4842396007961177`} *)

(* PDF takes forever *)
TimeConstrained[PDF[dist, 10. Degree], 5 60]

(* $Aborted *)
$\endgroup$
  • 2
    $\begingroup$ I think I'll need to toss my two responses so far because I have an exact answer. I should be done with that in a few minutes. $\endgroup$ – JimB Apr 2 '19 at 23:39
  • $\begingroup$ @JimB This is very useful. Thank you. $\endgroup$ – David Keith Apr 2 '19 at 23:52
4
$\begingroup$

Third times the charm. Here is what I think is the "exact" answer (and should work with other values of $\theta$ and $\alpha$.

The problem was getting an inverse function for determining $\phi$ from expr (which I'm now calling $X$ and putting $X$ in terms of degrees). The trick was to solve using the tangent of $X$ rather than just $X$ itself.

(* Define the transformation of ϕ *)
X = ArcTan[Cos[α] Cos[θ] - Cos[ϕ] Sin[α] Sin[θ],
    Sqrt[(Cos[θ] Sin[α] + Cos[α] Cos[ϕ] Sin[θ])^2 + Sin[θ]^2 Sin[ϕ]^2]]/Degree;

(* Set specific values for θ and α *)
parms = {θ -> 10 Degree, α -> 1 Degree};

(* X ranges from ... *)
xLower = X /. ϕ -> π /. parms // N
(* 9. *)
xUpper = X /. ϕ -> 0 /. parms // N
(* 11. *)

(* Solve for ϕ in terms of x *)
sol = Flatten[Solve[Tan[X Degree] == Tan[x Degree], ϕ]][[2]]

$$\phi \to \cos ^{-1}\left(\frac{1}{2} \csc (\alpha ) \csc (\theta ) (\cos (\alpha -\theta )+\cos (\alpha +\theta )-2 \cos (x {}^{\circ}))\right)$$

cdf = (π - ϕ /. sol /. parms)/π

$$\frac{\pi -\cos ^{-1}\left(\frac{1}{2} \csc ({1}^{\circ}) \csc (10 {}^{\circ}) (-2 \cos (x {}^{\circ})+\cos (9 {}^{\circ})+\cos (11 {}^{\circ}))\right)}{\pi }$$

pdf = D[cdf, x]

$$\frac{(\pi/180) \csc ({1}^{\circ}) \csc (10 {}^{\circ}) \sin (x {}^{\circ})}{\pi \sqrt{1-\frac{1}{4} \csc ^2({1}^{\circ}) \csc ^2(10 {}^{\circ}) (-2 \cos (x {}^{\circ})+\cos (9 {}^{\circ})+\cos (11 {}^{\circ}))^2}}$$

Plot[cdf, {x, xLower, xUpper}, PlotLabel -> "CDF"]
Plot[pdf, {x, xLower, xUpper}, PlotLabel -> "PDF"]

CDF

PDF

| improve this answer | |
$\endgroup$
  • $\begingroup$ This is exactly what I needed. In fact, I tried to solve the expression first, but did to think of the Tan[x] method you used. Thank you. $\endgroup$ – David Keith Apr 3 '19 at 19:40
  • $\begingroup$ That's nice work! +1 $\endgroup$ – ciao Apr 3 '19 at 22:09
  • 2
    $\begingroup$ @ciao Thanks. I learned a lot from answers David Keith gave on Wolfram Community so I'm happy I could give back a little. $\endgroup$ – JimB Apr 4 '19 at 2:39
3
$\begingroup$

I think you have to do it numerically. Here's one way using Interpolation. (And I had to use Quiet on the plot of the PDF because of issues I don't understand.) I also haven't included (yet) a rationale for construction of the cdf. Maybe later today.

n = 100;
t = Flatten[{{{11, 0}}, Table[{(expr /. {θ -> 10 Degree, α -> 1 Degree})/Degree, ϕ},
    {ϕ, π/n, π - π/n, π/n}], {{9, π}}}, 1];
f = Interpolation[t];
cdf[x_] := Piecewise[{{0, x <= 9}, {(π - f[x])/π, 9 < x < 11}, {1, x >= 11}}]
pdf[x_] := Piecewise[{{0, x <= 9}, {Max[0, cdf'[x]], 9 < x < 11}, {0, x >= 11}}]
Plot[cdf[x], {x, 8, 12}, PlotLabel -> "CDF"]
Quiet[Plot[pdf[x], {x, 8, 12}, PlotLabel -> "PDF"]]

CDF

PDF

| improve this answer | |
$\endgroup$
3
$\begingroup$

This is an extended comment rather than an answer.

I wonder if the underlying distribution is either a linearly transformed beta distribution or could be very well approximated by such a distribution. Here's why:

(* Define original distribution *)
expr = ArcTan[Cos[α] Cos[θ] - Cos[ϕ] Sin[α] Sin[θ], 
  Sqrt[(Cos[θ] Sin[α] + Cos[α] Cos[ϕ] Sin[θ])^2 + Sin[θ]^2 Sin[ϕ]^2]]
dist = TransformedDistribution[expr /. {θ -> 10 Degree, α -> 1 Degree}, ϕ \[Distributed] UniformDistribution[{-Pi, Pi}]];

(* Define a transformed beta distribution *)
upper = (expr /. {θ -> 10 Degree, α -> 1 Degree})/Degree /. ϕ -> 0;
lower = (expr /. {θ -> 10 Degree, α -> 1 Degree})/Degree /. ϕ -> π;
dist2 = TransformedDistribution[lower + (upper - lower) x, x \[Distributed] BetaDistribution[a, b]];

(* Take a large random sample from the original distribution *)
SeedRandom[12345];
y = RandomVariate[dist, 10000000]/Degree;

(* Determine values of a and b the result in the same mean and variance 
   as the sampled distribution, i.e., method of moments *)
sol = Solve[{Mean[y] == Mean[dist2], Variance[y] == Variance[dist2]}, {a, b}][[1]]
(* {a -> 0.5141003518067963`,b -> 0.48895483978530496`} *)

(* Show results *)
Show[Histogram[y, {0.01}, "PDF", PlotRange -> {{lower, upper}, {0, 2.5}}, ImageSize -> Large],
 Plot[PDF[dist2, x] /. sol, {x, lower, upper}, PlotRange -> {Automatic, {0, 2.5}}]]

histogram and fitted pdf

| improve this answer | |
$\endgroup$
  • $\begingroup$ Now that you mention it, it does resemble ArcSinDistribution[] somewhat. $\endgroup$ – J. M.'s technical difficulties Apr 2 '19 at 23:38
  • $\begingroup$ @J.M.isslightlypensive Yes, but from the random samples it doesn't seem quite symmetric. An ArcSinDistribution is a BetaDistribution[a,1-a] (I think). $\endgroup$ – JimB Apr 2 '19 at 23:42
  • 1
    $\begingroup$ @J.M.isslightlypensive Ooops! My memory is bad. I now looked it up and an arcsine distribution is a beta(1/2,1/2) and *generalized standard arcsine" is a beta(1-a,a). $\endgroup$ – JimB Apr 3 '19 at 2:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.