# Random Variate inside a function causing the function to behave unexpectedly

I define a Random Variable drawn from a custom probability density.

u[xi_, yi_, zi_, τ_, δt_, θ0_, k_] :=
RandomVariate[ProbabilityDistribution[(Sqrt[δt/τ]*(Cosh[(u*δt)/τ] +
(zi*Cos[θ0] + Sin[θ0]*(xi*Cos[k*δt] + yi*Sin[k*δt]))*Sinh[(u*δt)/τ]))/
(E^(((1 + u^2)*δt)/(2*τ))*Sqrt[2*Pi]), {u, -Infinity, Infinity}]]


I want to use this random variable inside the delayed function xfun:

xfun[xi_, yi_, zi_, θ0_, k_, δt_, τ_] := (xi - xi*Cos[2*θ0] - 4*Cos[2*k*Pi*δt]*Sin[θ0] +
2*xi*Cos[4*k*Pi*δt]*Sin[θ0]^2 + 2*yi*Sin[4*k*Pi*δt]*Sin[θ0]^2 +
4*E^((2*δt*u[xi, yi, zi, τ, δt, θ0, k])/τ)*Cos[2*k*Pi*δt]*Sin[θ0]*
(1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0] + 2*E^((δt*u[xi, yi, zi, τ, δt, θ0, k])/τ)*
(3*xi + xi*Cos[2*θ0] - 2*(xi*Cos[4*k*Pi*δt] + yi*Sin[4*k*Pi*δt])*
Sin[θ0]^2 - 2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0]))/(-4*(-1 + zi*Cos[θ0] +
xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
4*E^((2*δt*u[xi, yi, zi, τ, δt, θ0, k])/τ)*(1 + zi*Cos[θ0] +
xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]))


Now checking the special case (specified by the arguments of xfun) I find that I will sometimes get values greater than 1 for my function

xfun[1, 0, 0, Pi/2, 1, 0.01, 1]


This should be impossible since I can check, replacing the function u[,,,,] with the symbol u in the RHS of the function xfun, that the maximum value xfun should take for any u in this special case should be 1.

Maximize[(xi - xi*Cos[2*θ0] - 4*Cos[2*k*Pi*δt]*Sin[θ0] +
2*xi*Cos[4*k*Pi*δt]*Sin[θ0]^2 +
2*yi*Sin[4*k*Pi*δt]*Sin[θ0]^2 + 4*E^((2*δt*u)/τ)*Cos[2*k*Pi*δt]*Sin[θ0]*
(1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0] + 2*E^((δt*u)/τ)*(3*xi + xi*Cos[2*θ0] -
2*(xi*Cos[4*k*Pi*δt] + yi*Sin[4*k*Pi*δt])*Sin[θ0]^2 -  2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0]))/
(-4*(-1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] +  yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
4*E^((2*δt*u)/τ)*(1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] +
yi*Sin[2*k*Pi*δt]*Sin[θ0])) /. {θ0 -> Pi/2, k -> 1, δt -> 0.01, τ -> 1,
xi -> 1, yi -> 0, zi -> 0}, u]


Does anybody know what the source of the problem could be? Is there something about the behavior of the RandomVariate function that I am not understanding? Thanks!

In the definition of xfun you called u three times. Each time generated a different random variable. Use Module to structure xfun to call u only once and use that value all three times.

Clear["Global*"]

u[xi_, yi_, zi_, τ_, δt_, θ0_, k_] :=
RandomVariate[
ProbabilityDistribution[(Sqrt[δt/τ]*(Cosh[(u*δt)/τ] \
+ (zi*Cos[θ0] +
Sin[θ0]*(xi*Cos[k*δt] + yi*Sin[k*δt]))*
Sinh[(u*δt)/τ]))/(E^(((1 + u^2)*δt)/(2*τ))*
Sqrt[2*Pi]), {u, -Infinity, Infinity}]]


Modified definition of xfun

xfun[xi_, yi_, zi_, θ0_, k_, δt_, τ_] :=
Module[{urv = u[xi, yi, zi, τ, δt, θ0, k]},
(xi - xi*Cos[2*θ0] - 4*Cos[2*k*Pi*δt]*Sin[θ0] +
2*xi*Cos[4*k*Pi*δt]*Sin[θ0]^2 +
2*yi*Sin[4*k*Pi*δt]*Sin[θ0]^2 +
4*E^((2*δt*urv)/τ)*Cos[2*k*Pi*δt]*
Sin[θ0]*(1 + zi*Cos[θ0] +
xi*Cos[2*k*Pi*δt]*Sin[θ0] +
yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0] +
2*E^((δt*urv)/τ)*(3*xi + xi*Cos[2*θ0] -
2*(xi*Cos[4*k*Pi*δt] + yi*Sin[4*k*Pi*δt])*
Sin[θ0]^2 -
2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0]))/(-4*(-1 +
zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] +
yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
4*E^((2*δt*urv)/τ)*(1 + zi*Cos[θ0] +
xi*Cos[2*k*Pi*δt]*Sin[θ0] +
yi*Sin[2*k*Pi*δt]*Sin[θ0]))]


Testing,

SeedRandom[1234];

Max@Table[xfun[1, 0, 0, Pi/2, 1, 0.01, 1], {50}]

(* 1. *)


This is not an answer but just an observation. Occasionally there are errors such as the following (with just a slight modification of your function u):

u2[xi_, yi_, zi_, τ_, δt_, θ0_, k_, n_] :=
RandomVariate[ProbabilityDistribution[(Sqrt[δt/τ]*(Cosh[(u*δt)/τ] +
(zi*Cos[θ0] + Sin[θ0]*(xi*Cos[k*δt] + yi*Sin[k*δt]))*
Sinh[(u*δt)/τ]))/(E^(((1 + u^2)*δt)/(2*τ))*Sqrt[2*Pi]), {u, -Infinity, Infinity}], n]

x = u2[1, 0, 0, 1, 0.01, π/2, 1, 1000000];
`